CC BY
20
УДК 618.5.015
РОБАСТНАЯ ПРОЦЕДУРА АКТИВНОЙ ПАРАМЕТРИЧЕСКОЙ
ИДЕНТИФИКАЦИИ СТОХАСТИЧЕСКИХ ЛИНЕЙНЫХ ДИСКРЕТНЫХ СИСТЕМ
В.М. Чубич, А.Э. Прокофьева
На основе фильтра Изанлу - Фейкуриана - Джазди - Саймона предложена робастная процедура активной параметрической идентификации стохастических линейных дискретных систем, включающая в себя робастное оценивание параметров и планирование входных сигналов. Разработано соответствующее программно-математическое обеспечение. Эффективность применения процедуры продемонстрирована на примере модели системы управления положением двигателя постоянного тока при группированном характере расположения аномальных наблюдений.
Ключевые слова: активная параметрическая идентификация, аномальные наблюдения, критерий квазиправдоподобия, информационная матрица Фишера, планирование эксперимента, робастное оценивание, фильтр Изанлу - Фейкуриана - Джазди -Саймона, стохастическая линейная дискретная система.
Проблема идентификации, связанная с построением математических моделей динамических систем по экспериментальным данным, относится к одной из основных проблем теории и практики автоматического управления. Ее качественное решение способствует эффективному применению на практике современных математических методов и наукоемких технологий, например, при расчете и проектировании систем управления подвижными (в том числе авиационно-космическими) и технологическими объектами, построении прогнозирующих моделей (например, в экономике и бизнес-процессах), конструировании следящих и измерительных систем.
Первоначально методология построения динамических моделей развивалась в рамках пассивного подхода, при котором идентификация проводится в режиме нормальной эксплуатации исследуемой системы. Современная теория включает в себя также методы активной идентификации, предполагающие подачу на вход исследуемой системы определенным образом синтезированных управляющих сигналов. При этом трудности, связанные с необходимостью нарушения технологического режима, окупаются за счет повышения эффективности и корректности проводимых исследований, что вытекает из самой идеологии активной идентификации, базирующейся на сочетании приемов параметрического оценивания и концепции планирования эксперимента.
При предварительно выбранной модельной структуре процедура активной идентификации систем предполагает выполнение следующих этапов:
1) вычисление оценок неизвестных параметров по экспериментальным данным, соответствующим пробному входному сигналу;
323
2) синтез оптимального входного сигнала на основе полученных оценок и по некоторому выбранному критерию;
3) пересчет оценок неизвестных параметров по экспериментальным данным, соответствующим синтезированному сигналу.
Применительно к моделям в пространстве состояний можно отметить, что среди публикаций по активной параметрической идентификации стохастических систем заметно превалируют работы, в которых используются детерминированные уравнения состояний (в качестве примера приведем [1 - 6]). Исследований, в которых используются стохастические модели состояний, выполнено заметно меньше (выделим [7 - 11]). По-видимому, данное обстоятельство связано с заметной сложностью вычисления информационной матрицы Фишера, лежащей в основе процедур численного построения оптимальных планов эксперимента. Косвенным подтверждением сказанному является то, что некоторые исследователи вместо непосредственного вычисления информационной матрицы предпочитают работать с ее приближенным аналогом [12 - 15]. В настоящий момент сложилась определенная диспропорция между относительно недавно полученными теоретическими результатами в области планирования эксперимента (имеются в виду выражения для информационных матриц Фишера применительно к различным модельным структурам из [16 - 19]) и их прикладным применением.
При решении практических задач довольно часто возникают ситуации, связанные с наличием аномальных наблюдений в измерительных данных, что может быть обусловлено отклонением распределения шума измерений от гауссовского. Использование алгоритмов, не позволяющих учесть появление таких наблюдений, может привести к смещению оценок и/или некорректному решению задачи оценивания. В связи с этим становится актуальным разработка робастных процедур активной параметрической идентификации стохастических динамических систем, в которых применяется тот или иной метод робастного оценивания неизвестных параметров. В настоящей статье авторы предпринимают попытку разработать подобную процедуру применительно к моделям стохастических линейных дискретных систем, взяв за основу коррентропийный фильтр Изанлу - Фейкуриана - Джазди - Саймона (один из наиболее эффективных робастных фильтров в настоящее время [20]), материалы по планированию эксперимента из [10, 11] и по робастному оцениванию из [21].
Рассмотрим следующую управляемую, наблюдаемую и идентифицируемую модель стохастической линейной дискретной системы в пространстве состояний:
х(гк+1) = Ф (гк) х(гк) + V (гк )и (гк) + Т(гк) ъ(гк), (1)
У(*к+1) = Щ*к+1)х«к+1) + ^к+д, к = 0,1,...,К-1, (2)
324
где x(tk) - п -вектор состояния; u(tk) - детерминированный г -вектор управления (входа); w(tk) - p-вектор шума системы; у (tk+1) - ш-вектор измерения (выхода); v(tk+1) - ш-вектор шума измерений.
Предположим, что:
- случайные векторы w(tk), v(tk+1) образуют белые гауссовские последовательности, причем
Е [ w(tk)] = 0, Е [ ) М,Т (Ц) ] = )8Ы,
Е ^+1)] = 0, Е\v(tk+^Т (tг+1)] = Я^,
Е
Т
У )
0, k,I = 0,...Л-1
(здесь и далее Е[•] - оператор математического ожидания, 8ц - символ Кронекера);
- начальное состояние химеет нормальное распределение с параметрами
Е [ х(*0)] = X (¿0),
Е
(x(to) - х ад)(- х )
Т
=рад
и не коррелирует с w(tk), v(tk+1) при любых значениях переменной £ ;
- в выходных данных могут присутствовать аномальные наблюдения;
- неизвестные параметры сведены в вектор © = ($1,$2,...,#5 )е и могут входить в элементы матриц ), ), ), H(tk+l), Q(tk), Я(tk+1), P(to) и в вектор Xад в различных комбинациях.
Необходимо для математической модели (1), (2) с учетом высказанных априорных предположений разработать на базе робастного фильтра Изанлу - Фейкуриана - Джазди - Саймона процедуру активной параметрической идентификации стохастических динамических системы и провести численное исследование эффективности ее применения.
Оценивание неизвестных параметров математической модели, описываемой уравнениями (1), (2), будем осуществлять по данным наблюдений х в соответствии с некоторым критерием идентификации с Сбор данных происходит в процессе проведения идентификационных экспериментов, выполняющихся по некоторому плану .
Предположим, что можно произвести п запусков системы, причем сигнал подается на вход системы dl раз, сигнал и2 - d2 раз и т. д., наконец, сигнал и(1 - dq раз. В этом случае дискретный (точный) нормированный план эксперимента имеет вид
U1 U2 ■ • Uq
= d1 di . . d±
. V V V ,
, Uj е Wu, i = 1,2,...,q,
Nr
где Wu c ^ - область допустимых входных сигналов.
Обозначим через Y? =
ylJ (ti)
T
ylJ (tN )
T
j-ю реализацию
выходного сигнала (J = 1,2,..., dj), соответствующую l-му входному сигналу
UT =
ul (to)
T
ul (tN-i)
T
l = 1,2,...,q. Тогда в результате проведения
экспериментов по плану XV будет сформировано множество
5 = {(U
\ т q
Ц,Ylj), J = 1,2,...,dt, l = 1,2,...,q}, X dT
l=1
V
В качестве метода оценивания выберем метод максимального правдоподобия, обладающий такими важными для практики асимптотическими свойствами как несмещенность, состоятельность и эффективность [22, 23]. В силу того, что измерительные данные содержат аномальные наблюдения (распределение шума измерений может отличаться от нормального), будем вычислять квазиправдоподобные оценки [24], т.е. найдем такие значения 0 из области допустимых значений Wq , для которых
0 = arg min [^(0; 5)1 = arg min Г - ln L (0; 5)1. (3)
0eW@ QGWq
Здесь
/^„ч Ищу „ 1 ^¿И-1 J п]. с(0;х) = — 1п2р + — ЕЕ X В (Гк+1) +
2 21=1]=1 к=0 1 ^ ^ N—1Г .. -Гг .. -1—1
+2ЕЕ Е е(1к+1) В&+1) е]&+1), 2/=1 ]=1 к=0 -1 -1
где е ^к+1) и Би +1) определим по следующим рекуррентным уравнениям фильтра Изанлу - Фейкуриана - Джазди - Саймона из [25]:
(^+11*к) = ф(Ь)^^ I tk) + ¥(1к У ^к); Р] (<к+11 к) = ф(Ъ )р] (к I *к )Фг &)+т №к )гг (Гк); ее (Гк+1) = / (Гк+1)— Н(Гк+1)Х (Гк+11 Гк);
Ll (tk+1) = exp
-\T
£lJ (tk+1) " R~\tk+1)£lJ (tk+1)
2s
B1J (tk+1) = H(tk+1 (tk+11 tk)Lj (tk+1 )HT (tk+1) + R(tk+1);
T г .. -|-1
(tk+1) = PZ7 (tk+11 tk)L (tk+1)HT (tk+1) BlJ (tk+1)
JC
ij
(tk+11 tk+i) = xj (tk+11 tk) + Kl (tk+1)e' (tk+i); R(tk+11 tk+1)=[1 - Kj (tk+i)H(tk+i)] R(tk+11 tk).
Начальные условия для приведенных уравнений выбираются с учетом высказанных при постановке задачи априорных предположений
(xlJ (t0 110) = x (t0), Pj (t0 110) = P(to)), а параметр s - в соответствии с материалами [20].
При решении задачи нелинейного программирования (3) воспользуемся методом последовательного квадратичного программирования [26, 27]. Данный метод относится к методам первого порядка, реализован в рамках пакета Optimization Toolbox программной системы MATLAB, и в настоящее время является одним из самых передовых и эффективных.
Планирование входных сигналов. Применение теории оптимального управления при оценивании неизвестных параметров позволяет повысить точность оценивания параметров за счет сбора наиболее информативных данных.
Под непрерывным нормированным планом X будем понимать совокупность величин
U1 Р1
U 2 Р2
U
q
q
1,2,..., q. (4)
Рг ^ 0,1 Рг = 1, иг еО^, г РЧ \ г=1
В отличие от дискретного плана в непрерывном плане X веса Рг могут принимать любые значения, в том числе и иррациональные.
Информационная матрица М (X) плана (4) определяется соотношением
q
m (Х) = I PiM (и, (),
i=1
в котором информационные матрицы Фишера одноточечных планов M (Ui, 0) зависят от вектора оценок неизвестных параметров 0 (что позволяет говорить о локально-оптимальном планировании) и вычисляются в соответствии с [18].
Оптимальный план эксперимента находится для некоторого выпуклого функционала X от информационной матрицы M (X) путем решения следующей экстремальной задачи:
X* = arg min X[m(X)]. (5)
Воспользуемся критерием D-оптимальности, для которого X [M (X) ] = - lndet M (X) . Применение этого критерия способствует минимизации объема эллипсоида рассеяния оценок неизвестных параметров.
Возможны два подхода к решению экстремальной задачи (5): прямой и двойственный [1, 28-30]. Первый из них предполагает поиск минимума функционала X [M (X)] непосредственно (при этом рекомендуется s(s +1)
выбирать q
2
+1), второй - решение двойственной задачи. При
решении практических задач целесообразно применение так называемого комбинированного подхода, при котором для улучшения начального приближения плана сначала используют прямую процедуру, а затем - для получения окончательного результата - двойственную.
Представим, основываясь на методе последовательного квадратичного программирования, комбинированную процедуру построения непрерывных оптимальных планов при предварительно вычисленных оценках
параметров 0.
Шаг 1. Зададим начальный невырожденный план
x
U0
Pl
U2 Po
Uq Pq0
, U0 EÜU, P2
—, l
1,..., q.
Вычислим информационные матрицы одноточечных планов m(U°,0)для l = 1,...,q и положим k = 0.
Шаг 2. Считая веса плана р1,..., фиксированными, найдем план
Xk+1 = arg , min Х[M(Xk)].
min
Ulk ,..,Uk îWu
M (U
Вычислим k+1
информационный матрицы одноточечных планов 0), I = 1,...,Я.
Шаг 3. Зафиксировав точки спектра плана Хк+1, найдем план
x+1 = arg
min
P\,P2 ,...,Pkq
X
M (X+1)
Pk ^0, qpk = 1, l = 1,...,q. i=1
Шаг 4. Если для малого положительного d выполняется неравен-
ство
q
I
i=1
uk+1 - uk
+ ( pk+1 - pk )2
<
то положим Хо = Хк+1 (выполнение прямой процедуры закончено), к = 0 и перейдем на шаг 5. В противном случае примем к = к +1 и перейдем на шаг 2.
2
Шаг 5. Вычислим информационную матрицу М (Хк). Шаг 6. Найдем локальный максимум
ик = а^ тах ф иеаи
м-1(Хк) М (и)'
где под Sp[•] понимается след соответствующей матрицы. Если выполня-
ется условие
Яр
м~
1 (4) м [ик)'
<¿2, то закончим процесс. Если
Яр
мг
1 (X) м (ик)'
> £, перейдем на шаг 7, в противном случае будем
искать новый локальный максимум. Шаг 7. Найдем ц
Ч = а^ тт X М (Х^+1) 0<ч<1 L
Здесь Хк+1 =(1 ~ч)Хк +Ч%(ик), где Х(ик) - одноточечный план, разме-
щенный в точке и .
Шаг 8. Составим план
Хк+1 =(1 -Чк )Хк +ЧкХ(ик)
произведем его «очистку», следуя [28], положим к = к +1 и перейдем на шаг 5.
Приведенный алгоритм требует вычисления производных от информационной матрицы одноточечных планов по компонентам входного сигнала (аналитическое выражение производной и алгоритм ее вычисления представлен в [11]).
Практическое применение синтезированного оптимального плана затруднительно, поскольку веса представляют собой произвольные вещественные числа, заключенные в интервале от нуля до единицы. В случае заданного числа п возможных запусков системы необходимо провести «округление» непрерывного плана до дискретного [29], провести идентификационные эксперименты и пересчитать оценки неизвестных параметров.
Реализуем описанную процедуру активной параметрической идентификации в рамках программной среды МЛТЬЛВ и апробируем ее на модели системы управления положением, состоящей из антенны и двигателя постоянного тока [31].
£
Пусть первая компонента вектора состояния отвечает за угловое положение антенны, вторая - за ее угловую скорость; входным сигналом является напряжение на входе усилителя постоянного тока, управляющего двигателем; с помощью потенциометра измеряется угловое положение и выполняются все априорные предположения, высказанные при постановке задачи. Тогда модели состояния и наблюдения в дискретные моменты времени определяются соотношениями
х( /к+1)
1 —(1 - в4 ) 4
0 в~4Т
х(/к ) +
—(Т -1 +1 в~4Т ) 4 4 4
— (1 - в"4!? ) 4
и<4к ) +
у(/к+1) = [1 0]х(*к+1) + ^+1), к = 0,1,...,N -1,
где 4,4 - неизвестные параметры и ^0 ={1 <4 < 6, 0 <42 < 1] .
Положим Т = 0.1, N = 30, ((/к) = ( = 0.01, Я(/к+0 = Я = 0.1,
х(/0) =
"0" "0.01 0 "
0 , Р(/0) = 0 0.01
и зададим значение параметра фильтра
Изанлу - Фейкуриана - Джазди - Саймона а = 1.
В качестве области допустимых входных сигналов выберем
% ={ие □ N12 < и(/к) < 20,к = 0,1,...,N -1].
Для ослабления зависимости результатов оценивания от выборочных данных проведем пять независимых запусков системы, а полученные оценки неизвестных параметров усредним. Реализации выходных сигналов
получим компьютерным моделированием при истинных значениях пара* *
метров 4 = 4.600,42 = 0.787 [31], расположив подряд три аномальных
наблюдения с дисперсией шума Яд =100Я в середине выборки.
О качестве параметрической идентификации будем судить по значению относительной ошибки оценивания 8д, вычисляемой по формуле
4
1
>к Л О * Л О
(4 -4)2 + (4* -4*)2
* о * о
(4* )2 + (4*)2
где 00 = ) - вектор оценок неизвестных параметров.
Результаты выполнения робастной процедуры активной параметрической идентификации представлены в таблице (оптимальный план получился одноточечным).
Результаты выполнения робастной процедуры активной параметрической идентификации
Исходный и синтезированный входные сигналы вместе со значениями критерия оптимальности
Номер
запуска
системы
Оценки и ошибка оценивания
в
1
А
в
1
4,2188
0,6707
2
4,2263
0,7204
4,0909
0,6724
4
4,9774
0,9003
5
4,2188
0,6707
Х[ М (Х)] =-7.6551
Среднее значение
по запускам
4,3464
0,7269
0,084
0,081
0,112
0,084
0,085
0,089
1
4,7761
0,8134
2
4,3398
0,7388
4,7784
0,8330
4
4,8019
0,8661
5
4,4058
0,7642
Х[ М (X)] =-10.1402
Среднее значение
по запускам
4,6204
0,8031
0,038
0,057
0,039
0,046
0,042
0,045
3
3
Анализ таблицы показывает, что планирование входных сигналов с использованием комбинированной процедуры позволяет повысить качество оценивания на 4,4 % процента.
Заключение. Для моделей стохастических линейных дискретных систем на основе фильтра Изанлу - Фейкуриана - Джазди - Саймона разработана робастная процедура активной параметрической идентификации, включающая в себя робастное оценивание параметров и планирование входных сигналов. Рассмотрен случай вхождения неизвестных параметров в уравнения состояния, наблюдения, начальное условие и ковариационные матрицы шума системы и измерений. Эффективность разработанной роба-стной процедуры активной параметрической идентификации продемонстрирована на примере модели системы управления положением двигателя постоянного тока при группированном характере расположения аномальных наблюдений.
Работа выполнена при финансовой поддержке Министерства образования и науки РФ по государственному заданию (проект № 2.7996.2017/8.9).
Список литературы
1. Денисов В.И., Попов А.А. Пакет программ оптимального планирования эксперимента. М.: Финансы и статистика, 1986. 159 с.
2. Morelli E.A. Flight test of optimal inputs and comparison with Conventional Inputs // Journal of aircraft. 1999. Vol. 36. №2. P. 389 - 397.
3. Овчаренко В.Н. Планирование идентифицирующих входных сигналов в линейных динамических системах //Автоматика и телемеханика. 2001. № 2. С. 75 - 87.
4. Jauberthie С, Bournonville F., Coton P., Rendell F. Optimal input design for aircraft parameter // Aerospace science and technology. 2006. №10. P.331 - 337.
5. Jakowluk W. Design of an optimal input signal for plant-friendly identification of inertial systems // Przeglad Elektrotechniczny. 2009. R. 85 (10). P. 123 - 128.
6. Galvanin F., Marchesini R., Barolo M., Bezzo F., Fidaleo Optimal design of experiments for parameter identification in electrodialysis models // Chemical engineering research and design. 2016. № 105. P. 107 - 119.
7. Mehra R.K. Synthesis of optimal inputs for multiinput - multioutput (MIMO) systems with process noise. Part II. Time domain synthesis // System identification - advances and case studies. New York: Academic press, 1976. P. 230 - 249.
8. Goodwin G.C., Payne R.L Dynamic system identification: experiment design and data analysis. New York: Academic Press, 1977. 302 P.
9. Активная параметрическая идентификация стохастических линейных систем: монография / В.И. Денисов, В.М. Чубич, О.С. Черникова, Д.И. Бобылева. Новосибирск: Изд-во НГТУ, 2009. 192 с.
10. Чубич В.М. Оптимальная идентификация дискретных систем на основе метода статистической линеаризации // Информационные технологии и вычислительные системы. 2010. № 4. С. 47 - 56.
11. Чубич В.М. Информационная технология активной параметрической идентификации стохастических квазилинейных дискретных систем // Информатика и ее применения. 2011. Т. 5. №1. С.46 - 57.
12. Gupta N.K., Mehra R.K. Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations // IEEE Transactions on automatic control. 1974. Vol. 19. № 6. P. 774 - 783.
13. Zadrozny P.A., Mittnik S. Kalman-Filtering methods for computing information matrices for Time-Invariant, periodic, and generally time-varying VARMA models and samples // Computer Math. Applic. 1994. Vol. 28. P. 107 - 119.
14. Еланцева И.Л., Денисов В.И., Полетаева И.А. Активная идентификация авторегрессионных моделей со скользящим средним // Научный вестник НГТУ. 2006. №3(24). С. 27 - 34.
15. Kok M., Dahlin J., Schon T.B., Wills A. Newton-based maximum likelihood estimation in nonlinear state models // IFAC-PapersOnLine. 2015. Vol. 48. P. 398 - 403.
16. Klein A., Melard G., Zahaf T. Construction of the exact Fisher information matrix of Gaussian time series models by means of matrix differential rules // Linear Algebra and its Applications. 2000. Vol. 321. P. 209 - 232.
17. Klein A., Melard G., An algorithm for the exact Fisher information matrix of vector ARMAX time series // Linear Algebra and its Applications. 2014. Vol. 446. P. 1 - 24.
18. Чубич В.М. Вычисление информационной матрицы Фишера в задаче активной параметрической идентификации стохастических нелинейных дискретных систем // Научный вестник НГТУ. 2009. №1. С. 23 - 40.
19. Чубич В.М. Особенности вычисления информационной матрицы Фишера в задаче активной параметрической идентификации стохастических нелинейных непрерывно-дискретных систем // Научный вестник НГТУ. 2009. №1. С.41 - 54.
20. Чубич В.М., Прокофьева А.Э. Сравнительный анализ некоторых робастных фильтров для нестационарных линейных дискретных систем // Вестник Иркутского государственного технического университета. 2017. Т. 21. № 12. С. 123 - 137.
21. Чубич В.М., Прокофьева А.Э. Параметрическая идентификация стохастических линейных дискретных систем на основе робастной фильтрации // Вестник Иркутского государственного технического университета. 2018. Т. 22. № 2. С. 84 - 94.
22. Боровков А. А. Математическая статистика. Новосибирск: Наука, 1997. 772 с.
23. Ивченко Г.И., Медведев Ю.И. Введение в математическую статистику. М.: Изд-во ЛКИ, 2010. 600 с.
24. Мудров В.И., Кушко В. Л. Методы обработки измерений. Квазиправдоподобные оценки. М.: Радио и связь, 1983. 304 с.
25. Kulikova M.V. Square-root algorithms for maximum correntropy estimation of linear discrete-time systems in presence of non-Gaussian noise // Systems & Control Letters. 2017. Vol. 108. P. 8 - 15.
26. Antoniou A., Lu W.-S. Practical optimization: algorithms and engineering applications. New York: Springer, 2007. 669 p.
27. Измаилов А.Ф., Солодов М.В. Численные методы оптимизации. М.: ФИЗМАТЛИТ, 2008. 320 с.
28. Федоров В.В. Теория оптимального эксперимента. М.: Наука, 1971. 312 с.
29. Ермаков С.М., Жиглявский А.А. Математическая теория оптимального эксперимента. М.: Наука, 1987. 320 с.
30. Pronzato L., Pazman A. Design of experiments in nonlinear models. Asymptotic normality, optimality criteria and small-sample properties. New York: Springer, 2013. 399 p.
31. Квакернаак Х., Сиван Р. Линейные оптимальные системы управления / пер. с англ. М.: Мир, 1977. 650 с.
Чубич Владимир Михайлович, д-р техн. наук, доцент, зав. кафедрой, chuhich a ami.nstu.ru, Россия, Новосибирск, Новосибирский государственный технический университет,
Прокофьева Алина Эдуардовна, магистрант, pae. I994'a yandex. ru, Россия, Новосибирск, Новосибирский государственный технический университет
ROBUST PROCEDURE OF ACTIVE PARAMETER IDENTIFICATION OF STOCHASTIC LINEAR DISCRETE SYSTEMS
V.M. Chuhich, A.E. Prokofeva
On the basis of Izanloo - Fakoorian - Yazdi - Simon filter we proposed a robust procedure for the active parametric identification of stochastic linear discrete systems, including robust estimation of parameters and input signal design. We have developed the relevant software. The effectiveness of the procedure is demonstrated hy the example of a DC motor control system for the grouped location of outliers.
Key words: active parametric identification, outliers, quasi- likelihood criterion, Fisher information matrix, design of experiments, robust estimation, Izanloo - Fakoorian -Yazdi - Simon filter, stochastic linear discrete system.
Chubich Vladimir Mikhailovich, doctor of technical sciences, docent, head of chair, chubichaami. nstu. ru, Russia, Novosibirsk, Novosibirsk State Technical University,
Prokofeva Alina Eduardovna, master's, pae. [email protected], Russia, Novosibirsk, Novosibirsk State Technical University
