CC BY
40
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОПТИМАЛЬНОЕ УПРАВЛЕНИЕ
УДК 62-50
К ТЕОРИИ ОЦЕНИВАНИЯ СОСТОЯНИЙ ДИНАМИЧЕСКОГО ОБЪЕКТА © 2010 г. Н.Я. Коган, Ю.И. Неймарк, Л.В. Коган
Нижегородский госуниверситет им. Н.И. Лобачевского [email protected]
Поступила в редакцию 16.02.2010
На основе байесова подхода получены рекуррентные формулы оптимальной оценки состояния линейного нестационарного дискретного динамического объекта при аппроксимации в виде гауссовых сумм законов распределения начального состояния объекта и действующих на него возмущений.
Ключевые слова: состояние динамического объекта, фильтр Калмана-Бьюси, гауссовы суммы.
Рассматриваются динамический объект и измерительная система, описываемые разностными уравнениями:
= А х, + ВЛ*
/к (х — ц), где ц - математическое ожидание случайного вектора х, К - корреляционная
матрица, то есть
, ^ = 0,1,2,..., (1)
У,+1 = с*+1 х*+1 +п где х^ - п -мерный вектор состояния объекта в 5 -й момент времени, у, - т -мерный вектор измерения выхода объекта, 2^ и п5 - п -мерный и т -мерный векторы внешнего воздействия на объект и ошибок в измерителе соответственно. Матрицы А, и Б,, - квадратные размерности п х п . Матрица С имеет размерность т х п и вместе с Ах удовлетворяет условию наблюдаемости. Вектор п, ошибок измерителя будем считать гауссовым с нулевым математическим ожиданием Мп =0 и корреляционной матрицей М(пГпГ ) = ^ • Векторы П, ^ = 0,1,2, к, для разных моментов времени некоррелированы, то есть М(пП) = 0, 5 = I. Плотность вероятностей случайного вектора 2 ^ внешних воздействий на объект будем аппроксимировать взвешенной суммой нормальных законов в соответствии с формулой
^) = X*/к . * * -I/, *) • (2)
7 = 1 ^
Здесь и в дальнейшем будем обозначать плотность гауссова распределения в виде
= (2п) 2 |К7>| 2 ехр^- -2 (^ -%^)Т :
х К
%
Весовые коэффициенты qj ^ в (2) неотрица-
ь
тельны qJ 2 > 0 и нормированы л =1 •
У=1
Внешние воздействия на объект в разные моменты времени, то есть £,^ и 2,/, 5 Ф I, некоррелированы М{$,^ 2,^ ) = 0 при 5 Ф I •
Начальное состояние объекта будем считать случайным с распределением, плотность которого аппроксимируется взвешенной суммой нормальных законов
а / \
ё х0 (х0 ) = Хр, 0 /Кг.0 {х0 - '€г,0 ), (3)
г=1
а
где р10 > 0 и ^РI 0 =1. Кроме того, будем
I=1
считать, что случайные воздействия на объект и ошибки измерителя некоррелированы и не зависят от начального состояния объекта.
Частный случай такой задачи для одномерного объекта при аппроксимации гауссовой суммой закона распределения начального состояния исследован в работе [1].
Вывод рекуррентных формул оценивания аналогичен приведенному в [2] и состоит в отыскании полного статистического описания состояния динамического объекта как случайного вектора при учёте известных к данному моменту измерений выхода.
Найденная плотность условного распределения вероятности вектора х,+1 (я = 0,1, 2, к)
представляет собой взвешенную сумму гауссовых функций
ёхэ +1/7.5+1 (Х+^ _
^+1 Л (4)
_ X рк, 51+1/Кт. _+1 (х+1 ~^к, 51+1) ’ к=1
где У,+1 = (у1,у2у,+1) - вектор значений выхода, измеренных к моменту 5 + 1. Коэффи-
где
-^+1
циенты р ^ > 0 и £ рк ^+1 =1. Число Л.
к=1
ііІ+1 = аЬї+ 5 = 0,1,2,
Основной результат работы состоит в приводимых ниже формулах определения оптимальной оценки состояния объекта
■^+1
к=1
Кк „5+1 _ Рк „5+1СГ+А+1,
гк,з = Рі,зЧ],з ,
( - \
Г+1 = X Гк/ъ, Iу,5+1 - Сэ+12к к=1 М
•/йі. . +1 Cs+1^zk,,5 ) _
= (2п)-т2 К У2
к , і 1 1
Л+1 -Сі+12к У°к 1 (у^+1 -Сі +12к )
®к, s - Rs+1 +
гк,s = Аs + ,
С+1 (лхр. , АА + В^К7. ,В*С^). В формулах (7)-(10) индексы г, у,к связаны соотношением
к = (г - 1)Ь + у, г = 1,(!,,+1, у = 1,Ь .
По заданному к значения г и у определяются следующим образом: к
і = і
гауссовых составляющих в (4) при условии Ь >1, где Ь - число гауссовых составляющих в (2), увеличивается со временем по показательному закону
] =
(5)
к -
Ь,
+1, если к = аЬ если к = аЬ а =1, 2,к,аЪ8; к Ь, если к = аЬ
если к = а Ь
(6)
где векторы Хк ^+1 находятся из рекуррентных уравнений
4 +1 = А Ъл + в* \] +
+ Кк^+1[у+1 - Сэ+1(Д?Хг> + В)], (7)
5 = 0,1,2,к;
(8)
Рк ,5+1 0*к ,5
- бкХ+1 (С.^А + *5+1 )-1С5+1бк,5 , (9)
^ = 0,1,2,...;
бк,* = АР,* А + BsК],sВI,
а весовые коэффициенты в (6) определяются рекуррентными формулами
Рк,5+1 = {у+1 - 1 гк^), (10)
^ ,5+1
а = 1,2,к, а6“ (11)
где знак [ ] означает целую часть числа, стоящего в скобках.
Формулы (6)-( 11) и есть оптимальный фильтр состояния динамического объекта (1). Его можно трактовать как совокупность независимых фильтров Калмана-Бьюси, формулы (7)-(9), суперпозиция выходов которых, формула (6), с коэффициентами, определяемыми формулами (10), представляет выход искомого фильтра.
Каждый фильтр Калмана-Бьюси определяет оптимальную оценку .€^^+1 состояния объекта
(1) по отношению к I -й (/ =1,2,...,) гауссовой составляющей функции плотности вероятности состояния х как случайного вектора и у -й (у = 1,2, к, Ь) гауссовой составляющей функции плотности вероятности случайного вектора ^^ внешнего воздействия на объект. Связь между индексами г, у,к задаётся формулами (11).
Коэффициенты рк 5+1, определяемые по формуле (10), устанавливают вес выхода каж-
дого фильтра Калмана-Бьюси в формировании оптимальной оценки £^+1 состояния объекта.
Особенностью рассматриваемого алгоритма является то, что согласно (5) при Ь > 1 увеличивается число фильтров Калмана-Бьюси из-за роста числа гауссовых составляющих в функции плотности вероятности состояния как случайного вектора. Это обстоятельство является принципиальным препятствием для вычислительной реализации алгоритма по формулам (5)-(11). Свободным от этого препятствия является важный для приложений частный случай Ь = 1. При Ь = 1 внешнее воздействие на объект является гауссовым белым шумом. Благодаря этому не изменяется число гауссовых составляющих в функции плотности вероятности состояния как случайного вектора и, следовательно, не изменяется число фильтров Калма-на-Бьюси в алгоритме (5)-(10). В соответствии
с (5) и (11) у = 1, к = I, г = 1, а. Каждый из независимых фильтров Калмана-Бьюси отрабатывает оптимальную оценку X ^ состояния объекта по отношению к одной из гауссовых составляющих плотности вероятности вектора начального состояния.
В общем случае Ь > 1 число фильтров Кал-мана-Бьюси увеличивается столь сильно, что вычислительная реализация алгоритма оптимальной фильтрации по формулам (5)-(11) становится практически невозможной. Так, при Ъ = 2 за десять временных интервалов первоначальное число фильтров увеличивается примерно в 103 раз, а при Ь = 3 в 6 • 104 раз. Естественным выходом из этой ситуации прдставля-ется ограничение числа гауссовых составляющих в функции плотности вероятности текущего состояния как случайного вектора. Если, например, после одного или нескольких шагов по времени функция плотности вероятности «обостряется» за счет небольшого числа гауссовых составляющих, сумма весовых коэффициентов которых близка к единице, рк13 + рк2 8 +
+---+ Рк1 х >1 - е, где 8 - малое положительное
число, то естественно на следующем шаге сохранить только к1 фильтров, отвечающих этим
составляющим. Если «обострение» плотности вероятности происходит за счет большого числа гауссовых составляющих с близкими математическими ожиданиями, то можно использовать подход, состоящий в аппроксимации функции плотности вероятности текущего состояния как случайного вектора гауссовой суммой меньше-
го числа слагаемых. Один из возможных способов такой аппроксимации состоит в следующем. Пусть g1 (х) - исходная плотность вероятности. Суперпозицию гауссовых функций
I
XРк /ак (х - тк ) , аппр°ксимирукщих £1 (х) ,
к=1
ищем по следующей схеме:
* I * \
ёк+1 (х) = ёк (х) - Як/а* Iх - тк),
где
(к > > т*) = аг§ т 1п | (ёк (х) -
% >°к >тк -» .
- Чк /ок (х - тк ))2 А •
Процесс останавливается, когда
[ (ё1 (х) - Ч*/ * (х - Щ ))2 < 8 ,
-да ‘
где 8 - малое положительное число. Весовые коэффициенты суперпозиции нормируются
I
Рк = Чк 1X Чs.
5 = 0
Теперь перейдем к выводу представленного в формулах (6)-(11) алгоритма оптимальной фильтрации. Сначала сформулируем два утверждения, справедливость которых может быть легко проверена.
Рассмотрим два случайных вектора х и £,, плотности вероятности которых представляют суперпозицию гауссовых функций
а
ёX (х) = XР /Р. (Х - Х ) ,
1=1
ё* ©=к ч . /к. (ъ-1 .).
}=1 1
Утверждение 1. Случайные векторы 2 = Ах и ^2 = В2, где А и В - матрицы, распределены с плотностями вероятностей
а
ёz1 (21) = Хр/ар,аТ (21 - АХ1),
. =1 .
Ег2 (^2 )= XЧ] 1ВК ВТ (^2 - В\] ).
]=1 ]
Утверждение 2. Случайный вектор ю = 2-1 + ^2 распределен с плотностью вероятности
аЬ
ёш (®)= Хгк/зк (®-®к ), (12)
К=1
где гк = Рг Ч] , ®к = Лхг + В^] , бк = Ар А + + ВКВТ , к = (г - 1)Ь + у , г = 1, а, у = 1, Ь . Для заданного К значения 1 и у определяются следующими формулами: к
і = і
+ 1, если к = аЪ если к = аЪ а = 1,2, к, аЬ8;
7 =
к -
Ь,
Ь, если к Фа Ь
если к ф<хЬ
где
Рассмотрим еще один случайный вектор
у = Сх+1 - и • (20)
В соответствии с утверждением 2 и (15)-(18)
он распределен с плотностью вероятности
^+1
где
ук, і
_ +12к, і - у,
(21)
(22)
а = 1,2,к,аЪ.
Докажем справедливость формулы (4). Доказательство проведем по индукции путем непосредственного представления gx ^у 1 (х,,+1)
через функцию gх /у (х^) аналогичного вида,
используя для этого формулу Байеса.
Рассмотрим два независимых случайных вектора
2 = 4X + » и = у-п^ +1. (13)
Случайный вектор х5 распределён с плотностью вероятности
£** и* (Х) = Тї1*Ри /р* (Х - ), (14)
І+1
где коэффициенты р 2 неотрицательны и нормированы - р- > 0, р,5 = 1. Вектор у -
1 = 1
неслучайный.
В соответствии с утверждениями 1 и 2 вектор 2 распределен с плотностью вероятности вида (12):
^+1 , .
ёг (2)= Е гк,з/дк * (2 - 2к,* ) , (15)
к = 1
Вк = Сэ+хбк С я+1 + ^+1 • (23)
Согласно (1), (13) вектор 2 представляет собой состояние объекта х^+1 как случайный вектор, а функция (15) определяет плотность вероятностей gх ^у (х^+1) его априорного распределения, не учитывающего у5+1 - измерения выхода в (я +1) -й момент времени. Учёт этого измерения изменяет условное распределение вектора х ^+1.
Запишем формулу Байеса для случайных векторов г и V:
„ ёу/х (V)ёг (Ю
ё х/у(2)---------------тт---
ёу (v)
(24)
В (24) g(у) - гауссова плотность вероятности:
ёук М = +1 [V - (СХ+1г - у)] • (25)
Если в (13) за у взять измерение уз+1 выхода объекта, то согласно (1), (20) вектор V обращается в нуль и в левой части формулы (24) стоит апостериорная плотность вероятности состояния Х^+1 при условии известных измерений +1 =(^1, у2у+1) выхода объекта.
Подставим в (24) выражения входящих в неё функций, учитывая условие V = 0 :
ё г/у=0 ( —) =
( \^+1 -
/кс+1 (у+1 — Сх+1-) Е гк^ (- — -К^ ) (26)
К=1
(16)
бк, ^ = А рі, ^ аТ + вх к ^ , (17)
гк,* = Рі,* Чі,ц, +1 = ds ь, (18)
при этом индексы і, у,к связаны соотношения-
“^+1
г УЬк (Ув+1 С^+1—к)
к=1 к’"
Сначала преобразуем аргумент показательной функции, стоящей в числителе выражения (26):
ми (11) при і = 1, ^ и у = 1, Ь.
Вектор и имеет гауссово распределение с плотностью
ёы (М) = /ях + 1 (М - у) . (19)
(2 - 2к ^ )Т О- (2 - 2к ,$) +
+(У1+1 Св+1г) ^1+-1(^1+1 Сі+12).
(27)
Прибавляя и вычитая в (27) выражение ктЄ~Хк , где
G = от1 + С Кх С венства по компонентам вектора х^+1. Исполь-
, ^ Г-1 (28) зуя (33) , запишем это условие в виде
к = б*,*2*,* + Ч+Л+1 >’s+l, +1
преобразуем (27) к виду IУ к ,*+1 = г*+1. (38)
(z - G-1k)G(z - G-1к) + H - kTG-1к , (29) *=i
где
H zk, sQk,s zk, s + У s+1Rs+1 У s+1 • (30)
Подставляя (29), (30) в (26), умножая и деля 1
Подставляя (37) в (38) и приравнивая коэффициенты при одинаковых показательных функциях, получим
{я*+11 \Qk.sl\оI \1\вк,, | =1.
числитель на | в2 | и учитывая, что выражение Учитывая эт° условие, получим рекуррентные уравнения для пересчёта весовых коэффи-
(26) представляет апостериорную плотность веро-
циентов в (31) в виде (10):
ятности Х5+1 при условии известных измерений
выхода объекта, запишем (26) в виде (4): р,5+1 =-^- /в, (У5+1 - сх+\*к,5).
^ (X ) — 5+1 ’5
^+1^^+1 ^+1 Теперь преобразуем выражение (32) матема-
ds+1 , ^ ч (31) тических ожиданий гауссовых составляющих в
= X Рк+1^Рк 5+1+1 - хк,s+и ’ (31) к виду (7), (8). Используя (28), (32) и (16),
получим
к=1
где
s-i-lj п /'1—1 Xk ,s+1 Pk ,s+1(Qk,s zk, s + Cs+1Rs+1ys+1)
= G h, Pk +i =G , (32)
xk ,s+1~^ 4 ,s+1 - v , _ T
Y = Pk, s+1 (Qk ,s zk ,s + Cs+1Rs+1Cs+1zk, s —
_ I k,s +1
Рк,s+1 ~ ? — ct R-1 C z + CT R-1 y ) —
1 s+1 s+1fis +1^ s +1zk ,s ^s+1fis+1У s +1)
2 r | r і 2 v = zk, s + Pk, s+1Cs+1Rs+1(.ys+1 — Cs+1zk ,s ) —
Yк,s+1 = (2n) 2 RK,s 1 Rs+1 1 2 X
(34)
*IQ*,J2|G|2 Е 2 , + К г C (B - )1
d + Kk,s+1ІЛ+1 - Cs+1( Asxi,s + Bs - j,s ^
-і -i -i(H-hTG-1h) = As xi, s + Bs-j, s +
ds+1
Г„1 = E n,sfDk (>■ s+1 - Cs+1 zt,s). (35) Д"« получения уравне^ния пересчета корре-
k=i k’s ляционных матриц гауссовых составляющих в
Используя (28) и (30), преобразуем аргумент формуле (31) воспользуемся известным мат-
показательной функции в (34) к виду рттаьш тождеством
H - hTG-1h = (y+1 - Cs+1zk)T x (A-1 + BTC-1B)-1 = A - ABT (BABT + C)-1 BA
D~\ ( с z ) (36) для обращения матрицы G . В соответствии с
k ys+1 s+1 k ,s )■ (28), (32) получим рекуррентное уравнение для
Тогда выражение для уk,s+1 может быть за- отыскания Pk,s+1 в виде (9).
писано в виде:
У к ,s+l = Гк ,s I Rs+1 I 2 I Qk ,s I 2 X
i 1 (37) 1. Sorenson H.W. and Alspach D.L. Recursive
X | g I 21 d 12 f (y — c z ) Bayesian estimation using Gaussian sums // Automatica.
l^k,sl JDk /Л+1 ^s+lZk,s)- 1971. у. 7, N. 4. p. 465-479.
Теперь воспользуемся условием нормировки 2. Неймарк Ю.И., Коган Н.Я., Савельев В.П.
коэффициентов в (31), которое следует из ин- Динамические модели теории управления. М.: Наука,
тегрирования правой и левой частей этого ра- 1985.
Список литературы
A CONTRIBUTION TO THE STATE ESTIMATION THEORY OF A DYNAMIC OBJECT
N.Ya. Kogan, Yu.I. Neimark, L.V. Kogan
Recurrent formulas for optimal state estimation of a linear nonstationary discrete dynamic object have been obtained on the basis of the Bayesian approach. Distribution laws of the initial object state and perturbations acting upon it have been approximated by Gaussian sums.
Keywords: dynamic object state, Kalman-Bucy filter, Gaussian sums.
