Уравнения методов идентификации линейных дифференциальных уравнений

Математические заметки СВФУ Июль—сентябрь, 2014. Том 21, № 3

УДК 517.925.54:517.962.27/.8

УРАВНЕНИЯ МЕТОДОВ ИДЕНТИФИКАЦИИ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ А. О. Егоршин

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

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

1. Введение

Идентификацией называют оценивание параметров уравнений (дифференциальных или разностных) некоторого (проектируемого, изучаемого, управляемого, прогнозируемого) динамического процесса по экспериментальным или синтезированным (исходным) данным. Иногда термином «идентификация» называют оценивание характеристик интегрального описания динамического объекта, например, его импульсной функции, функции Грина.

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

Задачи идентификации — построение (оценивание характеристик) математических моделей по исходным данным, связанным с исследуемыми объектами — относятся к обратным (inverse) задачам моделирования в отличие от прямых задач моделирования подобия (simulating) — воспроизведения поведения моделируемого процесса по известным (заданным) его описаниям (часто

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта № 13—01—00329) и Сибирского отделения Российской академии наук (междисциплинарный проект №80).

© 2014 Егоршин А. О.

приближенным). Обратные задачи моделирования мы называем математическим моделированием. Предварительный выбор класса моделей при постановке задач математического моделирования есть один из способов регуляризации некорректных обратных задач восстановления операторов реальных физических объектов.

Идентификацией иногда называют специальный класс задач: математические задачи оценивания характеристик динамических объектов при стохастических исходных данных [1, с. 242]. В частности, это задачи математического моделирования, в которых заданы и модели измерительных ошибок как случайных переменных или процессов. В таких задачах структура (вид уравнения) объекта считается известной и совпадающей со структурой модели. В подобных задачах оценивания возможно, по крайней мере в принципе, строгое исследование статистических свойств получаемых оценок. Мы используем термин «идентификация» в более широком смысле: как математическое моделирование динамических объектов или процессов, имеющих неизвестное и, возможно, более сложное описание, чем принятая модель. Последнее, очевидно, типичная ситуация для большинства приложений.

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

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

Чаще всего (в частности, и по аналитическим соображениям) используются среднеквадратические критерии качества аппроксимации объекта более простой моделью [2, с. 10,16, 201; 3, 4]. При этом задачи математического моделирования и идентификации становятся задачами проектирования (поиска ближайшего элемента) на некоторые допустимые множества модели в соответствующем евклидовом или гильбертовом пространстве исходных данных [5].

2. Обозначения

Пусть на интервале наблюдения I = [0,4], задана равномерная Л-сетка Д из Ь + 1 точек Тг = гК, г = 0, Ь (т. е. Ь = 1/1г) и на ней заданы (возможно, с неконтролируемыми ошибками, измерительными и/или структурными, см. выше) значения (отсчеты) у = у(т^) некоторого решения у(т), т € I, изучаемого дифференциального уравнения или некоторого исследуемого динамического

или

процесса. Отсчеты образуют вектор y = yo = {y^lo G E = El = Eo+1 финитную последовательность y G l2[0, L] = l2. При квадратичной норме в E и

o

l2 (l|y||2 = 2 | yi |2) это эквивалентные объекты. о

Скалярное произведение элементов x, y в E и I2 обозначаем через (y, x) = (y, x)r = x* Ry, где самосопряженная матрица R > 0, что обозначает положительную определенность. Инволюция x* будет обозначать композицию двух инволюций: транспонирования и комплексного сопряжения. Символ x' будет обозначать элемент сопряженного пространства E'. Тогда x' = x*R = (•, x). Для операторов K : E ^ E имеет место равенство K' = R-1K*R. Считаем далее для простоты, что R = lo = I — единичная матрица. Тогда (y, x) = y'x = x*y, x' = (•, x) = x*, K' = K*. Полагаем, что ||y||2 = ||y||| = (y, y) = y'y = y*y.

Функционалы в E определяем через скалярное произведение. Орты в E обозначаем через e¿, а орты в Е' — через = (-,e¿) = е*, i = 0,L. Тогда координаты y i векторов y = {уг}° G E определяются как функционалы: y = (y, ei) = eiy = e*y.

Условимся об обозначениях для векторов и матриц как множеств. Конструкции вида {хг|т будут обозначать вектор-столбец, |xj |™ или |x;,..., xn| — вектор-строку, а {xj |mn — матрицу (i — номера строк) множества элементов вида x. Если элементы xi в {x¿Ц™ — (n — l + 1)-векторы-строки, а xj в |xj |™

или в |x;,...,xn| — (m — k + 1)-векторы-столбцы, то эти конструкции — матрица

rx. .lmn

{xij .

Через Е^гу = ..., e¡\ = \ei\lk обозначаем (L+1) х (I — к+1)-матрицу — совокупность соответствующих ортов в Е. Это также и оператор E^j : El~k+1 —> Е. Совокупность соответствующих ортов из E' есть (l — k + 1) х (L + 1)-матрица {e'i\k- Это также оператор Е^-у : Е —> El~k+1. Будем считать, что e¿ — единичные орты с одной единицей на i-й позиции: ei = {Si¡j(áijj — символ Кронекера). Тогда можно писать: E¿j = Ещ.

Поскольку мы приняли, что R = I, то Eq-j- = Е = Е = I¿, = I — единичная L + 1-матрица, а

Щц = Щц = lOi-fe+i.fe'^-fe+bOi-fe+i.L-il : Е

— так называемая «вырезающая» матрица. Используем обозначение 0kti для нулевых (k х 1)-матриц и 0k для k-векторов. Введем также оператор сдвига «вниз» такой, что 11e^ = e^+1. Его матрица есть 11 = {Si-1¡jпричем I* = I-1. Оператор сдвига определяют еще и так: sy& = yk+1. Легко видеть, что

I-1 = s.

3. Вариационная задача аппроксимации

Задача математического моделирования процесса у ставится как вариационная задача его кусочно-линейной динамической аппроксимации. Так названа задача наилучшего среднеквадратического приближения функции у в конечном интервале функцией у, удовлетворяющей в этом интервале линейному разност-

ному уравнению (РУ) с постоянными коэффициентами:

Оу = £>*?=/ Ук+гаЛ =0=| ]Г

I г=0 ) 0 I г=0 ) 0

По некоторым п локальным условиям (начальным [у]0 = (Уг|о-1, конечным [у]№+1 = {Уг}^+1 или иным координатам в кег О), определяющим проекцию уа € кег Оа С Е, а также по коэффициентам а* = |а*|0 = |а0,...,а0| уравнения (1) процесса у минимизируется расстояние

0

Задача (1), (2) обобщает классическую задачу п-полиномиальной аппроксимации. В последней задаче используется также критерий (2). В качестве условий минимизации вместо РУ (1) используется уравнение

о

дпу = 0 = £ Ук+г(-1)°-гсо. 0

Это уравнение с известными коэффициентами есть частный случай уравнения (1), в котором а! = |(-1)о-гС010.

Вектор а € ш (ш — допустимое множество для значений а) оценивается в

(1), (2) как направление в Ео+1. Если предполагать, что ао = 0, то можно считать, что ао = 1. Тогда ш С С = {$ € Е"+1 : = 1} и оцениваемыми параметрами уравнения в задаче (1), (2) будут коэффициенты а* = |а*10 1 = |а*,..., а0-11. Более надежна (не требует подобных априорных предположений) нормировка к единичной длине вектора коэффициентов: ||а|| = 1. Тогда ш С 8 = {^ € Еп+1 : ||^|| = 1}.

Минимизация функционала (2) при заданном векторе а коэффициентов уравнения (1) называется задачей сглаживания. Минимизация функционала

(2) и по вектору а коэффициентов уравнения (1) называется задачей идентификации.

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

о

некоторого решения дифференциального уравнения ^ У(г)а* = 0 в интервале

г=0

74, то ,7 = 0, у = у и задача (1), (2) представляет собой вариационный (в отличие от аналитического, основанного на теореме Гамильтона — Кели) метод равномерной (в интервале 74) дискретизации этого дифференциального уравнения [5, 6].

Пространство Е Э у называется пространством исходных данных, подпространство Ф = Фа = кег Оа С Е размерности п — ядро оператора О — называется подпространством модели. Вектор у называется исходной реализацией, вектор у — сглаженной, а вектор уа — оптимальной сглаженной реализацией. Она дает минимум функционалу (2) для заданного значения коэффициентов а в уравнении (1).

Введем в задачу (1), (2) параметр М < Ь, обозначающий длину М + 1 начального отрезка ум = {Уг}о1, М = Ьо,Ь, исходной реализации у. Здесь Ь0 + 1 > 2п — длина некоторого начального ее отрезка, для которого решение задачи (1), (2) может существовать и быть единственным. Соответствующие подзадачи (1), (2) с параметром М называем далее М-задачами, а их решения обозначаем через уа (М) и а(М) и называем (частичными) М-решениями.

Одна из целей статьи — получить приближенные разностные уравнения для вектор-функции а(М), М = Ьо,Ь (разд. 8). О минимальном значении 2п величины Ь0 + 1 см. ниже замечание 2.

Лемма 1. Для получения частичных М-решений (М < Ь) в задаче (1), (2) необходимо и достаточно в равенствах (1) положить Ь = М, М = Ьо, Ь.

Доказательство. Если Ь = М в (1), то на отсчеты уг, I = М + 2, Ь, в задаче (1), (2) ограничения (1) отсутствуют. Поэтому для г = М + 2,Ь имеют место равенства Уг = Уг. Следовательно, || • ||| = || • |||м+1.

Этот факт позволяет при постановке и решении М-задач (1), (2) не переходить из Е = Еь в Ем при М < Ь, а решать их, оставаясь в пространстве Е, положив в (1) Ь = М. Будем обозначать число условий (1) в этом случае через К = М - п, где М = Ь0,Ь, К = М0,М, М0 = Ь0 - п > п - 1, Ь0 > 2п - 1.

В статье речь идет о скалярных действительных уравнениях: коэффициенты сц, г = 0, п, и отсчеты у;, у; I = О, Ь, в (1), (2) считаются действительными скалярами.

4. Задача ортогонального проектирования

Представим вариационную задачу (1), (2) как задачу проектирования в евклидовом пространстве Е. Для этого запишем условия (1) минимизации функционала (2) в матричном виде. Это можно сделать двумя принципиально (как будет видно из дальнейшего) различными способами. Они существенны для решения задачи (1), (2).

Используем специальные сдвиговые матрицы, образуемые коэффициента-

ми и отсчетами:

А = Ак =

а0 0 0

0

ао а0

0

0 0 ао

V = VN =

У0, У1, . . . , У о 1 , Уо

У1 , У2, . . . , Уо , Уп+1

У2, Уз, . . . , Уп+1, Уп+2

Уз, У4, . . . , Уп+2, Уп+3

Ум, У№+1,

Уь-1,

Уь

Ленточнаятеплицева ((Ь+1)х(Ж +1))-матрицаА = Ам = А(а) (А7 = А(са)). Матрицу А = А(а) = Аа = Ам называем далее матрицей скользящего вектора а (МСВ а):

А = Ам = Аа = Ам(а) = |п0,П1,... |

где Пк = 7кП0 € Еь+1, П0 = |а*, 0%|*.

МСВ а образуется (п + 1)-вектором, поэтому число ее столбцов N + 1 на п меньше, чем строк. Если образующий вектор МСВ, например Л, имеет размерность N +1, то МСВ Л будет ((Ь + 1) х (п + 1))-матрицей, т. е. число ее столбцов п +1 будет на N меньше, чем строк. МСВ а — А(а) — и МСВ Л — Л(Л) (как и матрицы Ак) — встретятся далее. Последняя есть МСВ вектора множителей Лагранжа. В теореме 2 он определяет реализацию «ошибок» Луа = У — Уа • В М-задачах при М < Ь соотношение числа строк и столбцов будет иным. Столбцы МСВ есть специальные базисы некоторых подпространств в Е = В силу леммы 2 в М-задачах нет необходимости переходить в

пространства реализаций меньшей размерности М + 1.

Вторая специальная матрица в (3) — ганкелева (^ + 1) х (п + 1))-матрица V = Уу = V(у). Столбцы матрицы V* суть (п + 1)-векторы

^ = Ыкк+п = \у}Нп+г), к = 1Щ,

отсчетов исходной реализации у. Столбцы щ, к = О, -/V, матрицы V* = У* (у) есть (п + 1)-выборки [у]к(п+1) = {Уг}к+п отсчетов из сглаженной реализации у. Они формируют модель (1).

Введем так называемый вектор невязок уравнения модели (1). Это ^ + 1)-вектор значений оператора О = Оа уравнения (1) на исходной реализации у:

т = та = °ау = {а*Ык(п+1)}у = {а*«к}у

= {Ка)Еп+1 }У = {

Лемма 2. Для образа т = Оау действительного оператора Оа уравнения (1) на действительной исходной реализации у имеют место следующие тождества [3, 5]:

т = А*у = Уа или т =

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

Очевидные теперь факты констатирует следующая теорема и ее следствия.

Теорема 1. Задача (1), (2) является параметрической задачей проектирования в Е:

минимизировать ||у — у||2 при условиях А*у = 0 = Уа,

где приведенные два условия эквивалентны в Е.

Следствие 1. Вариационная задача проектирования (1), (2) —условной минимизации функционала (2) — является задачей поиска ближайшего элемента в множестве

£ = {Ф С Е : Ф = Фа = кега е ш}

допустимых (определяемых множеством ш допустимых значений вектора коэффициентов а) подпространств Ф. Они суть ортогональные дополнения Б± е Е к замкнутым линейным оболочкам Б = Б(А) = Ба, где А = Аа — МСВ а коэффициентов РУ (1).

Следствие 2. Пусть М — п = К < N. Тогда в М-задаче А = Ак = |пг|К — ((Ь + 1) х (К + 1))-матрица, V — ((К + 1) х (п + 1))-матрица, а на реализацию у длины Ь + 1 равенствами (1) (или (7)) наложено К + 1 < N +1 условий.

Два вида условий (1) задачи (1), (2), заданные в форме (7), позволяют, во-первых, легко дифференцировать эти условия по двум видам переменных этой задачи — сглаживания у и идентификации а. Во-вторых, они позволяют сопоставить задачу (1), (2) с известной более простой задачей оценивания и идентификации в Еп+1. Это сделано в следующем разделе.

С помощью так называемых «полных» (п + 1)-выборок -к = [у]к(п+1), к = 0, N (столбцов матрицы V*), а также «коротких» п-выборок — п-векторов состояний [у]к(п), — уравнения (1) могут быть записаны следующим образом:

^аУ = {а*[у]к(п+1)}Г=о = {а*-к }£=0 = 0

или

где

Ук+п

0

Ук+гаг = -а* [у\к{п) = -а*ьк, к = 0, Ж,

[у]к(

к

{у}к+п

-к,

[у]к

шк+п—1

«к-

Короткие п-выборки \у\к(п) = ^к, к = 0, N + 1, суть состояния модели (1). С помощью рекуррентной записи (9) уравнения (1) они позволяют вычислять реализацию у = у N = {Уг}^ последовательно, начиная с начальных условий [У]о = [У]о(п), при условии, что ап = 1, т. е. а € С Э ш. Значению к = N + 1 соответствует граничное состояние [у]№ +1(п). Его прогноз на отсчет уь+1 по формуле (9) выходит за границу интервала наблюдения 74 = [0, ЬЛ].

Итак, ганкелева (^ + 1) х (п + 1))-матрица выборок V из (3) в условиях минимизации задачи проектирования (7) может быть записана в виде

V = ^^ = К}£=о = {[у]

)}к

]к(п+1) }к=0

V = {у*к,Ук+п}к=о

{[у]к(п) ,ук+п}Г=о

ПЛП

где

У = Ум = {[у]*к

Матрица V коротких выборок ук — это матрица состояний модели

в интервале наблюдения 74 на сетке 7^ (кроме последнего для к = N + 1). Эти состояния определяют в (9) компоненты реализации у — последнего столбца матрицы V: у = — Уа в (12).

п— 1

5. Альтернативные подходы

Укажем на два известных альтернативных подхода к оцениванию коэффициентов уравнения (1). Первый подход — алгебраическая идентификация (АИ) — предполагает минимизацию некоторого функционала невязки (6) т = V(у)а

уравнения (8). Алгебраический подход включает в себя многочисленные и разнообразные группы относительно простых методов разомкнутой идентификации [3,6]. Алгебраический подход является самым простым и известным подходом к идентификации объектов вида (1) с тех пор (полвека назад), когда возможности вычислительной техники были минимальны. Второй подход предполагает минимизацию рассогласований выборок (метод М. Левина [1, с. 294]). Подобные подходы к оцениванию (когда считается искаженной ошибками матрица системы уравнений) называют методом ортогональной регрессии (ОР) [2, с. 287], тотальным методом наименьших квадратов (МНК).

Пусть ||Vу2 = 8р(У*У) — квадратичная норма матрицы. Тогда упомянутые три подхода (ВИ, ОР, и АИ) к оценке коэффициентов уравнения (1) наглядно и кратко могут быть сопоставлены с помощью следующей простой таблицы.

Подход Критерий Условия

ОР Зот = 9а = 0

АИ Зы = |ЩУН12 У-У = 0

Из определений (13) сразу вытекают следующие утверждения.

Лемма 3. (а) В задаче ВИ проектируется одна реализация у £ Е = Еь+1 длины Ь + 1 на подпространство размерности п в Е. В задаче ОР проектируется N + 1 выборок Vдлины п +1 на подпространство а± размерности п в Еп+1.

(b) В задачах ОР и АИ строки матрицы V — полные (п + 1)-выборки содержащихся в у отсчетов — рассматриваются как независимые векторы из Е п+1.

(c) Задача ОР — сумма N + 1 независимых (п + 1)-задач ВИ для строк матрицы V с одним оцениваемым вектором а.

Доказательство. Первое утверждение очевидно из первых двух строк в (13). Второе следует из того, что строки матрицы V в задачах ОР и АИ в (13) могут быть переставлены произвольным образом без изменения результатов решения этих задач. В задаче ВИ ганкелева структура не только исходной матрицы выборок V, но и сглаженной V однозначно определена тем, что отсчеты сглаженной реализации у однозначно упорядочены критерием оптимизации (расстоянием в Е) и теплицевой матрицей А* в левом тождестве (13) ВИ (лемма 2). Третье утверждение следует из того, что в задаче ВИ проектируется на подпространство (А) С Е размерности п одна (Ь + 1)-реализация из Е = Еь+1, а в задаче ОР на подпространство (а) = а± С Еп+1 (тоже размерности п) проектируются, как независимые, «облако» [2, с. 20] из N + 1 векторов в Еп+1 (столбцов матрицы V*).

Следствие 1. В задаче ВИ для определения переходного процесса (сглаженной реализации) у длины Ь +1 необходима и достаточна оптимизация п переменных — п-вектора его начальных (или любых других локальных, если ао = 0) условий. В задаче ОР — проектирования N + 1 выборок из Еп+1 на а± — число оптимизируемых начальных условий равно (^ + 1)п.

Доказательство. Число оптимизируемых переменных при сглаживании ВИ и ОР легко подсчитывается из формул (9) и вида матрицы V.

Определение 1. Модель М(а, у) называется замкнутой (с обратной связью), если представление ею реализации у в пространстве Е исходных данных у есть решение задачи Коши.

Следствие 2. В задаче АИ (13) исходная реализация у считается «выходом»» модели. Сглаживание отсутствует. Оптимизируемыми параметрами являются коэффициенты а.

(Ь) Ошибки выхода в виде невязок т = Vа (5) перенесены на вход модели.

(^ В модели М(а, у) АИ (14) отсутствует обратная связь.

Доказательство. В методах АИ, как следует из (13), в качестве модели отсчетов реализации ук (см. (11) и (12)) используется алгебраическое соотношение (возмущенное скользящее среднее):

п — 1

Ук+П =-^Ук+гаг+тк+1, к = О, Ж.

о

Минимизируется по а «вход» модели (14) — величина ||т||Ем+1. Ошибки в исходной реализации у отсутствуют. Поэтому нет и задачи сглаживания — вычисления реализации у.

(Ь) Соотношение (14) внешне похоже на модель (1) в виде (9), но имеет иной физический и информационный смысл. Ошибки в исходной реализации у перенесены в этой модели на ее «вход» т, подходящий функционал которого минимизируется. Это следует из (13) и (14).

(^ Из постановки задачи АИ в (13) и из равенства (14) видно, что исходная реализация у, используемая как выходная последовательность модели (14), не есть решение задачи Коши. Обратная связь, с помощью которой в (9) решается задача Коши для уравнения (8), в формуле (14) «разрушается» произвольной последовательностью — невязкой т. Она используется как «вход» модели (14).

Перенос ошибок с «выхода» на «вход» некорректен для физически реализуемого оператора. Это приводит к малой устойчивости методов АИ к ошибкам в исходных данных. Особенностью методов АИ является их сравнительная простота в использовании и в анализе. Например, из (14), используя определения (11), (12), получаем несовместную систему уравнений Vа « 0 и соответствующие оценки МНК, минимизирующие невязку т = V(у)а (6), (13):

Уа = т^ 0—>у&-Уа—> а =-

Замечание 1. Известно, что эти оценки смещены из-за ошибок в V [1, с. 283]. Поэтому методы АИ обычно предваряет обработка исходной реализации у для уменьшения ошибок в ней. Наилучший способ такой обработки — применение модели, адекватной исследуемому процессу. В нашем случае это модель (1). Именно совместные сглаживания у и оценивание а осуществляется в задаче ВИ. Это делается и в задаче ОР, но для строк {V*матрицы V.

Замечание 2. Из определений (11), (12) и формул (15) следует, что в лемме 1 Ьо + 1 > 2п. Здесь 2п — минимальное число отсчетов, необходимое для вычисления некоторой оценки вектора коэффициентов а по формуле (15). Эти отсчеты должны быть таковы, что квадратная (пхп)-матрица Vп-\ (см. (12)) — невырожденна. При этом минимальном числе таких отсчетов три оценки из таблицы (13) совпадают. Они вычисляются по формуле: а = — ^п_1[у]п(п).

Оценки а(М) всех частичных М-задач (1), (2) ВИ, ОР и АИ совпадают не только при Ь = 2п (и невырожденности матрицы Уп—{), но и еще в одном случае.

Лемма 4. Пусть реализация у — точное решение РУ (1) с Ь = М и некоторыми коэффициентами а = |ао, - - -, ап—1, 1|. Кроме того, пусть это решение таково, что начальная матрица выборок Уп_1 невырожденна. Тогда решения задач идентификации ВИ, ОР, АИ (13) для М = Ьо, Ь совпадают и могут быть получены по формуле (15).

Доказательство. При условии леммы Vа = Vа = 0, так как у = у. Условия минимизации задач и ВИ, и ОР, как видно из табл. (13), выполнены, все функционалы ^ минимальны и равны нулю на решении а. В условии

леммы для всех М = Ьо, Ь матрицы выборок Ук, К = М — п, имеют ранг п. Поэтому решение (15) единственно для всех задач из (13).

6. Аналитическое решение задачи ВИ

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

Формулы ортогонального проектирования на подпространство при заданном в нем или в его ортогональном дополнении базисе хорошо известны. Они легко могут быть получены из формулировки задачи (1), (2) в виде (7). Обозначим через С = С^ = (А, А) матрицу Грама системы векторов А = AN.

Теорема 2. (а) Формулы ортогонального проектирования в задаче (1), (2) могут быть получены с помощью метода множителей Лагранжа (МЛ); вектор МЛ равен А = С—1(у, А).

(Ь) Пусть П — проектор на кег = Фа = ^ = Б(А^), а Р = Р(А) = 7 — П — проектор на Б = Б(А) = Е 0 Фа, N = Ь — п. Тогда сглаженная реализация уа = уа(Ь) и перпендикуляр Дуа = у — уа могут быть записаны в виде

у а = Пу, П = Па = П№ = 7 — Р, Р = Ра = Р№, Дуа = у — уа = АА = Ру,

где

P = P

(c) Для частичных M-решений у2(М) и Ду2(М) M-задач (1), (2) имеем в (16), (17): А = Ак, П = Як, Р = Рк, К = М — п = Ь0 - п, N.

Важные формулы для квадрата длины перпендикуляра Ду2 = y — у2 дает

Теорема 3. Значение J функционала J (2) на проекции у2 (16) исходной реализации у длины L + 1 на подпространство Ф2 = ker Da = S± С EL+1 — ядро разностного оператора Da модели (1) — определяется равенствами

J = Ja = p2

Следствие 1. Длина p = ||Ду2 У = p(a) перпендикуляра Дуа = у — у2 не зависит от параметров сглаживания задачи (1), (2), а зависит только от вектора коэффициентов а.

Функцию p2(a) = J(a) от а (18) назовем функционалом идентификации.

Следствие 2. (а) Задача оптимизации модели (14) по вектору ее коэффициентов а и эквивалентная ей задача АИ из (13) становятся эквивалентными задаче ВИ (1), (2) из (13), если нормы в пространстве EN+1 невязок Va = m G EN + 1 определены как || * ||c-i, где C = C(а) — теплицева матрица Грама системы векторов A^) (3) из формул проектирования (17), (18) теоремы 2.

(b) Тогда критерии минимизации невязки m G EN+1 задачи АИ в (13) оптимизации коэффициентов а модели (14) суть функционалы идентификации в En+1: |^а||£-1 = ||m||C-i.

Доказательство . Результат теоремы 3 сводит вариационную задачу идентификации из (13) к задаче безусловной минимизации функционала идентификации (18). Вид функционала (18) доказывает следствие 2 теоремы 3.

Теорема 4. Задача оптимизации неизвестных коэффициентов РУ в вариационной задаче аппроксимации (1), (2) сводится к задаче безусловной минимизации следующего функционала идентификации в E = EL+1:

p2^) = J

Равенства теоремы 4 и тождества леммы 2 приводят к следующему результату.

Теорема 5. (а) Функционал идентификации, определенный в следствии 2 теоремы 3 как функционал в E N +1, в теореме 4 как функционал в E = EL+1, может быть представлен также как функционал в En+1 на сфере.

(b) Он может быть приведен к виду квадратичной формы с неконстантной по а на сфере S(c) = {а G En+1 : ||а|| = с} идентифицирующей матрицей Q(^):

J = Jvi = Ja = p2a = a*Qa = Ja

где

Q = Qvi = Q(^) = V *(A*A)-1V = V*C-1V.

Вектор а = а(М) = arg min Ja(M) = Jvi, на котором достигается глобальный минимум функционала идентификации, назовем решением частичной вариационной M-задачи идентификации (1), (2) для реализации длины M + 1,

Замечание 3. В задаче ОР из (13) имеем такую аналогию с формулами (20) [3,5,7]:

Jor = а* Qora, Qor = V *V/||

а

I 2

7. Уравнения сглаживания и фильтрации

При реализации полученных решений возникают следующие вопросы: нахождение минимума функционалов идентификации Ja(M) и, если велика длина реализации L (число N = L — n), обращение матрицы Грама C. Основу соответствующих рекуррентных (по M) алгоритмов составляют уравнения двусторонней — иначе встречной (counter) — ортогонализации Грама — Шмидта [5, 7-9]. Ниже даны уравнения последовательного решения M-задач (1), (2) для M = L^L.

Введем обозначения для проекторов P^j. Так будут обозначаться проекторы на подпространство Syj = S(Ayj) — линейную замкнутую оболочку указанных столбцов матрицы A: A-^j = \щ,..., щ |. Обозначаем Р-^ = Pk, А= Ак, Sq-^ = Sk Аналогичными индексами снабжаем проекторы П^гу = 7 — P~kl, IIq-д: = Щ на ортогональные дополнения Ф^гу к подпространствам Syj- Очевидно, что П_1 = I, P-i = 0.

Определим ортогонализирующие векторы fk, и fk, к = 0, N, встречных прямого (forward) и возвратного (backward) процессов ортогонализации Грама — Шмидта системы векторов Ак- В М-задачах векторы fk, где к = 0, К, К = М — п, определяют цепочку проекторов П^, К = 0, N, М-задач (1), (2). Здесь К+ п = М = L0,L.

По определению процессов последовательной ортогонализации [7]

к

/к=Як-1т^Як=1-^М/г,/г)-1

0

Теорема 6. (а) Пусть ак = \\fk\\~2, а-к = ||Д|Г2, а /о = /о = По- Ортогонализирующие векторы fk+i, fk+i Для k = — 1, N — 1 вычисляются с помощью следующих нелинейных уравнений процесса двусторонней (встречной) ортогонализации:

где

fk+i = 7- Л0*+1, ffc+1 = ffc - 7

ök+1 = a/Mfc+1, M/+1 = (//, 11/k), a/+1 = a/+1 = (7 - 0fe+10*+1)-1afe, ao = 11 no N 1 = ||a|-2.

(b) Обращаемые числа 1 — \dk+i |2 > 0 не равны нулю для всех к = 0, N — 1, если вектор а коэффициентов РУ (1), имеет хотя бы одну ненулевую компоненту, т. е. rank An = N + 1.

Следствие 1. Проекция ук+1 = у есть решение М-задачи сглаживания (1), (2) для М = к + п + 1 = К + п, К = М-п = к + 1, к= -1 ,М- 1.

Следствие 2. Частично М-сглаженная (М = к+п +1) реализация у^+1 = П^+1у = у(М) описывается следующими рекуррентными уравнениями:

ук+1 = ук - !к+1Лк+1 пк+1, где пк+1 =

— процесс обновления: ошибка прогноза Ук+1+п/к модели (9) на значение отсчета Ук+1+п.

Доказательство. Первое следствие вытекает из леммы 1 и следствия 2 теоремы 1. Второе — из структуры базиса А (3), составляющих его векторов Пк (4), формул (21) и уравнений (22).

Лемма 5. Пусть М = к + п + 1 = п, Ь, а к = -1, N — 1. Для вычисления в (23) прогнозов Ук+1+п/к7 к = О, N — 1, достаточно иметь решение (М — 1)-задачи фильтрации, т. е. иметь (М — 1)-оценку предыдущего состояния модели:

[у]к+1(п)/к = {уг/к }к=кп+1 (10).

Следствие 1. В уравнении (23) для продолжения рекурсии достаточно вычислять только п-вектор последних отсчетов {Уг/к+1}к+п+1 М-сглаженной реализации ук+1.

Следствие 2. Для прогнозирования и продолжения рекурсии достаточно вычислять в уравнениях (22) только п последних компонент [/]к = {/¿к}к+п и [/Ь = ПРЯМЫХ и возвратных векторов и fk, к = О, N, процесса

двусторонней ортогонализации.

Теорема 7. (а) Вычисление в уравнениях (22) п последних компонент [/]к = {/¿к}к+п и [/]к = {/¿к}к+п ортогонализирующих векторов /к и /к необходимо и достаточно для решения М-задач (1), (2) фильтрации и прогнозирования для всех М = К + п = п, Ь, К = к = О, N.

(Ь) Пусть ао = 0. Тогда вычисление п последних компонент векторов /к, /к в уравнениях (22) необходимо и достаточно и для решения задач сглаживания — полной и всех частичных, — а именно, для вычисления сглаженных реализаций уК = у(М) всех М = К + п-задач (1), (2) при М = пД. Если М = Ь, то получаем решение полной задачи сглаживания (1), (2).

Доказательство. Утверждение (а) представляется очевидным в свете леммы 5 и ее следствий. В обоснование утверждения (Ь) отметим, что если ао = 0, то РУ (9) разрешимо относительно младшего отсчета Ук. Тем самым полу-

п

чаем возможность решать это РУ в обратном времени: Ук/к = — Ук+г/каг/ао,

_ 1

к = N,0. Следовательно, могут быть вычислены все (К + п)-сглаженные реализации у/с = Ук = у(М), М = п,Ь, с граничными условиями {yi/k}k+l ^

Утверждение (а) теоремы 7 справедливо и для уравнений вариационной идентификации в реальном времени на основе теорем 5 и 6. Эти уравнения будут получены далее.

реализаций ] = 0, п, таких, что их невязки А* — матрица выборок V, т.е. (К-Щ,А} = V.

8. Уравнения для идентифицирующей матрицы

Лемма 6. (а) Пусть (Ь +1) х (п + 1)-матрица Кз- Щ — совокупность п + 1 j = О, п, таких, что их невязки А* 1

з 1о> А

(Ь) Пусть матрица Ш такая из матриц |"К Щ, что квадрат ее нормы \\ ШУ2 = 8р(Ш*Ш) = Бр(Ш, Ш} минимален.

Тогда идентифицирующая матрица Я (20) задачи ВИ (1), (2) есть Я = Ят = (Ш, Ш} — матрица Грама такой «минимальной»» матрицы реализаций Ш.

Доказательство. Утверждение леммы следует из того, что «минимальное» решение системы А*Ш = V есть Ш = Ша = А(-1)*У = АС-1 V [10, с. 37].

Лемма 7. Пусть с^ = (щ,Ак) = \с^,скк\* —последний столбец самосопряженной (к + I)-матрицы Ск = (Ак,Ак). Здесь ск = {щ,Ак^1), скк = {щ,^) = ||а||2. Пусть РкЬк = |Рк , 1|*Ьк — последний столбец обратной матрицы Скх, Ьк — диагональный правый нижний элемент, С--^ — (к + 1)-матрица С--1, окаймленная нулевыми последними строкой и столбцом, а Ак = С-1 — С-—1Л/о-Тогда для к = 0, N

Рк = Ак ск,

Р*ск = Ь-1,

А Р = / ,

Ьк = ак = \\/к\\-2.

Доказательство. Формулы (24.1), 24.2) следуют из формулы Фробени-уса (ФФ) обращения матриц, окаймляемых снизу и справа [10, с. 60]. Для матрицы С запишем ФФ в следующем виде:

С-1 = С--1/0 + Рк (Р*ск )-1Р* —^ Ак = А кск (с* АкСк )-1ск А к.

Здесь Е = |Ё*,1|* = Акск, Тк = -Ск_1/0ск. Формулы (24.3), (24.4) получаем из цепочки равенств:

АкРк = АкАкск = Ак(С-1 — С--1/0) (пк,А к} = РкПк — Рк-1Пк

= Пк — Рк-1 Пк = Пк_ 1Пк = /к — а-1 = (АкР, АкР}

= Р**СкРк = Ь-1|0 к, 1|Рк = Ь-1.

Обозначим через И1к = АкСк1Ук, к = 0,Ж, минимальные решения частичных М -систем уравнений А к Шк = Vk для М = к + п. Здесь С к = (А к, А к}, а

Ак = \щ\*,щ&Е, г = 0Д.

Теорема 8. Минимальные решения Шк +1 системы уравнений А к +1 Шк+1 = 14+1 как функции параметра к = — 1, N — 1 описываются следующими разпост-ными уравнениями:

Шк +1 = Шк + /к+1ак +1К+1 — Я*+1/к )>

-*

Здесь «к+1/к —прогноз системы Ак+1Шк = Ук+^к на к + 1 -ю строку г>к+1 матрицы Ук+1, — «к+1/к — ошибка этого прогноза, /к+1 = ЩПк+1 — завершающий вектор процесса прямой ортогонализации Грама — Шмидта для подсистемы

Ак+1, ак+1 = || /к+1 N—2.

Доказательство. Ввиду ФФ и формулы из леммы 6 для к = —1, N — 1 и = 0 имеем

Шк+1 = Ак+1( С,.,/1! + ^й+1 Ук+1 = Шк + /к+1ак+1ск+1Дк+1Ук+1.

Учтем, что Дк+1 = С-+1 — С—/0. Получим то, что требуется доказать: Ск + 1С-/0Ук +1 = (А к +1 +1 )С-/1оУк +1 = (Шк,Пк +1 )>

ск+1С-+1 Н+1 = |0т+1 > +1 = «к+1

-> Рк+1 ^к+1 = Ск+1Д к+ 1^к +1 = ^+1 — %+1/ к ,

где _

К+1/к = ~Рк+1Ук = (\Ук+1,г]к+1).

Следствие 1. Для идентифицирующей матрицы Грама Q = (Ж, Ж) имеют место следующие рекуррентные уравнения при значениях к = — I, N — 1 и нулевых начальных условиях Q-1 = 0:

Q к +1 = Qk + 9к +1йк+19&+1; где ^к+1 = «к+1 — )к +1/ к.

Доказательство. Два слагаемых в уравнении (26) ортогональны друг другу, поэтому

Q к +1 = (Шк + /к +1йк +1 Ч*к +1 ,Шк + /к+1йк +1 Ч*к +1) = (Шк,Шк) + дк +10,кЧ*к+1. П

Следствие 2. Пусть матрица А есть МСВ (3), (4). Тогда векторы )к+1/к прогнозов (Шк, Пк +1) в уравнениях (26) и (27) могут быть вычислены в режиме фильтрации, описанном в теореме 7, лемме 5 и ее следствии 2. Это означает, что вектор прогнозов )к+1/ к определяется лишь последними п строками матрицы реализаций Шк. Следовательно, для их вычисления необходимо и достаточно знать лишь п последних компонент ортогонализирующих векторов / и / .

Эти факты задачи идентификации (они анонсированы после доказательства теоремы 7) дополняют утверждения теоремы 7, леммы 5 и ее следствия 2. В этих утверждениях аналогичные факты касались задач фильтрации, прогнозирования в соответствии с общими уравнениями сглаживания (23).

Доказательство. Режим фильтрации, напомним, означает, что в уравнениях (22) теоремы 6 вычисляются только п указанных в теореме 7 последних компонент ортогонализирующих векторов /к и /к. Этого достаточно и для решения задачи фильтрации, и для вычисления прогнозов в соответствии с общими уравнениями сглаживания (23). Этих компонент достаточно и для вычисления указаных в следствии 2 прогнозов )к+1/к = (Шк,Пк +1).

Действительно, последняя формула — из второго уравнения в (26) — и структура вектора Пк +1 — она определена в формулах (3), (4) — показывают,

что для вычисления векторов прогнозов ) к+1/ к достаточно знать лишь последние п строк матрицы реализаций Шк. Как следует из первого уравнения в (26), для вычисления этих строк досточно знать лишь п последних компонент ор-тогонализирующих векторов / и / . Эти компоненты — они основа режима фильтрации — определены в теореме 7. Согласно этой теореме, а также лемме 5 и ее следствию 2 знание этих компонент необходимо и достаточно для решения задачи фильтрации в соответствии с общим уравнением сглаживания (23). Согласно доказываемому утверждению этих компонент необходимо и достаточно также и для решения задачи идентификации.

Уравнения (22), (26), (27) составляют систему нелинейных разностных уравнений для идентифицирующей матрицы (¿к = (3(М), М = к + п = Ьо, Ь (см. лемму 1) частичных М-задач (1), (2).

9. Об уравнениях идентификации в реальном времени

Так названы уравнения для оценок )(М) как функции показателей длины реализации — чисел М или К = М — п. Идентификация в реальном времени актуальна, в частности, для ряда задач автоматического управления с идентификатором (адаптивные и самонастраивающиеся системы управления).

В задачах АИ в (13) возможно получение точных уравнений для оценок (15). Они основаны на матричных уравнениях Риккати обращения матриц с аддитивными факторизованными приращениями [11, с. 20]. Система уравнений для вычисления матричной функции <3^+1 > к = Ьо — п, N — 1 (Ьо > 2п), если матрица QL0 (27) Ьо-задачи (1), (2) обратима (замечание 2), может быть записана в виде

^+1 = № к + 9к + 1ак+1^к +1) 1 = Qk 1 — ^Чк +1 (а 1 + ?к+^-1?к +1) ^/О+^АД

Используя эти общие уравнения в формуле (15) для обращения матриц

Ук+1Ук+1 = УкУк +^+1^+!,

получим известные уравнения рекуррентного метода наименьших квадратов [11, с. 279].

Точные уравнения для оценок в задачах ВИ и ОР получить не удается. Ниже получим уравнения для одного вида приближенных оценок коэффициентов а в задачах ВИ и ОР.

В задачах ВИ и ОР есть общие черты не только в постановке (13) (замечание 3). Для М-задач (1), (2) введем К + 1-векторы А„ и Аог множителей Лагранжа, функционалы

Лгм = Лг + А^гА* у = + А^гА*Уа, ,/огь = ,/ог + А^г уа, а для задачи ВИ еще и ((М +1) х (п + 1))-матрицу вида (4) скользящего вектора МЛ А„г. Заметим, что функционалы идентификации задач ВИ (19) и ОР (теорема 9) инвариантны к длине вектора коэффициентов а, поэтому их градиенты ортогональны этому вектору.

Для функционалов идентификации ./„¿м, <Л>гМ и их производных по вектору коэффициентов а можно получить следующие формулы, используя результаты работ [3, 5,12].

Теорема 9. (а) МЛ и функционалы идентификации М-задач ВИ и ОР имеют вид

Xvl = СХог = (1/\\а\\2^а, •Ля = а* А*Аа = а*QViа, Зот = а*аХ*г Хог = а*Яот а.

(Ь) Градиенты 3'от и 3.Vг — векторы производных функционалов 3ог и З.г по коэффициентам а — могут быть представлены также в сопоставимых видах:

¿'т =

З'ог = (Яог — 1п+1 • Х*тХог)а = (Яог — 1п+1 • Зот/\\а\\ )а.

Используем для минимизации функционалов З.г и 3от градиентные итерационные процедуры (ИП): а[з+1] = а[з] — Т-1^]. Здесь Т — положительно определенная матрица. Подставляя вместо Т матрицы = V*С_1V или Яот = V*V/\\а\\2, получим такие ИП для задач ВИ и ОР, которые сходятся к минимумам функционалов З.г 3от при определенных условиях:

Для ОР в (30(Ь)) получили ИП, применяемую для минимизации подобных функционалов: поиск собственного вектора (СВ), соответствующего минимальному характеристическому числу (МХЧ), с множителем Х*Х, ограничивающим рост длины вектора а при этих итерациях.

Итерации вида (30(а)), (30(Ь)) имеют, как показывает опыт, высокую скорость и широкую область сходимости. Для однозначной сходимости итераций в задаче ОР (30(Ь)) необходимо и достаточно выполнения двух условий: изолированность (не кратность) МХЧ и ненулевая проекция начального вектора а[о] на соответствующий СВ.

Теоретический анализ функционала (20) и итераций (30(а)) задачи ВИ весьма сложен. Методы такого анализа отсутствуют. Пока имеются лишь обнадеживающие экспериментальные результаты. Они показывают, что итерации (30а) задачи ВИ имеют еще большую скорость и более широкую область сходимости (при одинаковом относительном уровне ошибок в исходной реализации у, а именно он определяет остроту экстремума и скорость сходимости), чем итерации в задаче ОР. Это можно объяснить меньшим числом независимых переменных в задаче аппроксимации ВИ по сравнению с задачей ОР в (13) (следствие 1 леммы 3).

Пусть уровень ошибок невелик. Тогда получение приближения а(М) к точному решению а(М) с достаточной для практики точностью возможно за одну итерацию вида (30(а)):

а(М )= а [1](к + п), а = а/\\а\\,

где ащ{к + п) = <Э/71(о?[о]) • а: [о], к = Ко, Ж, Ко > п. Матрицы (З^1 вычисляются по формулам (28), благодаря разложению (27) матрицы Я.

Пусть ошибки в исходной реализации у отсутствуют или их уровень очень мал, т. е. равно нулю или очень мало МХЧ •7т1П = тта 3. Тогда скорость

сходимости ИП (30) очень высока, но матрицы Q вырожденны или плохо обусловлены. В этих случаях необходимо обращать матрицы Q + Те, где е > 0. Это означает сдвиг на е характеристических чисел матрицы Q или, иначе, прибавление к функционалам аппроксимации J (2) и идентификации J (20) регуляри-зующего слагаемого е||а||2.

ЛИТЕРАТУРА

1. Эйкхофф П. Основы идентификации систем управления. М.: Мир, 1975.

2. Линник Ю. В. Метод наименьших квадратов и основы теории обработки результатов наблюдений. М.: Физматгиз, 1961.

3. Егоршин А. О. Вычислительные замкнутые методы идентификации линейных объектов // Оптимальные и самонастраивающиеся системы (ред. В. М. Александров). Новосибирск: Изд-во ИАиЭ СО АН СССР, 1971. С. 40-53.

4. Aoki M., Yue P. C. On priory error estimates of some identification methods // IEEE Trans. Automat. Control. 1970. V. 15, N 5. P. 541-548.

5. Егоршин А. О. Идентификация стационарных моделей в унитарном пространстве // Автоматика и телемеханика. 2004. Т. 65, № 12. С. 29-48.

6. Егоршин А. О. Идентификации и дискретизация линейных дифференциальных уравнений с постоянными коэффициентами // Вестн. НГУ. Сер. Математика, механика, информатика. 2014. Т. 14, № 3. С. 29-42.

7. Егоршин А. О. Об одном способе оценки коэффициентов моделирующих коэффициентов для последовательностей // Сиб. журн. индустр. математики. 2001. Т. 3, № 2. С. 78-96.

8. Егоршин А. О. Об одной вариационной задаче сглаживания // Вестник Удмурт. ун-та. Математика. Механика. Компьютерные науки. 2011. № 4. С. 9-22.

9. Егоршин А. О. Об одной вариационной задаче кусочно-линейной динамической аппроксимации // Вестн. Удмурт. ун-та. Математика. Механика. Компьютерные науки. 2012. № 4. С. 30-45.

10. Гантмахер Ф. Р. Теория матриц. М.: Наука, 1966.

11. Эльясберг П.Е. Определение движения по результатам измерений. М.: URSS, 2011.

12. Егоршин А. О. Об отслеживании параметров экстремума в вариационной задаче идентификации // Вестн. НГУ. Сер. Математика, механика, информатика. 2011. Т. 11, № 3. С. 95-114.

Статья поступила 25 ноября 2014 г. Егоршин Алексей Олегович

Институт математики им. С. Л. Соболева СО РАН, пр. Академика Коптюга, 4, Новосибирск 630090 egorshin@math.nsc.ru