РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ СТОХАСТИЧЕСКИХ ОБЪЕКТОВ С ИСПОЛЬЗОВАНИЕМ КРИТЕРИЯ МИНИМУМА ВЕРОЯТНОСТИ ОШИБКИ ОЦЕНИВАНИЯ
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ И СИСТЕМЫ
УДК 519.248: 681.5.001.3
С. В. СОКОЛОВ, П. А. КУЧЕРЕНКО
РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ
СТОХАСТИЧЕСКИХ ОБЪЕКТОВ С ИСПОЛЬЗОВАНИЕМ КРИТЕРИЯ МИНИМУМА ВЕРОЯТНОСТИ
ОШИБКИ ОЦЕНИВАНИЯ
Рассматривается проблема нелинейной параметрической идентификации стохастических объектов. Предлагается метод решения задачи идентификации параметра дискретного наблюдателя с использованием критерия минимума вероятности ошибки оценивания. Для иллюстрации эффективности предлагаемого подхода рассматривается численный пример.
Ключевые слова: параметрическая идентификация, ошибка оценивания, минимум критерия, фильтр Калмана.
Введение. Расширение области практического использования методов и алгоритмов стохастической параметрической идентификации обусловливает устойчивый рост интереса к развитию теории идентификации и разработке новых подходов к решению существующих в данной сфере проблем. Как показывает анализ публикаций, большинство работ посвящено вопросам, связанным с идентификацией параметров модели наблюдаемого стохастического объекта, в то время как вопросы, касающиеся определения параметров измерителя (наблюдателя) вектора состояния объекта, остаются практически не освещенными [1—9]. При этом задачи подобного рода часто возникают в различных областях телекоммуникации и связи, радионавигации, метрологии и др. К числу наиболее распространенных можно отнести задачи определения характеристик тракта передачи (времени распространения сигналов, параметров самой среды передачи и др.), а также определения коэффициентов усиления аппаратуры приемника.
В настоящей статье предлагается метод решения задачи параметрической идентификации, позволяющий, во-первых, освободиться от присущих известным методам ограничений (таких, как линейность модели измерителя относительно параметров, необходимость нормального вида распределения аддитивных шумов объекта и помех наблюдаемых сигналов и пр.), а во-вторых, повысить потенциальную точность процедуры идентификации за счет использования обобщенного вероятностного критерия, зависящего в общем случае нелинейно от апостериорной плотности распределения вероятности вектора состояния.
Следует отметить при этом, что вопросы идентификации по критериям, основанным на минимизации среднеквадратического отклонения ошибки оценивания, рассматривались ранее в работах [3—6]. Однако эти критерии в силу неравенства Чебышева являются (как это будет видно из последующих рассуждений) лишь частными случаями обобщенного нелинейного критерия, составляющего основу предлагаемого в настоящей статье подхода.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
6 С. В. Соколов, П. А. Кучеренко
Для упрощения изложения остановимся подробнее на скалярных уравнениях, определяющих нелинейную модель наблюдаемого объекта и структуру его измерителя.
Постановка задачи параметрической идентификации. Пусть дискретный объект задан нелинейным разностным уравнением
xk = f (xk−1) + n , x1 = x(1) ,
(1)
где n — возмущающее воздействие (шум) с известной функцией плотности распределения
вероятности q(n) ; xk — переменная состояния объекта в k-й момент времени; f — извест-
ная нелинейная функция, x1 — значение переменной состояния объекта в начальный момент времени.
Наблюдение за переменными состояния в дискретном времени осуществляется измерителем, описываемым также нелинейным (как относительно параметра наблюдателя, так и относительно переменной состояния объекта) уравнением следующего вида:
zk = χ(c, xk ) + w ,
(2)
где zk — дискретный отсчет сигнала измерителя; c — искомый параметр измерителя; w — шум измерителя с известной функцией плотности распределения вероятности g(w) ; χ — из-
вестная нелинейная функция наблюдения. При этом отметим, что идентифицируемость параметров измерителя обеспечивается
[10] единственностью апостериорной плотности распределения вероятности (АПРВ) вектора состояния наблюдаемого объекта, содержащей максимально полную информацию о совокупности переменных его состояния и параметрах измерителя.
Для упрощения совокупности дискретных отсчетов сигнала измерителя zi , i=1…k, обо-
значим через z1k . В рассматриваемом общем нелинейном стохастическом случае задача идентификации
неизвестного параметра c может быть сформулирована как задача нахождения его значения, доставляющего оптимум некоторому обобщенному вероятностному критерию J , зависящему от АПРВ вектора состояния. В качестве условия оптимизации (минимизации) критерия J используем далее условие минимума апостериорной плотности распределения вероятности текущей ошибки оценивания σ переменных состояния объекта на выбранном интервале ее
предельно допустимого изменения — от σmin до σmax , т.е.
σmax
∫min J = min cc
ρ(σk | z1k )dσk ,
σmin
где σk = xk − xˆk — текущая ошибка оценивания, xˆk — текущая оценка переменной состояния объекта; ρ(σk | z1k ) — АПРВ ошибки оценивания.
Учитывая линейную зависимость значений ошибки σk и переменной состояния xk , выразим АПРВ ошибки оценивания ρ(σk | z1k ) через АПРВ переменной состояния p(xk | z1k ) (выражение для которой будет получено ниже):
ρ(σk | z1k ) = p(σk + xˆk | z1k ) .
В этом случае минимизация критерия может быть представлена следующим образом:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
Решение задачи нелинейной параметрической идентификации
7
σmax
σmax
∫ ∫min J = min cc
ρ(σk
| z1k )dσk
= min
c
p(σk + xˆk | z1k )d σk .
σmin
σmin
(3)
В результате поставленная задача сводится к нахождению АПРВ p(σk + xˆk | z1k ) и по-
следующему определению значения искомого параметра из условия минимума критериально-
го выражения в формуле (3).
Синтез алгоритма нелинейной параметрической идентификации. Для определения
АПРВ p(σk + xˆk | z1k ) предварительно используем выражение для АПРВ p(xk | z1k ) с после-
дующей соответствующей заменой переменных. Известно [11], что АПРВ информационного параметра x для k-го момента времени оп-
ределяется выражением
∞
∫ p(xk−1 | z1k−1) p(xk | xk−1)dxk−1 p(zk | xk )
p(xk | z1k ) = −∞
h
, (4)
где
∞∞
∫ ∫h = p(xk−1 | z1k−1) p(xk | xk−1)dxk−1 p(zk | xk )dxk .
−∞ −∞
Условная плотность вероятности p(xk | xk−1) в формуле (4) может быть определена из
уравнения (1) при известном виде плотности распределения вероятности значений шума n (в предположении их взаимной статистической независимости):
p(xk | xk−1) = q(xk − f (xk−1)) .
Аналогичным образом из уравнения (2) можно определить и входящую в формулу (4) функцию правдоподобия:
p(zk | xk ) = g(zk − χ(c, xk )) .
Так как АПРВ p(xk−1 | z1k−1) в равенстве (4) является известной функцией, рекуррентный алгоритм определения АПРВ переменной состояния для k-го момента времени при наличии дискретных отсчетов сигнала z1k принимает следующий вид:
∞
∫ p(xk−1 | z1k−1)q(xk − f (xk−1))dxk−1g(zk − χ(c, xk ))
p(xk | z1k ) = −∞
h* (c)
,
(5)
где
∞∞
∫ ∫h*(c) =
p(xk−1 | z1k−1)q(xk − f (xk−1))dxk−1g(zk − χ(c, xk ))dxk .
−∞ −∞
Произведя соответствующую замену переменных в уравнении (5) и обозначив критериальное выражение через Ω(c) , задачу поиска минимума критерия (3) можно представить в
следующем виде:
min J = min Ω(c) ,
cc
(6)
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
8 С. В. Соколов, П. А. Кучеренко
где
σmax
∫Ω(c) = p(σk + xˆk | z1k )dσk = σmin
∫ ∫=
σmax σmin
⎛ ⎜ ⎜ ⎜ ⎜
∞ −∞
p( xk −1
|
z1k −1 )q((σk
+
xˆk
)
− f (xk−1))dxk−1g(zk h* (c)
− χ(c, σk
+
⎞
xˆk )) ⎟
⎟ ⎟
dσk
⎟
.
⎜⎟
⎝⎠
Здесь важно отметить, что в общем случае решения поставленной задачи оценка пере-
менной состояния xˆk , входящая в формулу (6), представляет собой некоторый функционал
(оператор) L от апостериорной плотности распределения вероятности переменной состояния,
( )т.е. xˆk = L p(xk | z1k ) , и, следовательно, в силу выражения (1) является нелинейной функци-
ей от искомого параметра c : xˆk = U (c) .
Тогда критериальное выражение в уравнении (6) окончательно можно представить в следующем обобщенном виде:
∫ ∫ ( ) ( )Ω(c)
=
σmax
⎛ ⎜ ⎜ ⎜
∞ −∞
p( xk −1
|
z1k −1)q
σmin ⎜
(σk + U (c)) − f (xk−1)
h* (c)
dxk −1 g
zk − χ(c,σk + U (с))
⎞
⎟
⎟ ⎟
dσk
.
⎟
(7)
⎜⎟
⎝⎠
Идентификация параметра, удовлетворяющего условию минимума критерия (6), предполагает минимизацию полученного критериального выражения (7). Для этой цели в зависимости от конкретного вида получаемой функции Ω(c) можно использовать известные и ши-
роко применяемые методы оптимизации: градиентный, метод Ньютона, метод сопряженных направлений, различные прямые методы и др. Выбор метода определяется особенностями наблюдаемого объекта и его измерителя.
Многомерный метод нелинейной параметрической идентификации. Минимизацию критерия (6) можно легко обобщить для многомерного случая, когда уравнения объекта и наблюдателя описываются следующими векторными уравнениями:
xk = f (xk−1) + n ; zk = χ(C, xk ) + w ,
где xk и xk−1 — n-мерные векторы переменных состояния в k-й и (k –1)-й моменты времени (шаги) соответственно; zk — m-мерный вектор сигналов измерителя; n — n-мерный вектор шума с известной n-мерной функцией плотности распределения вероятности q(n) ; C — вектор (или матрица) параметров наблюдателя; w — m-мерный вектор шума с известной m-мерной функцией плотности распределения вероятности g(w) .
Поскольку методика синтеза алгоритмов, соответствующих многомерному случаю, является полностью аналогичной изложенной выше для скалярного случая, приведем окончательную форму критериального выражения для рассматриваемого векторного случая:
σmax σmax
∫ ∫Ω(C) = ... p(σk + xˆ k | z1k )dσk = σmin σmin
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
Решение задачи нелинейной параметрической идентификации
9
∫ ∫ ∫ ∫=
⎛
σ
max
σmax
...
⎜ ⎜ ⎜
⎜σmin σmin
∞∞
...
−∞ −∞
p (x k −1
|
z1k −1 )q((σ k
+ xˆ k ) − f (xk−1))dxk−1g(zk h* (C)
− χ(C,σk
+
xˆ k
⎞ )) ⎟
⎟ ⎟
dσ
k
;
⎟
⎜⎟
⎝⎠
∫ ∫ ∫ ∫h*(C)
=
∞ ∞⎛∞ ∞
−∞
...
−∞
⎜⎜⎝
−∞
...
−∞
p(x k −1
| z1k−1)q((σk
+ xˆ k ) −
f (xk−1))dxk−1g(zk
− χ(С, σk
⎞ + xˆ k )) ⎟⎟⎠ dxk ,
где xˆ k — вектор оценок переменных состояния; z1k — совокупность векторов сигналов
zi , i=1...k; σk — вектор ошибки оценивания.
Эффективность использования предложенного подхода проиллюстрируем на следующем примере.
Пример. Рассмотрим стохастический дискретный объект, заданный нелинейным разностным уравнением
xk = 3xk−1 − (xk−1)2 + n , x1 = 1,
(8)
где n — белый гауссов шум с дисперсией Dn = 0, 02 и нулевым средним.
Наблюдение за переменными состояния объекта осуществляется измерителем, описы-
ваемым следующим нелинейным уравнением:
zk = c(2xk − 0,5(xk )2 ) + w ,
(9)
где c — искомый параметр наблюдателя (для рассматриваемого примера выберем исходное
значение этого параметра c = 1,5 ); w — белый гауссов шум с нулевым средним и дисперсией
Dw = 0, 25 .
Определение оценок произведем с использованием рекуррентного алгоритма калмановской фильтрации. Линеаризованные уравнения (8) и (9) принимают при этом следующий вид:
xk = Ak xk−1 + Bk + n , x1 = 1; zk = c(Ek xk + Pk ) + w ,
где Ak = 3 − 2xˆk−1 , Bk = (xˆk−1)2 , Ek = 2 − xˆk−1 , Pk = 0,5(xˆk−1)2 — коэффициенты, полученные
в результате линеаризации функций f (xk−1) и χ(c, xk ) в окрестностях оценок переменной состояния объекта для (k –1)-го шага.
Для определения текущего значения оценки переменной состояния объекта использовался оптимальный фильтр Калмана [11], который в рассматриваемом случае определяется как
xˆk = U (c) = Ak xˆk−1 + Bk + Kk ( zk − c(Ek ( Ak xˆk−1 + Bk ) + Pk )) , xˆ1 = 0,5 ;
Kk
=
c
Rk Dw
,
Rk
⎛1 = ⎜⎝⎜ Rk−1 + Dn
+
c2 Dw
⎞−1 ⎟⎠⎟
,
R1
= 0,5 .
Критериальное выражение Ω(c) на k-м шаге алгоритма для данного примера примет
вид
∫σmax
Ω(c) =
∫ ∫σmin
p(σk
+ xˆk
|
z1k )d σk
⎛
⎜
σmax ⎜
=
σmin
⎜ ⎜ ⎜
Dn
1 Dw
2π
∞ −∞
p( xk −1
|
z1k
−1
)e
−
τ2 2 Dn
−
dxk −1e
υ2 2Dw
h* (c)
⎞
⎟
⎟
⎟ ⎟
dσk
,
⎟
(10)
⎜⎟
⎝⎠
τ = σk + xˆk − 3xk−1 + (xk−1)2 , υ = zk − c(2(σk + xˆk ) − 0,5(σk + xˆk )2 ) .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
10 С. В. Соколов, П. А. Кучеренко
АПРВ для первой итерации алгоритма выберем нормальной с дисперсией D0 = 0,5 и нулевым математическим ожиданием. При этом отклонение среднего значения апостериорной плотности от начального значения переменной состояния не оказывает в дальнейшем существенного влияния на качество процедуры идентификации.
На рис. 1 представлен полученный в результате моделирования график входящей в формулу (10) АПРВ текущей ошибки оценивания для k-го шага алгоритма (k=50), которая, являясь функцией текущей ошибки σk , зависит также и от значений искомого параметра c :
( p(σk + xˆk | z1k ) = p(σk + U (c) | z1k ) = V (c, σk ) ).
0
1
V(c, σk)
2
3 с, о.е. k
40
20
0 0 –0,5
–1 σk
Рис. 1
Определение интегралов в выражении (10) производилось численно с использованием квадратурных формул с шагом ∆ = 0,03 . Бесконечные пределы интегрирования по переменной состояния x были заменены на конечные значения, удовлетворяющие точностным требованиям к алгоритму оценки ( xmin = 0 , xmax = 4 ).
На рис. 2 приведен график зависимости функции критериального выражения Ω(c) от искомого параметра при k=50. Границы интервала интегрирования по текущей ошибке оценивания также выбирались исходя из требований, обеспечивающих необходимую точность алгоритма оценки ( σmin = −1,5 , σmax = 0, 2 ).
Ω(с), о.е.
0,4
0,3
0,2 0 1 2 3 4 с, о.е. Рис. 2.
Как показали результаты моделирования, вид приведенной на рис. 2 зависимости является характерным для критериальных выражений, получаемых на различных итерациях алгоритма (7).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
Решение задачи нелинейной параметрической идентификации
11
Здесь важно отметить, что, являясь многоэкстремальными, критериальные выражения на различных шагах алгоритма (см. рис. 2) принимают наименьшие значения при c = 2 .
Для минимизации функции критериального выражения на очередном шаге алгоритма, т.е. для однозначного определения численного значения текущей оценки искомого параметра, целесообразно, задав некоторый интервал возможных значений параметра c (в рассмотренном примере 0 ≤ c ≤ 5 ), применить один из методов прямой минимизации. В данном случае использовался модифицированный симплексный метод прямой минимизации Нелдера — Мида — метод Бокса, обладающий достаточной вычислительной эффективностью и удобной программной реализацией [12].
Результаты компьютерного моделирования процедуры нелинейной параметрической идентификации показали, что если количество дискретных значений сигнала измерителя больше 250, то отклонение оценки параметра наблюдателя от его истинного значения c = 2 не превышает 9,7 %.
Таким образом, результаты проведенных исследований подтверждают принципиальную возможность эффективной реализации метода нелинейной параметрической идентификации с использованием критерия минимума АПРВ текущей ошибки оценивания. При этом важно отметить, что упомянутые выше методы параметрической идентификации [3—6] в изложенной постановке рассмотренную задачу решить не позволяют.
Заключение. Полученное выражение (7) определяет самый общий вид алгоритма нелинейной параметрической идентификации, обладающего рядом принципиально новых свойств. К их числу следует отнести:
— более высокий по сравнению с традиционными методами уровень потенциальной точности процесса идентификации благодаря использованию обобщенных вероятностных критериев, зависящих в общем случае нелинейно от апостериорной плотности распределения вероятности вектора состояния объекта;
— инвариантность к виду плотности распределения вероятности шума как объекта, так и измерителя;
— возможность применения метода для нелинейных объектов и наблюдателей, в том числе, при нелинейной зависимости функции наблюдения от искомого параметра.
Таким образом, предложенный метод нелинейной параметрической идентификации на основе обобщенных вероятностных критериев может быть весьма эффективно использован в различных областях связи, управления, метрологии и т.п.
СПИСОК ЛИТЕРАТУРЫ
1. Гроп Д. Методы идентификации систем. М.: Мир, 1979.
2. Льюнг Л. Идентификация систем. Теория для пользователя. М.: Наука, 1991.
3. Эйкхофф П. Основы идентификации систем управления. М.: Мир, 1975.
4. Сейдж Э., Мелса Дж. Идентификация систем управления. М.: Мир, 1974.
5. Пугачев В. С., Казаков И. Е., Евланов Л. Г. Основы статистической теории автоматических систем. М.: Машиностроение, 1974.
6. Пащенко Ф. Ф. Введение в состоятельные методы моделирования систем. Идентификация нелинейных систем. М.: Финансы и статистика, 2007.
7. Петров Б. Н., Уланов Г. М., Гольденблат И. И., Ульянов С. В. Теория моделей в процессах управления. М.: Наука, 1978.
8. Справочник по теории автоматического управления / Под ред. А. А. Красовского. М.: Наука, 1987.
9. Штейнберг Ш. Е. Идентификация в системах управления. М.: Энергоатмоиздат, 1987.
10. Пугачев В. С., Синицын И. Н. Стохастические дифференциальные системы. М.: Наука, 1985.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
12 Ю. В. Лужков
11. Тихонов В. И., Харисов В. Н. Статистический анализ и синтез радиотехнических устройств и систем. М.: Радио и связь, 1991.
12. Банди Б. Методы оптимизации: Вводный курс. М.: Радио и связь, 1988.
Сергей Викторович Соколов Павел Александрович Кучеренко
Сведения об авторах — д-р техн. наук, профессор; Ростовский государственный университет
путей сообщения, кафедра автоматики и телемеханики на ж.-д. транспорте, Ростов-на-Дону — аспирант; Ростовский государственный университет путей сообщения, кафедра автоматики и телемеханики на ж.-д. транспорте, Ростов-наДону; E-mail: pavelpost83@mail.ru
Рекомендована кафедрой автоматики и телемеханики на ж.-д. транспорте
Поступила в редакцию 21.05.08 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
УДК 519.248: 681.5.001.3
С. В. СОКОЛОВ, П. А. КУЧЕРЕНКО
РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ
СТОХАСТИЧЕСКИХ ОБЪЕКТОВ С ИСПОЛЬЗОВАНИЕМ КРИТЕРИЯ МИНИМУМА ВЕРОЯТНОСТИ
ОШИБКИ ОЦЕНИВАНИЯ
Рассматривается проблема нелинейной параметрической идентификации стохастических объектов. Предлагается метод решения задачи идентификации параметра дискретного наблюдателя с использованием критерия минимума вероятности ошибки оценивания. Для иллюстрации эффективности предлагаемого подхода рассматривается численный пример.
Ключевые слова: параметрическая идентификация, ошибка оценивания, минимум критерия, фильтр Калмана.
Введение. Расширение области практического использования методов и алгоритмов стохастической параметрической идентификации обусловливает устойчивый рост интереса к развитию теории идентификации и разработке новых подходов к решению существующих в данной сфере проблем. Как показывает анализ публикаций, большинство работ посвящено вопросам, связанным с идентификацией параметров модели наблюдаемого стохастического объекта, в то время как вопросы, касающиеся определения параметров измерителя (наблюдателя) вектора состояния объекта, остаются практически не освещенными [1—9]. При этом задачи подобного рода часто возникают в различных областях телекоммуникации и связи, радионавигации, метрологии и др. К числу наиболее распространенных можно отнести задачи определения характеристик тракта передачи (времени распространения сигналов, параметров самой среды передачи и др.), а также определения коэффициентов усиления аппаратуры приемника.
В настоящей статье предлагается метод решения задачи параметрической идентификации, позволяющий, во-первых, освободиться от присущих известным методам ограничений (таких, как линейность модели измерителя относительно параметров, необходимость нормального вида распределения аддитивных шумов объекта и помех наблюдаемых сигналов и пр.), а во-вторых, повысить потенциальную точность процедуры идентификации за счет использования обобщенного вероятностного критерия, зависящего в общем случае нелинейно от апостериорной плотности распределения вероятности вектора состояния.
Следует отметить при этом, что вопросы идентификации по критериям, основанным на минимизации среднеквадратического отклонения ошибки оценивания, рассматривались ранее в работах [3—6]. Однако эти критерии в силу неравенства Чебышева являются (как это будет видно из последующих рассуждений) лишь частными случаями обобщенного нелинейного критерия, составляющего основу предлагаемого в настоящей статье подхода.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
6 С. В. Соколов, П. А. Кучеренко
Для упрощения изложения остановимся подробнее на скалярных уравнениях, определяющих нелинейную модель наблюдаемого объекта и структуру его измерителя.
Постановка задачи параметрической идентификации. Пусть дискретный объект задан нелинейным разностным уравнением
xk = f (xk−1) + n , x1 = x(1) ,
(1)
где n — возмущающее воздействие (шум) с известной функцией плотности распределения
вероятности q(n) ; xk — переменная состояния объекта в k-й момент времени; f — извест-
ная нелинейная функция, x1 — значение переменной состояния объекта в начальный момент времени.
Наблюдение за переменными состояния в дискретном времени осуществляется измерителем, описываемым также нелинейным (как относительно параметра наблюдателя, так и относительно переменной состояния объекта) уравнением следующего вида:
zk = χ(c, xk ) + w ,
(2)
где zk — дискретный отсчет сигнала измерителя; c — искомый параметр измерителя; w — шум измерителя с известной функцией плотности распределения вероятности g(w) ; χ — из-
вестная нелинейная функция наблюдения. При этом отметим, что идентифицируемость параметров измерителя обеспечивается
[10] единственностью апостериорной плотности распределения вероятности (АПРВ) вектора состояния наблюдаемого объекта, содержащей максимально полную информацию о совокупности переменных его состояния и параметрах измерителя.
Для упрощения совокупности дискретных отсчетов сигнала измерителя zi , i=1…k, обо-
значим через z1k . В рассматриваемом общем нелинейном стохастическом случае задача идентификации
неизвестного параметра c может быть сформулирована как задача нахождения его значения, доставляющего оптимум некоторому обобщенному вероятностному критерию J , зависящему от АПРВ вектора состояния. В качестве условия оптимизации (минимизации) критерия J используем далее условие минимума апостериорной плотности распределения вероятности текущей ошибки оценивания σ переменных состояния объекта на выбранном интервале ее
предельно допустимого изменения — от σmin до σmax , т.е.
σmax
∫min J = min cc
ρ(σk | z1k )dσk ,
σmin
где σk = xk − xˆk — текущая ошибка оценивания, xˆk — текущая оценка переменной состояния объекта; ρ(σk | z1k ) — АПРВ ошибки оценивания.
Учитывая линейную зависимость значений ошибки σk и переменной состояния xk , выразим АПРВ ошибки оценивания ρ(σk | z1k ) через АПРВ переменной состояния p(xk | z1k ) (выражение для которой будет получено ниже):
ρ(σk | z1k ) = p(σk + xˆk | z1k ) .
В этом случае минимизация критерия может быть представлена следующим образом:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
Решение задачи нелинейной параметрической идентификации
7
σmax
σmax
∫ ∫min J = min cc
ρ(σk
| z1k )dσk
= min
c
p(σk + xˆk | z1k )d σk .
σmin
σmin
(3)
В результате поставленная задача сводится к нахождению АПРВ p(σk + xˆk | z1k ) и по-
следующему определению значения искомого параметра из условия минимума критериально-
го выражения в формуле (3).
Синтез алгоритма нелинейной параметрической идентификации. Для определения
АПРВ p(σk + xˆk | z1k ) предварительно используем выражение для АПРВ p(xk | z1k ) с после-
дующей соответствующей заменой переменных. Известно [11], что АПРВ информационного параметра x для k-го момента времени оп-
ределяется выражением
∞
∫ p(xk−1 | z1k−1) p(xk | xk−1)dxk−1 p(zk | xk )
p(xk | z1k ) = −∞
h
, (4)
где
∞∞
∫ ∫h = p(xk−1 | z1k−1) p(xk | xk−1)dxk−1 p(zk | xk )dxk .
−∞ −∞
Условная плотность вероятности p(xk | xk−1) в формуле (4) может быть определена из
уравнения (1) при известном виде плотности распределения вероятности значений шума n (в предположении их взаимной статистической независимости):
p(xk | xk−1) = q(xk − f (xk−1)) .
Аналогичным образом из уравнения (2) можно определить и входящую в формулу (4) функцию правдоподобия:
p(zk | xk ) = g(zk − χ(c, xk )) .
Так как АПРВ p(xk−1 | z1k−1) в равенстве (4) является известной функцией, рекуррентный алгоритм определения АПРВ переменной состояния для k-го момента времени при наличии дискретных отсчетов сигнала z1k принимает следующий вид:
∞
∫ p(xk−1 | z1k−1)q(xk − f (xk−1))dxk−1g(zk − χ(c, xk ))
p(xk | z1k ) = −∞
h* (c)
,
(5)
где
∞∞
∫ ∫h*(c) =
p(xk−1 | z1k−1)q(xk − f (xk−1))dxk−1g(zk − χ(c, xk ))dxk .
−∞ −∞
Произведя соответствующую замену переменных в уравнении (5) и обозначив критериальное выражение через Ω(c) , задачу поиска минимума критерия (3) можно представить в
следующем виде:
min J = min Ω(c) ,
cc
(6)
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
8 С. В. Соколов, П. А. Кучеренко
где
σmax
∫Ω(c) = p(σk + xˆk | z1k )dσk = σmin
∫ ∫=
σmax σmin
⎛ ⎜ ⎜ ⎜ ⎜
∞ −∞
p( xk −1
|
z1k −1 )q((σk
+
xˆk
)
− f (xk−1))dxk−1g(zk h* (c)
− χ(c, σk
+
⎞
xˆk )) ⎟
⎟ ⎟
dσk
⎟
.
⎜⎟
⎝⎠
Здесь важно отметить, что в общем случае решения поставленной задачи оценка пере-
менной состояния xˆk , входящая в формулу (6), представляет собой некоторый функционал
(оператор) L от апостериорной плотности распределения вероятности переменной состояния,
( )т.е. xˆk = L p(xk | z1k ) , и, следовательно, в силу выражения (1) является нелинейной функци-
ей от искомого параметра c : xˆk = U (c) .
Тогда критериальное выражение в уравнении (6) окончательно можно представить в следующем обобщенном виде:
∫ ∫ ( ) ( )Ω(c)
=
σmax
⎛ ⎜ ⎜ ⎜
∞ −∞
p( xk −1
|
z1k −1)q
σmin ⎜
(σk + U (c)) − f (xk−1)
h* (c)
dxk −1 g
zk − χ(c,σk + U (с))
⎞
⎟
⎟ ⎟
dσk
.
⎟
(7)
⎜⎟
⎝⎠
Идентификация параметра, удовлетворяющего условию минимума критерия (6), предполагает минимизацию полученного критериального выражения (7). Для этой цели в зависимости от конкретного вида получаемой функции Ω(c) можно использовать известные и ши-
роко применяемые методы оптимизации: градиентный, метод Ньютона, метод сопряженных направлений, различные прямые методы и др. Выбор метода определяется особенностями наблюдаемого объекта и его измерителя.
Многомерный метод нелинейной параметрической идентификации. Минимизацию критерия (6) можно легко обобщить для многомерного случая, когда уравнения объекта и наблюдателя описываются следующими векторными уравнениями:
xk = f (xk−1) + n ; zk = χ(C, xk ) + w ,
где xk и xk−1 — n-мерные векторы переменных состояния в k-й и (k –1)-й моменты времени (шаги) соответственно; zk — m-мерный вектор сигналов измерителя; n — n-мерный вектор шума с известной n-мерной функцией плотности распределения вероятности q(n) ; C — вектор (или матрица) параметров наблюдателя; w — m-мерный вектор шума с известной m-мерной функцией плотности распределения вероятности g(w) .
Поскольку методика синтеза алгоритмов, соответствующих многомерному случаю, является полностью аналогичной изложенной выше для скалярного случая, приведем окончательную форму критериального выражения для рассматриваемого векторного случая:
σmax σmax
∫ ∫Ω(C) = ... p(σk + xˆ k | z1k )dσk = σmin σmin
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
Решение задачи нелинейной параметрической идентификации
9
∫ ∫ ∫ ∫=
⎛
σ
max
σmax
...
⎜ ⎜ ⎜
⎜σmin σmin
∞∞
...
−∞ −∞
p (x k −1
|
z1k −1 )q((σ k
+ xˆ k ) − f (xk−1))dxk−1g(zk h* (C)
− χ(C,σk
+
xˆ k
⎞ )) ⎟
⎟ ⎟
dσ
k
;
⎟
⎜⎟
⎝⎠
∫ ∫ ∫ ∫h*(C)
=
∞ ∞⎛∞ ∞
−∞
...
−∞
⎜⎜⎝
−∞
...
−∞
p(x k −1
| z1k−1)q((σk
+ xˆ k ) −
f (xk−1))dxk−1g(zk
− χ(С, σk
⎞ + xˆ k )) ⎟⎟⎠ dxk ,
где xˆ k — вектор оценок переменных состояния; z1k — совокупность векторов сигналов
zi , i=1...k; σk — вектор ошибки оценивания.
Эффективность использования предложенного подхода проиллюстрируем на следующем примере.
Пример. Рассмотрим стохастический дискретный объект, заданный нелинейным разностным уравнением
xk = 3xk−1 − (xk−1)2 + n , x1 = 1,
(8)
где n — белый гауссов шум с дисперсией Dn = 0, 02 и нулевым средним.
Наблюдение за переменными состояния объекта осуществляется измерителем, описы-
ваемым следующим нелинейным уравнением:
zk = c(2xk − 0,5(xk )2 ) + w ,
(9)
где c — искомый параметр наблюдателя (для рассматриваемого примера выберем исходное
значение этого параметра c = 1,5 ); w — белый гауссов шум с нулевым средним и дисперсией
Dw = 0, 25 .
Определение оценок произведем с использованием рекуррентного алгоритма калмановской фильтрации. Линеаризованные уравнения (8) и (9) принимают при этом следующий вид:
xk = Ak xk−1 + Bk + n , x1 = 1; zk = c(Ek xk + Pk ) + w ,
где Ak = 3 − 2xˆk−1 , Bk = (xˆk−1)2 , Ek = 2 − xˆk−1 , Pk = 0,5(xˆk−1)2 — коэффициенты, полученные
в результате линеаризации функций f (xk−1) и χ(c, xk ) в окрестностях оценок переменной состояния объекта для (k –1)-го шага.
Для определения текущего значения оценки переменной состояния объекта использовался оптимальный фильтр Калмана [11], который в рассматриваемом случае определяется как
xˆk = U (c) = Ak xˆk−1 + Bk + Kk ( zk − c(Ek ( Ak xˆk−1 + Bk ) + Pk )) , xˆ1 = 0,5 ;
Kk
=
c
Rk Dw
,
Rk
⎛1 = ⎜⎝⎜ Rk−1 + Dn
+
c2 Dw
⎞−1 ⎟⎠⎟
,
R1
= 0,5 .
Критериальное выражение Ω(c) на k-м шаге алгоритма для данного примера примет
вид
∫σmax
Ω(c) =
∫ ∫σmin
p(σk
+ xˆk
|
z1k )d σk
⎛
⎜
σmax ⎜
=
σmin
⎜ ⎜ ⎜
Dn
1 Dw
2π
∞ −∞
p( xk −1
|
z1k
−1
)e
−
τ2 2 Dn
−
dxk −1e
υ2 2Dw
h* (c)
⎞
⎟
⎟
⎟ ⎟
dσk
,
⎟
(10)
⎜⎟
⎝⎠
τ = σk + xˆk − 3xk−1 + (xk−1)2 , υ = zk − c(2(σk + xˆk ) − 0,5(σk + xˆk )2 ) .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
10 С. В. Соколов, П. А. Кучеренко
АПРВ для первой итерации алгоритма выберем нормальной с дисперсией D0 = 0,5 и нулевым математическим ожиданием. При этом отклонение среднего значения апостериорной плотности от начального значения переменной состояния не оказывает в дальнейшем существенного влияния на качество процедуры идентификации.
На рис. 1 представлен полученный в результате моделирования график входящей в формулу (10) АПРВ текущей ошибки оценивания для k-го шага алгоритма (k=50), которая, являясь функцией текущей ошибки σk , зависит также и от значений искомого параметра c :
( p(σk + xˆk | z1k ) = p(σk + U (c) | z1k ) = V (c, σk ) ).
0
1
V(c, σk)
2
3 с, о.е. k
40
20
0 0 –0,5
–1 σk
Рис. 1
Определение интегралов в выражении (10) производилось численно с использованием квадратурных формул с шагом ∆ = 0,03 . Бесконечные пределы интегрирования по переменной состояния x были заменены на конечные значения, удовлетворяющие точностным требованиям к алгоритму оценки ( xmin = 0 , xmax = 4 ).
На рис. 2 приведен график зависимости функции критериального выражения Ω(c) от искомого параметра при k=50. Границы интервала интегрирования по текущей ошибке оценивания также выбирались исходя из требований, обеспечивающих необходимую точность алгоритма оценки ( σmin = −1,5 , σmax = 0, 2 ).
Ω(с), о.е.
0,4
0,3
0,2 0 1 2 3 4 с, о.е. Рис. 2.
Как показали результаты моделирования, вид приведенной на рис. 2 зависимости является характерным для критериальных выражений, получаемых на различных итерациях алгоритма (7).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
Решение задачи нелинейной параметрической идентификации
11
Здесь важно отметить, что, являясь многоэкстремальными, критериальные выражения на различных шагах алгоритма (см. рис. 2) принимают наименьшие значения при c = 2 .
Для минимизации функции критериального выражения на очередном шаге алгоритма, т.е. для однозначного определения численного значения текущей оценки искомого параметра, целесообразно, задав некоторый интервал возможных значений параметра c (в рассмотренном примере 0 ≤ c ≤ 5 ), применить один из методов прямой минимизации. В данном случае использовался модифицированный симплексный метод прямой минимизации Нелдера — Мида — метод Бокса, обладающий достаточной вычислительной эффективностью и удобной программной реализацией [12].
Результаты компьютерного моделирования процедуры нелинейной параметрической идентификации показали, что если количество дискретных значений сигнала измерителя больше 250, то отклонение оценки параметра наблюдателя от его истинного значения c = 2 не превышает 9,7 %.
Таким образом, результаты проведенных исследований подтверждают принципиальную возможность эффективной реализации метода нелинейной параметрической идентификации с использованием критерия минимума АПРВ текущей ошибки оценивания. При этом важно отметить, что упомянутые выше методы параметрической идентификации [3—6] в изложенной постановке рассмотренную задачу решить не позволяют.
Заключение. Полученное выражение (7) определяет самый общий вид алгоритма нелинейной параметрической идентификации, обладающего рядом принципиально новых свойств. К их числу следует отнести:
— более высокий по сравнению с традиционными методами уровень потенциальной точности процесса идентификации благодаря использованию обобщенных вероятностных критериев, зависящих в общем случае нелинейно от апостериорной плотности распределения вероятности вектора состояния объекта;
— инвариантность к виду плотности распределения вероятности шума как объекта, так и измерителя;
— возможность применения метода для нелинейных объектов и наблюдателей, в том числе, при нелинейной зависимости функции наблюдения от искомого параметра.
Таким образом, предложенный метод нелинейной параметрической идентификации на основе обобщенных вероятностных критериев может быть весьма эффективно использован в различных областях связи, управления, метрологии и т.п.
СПИСОК ЛИТЕРАТУРЫ
1. Гроп Д. Методы идентификации систем. М.: Мир, 1979.
2. Льюнг Л. Идентификация систем. Теория для пользователя. М.: Наука, 1991.
3. Эйкхофф П. Основы идентификации систем управления. М.: Мир, 1975.
4. Сейдж Э., Мелса Дж. Идентификация систем управления. М.: Мир, 1974.
5. Пугачев В. С., Казаков И. Е., Евланов Л. Г. Основы статистической теории автоматических систем. М.: Машиностроение, 1974.
6. Пащенко Ф. Ф. Введение в состоятельные методы моделирования систем. Идентификация нелинейных систем. М.: Финансы и статистика, 2007.
7. Петров Б. Н., Уланов Г. М., Гольденблат И. И., Ульянов С. В. Теория моделей в процессах управления. М.: Наука, 1978.
8. Справочник по теории автоматического управления / Под ред. А. А. Красовского. М.: Наука, 1987.
9. Штейнберг Ш. Е. Идентификация в системах управления. М.: Энергоатмоиздат, 1987.
10. Пугачев В. С., Синицын И. Н. Стохастические дифференциальные системы. М.: Наука, 1985.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3
12 Ю. В. Лужков
11. Тихонов В. И., Харисов В. Н. Статистический анализ и синтез радиотехнических устройств и систем. М.: Радио и связь, 1991.
12. Банди Б. Методы оптимизации: Вводный курс. М.: Радио и связь, 1988.
Сергей Викторович Соколов Павел Александрович Кучеренко
Сведения об авторах — д-р техн. наук, профессор; Ростовский государственный университет
путей сообщения, кафедра автоматики и телемеханики на ж.-д. транспорте, Ростов-на-Дону — аспирант; Ростовский государственный университет путей сообщения, кафедра автоматики и телемеханики на ж.-д. транспорте, Ростов-наДону; E-mail: pavelpost83@mail.ru
Рекомендована кафедрой автоматики и телемеханики на ж.-д. транспорте
Поступила в редакцию 21.05.08 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 3