Активная параметрическая идентификация нелинейных дискретных систем на основе линеаризации во временной области и оптимального управления

УДК 681.5.015

АКТИВНАЯ ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ НЕЛИНЕЙНЫХ ДИСКРЕТНЫХ СИСТЕМ НА ОСНОВЕ ЛИНЕАРИЗАЦИИ ВО ВРЕМЕННОЙ ОБЛАСТИ

Н НПТНМАЛЬНОГО УПРАВЛЕННЯ1

В.М. Чубич, О.С. Черникова

Получены результаты для случая, когда неизвестные параметры содержатся в уравнениях состояния и наблюдения, начальных условиях и ковариационных матрицах помех динамики и ошибок измерений. Рассмотрен пример оптимального оценивания параметров одной модельной структуры.

Ключевые слова: линеаризация, оценивание параметров, метод максимального правдоподобия, планирование оптимальных входных сигналов, информационная матрица Фишера, метод Шатровского.

ВВЕДЕНИЕ

Проблема идентификации относится к одной из основных проблем теории и практики автоматического управления. Ее качественное решение способствует эффективному применению современных математических методов и сложных наукоемких технологий при проектировании различных систем управления подвижными и технологическими объектами, построении прогнозирующих моделей, конструировании следящих и измерительных систем и др.

По способу проведения эксперимента идентификацию можно разделить на пассивную и активную. При пассивной идентификации для построения математической модели используются реально действующие в системе сигналы и тем самым нормальный режим эксплуатации не нарушается. Методы пассивной идентификации достаточно полно описаны, например, в книгах [1, 2]. Активная идентификация, напротив, предполагает нарушение технологического режима и подачу на вход изучаемой системы специальным образом синтезированного сигнала. Его находят в результате решения экстремальной задачи для некоторого предварительно выбранного функционала от информа-

1 Работа выполнена при поддержке Федерального агентства по образованию в рамках ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009 — 2013 гг.

ционной (или дисперсионной) матрицы вектора оцениваемых параметров. Трудности, связанные с необходимостью нарушения технологического режима, должны окупаться повышением эффективности и корректности проводимых исследований, что обусловлено самой идеологией активной идентификации, базирующейся на сочетании приемов параметрического оценивания с теорией планирования эксперимента [3, 4].

Процедура активной параметрической идентификации (оптимального оценивания параметров) моделей динамических систем предполагает:

— вычисление оценок параметров по измеренным данным, соответствующим некоторому пробному сигналу;

— синтез на основе полученных оценок оптимального по некоторому выбранному критерию сигнала;

— пересчет оценок неизвестных параметров по измеренным данным, соответствующим синтезированному сигналу.

Целесообразность применения активной параметрической идентификации для построения математических моделей стохастических линейных систем показана в работах [5—7]. Тем не менее, данная область исследований остается еще недостаточно изученной, возможности применения в ней методов оптимального оценивания параметров выявлены далеко не полностью. В настоящей статье приведены результаты дальнейших исследований авторов в рамках указанной проблемы

применительно к гауссовским нелинейным дискретным системам, причем входной сигнал синтезируется путем решения соответствующей задачи дискретного оптимального управления методом Шатровского.

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим следующую модель управляемой, наблюдаемой, идентифицируемой динамической системы в пространстве состояний:

х

у(к + 1) = А[х(к + 1), к + 1] + у(к + 1),

к = 0, 1, ..., N - 1,

где f [х(к), и(к), к] — п-мерная векторная функция указанных аргументов; х(к) — п-мерный вектор состояния; и(к) — детерминированный г-мерный вектор управления (входа), Г(к) — матрица размера п х р; ^(к) — р-мерный вектор возмущения; к — дискретное время; у(к + 1) — т-мерный вектор измерения (выхода); А[х(к + 1), к + 1] — т-мерная векторная функция указанных аргументов; у(к + 1) — т-мерный вектор ошибки измерения.

Предположим, что f [х(к), м(к), к] и А[х(к + 1), к + 1] непрерывны и дифференцируемы по х(к), и(к) и х(к + 1) соответственно; случайные векторы ^(к) и у(к + 1) образуют взаимно независимые стационарные белые гауссовские последовательности с нулевыми средними и ковариационными матрицами О и Я соответственно, т. е.

Е [ЦА)] = 0,

Е [ЦА)^/)] = 06^

Е [у(к + 1)] = 0, Е [у(к + 1)уг(/ + 1)] = Я5Й,

Е[у(к + 1)^т(/)] = 0, к, I = 0, 1, ..., N - 1,

(здесь и далее Е [•] — оператор математического ожидания, 8Й — символ Кронекера); начальное состояние х(0) имеет нормальное распределение

со средним х (0) и ковариационной матрицей Р(0), т. е.

Е [х(0)] = X (0),

Е{[х(0) - X (0)][х(0) - X (0)]т} = Р(0),

и не коррелирует с м>(к) и у(к + 1) при любых значениях переменной к; неизвестные параметры сведены в вектор © = (9Х, 02, ..., 9^), включающий в себя элементы вектор-функций ^х(к), и(к), к], А[х(к + 1), к + 1], матриц г(к), О, Я, Р(0) и вектора

X (0) в различных комбинациях.

Необходимо для математической модели (1), (2) с учетом высказанных априорных предположений разработать процедуру активной параметри-

ческой идентификации на основе линеаризации и оптимального управления, исследовать ее эффективность и целесообразность применения.

2. ЛИНЕАРИЗАЦИЯ МОДЕЛИ

Считая значение вектора неизвестных параметров © фиксированным, выполним линеаризацию во временной области нелинейной модели (1), (2) относительно номинальной траектории {хн(к + 1), к = 0, 1, ..., N - 1}, для которой

Хн (А + 1) = /[ Хн( А), Ин( А), А], А = 0, 1, N - 1,

Хн(0) = Х(0).

Разложив для каждого к вектор-функции f[х(к), и(к), к] и А[х(к + 1), к + 1] в ряды Тейлора в окрестностях точек [хн(к), ин(к)] и хн(к + 1) соответственно и отбросив члены второго и более высоких порядков, запишем уравнения линеаризованной модели

х(А + 1) = / [хн(А), ин(А), А] +

+ - і-- н - - - - - 7 н V - - , >  - - [х(А) - хн(А)] +

3/ [Хн

+ 3/[ Хн

3и

у(А + 1) = А[хн(А + 1), А + 1] +

[х(А + 1) - хн(А + 1)] +

+ ЗА[хн(А+ 1) , А+ 1-

Зх(А+ 1) н

+ у(А + 1),

для которой и будем решать поставленную задачу. С учетом обозначений

где

я(А) = а(А) + ^(А)и(А);

а(А) = / [хн(А), ин(А), А]

3/[Хн(А), ин(А), А] _ 3/[Хн(А), ин(А), А]

Зх

Хн

3и

ин(А),

Ш(А) = 3/[Хн(ин(А]

( ) 3и(А) ,

= 3/[ Хн ( А) , ин ( А) , А] .

( ) Зх(А) ;

А(А + 1) = А[хн(А + 1), А + 1] ЗА[хн(А + 1), А + 1]

Зх

Хн(А + 1);

Н

Зх

соотношения (4) и (5) определяют модель гауссовской линейной нестационарной системы, описываемой уравнениями

х

у

Заметим, что изложенный способ линеаризации не применим к неоднозначным функциям и нелинейностям, имеющим угловые точки и разрывы. Для линеаризации таких нелинейностей можно воспользоваться методом статистической линеаризации (см., например, работы [8, 9]).

3. ОЦЕНИВАНИЕ НЕИЗВЕСТНЫХ ПАРАМЕТРОВ

Априорные предположения, высказанные в § 1, и выполненная в § 2 линеаризация во временной области моделей состояния и наблюдения относительно выбранной детерминированной номинальной траектории (3) позволяют воспользоваться для оценивания параметров методом максимального правдоподобия. В соответствии с этим методом

необходимо найти такие значения параметров О, для которых

О = аг§ тт [-1^(0; 5)] = аг§ тт [х(©; 2)],

0ЕЙ0 0ЕЙ0

где в соответствии с [2, 10] можно записать:

N - 1

22

х(©; 2) = ^ 1п2п + 1 £ [е(А + 1)]ТВ Х(А + 1) X

х

1

к = 0

N - 1

[е(А + 1)] + і £ 1ndetB(k + 1).

к = 0

Здесь О0 — область допустимых значений параметров; L(©; Е) — функция правдоподобия; Е — выборочные данные; %(©; Е) — взятая со знаком «—» логарифмическая функция правдоподобия;

є (А + 1) = у(А + 1)-у (А + 1|А); у (А + 1|А) и B(k + 1) определяются по соответствующим рекуррентным уравнениям дискретного фильтра Калмана [11]:

Х (А + 1|А) = Ф(А) Х (А|А) + я(А);

Р(А + 1|А) = Ф(А)Р(А |А)ФТ(А) + ГДООГ^А); у (А + 1|А) = Н(А + 1)Х(А + 1|А) + А(А + 1);

В(А + 1) = Н(А + 1)Р(А + 1|А) + НТ(А + 1) + Р; К(А + 1) = Р(А + 1^)#^ + 1) + В-1 (А + 1);

Х (А + 1|А + 1) = Х (А + 1|А) + К(А + 1)б(А + 1);

Р

для к = 0, 1, ..., N - 1 с начальными условиями

х (0|0) = х (0), Р(0|0) = Р(0).

4. СИНТЕЗ ВХОДНОГО СИГНАЛА

Точность оценивания неизвестных параметров можно повысить, оптимизируя определенным образом выбранный критерий от информационной матрицы Фишера

М(и; ©) = [^(«(0), и(1), ..., - 1); ©),

/, у = 1, 2, ..., 5,

которая связана с логарифмической функцией правдоподобия соотношением [2]

М(и ©) = - Е

'321п L ( © ; Е )' 39/30/

Здесь усреднение берется по выборочному пространству. Зависимость ИМФ от неизвестных параметров обусловливает локально-оптимальное планирование.

В теории оптимального эксперимента наиболее часто используются критерии ^-, А-, Е-оптималь-ности с вполне определенной статистической интерпретацией [3, 4]. Следуя работе [12], будем использовать в качестве критерия оптимальности след ИМФ:

и * = а^шах $>М(и; ©).

и ЕЙ,

Для гауссовской линейной нестационарной модели (10), (11) в работе [13] получена формула для ИМФ, допускающая ее разложение на сумму двух слагаемых, одно из которых зависит от входного сигнала, а другое — нет:

М (и; ©) = ИУ (и; ©) + Уу(©).

Это позволяет задачу (12) представить в виде и* = ащшах ДрИЩ ©).

ие

Воспользовавшись материалами работы [13], запишем

/1 = ^Ж(и; 9) = £ И^и; 9) =

і = 1

N -1

-Т,

£ [ ХА (А + 1)2А1(А + 1) ХА (А + 1) +

к = 0

+ 2ХА

где

где

ха (А + 1) = Е

х

301 Зх

39,

^(к + 1) =

^Аі'сА + 1) ^Аі2

А + 1) ^2’12(А + 1) ...

0

^ + 1 1( А + 1) 0

причем

= £ 3ЯГ( А + 1 ) В-1 (А + 1) 3Н( А + 1 ) ;

А1 £ 39. 4 7 39;

і= 1 і і

Т

гЦ = 3Н ( А + 1 ) В-1(А + 1)Н(к + 1), у = 2,7+1

39/ -1

= НТ(к + 1)В-1(А + 1)3Н ( А + 1 ), / = 2, 5 + 1;

39; -1

= НТ(А + 1)В“\А + 1)Н(к + 1), / = 2, 5 + 1;

^(к + 1) =

0 дНТ{к + 1 )в-1

-1

00

0 НТ(к + 1)в-1(к + 1) ...

0

0 0 ... НТ

А

А(А + 1) =

В силу симметричности матрицы 2А 1(А + 1) и того, что

Ха (А + 1) = Фа( А) Ха (А) + ^а( А) и (А) + Ьа (А),

А = 0, 1...N - 1;

ха

ХТ(0) ЗхЮ) ... 3ХІ0

39,

39"

Ф

8Ф

30! 1 ’ 30!

дФ

30, 1 ; 30,

Фа (А) =

0 ...

Ф(к) -К(к)И(к) ...

0 ...

0

0

Ф

^А (А) =

У( А) 3у ( А) 3 9 1

3у

3 9„

Ьа (А) =

а

3а ( А- - ^^( А)3А( А) 391 ^ 7 391

3а

39"

39"

К (к) = Ф(к)К(к), критерий (13) можно привести к виду

N - 1

/1 = £ [^0(к) + ха (к)^(к) + иТ(к)^2(к) +

к = 0

+ ХІ (к)^1(к)и(к) + ХІ (к)^2(к) ХА (к) +

+ и (к)Л(к)и(к)],

в котором

^(к) = ЬА (А)^а 1(к + 1)Ьа (к) +

+ 2 ^ (к^А2(к + 1)АА(к + 1);

^(к) = 2 ФІ (к^А2(к + 1)АА(к + 1) +

+ 2 ФІ (к)^11(к + 1)ЬА(к);

^(к) = 2 (к)^А2(к + 1)АА (к + 1) +

+ 2 < (к)^11(к + 1)ЬА(к);

^(к) = 2 ФІ (А)^а 1(к + 1) (к);

ВД = ФІ (А)^а 1(к + 1)Фа (к);

^э(к) = < (А)^а 1(к + 1)^(к).

Поскольку ^0(к) в соотношении (15) не зависит от и, при максимизации критерия /1 это слагаемое можно не учитывать. В результате приходим к критерию

N -1 Т Т

/2 = £ [ХА (к)^(к) + иТ(к)^2(к) + ХА (к)^1(к)и(к) +

к = 0

+ ХІ (к)^2(к) ХА (к) + иТ(к)^3(к)и(к)],

который вместе с системой (14) определяет задачу синтеза оптимальных входных сигналов как задачу дискретного оптимального управления с суммарным показателем качества [14]. Для ее решения воспользуемся симплексным методом Нелдера— Мида.

Сформулированную задачу можно свести к следующей задаче оптимизации конечного состояния [14]:

Будем считать, что выполнены все априорные предположения (см. § 1), причем Е [м(к)м(')] = 0,68к = 08к; Е [у(к + 1)у(г + 1)] = 0,68К = Я8К; х(0) е N(0; 0,01).

Выполнив линеаризацию модели (16), (17) относительно номинальной траектории

Хн( к + 1) = [1 - т2) хн (к) + -0— ехр {0,25 [ мн( к) - хн( к)]},

V 01

к = 0, 1,..., N - 1;

Хн(0) = 0,

получим линеаризованную модель вида

а(к) = 0^ехр{0,25[мн(к) - хн(к)]} х

01

X {1 - 0,25[ин(к) - хн(к)] + 0,25и(к)};

Ф(к) = 1 - 02 - 0025 ехр{0,25[ин(к) - хн(к)]};

0

0

*а(к+ 1) ,Х(к+1).

= Х^) ^ шах

и Є О тт

Фа

Х( к) + Х^( к )^і( к )и (к) +

+ Х^( к) _02 (к) хА (к) + ит( к) _03( к) и( к) +

- т т

. + хА

к = 0, 1,..., N - 1;

Г-Т л 1

ха

-Х(0)-

дх (0 -д 0 1

0

дх (0 -

д 0

Для решения этой задачи воспользуемся методом последовательного улучшения управлений Л.И. Шатровского [15, 16], адаптировав его к дискретной задаче.

5. ОПТИМАЛЬНОЕ ОЦЕНИВАНИЕ ПАРАМЕТРОВ ОДНОЙ МОДЕЛЬНОЙ СТРУКТУРЫ

Рассмотрим следующую модель нелинейной дискретной системы:

х(к + 1) = Г1 - х(к) + 0^1ехр{0,25[м(к) - х(к)]} +

Г 01 01

+ 0г м

01

у(к + 1) = х(к + 1) + у(к + 1),

к = 0, 1, ..., N - 1,

где 01 и 02 — неизвестные параметры системы, причем

3 < 01 < 10; 0,05 < 02 < 1,25.

А(к + 1) = 0; Н(к + 1) = 1.

Таким образом, необходимо оценить параметры 01 и

02, входящие в выражения для а(к), Ф(к) и Г(к).

Считая, что для номинальной траектории (18) ин(к) = и(к), к = 0, 1, ..., N - 1, обеспечим наилучшее приближение построенной линеаризованной модели к своему нелинейному аналогу.

Определим область допустимых входных сигналов

= {и е < и(к) < 20, к = 0, 1, ..., N - 1}.

Для ослабления зависимости результатов оценивания от выборочных данных произведем пять независимых запусков системы, и усредним полученные оценки неизвестных параметров. Реализации входных и выходных сигналов получим компьютерным моделированием при истинных значениях параметров 01 = 4, 02 = 0,5 и N = 31. Результаты выполнения процедуры активной параметрической идентификации представлены в таблице.

При решении реальных задач истинные значения параметров неизвестны. Будем судить о качестве оценивания в пространстве откликов по значению коэффициента

к - II ^ - ^р || _ кУ ІІ

К - м

г-1 - 2

2

_ |к=_о___________________________________

2

I к = 0

где 7Ср _ {уСр(к + 1), к _ 0, 1, ..., N - 1}, 7ср _ {Уср (к +

+ 1|к + 1), к _ 0, 1, ..., N - 1}, 7с*р _ {У*р (к + 1|к + 1), к _ 0, 1, ..., N — 1} — усредненные по всем запускам последовательности измерений для вектора 0, равного 0*,

1

Результаты активной параметрической идентификации

Входной сигнал

Номер запуска Оценки параметров

01 0 2

1 4,686 0,432

2 3,277 0,553

3 5,067 0,423

4 6,482 0,609

5 7,424 0,282

Средние значения 5,387 0,460

1 4,190 0,561

2 4,043 0,502

3 3,785 0,461

4 3,803 0,498

5 4,495 0,567

Средние значения 4,063 0,518

1 4,695 0,492

2 4,231 0,506

3 3,139 0,497

4 3,625 0,522

5 4,723 0,538

Средние значения 4,083 0,511

Исходный

Синтезированный методом Нелдера — Мида

Синтезированный методом Шатровского

Рис. 1. Тестовый сигнал и для анализа качества прогнозирования на основе результатов, приведенных в таблице

0 ср, 0 ср соответственно, при входном сигнале и = {и(к), к = 0, 1, ..., N - 1}, представленном на рис. 1.

Определим при помощи равенства у (к + 1|к + 1) =

= А(к + 1) + Н(к + 1) х (к + 1|к + 1) величины {у (к + + 1|к + 1), к = 0, 1, ..., N - 1} для вектора 0, равного 0*,

0ср, 0Ср . Полученные результаты представлены на рис. 2.

Сопоставление рис. 2, б и в показывает практически одинаковое качество прогнозирования при оценках параметров, найденных по синтезированному входному сигналу методами Нелдера—Мида и Шатровского. В обоих случаях коэффициент к¥« 1,36. Учитывая, что метод Шатровского позволяет получать, вообще говоря, не оптимальное, а достаточно хорошее допустимое управле-

6

5

4

3

YCp

^ср

О Н—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і

1 5 9 13 17 21 25 29 к

а

о н—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і—і

1 5 9 13 17 21 25 29 к

б

1 5 9 13 17 21 25 29 к

6

к

Рис. 2. Оценки Ycp и Ycp, соответствующие исходному

ние, применение методов нелинейного программирования должно обеспечивать более качественные результаты оценивания.

ЗАКЛЮЧЕНИЕ

Для гауссовских нелинейных дискретных систем рассмотрена и решена задача активной параметрической идентификации с применением методов теории оптимального управления. Рассмотрен случай вхождения неизвестных параметров в уравнения состояния и наблюдения, начальные условия и ковариационные матрицы помех динамики и ошибок измерений. Разработана процедура оптимального оценивания параметров, целесообразность применения которой продемонстрирована на модельном примере. Выполненные исследо-

вания показали, что применение методов оптимального управления при решении задач активной параметрической идентификации стохастических динамических систем не только принципиально возможно, но и дает положительный эффект.

ЛИТЕРАТУРА

1. Цыпкин Я.З. Информационная теория идентификации. — М.: Наука, 1995. — 336 с.

2. Льюнг Л. Идентификация систем: Теория для пользователя. — М.: Наука, 1991. — 432 с.

3. Федоров В.В. Теория оптимального эксперимента (планирование регрессионных экспериментов). — М.: Наука, 1971. — 312 с.

4. Ермаков С.М., Жиглявский А.А. Математическая теория оптимального эксперимента. — М.: Наука, 1987. — 320 с.

5. Денисов В.И., Еланцева И.Л., Чубич В.М. Активная идентификация стохастических линейных дискретных систем, описываемых моделями в пространстве состояний и ARMAX — моделями // Сибирский журнал индустриальной математики. — 2о0о. — Т. 3, № 1(5). — С. 87—100.

6. Денисов В.И., Чубич В.М., Черникова О.С. Активная параметрическая идентификация стохастических линейных дискретных систем во временной области // Там же. — 2003. — Т. 6, № 3(15). — С. 70—87.

7. Активная параметрическая идентификация стохастических линейных систем: монография / В.И. Денисов и др. — Новосибирск: Изд-во НГТУ, 2009. — 192 с.

8. Казаков И.Е., Доступов Б.Г. Статистическая динамика нелинейных автоматических систем. — М.: Физматгиз, 1962. — 332 с.

9. Синицын И.Н. Методы статистической линеаризации // Автоматика и телемеханика. — 1974. — № 5. — С. 82—94.

10. Astrom K.J. Maximum likelihood and prediction errors methods // Automatica. — 1980. — Vol. 16. — P. 551—574.

11. Огарков М. А. Методы статистического оценивания параметров случайных процессов. — М.: Энергоатомиздат, 1980. — 208 с.

12. Mehra R.K. Optimal inputs for linear system identification // IEEE Trans. Automatic Control. — 1974. — Vol. AC-19, N 3. — P. 192—200.

13. Чубич В.М. Вычисление информационной матрицы Фишера в задаче активной параметрической идентификации стохастических нелинейных дискретных систем // Научный вестник НГТУ. — 2009. — № 1 (34). — С. 23—40.

14. Пропой А.И. Элементы теории оптимальных дискретных процессов. — М.: Наука, 1973. — 256 с.

15. Шатровский Л.И. Об одном численном методе решения задач оптимального управления // Журнал вычислительной математики и математической физики. — 1962. — № 2. — С. 488—491.

16. Афанасьев В.Н., Колмановский В.Б., Носов В.Р. Математическая теория конструирования систем управления. — М.: Высшая школа, 1998. — 574 с.

Статья представлена к публикации членом редколлегии

В.Н. Афанасьевым.

Чубич Владимир Михайлович — канд. техн. наук, доцент,

И сhubich_62@ngs.ru,

Черникова Оксана Сергеевна — канд. техн. наук, доцент,

И chernicova@ngs.ru,

Новосибирский государственный технический университет,

Ш (383) 346-27-76.