Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков)

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

Для оценки параметров непрерывных по времени динамических систем существуют два подхода. Первый состоит в переходе от непрерывной ММ к дискретной и оценивании параметров дискретной ММ. Способы дискретизации приведены в п. 3.3.3. Затем полученные оценки в параметров дискретной ММ пересчитываются в оценки исходной ММ путем континуализации (см. п. 3.3.3).

Второй подход состоит в сведении исходной ММ системы к некоторой вспомогательной статической ММ. Лля этого вводится функционал качества 1 £?(•), характеризующий погрешность приближения ММ к исходной системе, вычисляемую по траектории системы на заданном промежутке [О,Т]. Примером ттакого функционала является квадратичная ошиб-

ка Qi = fQ e2(t)dt, где = y(t) - г/(*), y{t) - реальный выход системы при входном воздействии ix(<)» £(0 ~ процесс в ММ. Поскольку ММ задается вектором параметров 0, т.е. y(t) = y(t} в) (начальные условия считаем фиксированными и одинаковыми в системе и в ММ), то функционал Q оказыва-

ется функцией конечномерного вектора параметров Q =	Q(e).
Теперь задачу оценивания можно поставить как поиск
в = argminQ(0),	(4.13)
6

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

[80, 22].

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

оналов качества:
-	интегральная абсолютная ошибка Q2 = /J
-	максимальная ошибка Q3 = sup 0 < K T \e(t)\\
-	терминальная ошибка на промежутке [0,7] Q4 = \s(T)\.
		j	*
Если ММ системы имеет вид уравнения состояния ^			=
= F(x, и,		0), где х G Rn - вектор состояния, то можно	сфор-
мировать		показатель качества из ошибки по состоянию. На-

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

151

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

	Qt=	J[О | И 0 - « ( 0 1 Г Л		(4.14)
или в более общей форме:
Q6=	[ (x{t)-x{t))R(x{t)-x{t))Л,			(4.15)
	Jo
т	0 - некоторая		симметричная	положительно
где R = R >
определенная матрица весовых			коэффициентов.

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

1. Значение функционала качества зависит от длины про-

межутка Г и от начальных условий х(0).		Вследствие	этого
величину	Т приходится брать достаточно большой, чтобы
за время	[О,Т] переходный процесс успел	закончиться.	Ми-

нимизация должна проводиться либо при одном, фиксированном значении х(0), либо при "наихудшем" по некоторому ограниченному множеству значении х(0) (т.е. при минимаксном критерии качества), либо по усредненному показателю Q(0) = MQ(6) (при заданном "априорном" распределении вектора х(0)).

2. Чтобы вычислить Q(6), требуется определить траекторию системы на всем промежутке [О, Т], т.е. оценка точности и коррекция ММ могут проводиться только после окончания переходного процесса в системе.

3. Зависимость Q(0) обычно весьма сложна (в частности, неквадратична), вследствие чего условие оптимально-

сти dQ/dO =	0 представляет собой нелинейное	относительно
в уравнение.	Это значит, что уравнение ММ	не сводится к

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

Последний недостаток особенно снижает привлекательность подхода, поскольку существуют другие варианты, позволяющие сводить задачу к линейной. Например, можно перейти к дискретной ММ (4.11) и задать показатель качества

152

в виде	N
	QV) = Y,(Vk-0Txk)\	(4.16)
	к—т

Избавиться от последнего недостатка можно также, перейдя от функционалов ошибки, характеризующих разницу в выходах (или в состояниях) исходной системы и модели, к функционалам невязки, зависящим от погрешности уравнения ММ. Например, заменив в (4.14) векторы состояния я(£), x(t) их производными, приходим к показателю

Q7=	[	| | i ( < ) - F ( z ( 0 , u ( 0 , 6 0 H 2 ^	(4.17)
	JО

который в случае линейности F(x(t), u(t)} в) по в является квадратичным по в. Надо сказать, что неквадратичность задачи не мешает применению целого ряда методов оптимизации [80]. По поводу оценивания параметров нелинейных систем см. также [27].

Продемонстрируем различные подходы на примере.

П р и м е р 4.1.6. Пусть структура ММ задана дифференци-

альным уравнением первого	порядка
у =	ау + Ьи,	(4.18)

параметры которого a, b подлежат оцениванию по результатам измерений ук = y{tk), ик = u{tk), при tk = kh, к = 0, 1,... , N.

Перечислим возможные способы формализации задачи оценивания параметров в (4.18).

	а) Лискретизируя (4.18) по методу Эйлера (см.			п. 2.2.3),
получаем дискретную ММ в виде:
			yk+i = {l + ah)yk + bhuk.	(4.19)
Лля		сведения	(4.18) к (4.1) положим 01 = 1-аЛ, e2 = bh}	у = ук+х,
х\	=	Ук > х2 =	ик. После расчета оценок параметров	(4.19) вХ}
02	необходимо пересчитать их в оценки параметров исходной
системы (4.18) по обращенным формулам
			а =	(4.20)

б) Заметим, что в формулах (4.20) присутствует деление на величину шага Л, которая может быть малой величиной.

153

Чтобы избежать потери точности, можно попробовать взять дискретную ММ в виде

	=аук	+	Ьчк.	(4.21)
Лля сведения к (4.1)	полагаем
У = Vk+lh Ук »	^1=^1» 02 =	6;	аг! =	х2 = ик.

Отметим также, что формы (4.19), (4.21) не изменяются, если

"продифференцировать" их по времени, перейдя	к прираще-
ниям Ук Ук — Ук-\)^к —• Uk — хсjb х- Этот прием	используется

для повышения точности оценок за счет увеличения разнообразия векторов х*.

в) Можно дискретизировать (4.18) по методу Тастина (см.

п. 3.3.3). Лля этого в передаточной функции системы (4.18)
L	С\ I	—1
W(p) = у _ а	делается подстановка р —• j r | ^	и результат
умножается	на ^ ^ z . Получим дискретную	передаточную

функцию
W ° ( z ) = (2 —	2 + ah'
соответствующую разностному	уравнению
=	+ А*'-	(4-22)

Заметим, что если вычесть из обеих частей (4.22) ук и разделить на Л, то получим модель

структура которой совпадает с (4.21). Отличие от варианта б) есть лишь в формулах пересчета оценок - прямых:

д	а	д	Ь
~	1 — ah/2'		" 1 - ah/2

и обратных:

а =	, 6 = ^ ( 1 + 0^/2) .

1+0 х Л/2

154

г) Проинтегрировав обе части (4.18) по отрезку [0,**], получим ММ (4.1) при у = укУо'

	rtk		Г*к
х\ =	Jо	y(t)dt;x 2	= / u(t)dt;6i = а;02 = Ь.
			Jо

Приближенно МОЖНО ВЗЯТЬ Хх = Л^*=12/лЯ2 = д) Наконец, можно поставить задачу оценивания параме-

тров (4.18) как минимизацию функции интегральной ошибки:

Q i M ) = /	Ы О - у ( М . ь ) ) 2 ^
Jо
где у(1,а,6) - решение уравнения ^ =		ау + Ьи с начальным
условием у(0, а, 6) = у(0), или интегральной невязки
т
Q*(a>b) = l	- Ы*) -	Ы*))3л

4.2.Регрессионный анализ и метод наименьших квадратов

После приведения ММ к стандартной форме (4.1) вычисляются оценки параметров ММ по результатам ряда наблюдений за входами и выходами системы. Пусть имеются результаты N наблюдений хц, x,n, yt, i = которые объединены в матрицу из N строк и п столбцов (N х n-матрицу) X = {^«j} и N-вектор У = col(3/1, ...,Удг).

Начинал с классической работы К. Гаусса (1809 г.), популярен способ вычисления оценок, называемый методом наименьших квадратов (МНК) и состоящий в выборе оценки |0|, минимизирующей выражение (4.16) - среднеквадратическую ошибку . МНК-оценка удовлетворяет системе так называе-

мых "нормальных	уравнений"
	ХТХв = XTY	(4.24)
и имеет вид	6 = CXTY,	(4.25)
где С = {XX)-1	- так называемая дисперсионная	матри-

ца. Для вычисления в не требуется никаких предположений

155

о свойствах измерений, однако чтобы сделать выводы о достоверности и точности ММ, такие предположения нужны. Часто принимается следующая гипотеза: если представить модель измерений в виде

У = xe + ip,

(4.26)

где <р = со1((^!, ...,<^п) ~ вектор погрешностей (невязок) ММ, то <Pi случайны, независимы, центрированы (M<pi = 0) и равноточны (имеют одинаковую дисперсию а2 = М<р2). Эти предположения традиционны для раздела теоретической и прикладной статистики, называемого регрессионным анализом (РА). При их выполнении МНК-оценки 0 обладают следующими свойствами:

1)МО = 0 (несмещенность);

2)cov(0) = а2 С, где cov(0) = М{в - 0)(0 - 0)т - ковариационная матрица оценок, в частности

М(0, - 6 { ) 2 = <тсц;

(4.27)

3) МНК-оценки эффективны: для любой несмещенной оценки 0« выполняется неравенство М||0<, - в\\2 > М\\9 - 0||2;

4) если обозначить через У = Хв N-вектор прогноза выходных величин У по модели (4.26), а через R = ||У — У||2 =

=YliLiiVi ~~ Уг)2 ~ среднеквадратическую ошибку прогноза

(остаточную сумму	квадратов), то MR = (N - г)а2,	где г -
ранг матрицы X X.	Как правило (если входы модели (4.1)
не связаны линейной зависимостью), г = п и, значит,		оценку

дисперсии погрешности модели (4.26) можно брать в виде

После построения ММ (4.26) немедленно возникает вопрос о

ееточности. Этот вопрос имеет две стороны:

а) нельзя ли без существенной потери точности заменить

модель (4.26) более простой, т.е. упростить?

б) является ли модель (4.26) адекватной исходной системе, т.е. не следует ли ее усложнить?

Ответ в обоих случаях требует дополнительных предположений. В первом варианте будем считать, что модель (4.26) адекватна и выполнены предположения РА. Тогда можно поставить задачу о проверке гипотезы 0, = 0 равенства нулю

156

некоторого коэффициента ММ. На языке математической статистики это задача проверки значимости коэффициента 0,. В силу свойств 2), 3) МНК-оценок и при дополнительном условии нормальности <р (которое не кажется слишком обременительным на фоне остальных допущений) для проверки гипотезы 0,- = 0 можно пользоваться следующей процедурой.

1. По заданному уровню надежности р и числу наблюдений N из таблиц распределения Стьюдента (см., например,

[1]) или с помощью программы tinv тулбокса				STATISTICS
системы	MATLAB выбирается	значение	а р ,	такое что
P{\t\ < ар}	= р, где t - случайная	величина,	распределенная

по закону Стьюдента с N степенями свободы. Обычно берут

р= 0.9,0.95,0.99.

2.Определяются границы доверительного интервала для параметра 0; при уровне надежности р:

	[6i ~	Ъ +	<храу/Е$.			(4.29)
3.	Если интервал (4.29)	содержит точку		0, =	0, т.е.	если
<	ckpijy/cii, то гипотеза	0 = 0	принимается,		т.е.	в ММ

можно положить 0, = 0. В противном случае коэффициент 0; с надежностью р считается значимым, т.е. гипотеза 0, = 0 отвергается.

Отметим, что при N > 20...30 вместо таблиц и программ распределения Стьюдента можно использовать таблицы или программы нормального распределения и брать, в частно-

сти, а0.9 =	1 -65, а0.95 =		1 96, а0.99 = 2.58. Впрочем, при меньших
N, как указывалось в п. 3.5, методы						математической стати-
стики вряд ли дадут надежные и нетривиальные выводы.
Лля ответа на второй вопрос требуется выполнение пред-
положений регрессионного				анализа		и возможность		проведе-
ния в каждой точке				нескольких независимых слу-
чайных измерений у ц , . ч т о б ы по ним построить								альтер-
нативную оценку дисперсии				ошибки
=	Tvirriy Е		Х >	-	*)'.	Й = i Е	УН •	(4-30)
		1=1 ; = 1				j'= l
Гипотеза	об	адекватности		ММ	(4.26) отвергается,			если ве-
личина F	=	- 2	называемая дробь Фишера) оказывается
		т-у (так

157

достаточно большой: F > Fp.	Порог Fp ищется	по задан-
ному уровню надежности р из таблиц распределения		Фишера
[1]. В тулбоксе статистических	вычислений STATISTICS па-

кета MATLAB имеется соответствующая функция flnv (см. п. D.6).

Возникает еще один практически важный вопрос: нельзя ли повысить точность и надежность результатов эксперимен-

та путем его планирования,	т.е. специального		выбора точек
х,п? Очевидно, для повышения точности			оценок нужно
стараться сделать дисперсионную матрицу С			"поменьше", а
матрицу X - "побольше".	О	том, как это сделать, написа-
но в книгах по планированию		регрессионных	экспериментов

(например, [1, 110]). Некоторые функции имеются в тулбоксе STATISTICS.

Ло сих пор мы придерживались мнения о стохастичности модели (4.26). Если же более оправдан детерминистский под-

ход и на возмущения в (4.26) накладывается	единственное
ограничение
Ы <	(4-31)

то точность оценивания, как, впрочем, и сами оценки, можно

определить из решения соответствующей минимаксной зада-
чи:		п
min	max	|г/, — >	0,
в	\<i<N	^	3	' Я

j = i

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

Во всех случаях выбор подхода, числа учитываемых параметров ММ, плана эксперимента является неформальным актом, основанным на опыте и интуиции исследователя. Здесь оказывается полезным следующий "принцип надежности": "чем проще модель, тем реже она обманет", или "чем лучше модель объясняет прошлое, тем хуже она прогнозирует будущее". В иной форме эту же мысль выразил Р. Акофф: "Степень понимания явления обратно пропорциональна числу переменных в его описании". Практически часто применяют следующее правило: число независимых наблюдений должно быть в 3-5 раз больше числа параметров модели. Наконец,

158

не следует забывать принцип "равнопрочности", который настойчиво пропагандировал знаменитый российский математик, механик и кораблестроитель А.Н. Крылов: "Точность результатов не может быть выше точности исходных данных; точности промежуточных вычислений должны быть согласованы".

4.3.Адаптивные модели и рекуррентные методы

Если наблюдений при построении ММ оказывается много и они поступают последовательно во времени, то часто оказывается удобным обрабатывать их по мере поступления. Иногда такой способ обработки единственно возможный, например, когда исходные данные из-за большого объема не помещаются в памяти компютера или когда свойства системы меняются во времени ("дрейфуют") и параметры ММ требуют постоянной коррекции. Модели, которые изменяются (подстраиваются, адаптируются) в процессе наблюдений за системой, называются адаптивными. Соответственно методы оценивания параметров таких моделей путем коррекции

по текущим наблюдениям называются адаптивными,				или ре-
куррентными [6, 59, 95].
Общий вид рекуррентного алгоритма следующий:
		0* = Ф(0*-1		(4 -3 2 )
где к = 1,2,... - номер шага наблюдения. Многие			алгоритмы
можно привести к рекуррентной форме.
Пример 4.3.1.		Пусть ММ системы имеет вид у, = 0 + y?i,
где	0 = Myi - оцениваемый параметр, ур, - помеха,			Му>{ = 0.
Обычная оценка среднего арифметического 0, =				У* мо~
жет быть преобразована к рекуррентной форме:
		вк = ek .x + l(yk -		(4.33)
(при 0О = 0). Обобщением (4.33) является метод			стохастиче-
ской	аппроксимации	[80, 58]:
		0* = 0*-i + lk(Vk - 0*-i),		(4.34)
где	> 0 - коэффициенты поправок (шаги) алгоритма, кото-
рые для обеспечения сходимости оценок должны			удовлетво-

159

рять условиям Роббинса-Монро:

оооо

i=l к=1

П р и м е р 4.3.2. Можно показать, что МНК-оценки параметров ММ (4.1) удовлетворяют рекуррентным соотношениям

Ok =	+	Ткх{к) (ук		- б1_1Х(к)),		(4.35)
Г, =	Г^	-	h f	m	^ ,	(4.36)
			1 + х(к)		Tk-\x(k)
где х(к) = c o l ( x j t i , х к п )		-	вектор	к-го	измерения	входов ММ,

Гь = Тк > 0 - симметричная, положительно определенная пхп -

матрица.	При этом Г* =	х (0а ? (, ')Т ) ~ дисперсионная
матрица,	и, значит, соотношения (4.35) соответствуют МНК

лишь при определенном выборе начальной матрицы Г0. Однако если к —• оо, то оценки, получаемые из (4.35), приближаются к МНК-оценкам при произвольных начальных условиях 0О, Г0. Практически часто берут Г0 = 7I, где I - единичная матрица, 7 > 0 - достаточно большое число.

Если (4.36) не учитывать, а брать Г* = 7*1? где 7* > О удовлетворяют условиям Роббинса-Монро, то получим обобщенный (многомерный) алгоритм стохастической аппроксимации:

0к = 0к.! + 7к{Ук - Ci х(к))х{к) .	(4.37)
Алгоритм (4.37), в свою очередь, является частным	случаем
градиентного алгоритма. 1
0k = 0k-i-lkV$Q{Ok-i,Vk,Xk)	(4-38)

и получается из (4.38) при выборе целевой функции текущей

ошибки Q(6,y,x)	в виде
	Q(o,v,x)_=1b(v-e ')•
1 Ч е р е з V e Q	обозначается градиент	функции векторного аргумента
в =	т.е. вектор из частных	производных

160