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

Бели погрешности модели (4.1) не случайны, а обладают лишь ограниченностью, т.е. стохастический подход неприменим, то алгоритмы (4.37), (4.38) сохраняют работоспособность при определенном выборе длин шагов 7*. Например, при выполнении (4.31) алгоритм (4.37) будет обеспечивать при достаточно большом числе наблюдений к неравенство \Ук — Окх(к)\ < А, если Д > Д^,, векторы х(к) ограничены ||х(А:)|| < const и коэффициенты шага выбраны в виде

		71 = ( РЙГ при I » > д .									(4.39)
				[	0	при	\ук - ек_хх{к)\		<	Д,
где 0		<	7	< 2 ^	д ^ .	Алгоритм		получил		известность под
названием				"полоска".		Свойства		подобных		и	более общих
алгоритмов				изучены		в рамках развитого В.А.					Якубовичем
и его учениками метода рекуррентных									целевых неравенств
[6,	95,	32,	28].
	Градиентные				адаптивные		алгоритмы		(4.38)		применяются

для идентификации статических и дискретных динамических систем. Что делать, если ММ системы непрерывна? Аналогом градиентных алгоритмов для непрерывных динамических систем являются так называемые алгоритмы скоростного градиента [6, 32, 28]. Поясним идею метода скоростного градиента для случал, когда ММ системы описывается уравнением состояния

i =	F(x,0,*)	(4.40)
и требуется достичь цели
Q{x{t)ye{0,0	— 0 при f — оо,	(4.41)

где Q(x,0 > 0 - некоторая гладкая целевая функция. Поскольку Q не зависит от 0, ее градиент равен нулю и градиентный алгоритм типа (4.38) неприменим. Однако вместо задачи уменьшения Q можно поставить вспомогательную задачу уменьшения величины

д<	VeQTF(x)e}t)
Q(x,0,O = ^f +

- скорости изменения Q вдоль траекторий (4.40). Градиентный алгоритм, построенный по целевой функции Q(x)6)t)) на-

6 Б. Р. Андриевский и др.	161

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

§ = " 7

(4-42)

где 7 > 0, или в конечной форме

в = -yVeQ.

(4.43)

Можно предположить, что если существует 0*, такое что Q^x.O^.t) убывает для всех x,t} то, начиная с некоторого момента времени, будет убывать и Q(x,t), что приведет к достижению исходной цели (4.41). Это действительно будет так при некоторых дополнительных предположениях [6, 63, 106].

Градиентные алгоритмы и алгоритмы скоростного градиента можно применять не только для оценивания параметров, но и для оценивания состояний динамических систем, для фильтрации, управления и т.д. Благодаря своей простоте они находят применение в бортовых вычислительных устройствах, при микропроцессорной реализации систем управления, контроля, диагностики в различных областях. В последние годы адаптивный подход применяется также в вычислительных процедурах [8].

4.4.Принципы выбора модели

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

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

162

ММ: у = Oq + 0\u. Если зависимость выхода от входа носит экстремальный характер, то берут квадратичную функцию у = 0О + + 02и2, а если есть основания думать, что зависимость выхода от входа имеет перегиб, то начинают с

кубической функции у = 0О + вхи + в2и2 + в3и3.

Аналогично действуют и в случае нескольких входов и выходов (см. п. 4.1 , а также [46], с. 39 - 45). Возможен различный выбор ММ в разных областях переменных. Например, если построение ММ выполняется с целью оптимизации, то вдали от экстремума можно ограничиться линейной моделью, а при приближении к экстремуму переходить на квадратичную. В любом случае предпочтительнее модели, в которые постоянные коэффициенты входят линейно (см. п. 4.1). В последние годы приобрели популярность модели, основанные на нейронных сетях [99, 100].

Если точность моделей с постоянными коэффициентами недостаточна, т.е. исследуемая система нестационарна, то в модель вводят зависимость коэффициентов от времени (дрейф) Дрейф может быть монотонным или периодическим, причем в большинстве случаев достаточно ограничиться простейшими моделями дрейфа - линейными 0; = а, + bit или гармоническими 0, = at cos(wtf) + biS\n(<Jit). Если неясно, как выбрать частоты и>{} или есть аргументы в пользу того, что параметры испытывают нерегулярные колебания, то стоит попробовать хаотические модели дрейфа.

Если возникает дилемма: "выбрать ММ детерминированную или стохастическую", то предпочтение следует отдать детерминированной ММ. И только если не удается обойтись без случайности, то вводят ее, причем сначала в наиболее простой форме, например приведенной к выходу: y = F(u}0)+

где <р - случайное возмущение. Точно так же переходить к нечеткой ММ следует только при неадекватности стохастических ММ (что в действительности бывает не так уж редко).

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

163

описываются моделями первого-второго порядков. Например, для описания монотонных (так называемых "апериодических" ) переходных процессов часто достаточно апериодического звена - линейной стационарной дифференциальной ММ первого порядка, имеющей передаточную функцию (см. п. 3.3.2) от входа к выходу W(p) = Tp^fl1 г д е ^ ~ коэффициент усиления статической ММ, Т - постоянная времени звена. 1

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

ющее передаточной функции W(p) =	1	р	$	+	%——		, где f
					р + 1
- показатель затухания колебаний (|£| <		1).			Если	форма ко-

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

порядка. Иначе говоря, "не р е ш а й	с л о ж н у ю з а д а ч у ,
не р е ш и в с н а ч а л а п р о с т у ю "	(см. п. 1.3).

Отметим, что если ММ системы есть передаточная функция в непрерывном времени, то степень ее числителя не должна превосходить степени знаменателя. Иначе ММ будет включать чистое дифференцирование и ее нельзя будет привести к уравнениям состояния. Кроме того, на входные сигналы конечной мощности такая ММ будет реагировать сигналами бесконечной мощности, что говорит о ее негрубости. Понятие грубости ММ было введено в конце 20- х годов А.А.Андроновым и означает, что малые изменения исходных данных (входных переменных параметров) приводят к малым изменениям результатов (выходов). Грубость - очень важное свойство модели, так как реальные исходные данные для ее построения и применения всегда содержат

1 Величина Т определяется тем свойством, что за время 3Т переходная составляющая процесса в апериодическом звене уменьшится до 5% от начального уровня.

164

ошибки. Если бы эти ошибки существенно сказывались на результатах, моделью просто нельзя было бы пользоваться. Резюмировать сказанное можно в виде принципа грубости: "без о ш и б к и нет м о д е л и , а п о э т о м у н е г р у б ы е м о д е л и п л о х и е " .

Перечисленные правила следует принимать не как законы, а как рекомендации. В мире моделей царствует плюрализм, и для достижения успеха нужно испытать несколько вариантов моделей. При этом самая полная модель необязательно самая точная, а самая точная необязательно самая хорошая. Здесь уместно перефразировать слова Э.Хемингуэя: 1 м о ж н о п р е н е б р е г а т ь ч е м у г о д н о , н у ж н о т о л ь к о т о ч н о з н а т ь , к а к э т о п о в л и я е т на р е з у л ь т а т .

4.5.Культура компьютерных вычислений

До сих пор мы говорили о постановках задач и об алгоритмах их решения. Однако успех решения задачи с применением компьютера определяется, очевидно, успехом каждого этапа цепочки "модель-алгоритм-программа". Это значит, что самая точная математическая модель, самый великолепный алгоритм решения задачи могут быть безнадежно испорчены при плохой программной реализации. Одной из возможных причин этого является несовпадение машинной арифметики с обычной из-за конечности разрядной сетки компьютера. Возникающие ошибки могут привести к большим неприятностям, если их не контролировать и не соблюдать элементарные правила организации вычислений. Правила эти неформальны и напоминают правила хорошего тона. Уровень их выполнения определяет уровень вычислительной культуры пользователя компьютера. Поясним на примерах основные из этих правил. Начнем с двух парадоксальных примеров. На первый взгляд, в них нарушаются законы классической математики.

П р и м е р 4.5.1. Выполняя вычисления на обычном быто-

вом микрокалькуляторе,2 легко убедиться, что 108+1—108 = 0, 1020+106-102° = 0, но в то же время 108-108+1 = 1, 1020-1020+106 =

1 Э.Хемингуэй давал следующий совет начинающим писателям: "Можно опускать что угодно, нужно только точно знать, что опускаешь".

2Действия выполняются слева направо.

165

= 106. Таким образом, в машинной арифметике	нарушают-
ся законы коммутативности и ассоциативности	действий. 1
Применимость основных законов элементарной	математики
ставится под сомнение.

Пример 4.5.2. Хорошо известно, что limn_>ooen			= е,	где
е„ =	+ ^	, е = 2.718281828... - основание натуральных
логарифмов.		Попробуем приближенно вычислить	этот	пре-

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

пе„

105	2.7182682
107	2.7182817

1.1-107 3.0041656

1.2• 107 3.3201163

2.0• 107 7.3890551

3.0 • 107

1

При вычислении по данному алгоритму с помощью пакета MATLAB 5 в среде PC Windows ошибка хоть и оказывается меньшей, но характер результата, как это видно из рис. 4.1, остается тем же.

	бе, %		n
140			s бе, %
120
100
80
60
40	£ ю
20		10"	10"

Рис. 4.1. Относительная ошибка вычисления е в обычном (а) и логарифмическом (6) масштабах.

Таким образом, в компьютере нарушается представление о пределе последовательности как о числе, к которому при-

1 Заметим, что при использовании пакета MATLAB 5 на вычислительной платформе PC Intel в среде Windows 95/98 указанный эффект наблюдается при действии с величинами порядка 1015.

166

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

О п р е д е л е н и е . Машинным эпсилоном называется наименьшее представимое в компьютере число е, удовлетво-

ряющее условию 1 + е > 1, т.е. € т				= min{e : 1 + е > 1}.
П	р а в и л о	1.	Величина	ет характеризует	наимень-
шую	относительную	погрешность		вычислений и зависит от
конкретного экземпляра			компьютера.
Величина ет может быть вычислена с помощью					простой

программы. Такая программа содержится в современных пакетах программ и производит настройку внутренних параметров пакета под данный компьютер. Например, в системе MATLAB, чтобы узнать величину ет Для данного компьютера, достаточно просто набрать команду eps). Лля версии MATLAB 5 PC Windows, £m « 2.22 • 10"16. В микрокалькуляторах и персональных компьютерах величина ет колеблется вблизи 10"7, в зависимости от конкретного компьютера

(при	обычной	точности	вычислений).	Очевидно, что если
€ т >	10""*, то	на данном	компьютере	нельзя гарантировать,

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

Лва следующих примера показывают, что эквивалентные математические формулы могут быть неэквивалентными по

точности вычислений.

П р и м е р 4.5.3. Пусть ет =

10~2

и требуется

решить

урав-

нение х2 + 9.9х - 1

= 0.

Очевидно, все

результаты

можно

округлять до двух значащих

цифр. Если

применять для ре-

шения уравнения х2

+рх

+ q =

0 формулу х1 2 =

^ ^

где

D2

= р2 - 4 т о

получим \[Ъ «

10, хх

« 0.05, х2 «

-9.95.

Однако верный ответ для хх будет хх

«

0.1.

Ошибка по-

лучена при вычитании близких чисел.

Если

же

вычислять

Xi по эквивалентной

формуле хх 2

= —2

^ ,—, то

получаем

о

р ± V D

x i ^

9 Q + Ю ~

0-1-

Правда, по этой формуле х2 будет вычи-

слено уже с двукратной

погрешностью.

167

Пример 4.5.4. Известно, что при статистических		рас-
четах оценку дисперсии случайной величины по	измерени-
2	г д е	®	=
ям можно вычислить по формуле <т =
2		Х
jfX^x®»» или по эквивалентной формуле <7 = ^		~~ У-

Попробуем оценить дисперсию выборки хх = 1234.1, х2 =

1234.2, х3 = 1234.3.

Вычисляя на

микрокалькуляторе,

име-

ем: х

=

1234.2,

х2 = 1.5232496, jEx?

= 15232496 и по пер-

вой формуле <72 = 0.

По второй формуле получаем

результат

<72 = |(0 . 1 2 +0 2 + 0.12) «

0.67 • 10"2, который

легко получить и

устно.

При

выборе

формулы

и порядка

вычи-

П р а в и л о

2.

слений

избегать

вычитания

близких

чисел

и деления

на

малые

величины.

Пример 4.5.5.

Найдем

решение двух

очень

похожих си-

стем из двух линейных уравнений:

Система 1:

Решение:

Г х + 5у

=

17,

Г х =

17,

\

1.5х + 7.501 у =25.5.

\ у = 0.

Система 2:

Решение:

Г х + 5у =

17,

Г х = 32,

\

1.5® + 7.499у

= 25.5.

\ у = - 3 .

Столь большие различия в ответах

возникли

из-за то-

го, что

матрицы

коэффициентов

систем

Л;,г

=

1,2

п л о х о

о б у с л о в л е н ы :

величины Д, = (letЛ; малы.

Действитель-

но, Дх =

0.001,

Д2

= -0.001 (Д = апа22

- ai2a2i).

П р а в и л о

3.

Избегать плохо

обусловленных

матриц.

Если

это не удается,

то использовать

для работы

с ними

спе-

циальные

методы

[112,

104, 20].

Пример 4.5.6.

Числа, представимые,

например,

в типич-

ном микрокалькуляторе, лежат в диапазоне 10~88 < |х| < 1086. При выходе результата действия за левую границу диапазона выдается сообщение "underflow" ("исчезновение порядка"), при выходе за правую границу - "overflow" ("переполнение"). 1 При переполнении обычно говорят, что плохи исход-

1 В пакете MATLAB при исчезновении порядка переменной присваивается нулевое значение, а при переполнении - значение Inf

("бесконечность").

168

ные данные, а при исчезновении порядка полагают результат равным нулю. Но не следует торопиться. Пусть, например, вычисляется величинах при а = Ю""30,Ь= Ю"6 0 ,с= Ю~40, d = Ю-50. Если выполнять действия в следующем порядке: х = а • b/c/d, то компьютер сообщит об исчезновении порядка; если вычислить х = \ jcjd • а • Ь, то получим переполнение. Если же вычислить, например, х = а/с - fc/d, то получим правильный ответ х = 1. Этот же ответ можно получить, если отмасштабировать переменные, например, умножив на Ю40.

П р а в и л о		4.	При переполнении				или исчезновении
порядка	следует	попытаться			изменить		последовательность
действий,	ввести	масштабные			множители		и т.д. При исчез-
новении порядка не всегда				следует		обнулять результат.
Пример 4.5.7.		Пусть		снова е т		= 10"2	и требуется вычи-
слить сумму
		5 =	100 + 0.1 +			... + 0.1
				2000 слагаемых

Если вести суммирование слева направо, то получим 5 = 100, поскольку с учетом округления до двух значащих цифр 100+ +0.1 = 100. Если вычислять справа налево, то после сложения тысячи слагаемых получим 100, дальнейшее прибавление по 0.1 ничего не изменит и результат окажется: 5 = 200, что, конечно, уже ближе к истине. Но правильный результат 5 = 300 можно получить, только если сложить 1000 чисел по 0.1, затем еще 1000 чисел по 0.1, а после этого сложить промежуточные суммы.

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

Аналогичное правило действует при перемножении большого числа сомножителей.

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

169

Пусть, например, требуется	вычислить	величину	S =
= YlkLi р" с точностью до Ю-3.	Если вести	вычисления до
тех пор, пока общии член ряда -К не станет меньше 10			, т.е.
к
до г = 32, то получим 5 = 1.610, в то время как на самом
деле S = -g2- = 1.650... Если же попытаться вычислить			таким
образом приближенную сумму ряда J^jfcLi р		т о погрешность

останется бесконечной, как бы мала ни становилась величина

р поскольку ряд расходится.
П р а в и л о			6. Нужно		помнить, что остановка				итераци-
онного процесса			хх, х2 ,...		по косвенному		критерию	(например,
по критерию		\хп	— xn-i|	< е или		|F^rn)|	< е - в задаче		решения
уравнения	F(x)		= 0, по критерию			| | ^ ( £ п)|| < s - в задаче опти-
мизации	f(x)	и т. д.)		не гарантирует			достижения		заданной
погрешности		\хп — limn^oo хп\ < е.
П р и м е р 4.5.9. [104]. Пусть требуется вычислить									интегра-
лы		^
Еп=		I	xnex'ldx)		п = 1,2,..., 10 при em =			10"7.
		Jо
Для величин Еп можно получить рекуррентное соотноше-
ние				Еп	= 1-пЕп-х.				(4.44)

Попытавшись вычислить интегралы рекуррентно, пользуясь (4.44), получим: Ех = 0.367879,..., Е7 = 0.110160, Es = 0.118720, Eq = —0.0684800, т.е. уже для Е9 получен бессмысленный результат. Причина в том, что начальная ошибка округления быстро накапливается: при вычислении Е2 она умножается на 2, затем наЗ, 4,... ,9. Попробуем переписать формулу (4.44) в другом виде:

Еп.г =

(4.45)

и будем вычислять в обратном порядке, начав с произвольно-

го начального приближения для Еп,			п > 10. Например, начав
с Е2о = 0, получим:
£,9 = 0.0500000,		Еи=	0.0627322,
Ец	= 0.0500000,	£ 1 3 =	0.0669477,
Е17	= 0.0527778,	Еи=	0.0717733,
Еи	= 0.0557190,	£ „ =	0.0773523,
Ехь	= 0.0590176,	£ 1 0 =	0.0838771,
		Е9=	0.0916123.

170