ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013
Управление, вычислительная техника и информатика
№ 4(25)
УДК 519.254.1
Д.В. Иванов, О.В. Усков
РЕКУРРЕНТНОЕ ОЦЕНИВАНИЕ БИЛИНЕЙНЫХ ЛЯХ-СИСТЕМ С ПОМЕХОЙ НАБЛЮДЕНИЯ ВО ВХОДНОМ СИГНАЛЕ
Предложен рекуррентный алгоритм для оценивания параметров билинейных АКХ с помехой наблюдения во входном сигнале. Доказана сильная состоятельность получаемых оценок параметров. Результаты моделирования подтвердили высокую эффективность предложенного алгоритма.
Ключевые слова: рекуррентное оценивание, стохастическая аппроксимация, помеха наблюдения, билинейные системы, метод наименьших квадратов.
Билинейные системы - это класс нелинейных систем с простой структурой. Билинейные системы являются простейшим обобщением линейных динамических систем: выходной сигнал зависит не только от входных и выходных сигналов, но и от произведения входного сигнала на выходной. Моделирование физических процессов с помощью билинейных систем находит применение во многих областях науки, таких, как ядерная физика, электрические сети, химическая кинетика, гидродинамика и т.д. [1].
Модели ошибки уравнения (АЯХ-модели) [2] - наиболее распространенный вид моделей параметризации шума. Идентификация моделей ошибки уравнения сводится к классической задаче регрессионного анализа и может быть решена методом наименьших квадратов. Однако во многих практических задачах помеха содержится также и во входном сигнале, в этом случае классический метод наименьших квадратов не позволяет получать состоятельные оценки.
В настоящее время активно развиваются методы идентификации билинейных динамических систем, такие, как инструментальные переменные [3], компенсирующий смещение метод наименьших квадратов [4], метод максимального правдоподобия [5] и методы на основе высших статистик [6]. Рекуррентные методы идентификации билинейных систем, которые могут быть получены из рекуррентных методов идентификации линейных систем, приведены в [7]. Некоторые предложенные методы используют подход на основе рекуррентных методов идентификации систем с ошибкой в уравнении, модифицируя при этом функцию ошибки, например улучшенный метод наименьших квадратов [8].
В статье предложен рекуррентный алгоритм оценивания параметров билинейных АКХ-систем с помехой во входном сигнале на основе стохастической аппроксимации.
Пусть билинейная динамическая система описывается стохастическими уравнениями с дискретным временем I =... -1,0,1 :
1. Постановка задачи
Г
т=0 т=0 к=1
где хг, ж г - ненаблюдаемая и наблюдаемая входные переменные; 2 ^ - наблюдаемая выходная переменная; §Д0 - помеха в уравнении; §2(0 - помеха наблюдения во входном сигнале;
Пусть выполняются следующие предположения:
10. Множество Ё, которому априорно принадлежат истинные значения параметров устойчивой, управляемой и идентифицируемой билинейной системы, является компактным.
20. Помехи {§Д0} и {§2(0} статистически не зависят между собой:
ЕЫ1)/ ^(1)} = 0, Е{§2(0/ ^(2)} = 0,
£{§2(0/^(1)} < щ(1) <«>, £{(§2(0/^(2)} = Щ(2) <”,
где /г®, - ст-алгебры, индуцированные семействами случайных величин
{§! У), t 6 Тг } и (§2 У), / 6 Тг } , Т = (^, ^ <г, ^ eZc}, Zc - множество целых чисел, Ш(), - случайные величины Е () < п§, Е (щ(2^) < п§2 , где Е - оператор
математического ожидания.
30. {§1 (-)}, {§2 (-)} статически не зависят от {хг } .
40. Последовательности {хг} - стационарные в узком смысле с дробно-рациональной плотностью случайные сигналы с £{(х1 )2} = ст2 >0. Для некоторых пх > 0 : |хг| < пх п.н.
50. Априорно известно отношение дисперсий помех у = ст2/ ст2.
2. Рекуррентный алгоритм идентификации
Уравнение (1) может быть представлено в форме линейной регрессии:
Уг =фГ 0 + ег,
(2)
где
Фг = ( (г) (г) (г) ) , Ф 2 (г) = (zг■-1,-2-г У , К (0 ^г-!---Ж-, )
К* (г) = (
2г-1 -, Жг2г-
00 = ( \ а
сь О N (с011)-.<
Бг II 1
Ж 2- ..................Ж 2- , Л
г-Г2 г-1’ > г-Г2 г-3)
)т
сТо ) , ь0 =(ь01)..-ь0Г)) , а0 =(«01)...«0Г1)У
)т
■(г21) с(г2г3(г2))
/ч ... С-/-Ч
г2 г3( т)
Бг = §1 (г) - £ а0И)§2 О' - т) - £ £ ^)2-к 0' - т).
т=0
т=0 к=1
Из предположений 10 и 20 следует, что обобщенная ошибка имеет нулевое среднее значение и ее локальная дисперсия с вероятностью 1 будет равна
1 Ы
= Цт77 £ Е ( (a0, C0, г) )2 )=ст12 + СТ2 аТа0 + СТ2 СТ2сТ с0 =
N Л/ . . ' '
= ст2(у + ат0 а0 + ст;;с’Т с0) = ст^ю(а0, с0).
Определим оценку 9(Ж) неизвестных параметров 9 из условия минимума
суммы взвешенных квадратов обобщённых ошибок ((а0,с0, I))2 с весом ю(а, О [9], т.е.
(у, -фГ9)
N
ШІП V -
. им (Ь, а, с)
■ = ІШІП- "
9єВ г=і У + а а + ст2с с Єє]В ю(а,с)
(3)
тогда оценки неизвестного вектора 9 можно получить с помощью стохастически градиентного алгоритма минимизации функции (3):
9(І +1) = 9(0 -агУ9
(+1 -фТ+19(І) )2
у+ат (г) а (г)+а12сг (г)с (і )
где аі - последовательность, удовлетворяющая условиям
ад ад
б0. Vа; = ^ аі ^аі+1 и Еаі <ад аі І > 1;
(4)
7°. Vа;^1(І) <^ Хаг^2(І) <ад п.н.
;=1 ;=1
Теорема 1. Пусть динамическая система описывается уравнениями (1) и выполняются предположения 1°-7°, тогда оценки, определяемые алгоритмом (4), ли-
бо 9(0-
->9° п.н., либо 9(;')-
-»ад.
Доказательство. Доказательство состоятельности получаемых с помощью (4) оценок основывается на методе непрерывных моделей [10, 11]. Построим асимптотическую непрерывную детерминированную модель алгоритма (4). Минимизируемую в (3) функцию можно представить в виде
где
3 (9) = ст1 +
Нф = Ііт Е
і + (9-9°) Нф(9-9°)
Т —2 Т
у + а а + о2с с
> °,
ф(°) = ( (і) і фТх (І)! ф^(І)), фг (І) = (
ф хг (І) =(xizi-1,^^^, Х; (°)
ХІ-17І-1, —, ХІ-12І-Гз(1)
'І-1>*** і-г )
ХІ-Гі 2І-1 , * * * , ХІ-Гі 2І-г3 (г2 ) 2
что следует из 10, 40.
В данном случае асимптотическая непрерывная детерминированная модель имеет вид
| = -У0 3 (9).
(5)
Здесь точка означает производную по времени. Связь между уравнениями (4) и (5)
к-1
устанавливается с помощью фиктивного времени !к = ^а, [1°, с. 89].
г =°
;=°
г =°
Пусть функция Ляпунова
V (9) = 3 (9), (6)
так как она непрерывно дифференцируема и
К (9) = V (9) 3 (9) = -||у9 3 (9)2,
У9 V(9)3(9) < 0, (7)
то множество В = {9 е Кг+ Г1+'з(0)+-+'з(г2)+ Г2+2 : V(9) = о} состоит из стационарных
точек 3(9) [10, с.114].
Для перехода от (4) к непрерывной модели необходимо показать, что для {§1(0}, {§2(0} и {а,} выполняется равенство [14, с.12]:
1
Ит Нт Бир
Т ^0 и^вд Т
а,6, (Х §2(0)
= 0, (8)
где для Т > 0 к(п,Т) = шахIк: ^а, < т!.
I , = И J
Для выполнения равенства (8) необходима ограниченность последовательности {9(0}, что подразумевает ограниченность роста функции У93 (9) при ||9| ^ вд. В нашем случае
Ит У93(9) = 0.
11е1Ьвд
Из ограниченности сумм в условии 70 и последовательности {9(0} следует
вд
ограниченность суммы ^а,6, (9(0, ^(0, §2(0) < вд, откуда следует выполнение
,=1
(8). 0 0 Из теорем, приведенных в [10, с. 12, с. 292], следует, что при выполнении 10-70 и (7), (8) последовательность {9(0} ограничена и при , ^ вд {9(0} стремится к
точкам множества В = {9 е КГ+Г1+Гз(0)+-+гз(г2)+г2+2 : V(9) = 0}.
Исследуем непрерывную модель (4): покажем, что
В* = {9 е Кг+г+Гз(0)+-+Г3(Г2)+Г2+2 : 9 = 90},
т.е. множество В* состоит из одной единственной точки 90.
Для этого рассмотрим функцию
ит Нфи 3' (и) =, и Би
где и = (и и )Т е иГ+Г1 +г3(0)+-+г3(г2)+г2 + 3
и иг+г1 + Г3(0)+-+ г3(г2)+г2 +3^
Нф = Ит Е
Б=
1 °1хг °1хг1+1 °1хГ3 (°)+- . .+Г3 (Г2 )+Г2+1
°гх1 °гхг °гхг1+1 °гхг3 (0)+_. -+Г3 (г^+Г +1
°,+м °г1+1хг 1/у °г1+1хг3 (°)+- . .+Г3 (Г2 )+Г2 +1
°Г3 (°)+- . .+Г3 (Г2 )+Г +1х1 °Г3 (°)+- . .+Г3 (Г2 )+Г2 +1хг °Г3 (°)+.. -+Г3 (Г2 )+Г2 +1х/|+1 - . +Г3 (г )+Г +11^
где 1г+1 и 1Гз(0)+ + Гз(Г2) - единичные матрицы размерностей г + 1 и
г3(0) + — + г3 (г2) соответственно.
Очевидно, что
пип3 (9) = пип3'(и ) = 3 (90 ) = лш1п , (9)
9 и
где Лш1п - минимальное собственное число регулярного пучка форм (так как Б -положительно определенная матрица), т.е. Лш1п - наименьший корень уравнения
ай( Н ф-лб) = 0.
Пусть Л -Д0) < <л(г+г1+г3(0)+-+г3(г2)+г2 +3)=Л и и и
Пусть Лш1п Л -----------Л Лшах и и1 ,...,иг+г1+г3(0)+—+г3(г2)+г2+3
какие-либо соответствующие им главные собственные векторы. Тогда Лк, где
к = 1, г + г1 + г3 (0) + — + г3 (г2) + г2 + 3 , являются стационарными значениями функции 3'(и), которые достигаются при и , равных и1,...,иг+Г1+г3(0)+ +г3(г2)+Г2+3
соот-
ветственно. Следовательно, стационарные значения функции 3(9); У93(9) = ° достигаются в точках
(м(2) и(Г+Г1 + Г3(°)+...+Гз(Г2)+Г2 +3) ^Т
u^,...,
,0)
( и(2)
Г+Г1 + Гз (°)+. - .+Г3 (г2 )+ Г2 +3
и(Г + Г1 + Гз(°) + . - .+Г3(Г1)+Г1 + 3) У Г+Г + Г3 (°)+ . -.+Г3 (Г2 )+ Г2 +3
и
(1)
г+г+Г3 (°)+. . .+Г3 (г2 )+Г2+3
г+ Г1 + г3(°)+ * - *+г3 (г2 )+ г2 + 3 (1)
ч иг-+Г1 + г3(°)+.-.+г3(г2)+ г2 +3
причем из (9) следует, что 91 =9 .
Остается показать, что
V2 3 (9) > ° (1°)
лишь в одной стационарной точке 9 = 91 =9°.
Задача определения минимума функции 3 (9) эквивалентна задаче на условный экстремум
" (11)
шіп ит Нфи, ит Би = 1.
Задача (11) может быть решена с помощью метода неопределенных множителей Лагранжа. Тогда необходимые условия запишутся в виде
(Нф - ХБ)и = °, и Би = 1,
(12)
где X - неопределенный множитель Лагранжа. Множеством решений системы
(12) являются Хе{Л1,...,Лг+(. + Гз(0)+ + Гз(Г2)+Г2+3} и соответствующие им главные
собственные вект°ры и1 ,..., иг+Г1+Г3(0)+ __ +Г3(Г2)+г2+3 .
Исследуем матрицу Нф -ХВ на положительную определенность. Из (12) следует, что
Л(1) Нф<Л(1) \Нф|,
где Л(1) |Нф| и Л(1) |Нф| - минимальные собственные числа матриц Нф и Нф соответственно.
В свою очередь, по теореме Штурма [13, с. 146]
Л(1) |Нф|<Л(2) |Нф| или Л(1) |Нф<Л(2) |Нф|. (13)
Из (13) следует, что матрица Нф -ХВ неотрицательно определена лишь при
Х = Лт1П и (10) выполняется в 91 =90, т.е для всех Х>Лтт матрица
Нф - ХВ имеет отрицательные собственные значения, откуда непосредственно следует (4).
В формуле (4) используется дисперсия выходного сигнала, которая обычно неизвестна. Согласно теореме Манна - Вольда [14]: если случайная величина с?2 сходится почти наверное соответственно к постоянной —2, то любая непрерывная функция 3 (—) сходится почти наверное к постоянной 3 (С2)
— п. >-2, з(-2) п.н. > з(с2). (14)
Следовательно, если заменить в (3) —2 оценками —2, оценки параметров 9 останутся сильно состоятельными. Состоятельная и несмещенная оценка дисперсии —2 может быть получена как
N N
—2 = ^ - !)-1 £ (2г - ^ 7.
г =1 г=1
Вычисление дисперсии (13) может быть представлено в виде рекуррентной процедуры:
Щг+1 = 7 +(Щ+1 - Ъ )/0' + 1) ,
<-2 (1 + 1) = —2(/ + 1) + ((Щ Щ )2 -—1(1 + 1)) /1.
3. Результаты моделирования
Предложенный алгоритм (4) был реализован в МаИаЪ и сравнен с рекуррент-
ным алгоритмом наименьших квадратов и рекуррентным методом расширенных инструментальных переменных. Динамическая система описывается уравнениями
Щ - 0,7 7-1 + 0,4Щ-2 = 0,3Х + 0,7 Х-1 + 0,2Х-2 + 0,2ХЩ-1 + §1(г■),
(15)
На вход подавался входной сигнал:
хг + 0,5 хг - =Сг + 0,8Сг - + 0,6Сг_2,
где (^- - белый шум.
Отношение «помеха - сигнал»: ст1 /стг и 0,2 , ст2/стг и 0,5.
Начальные значения параметров равны 0.
На рис. 1. представлены графики погрешности оценок параметров, определяемой по формуле
80г = ||® ,■ -^Ц/11(00 )-100%.
Рис. 1. График погрешности оценок параметров, %: 1 - рекуррентный метод наименьших квадратов; 2 - рекуррентный метод инструментальных переменных; 3 - алгоритм (4)
Заключение
В работе предложен рекуррентный алгоритм для оценивания параметров билинейной ARX-системы с помехой наблюдения во входном сигнале. Для получения сильно состоятельных оценок не требуется информация о законах распределения помех, достаточно знать отношение дисперсий помех. В среде Matlab создано программное обеспечение, результаты моделирования подтверждают эффективность работы предложенного алгоритма. Полученные результаты могут послужить основой для создания новых высокоэффективных автоматизированных систем управления технологическими процессами (АСУТП). Дальнейшие исследования могут быть направлены на построение алгоритмов идентификации при автокоррелированных помехах.
ЛИТЕРАТУРА
1. Mohler R.R. Bilinear Control Processes: with Applications to Engineering, Ecology, and
Medicine. New York: Academic Press, 1973.
2. ЛьюнгЛ. Идентификация систем. Теория для пользователя. М.: Наука, 1991. 432 с.
3. Ahmed M.S. Parameter estimation in bilinear systems by instrumental variable methods // Int.
J. Control. 1986. V. 44. No. 4. P. 1177-1183.
4. Ekman M. Modeling and Control of Bilinear Systems: Application to the Activated Sludge Process: PhD thesis, 2005.
5. Gabr M.M. Subba Rao T. On the identification of bilinear systems from operating records // Int. J. Control. 1984. V. 40(1). P. 121-128.
6. Tsoulkas V., Koukoulas P., Kalouptsidis N. Identification of input-output bilinear systems using cumulants // Proc. 6th IEEE Int. Conf. on Electronics, Circuits and Systems. Pafos, Greece, 1999. P. 1105-1108.
7. Fnaiech F., Ljung L. Recursive identification of bilinear systems // Int. J. Control. 1987. V. 45(2). P. 453-470.
8. Zhu Z., Leung H. Adaptive identification of bilinear systems // Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing, Phoenix, Arizona (March), 1999. P. 1289-1292.
9. Кацюба О.А. Теория идентификации стохастических динамических систем в условиях неопределенности: монография. Самара: СамГУПС, 2008. 119 с.
10. Chen H.F. Stochastic Approximation and Its Applications. Kluwer, Dordrecht, 2005.
11. Деревицкий Д.П., Фрадков А.Л. Прикладная теория дискретных адаптивных систем управления. М.: Наука, 1991. 215с.
12. Беллман Р. Введение в теорию матриц. М.: Наука, 1989. 376 с.
13. Mann H.B., Wald A. On stochastic limit and order relationship // The Annals of Mathematical Statistics. 1943. V.14(3). P. 217-226.
Иванов Дмитрий Владимирович Усков Олег Владимирович
Самарский государственный университет путей сообщения
E-mail: [email protected]; [email protected] Поступила в редакцию 3 мая 2012 г.
Ivanov Dmitriy V., Uskov Oleg V. (Samara State University of Transport). Recursive estimation of bilinear ARX systems with input-error.
Keywords: recursive estimation, stochastic approximation, error-in-variable, bilinear systems, least squares method
There is considered the problem of parameter estimation of bilinear ARX systems with noise in the input signals, described by the equations:
r r1 r2 r3(m)
zi - Z b0m) zi - m = Z a0m) xi-m + Z Z c0mk) xi - mzi-k + ^0'Х Wi = Xi + ^2(i), m=1 m=0 m=0 k=1
where xt, wt - unobserved and the observed input variables; zt - the observed output variable; ^1(i) - noise in the equation; ^2(i) - the noise in the input signal.
We propose a recursive algorithm for estimation of parameters, which is a generalization of the method of least squares. It is proved that that under non-restrictive conditions on the signals and noises, the proposed algorithm gives a strongly consistent estimators. The simulation results confirmed the high efficiency of the algorithm.