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

В этом фрагменте программы с помощью функции max определяется длина реализации Т, задается число точек N для вычисления спектров и находится период дискретности Т0. Поскольку заданное число точек не совпадает с полученным в результате моделирования, а также в виду неравномерности моментов времени t} далее производится линейная интерполяция процедурой interpl . Для этого вводится массив моментов отсчета t i с постоянным шагом ТО и формируются подлежащие последующей обработке массивы y l i , y2i: 1

ti=0:T0:(Nt-1)*T0;

y l i = i n t e r p l ( t , y ( : , 3 ) , t i ) ; y 2 i = i n t e r p l ( t , y ( : , 1 ) , t i ) ;

После этого выполняем ДПФ процедурой fft и находим энергетические спектры сигналов:

Y l = f f t ( y l i , N ) ;

Y2=fft(y2i,N);

Sl=abs(Yl).~2/N;

S2=abs(Y2).~2/N;

Далее вычисляем частоту Найквиста и>п (wn) и шаг по частоте 6ш (dw):

wn=pi/TO;

dw=2*wn/N;

В рассматриваемом примере задано Т = 5 с, N = 4096, откуда следует ~ 2570 с"1, 6ш и 1.26 с"1. Частота Найквиста значительно превышает диапазон частот исследуемых сигналов. Построим графики спектров процессов на входе и выходе фильтра в диапазоне частот [0,100] с"1.

отах=100;

w=0:dw:отах;

lw=length(w); figure subplot(211)

1 Можно использовать и процедуру i n t e r p l q . Тогда соответствующий оператор следует записать с транспонированием массива t i , а именно как

y l i e i n t e r p l ( t , у ( : , 3 ) , t i ' ) . Это связано с тем,	что процедура i n t e r p l q
"чувствительна" к совпадению числа столбцов	аргументов.

71

stem(w,Sl(1:lw),'к.'),grid t i t l e ( 'S_l(\omega) ') subplot(212) stem(w,S2(l:lw),, k.'), grid

xlabel('\omega,	1 / c ' )
title(, S_2(\omega)1 )
xlabel('\omega,	1 / c ' )

Результат выполнения программы показан на рис. 3.6.

Рис. 3.6. Спектры процессов на входе (Si) и выходе (52 ) ф и л ь т р а (3.10).

3.2.3.Исследование экологической системы

Выше приведена упрощенная математическая модель системы "хищник-жертва". Эта модель представляется уравнением Лотки-Вольтерра (3.7). Получим фазовый портрет системы (3.7) на плоскости (2/1,3/2) Для некоторых безразмерных значений параметров и переменных. 1

На рис 3.7 приведена 5-модель системы (3.7). Начальные условия интегрирующих блоков заданы переменными у1_0,

у2_0; параметры	модели (3.7) а,Ь, с, d заданы одноименны-
ми переменными.	Лля получения фазового портрета следу-

1 Аналитическое исследование задачи приведено в [52, 90].

72

			_	—1•	1
					— —
					S	— К Ю
Z		-4 ь > Г				у1
	*
		•	• —	>	-1	-
			»4.	~~	S	— К Ю
•	ж	• * [ c > J				У2
Рис.	3.7. S-модель системы			"хищник - жертва" .

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

Программа построения фазового портрета экологической системы

а=2; Ъ=1; с=1; d=5; у2.0=2;

Yl=[1.5:0.5:5];

f o r k=l:length(Yl) yl_0=Yl(k);

s i m ( ' s 2 . 2 ' )

p l o t ( y ( : , 1 ) , y ( : , 2 ) ) , h o l d on xlabel{'y _ 1'>,ylabel('y _ 2') axis( [0 12 0 i n f ] )

end hold off

В программе заданы фиксированное начальное значение 2/г(0) = 2 и массив Y1 начальных значений у\(0). Оператором for - end организован цикл по параметру к. Число итераций равно числу элементов массива Y1. После моделирования SIMULINK-программой s2_2 выводится фазовая кривая на плоскости (2/1,3/2)- Чтобы получить фазовый портрет, т.е. изобразить на одном графике множество фазовых кривых, после оператора plot выполняется оператор hold (фиксация) графика (устанавливается опция on). По окончании

73

цикла фиксация снимается (hold off). 1 Оператор axis служит для задания диапазона осей графика. Параметр inf задает автоматический выбор границы по соответствующему направлению. Результат моделирования представлен на рис. 3.8. Стрелочками помечено направление движения по фазовой траектории со временем. Довольно сложно (хотя и возможно) запрограммировать вывод этих стрелок, поэтому на данном графике они добавлены "вручную". Как видно из рисунка, в системе имеют место нелинейные колебания относительно некоторого состояния равновесия, амплитуда которых зависит от начальных условий. Часто эти кривые трактуются таким образом, что экологическая система имеет естественное состояние равновесия и уничтожение "хищников" (точка а на графике) приводит через некоторое время к резкому уменьшению числа "жертв" (точка 6). Если также учесть, что уменьшение численности популяции может привести к ее полному исчезновению, то следует сделать вывод о необходимости соблюдать осторожность при воздействии на экологические системы.

Рис. 3.8. Фазовый портрет системы "хищник - жертва" .

1 Заметим, что последний оператор имеет смысл только в том довольно редком случае, когда нужно выполнить вывод в то же графическое окно, стерев при этом находящееся в нем изображение, которое до этого "накапливалось". Обычно вывод происходит в следующее окно и оператор hold off не обязателен.

74

3.3.Модели состояния динамических систем

3.3.1.Модели общего вида

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

[6, 95]. 1	Например, для ЯС-цепочки переменная состояния
есть	поскольку значения E\(t) и входного воздействия

Eo(s) при s > t однозначно определяют (в силу (3.5)) значение Ex(s) при s = t. Для модели динамики популяций (3.7) состоянием является вектор х = со\(yi,y2 ).

В общем случае уравнения состояний - это системы дифференциальных или разностных уравнений первого порядка вместе с уравнениями для выходных величин. Начальное состояние представляет "память" системы о прошлом. Модель состояния непрерывной динамической системы записывается в виде [6, 95, 62]

(3.15)

Vi = Gi{xi,...,xn,ul,...,um,t)

У1 = Gi(xu...,xn,ui,...,um,t)

где иХ)...,ит - входные переменные; y\}...}yi - выходные переменные; хь ...,хп - переменные состояния. Вводя векторные обозначения, можно записать (3.15) в более компактном виде:

х = F{x,u,t),	y = G{x,u,t),	(3.16)
где х = со1(хх,х п ); и = col(iib ...,um); у = соl(j/i,		Для

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

1 Здесь и далее вектор (столбец) с компонентами ii,i2, ... ,i n обозна-

чается X = со\(хх,хп).

75

Рис. 3.9.

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

Ф(х,х,М) = 0,

и,0 = 0,

(3.17)

частным случаем которых являются неявные системы

0 = 0 .

(3.18)

3.3.2.Линейные модели и линеаризация

Часто вместо (3.15) используют упрощенные ММ, основанные на том, что процессы в системе протекают, мало отклоняясь от некоторой так называемой опорной траектории {x(f),£(f),y(f)}, удовлетворяющей уравнениям

	х = F(x,tz),			у = G(x,u).						(3.19)
Тогда можно записать приближенную линеаризованную										мо-
дель в отклонениях от этого режима:
	х = A{t)x		+	5(0^,
	y = C{t)x		+	D{t)u,						(3.20)
где х = х — х, у = у — у, й = и — й,
- d F u l
C(t) = f £ (x(«), € (0),				D(t)		=	fg(S(0,«(<))•
Пример 3.3.1.	x	= -x2 u;	x,u		6	Д1.		Линеаризуем		вбли-
зи траектории, соответствующей					u(t)		=	1.	Имеем х — —х2,
откуда либо x(t)	=	0 (при х(0) =			0),		либо		х(*) = 1 /(t	- а).
Рассмотрим второй		случай:
= ш ( - л ) „ 1				= - 2 х ( 0 = - г ^ ,

76

5 ( 0 = Ш		=	= " ( Г Г ^ -
В отклонениях х	= х — 1 /(t	— а)}й	= и — 1 линеаризованное
уравнение имеет	вид

Если расчетный режим является установившимся, т.е. не зависит от времени, то коэффициенты в (3.20) также не зависят от времени: A(t) = Л, B(t) = В и т.д. Такие системы называются стационарными. Особенно часто на практике встречаются стационарные линейные непрерывные системы, описываемые более простыми уравнениями

х = Ах + Ви, у = Сх.

(3.21)

Матрицы А, 5, С являются параметрами модели (3.21). Если линеаризация приводит к большим погрешностям, то стараются, по возможности, выбрать математическую модель, линейную по параметрам

х = Аф(х)и))
где Л - матрица параметров порядка nxN; V>(v) G	- нели-

нейная вектор-функция. К этому классу относятся, в частности, билинейные объекты, например х = ахх + а2хи + а3и, где

А = [а!,а2,а3];	= col(x,хи,и).
Сказанное выше относится и к уравнениям дискретных по
времени систем.	Уравнения	дискретной системы	в общем
случае имеют вид
Xfc+i = F(xkluk)}		ук = G{xk}uk).	(3.22)

Дискретным аналогом уравнений линейной стационарной системы (3.21) являются уравнения

Zjt+i = Рхк + Quik, ук = Rxk.

(3.23)

Наряду с уравнениями состояния широкое применение находят также модели в переменных "вход-выход" и модели, описываемые передаточными функциями. Для непрерывного времени уравнение "вход - выход" имеет вид

A(p)y(t) = B(p)u(t),

(3.24)

77

где р = d/dt - символ дифференцирования по времени; Л(А) =

=	+ an-iAn"1 + ... + а0;	5(A) = Ьт Хт +	..-.+ 60,	при-
чем в	(3.24) всегда т < п.	Дробно-рациональная функция
W(А)	= Б(А)/Л(А) называется передаточной функцией			систе-
мы (3.24), а полином - ее характеристическим			полиномом.1
Если уравнение (3.24) получено из (3.21), то
	W(X) = C(AIn - А)'1 В}			(3.25)
	А(А) =	det(AIn - А).		(3.26)

Они справедливы и в случае, когда вход и выход системы (3.21) являются векторами, при этом W(А) = Л"1(А)В(А) - матрица. Пользуясь (3.25), можно показать, что замена переменных состояния в (3.21) по формуле х' = Тх , где Т - неособая п х n-матрица (detT = 0), не приводит к изменению передаточной функции (3.26). Это значит, что обратный переход от описания "вход - выход" к уравнениям состояния (3.21) неоднозначен: при сохранении передаточной функции базис в пространстве состояний можно выбирать по-разному. На практике применяются несколько типовых способов перехода от передаточной функции к уравнениям состояния. Эти

способы соответствуют так называемым каноническим					пред-
ставлениям системы [6, 95].		Опишем	один из них, приводя-
щий к управляемому	каноническому		представлению.		Вместо
(3.23) вводятся два	уравнения:
	А(р) п = и,				(3.27)
	у = В(р) 7/,				(3.28)
где 7} - вспомогательная переменная. Очевидно, что					переда-
точные функции (3.24) и (3.27), (3.28)			совпадают. В качестве
вектора состояния в (3.23) берется			х = со1(7у,т),...
так что Х{ = rfx~~x\	Из (3.27)	и соотношений		х; = rj(i) = х,-+1|
г = 1,2,... ,п - 1, выводится		форма матрицы А и вектора В в
(3.21), а из (3.28), записанного в виде у = b0Xi + biX2-\					\-Ьтхт,
1 Л - комплексное число, аргумент полиномов А(Х)}				В(\).

78

получаем строку С:

0

1

0

. .

0

0

' 0 '

					, В =
0	0	1 . .	0	0	0
•	•	•	•	•
0	0	0 . .	0	1	0
—а0	— a i		• - f l n - 2		1
		С	= [Ь0.. .Ьт	0... 0].

Если для системы (3.24) наблюдению доступна производная уМ от величины у при г < п - т - то она может быть получена, если в найденных уравнениях сохранить Л, В в форме

(3.29) и взять С = [0,... ,60 ,... ,6т,... ,0].
Если в (3.24) 7П = п (такие передаточные функции			называ-
ются несобственными),	то систему (3.24) нельзя привести к
виду (3.21), но можно привести к виду
х = Ах + Ви,		y = Cx + Du,	(3.30)
где Л, В имеют вид	(3.29); С	= [Ь0 - во^щ • • • >br»-i	~ an-i^n];
D = 6П .

3.3.3.Дискретизация и континуализация

Дискретизация. Если исходное описание линейной системы непрерывно, можно перейти к дискретному описанию с помощью следующей процедуры.

Пусть состояние x(t)		системы (3.21) доступно измерению в
дискретные	моменты	времени t^ = kh,			к = 0,1,...,	где	h >
> 0 - шаг дискретности. Пусть u(t)				постоянно на		промежут-
ках между	моментами	коррекции f*.			Тогда динамику		век-
торов Xk =	x(tk) можно описать разностными уравнениями
(3.23), в которых матрицы			Р и Q определяются соотношени-
ями	Р = е д \		Q = A'l{P-In)B.				(3.31)
Здесь eAh -	экспоненциал матрицы Л, определяемый форму-
лой
	еАА = In + Ah + \A2h2...			=	£		(3.32)
					Jb=0

Если предположение о кусочном постоянстве u(t) не выполняется, то переход от (3.21) к (3.23) является приближенным, но

79

его точность растет по мере уменьшения шага h , если скорость изменения входа (величина u(t)) ограниченная. При достаточно малых h для вычисления еАН можно удерживать лишь первые несколько членов ряда (3.32) или аппроксимировать сумму (3.32) каким-либо способом.

Например, при переходе от (3.21) к (3.23) можно пользоваться формулой eAh « In + Ah , соответствующей численному интегрированию (3.21) методом Эйлера. При такой аппроксимации передаточные функции дискретной и непрерывной систем будут связаны соотношением

И^д(А) = Wn

,

(3.33)

т.е. при переходе к дискретному времени в передаточной функции W(p) системы (3.24) нужно заменить р на (1 - z)/h. Если матрица А - гурвицева, т.е. ЯеАДЛ) < 0, то метод Эйлера дает устойчивую аппроксимацию лишь при

h < mm (2КеАДЛ))/|А,(Л)|2,

(3.34)

где А,(А) - собственные числа матрицы А (корни полинома А(А)). Целый ряд способов перехода от (3.21) к (3.24) основан на аппроксимации матрицы еАН матричными дробями Паде (дробями, "числителем" и "знаменателем" которых являются матричные многочлены). Частными случаями этих способов является метод Тастина (формула Паде порядка (1,1)):

	( i - ^ ) " 1 ,		(3.35)
приводящий к соотношению между передаточными			функция-
ми
= гЬ^-(Ятд^) •			(336)
а также метод Дэвисона	(формула Паде порядка (2,2)):
е " « (I + 4А +	(I - 4А +	.	(3.37)

Отметим, что формулы (3.35) и (3.37) дают устойчивые аппроксимации при h > О (разумеется, если А - гурвицева).

Заметим, что формула (3.31) для вычисления матрицы Q применима, если det А ф 0. Трудностей, связанных с вычислением Q при вырожденной матрице Л, можно избежать, если

80