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

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

За последние годы появился целый ряд книг по математическому моделированию, написанных специалистами в различных областях знаний: прикладной математике [90, 13], по автоматизированному проектированию [72, 97], кибернетике и теории управления [68, 74]» Хотя отдельные вопросы трактуются в них по-разному, в целом они подтверждают, что математическое моделирование уже сформировалось как самостоятельная научная дисциплина.

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

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

1.Не решай сложную задачу, не решив простую (принцип простоты).

2.Без ошибки нет модели, а потому негрубые модели - плохие (принцип А.А. Андронова).

3.Можно пренебрегать чем угодно, нужно только точно знать, как это повлияет на результат (принцип Э.Хемингуэя).

201

4.Чем проще модель, тем реже она обманет (принцип надежности).

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

6.Цель расчетов - не числа, а понимание. Прежде чем решать задачу, подумай, что делать с ее решением (законы Р.Хемминга).

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

То, чего нету, умножь на два:

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

ввоздухе тающего без следа, словно кусочек льда.

Приложение А. Использование системы Scilab

А.1. Принципы построения системы

В конце 1990-х годов всем, интересующимся вычислительными аспектами и автоматизированным проектированием систем управления, стал доступен пакет Scilab, аналогичный по своим возможностям пакету MATLAB, но имеющий перед ним одно бесспорное преимущество: Scilab распространяется абсолютно свободно, т.е. бесплатно. Система Scilab разработана специалистами из французского Научно-исследова- тельского института по информатике и автоматике (INRIA) и может быть скопирована (вместе с документацией) из сети Интернет: http://www-rocq.inria.fr/scilab/.

Ниже приводятся основные сведения о системе и демонстрационные примеры конструкций языка и возможностей системы Scilab [128].

Scilab состоит из трех отдельных частей: интерпретатора, библиотек функций (процедуры Scilab) и библиотек подпрограмм на языках Фортран и С. Подпрограммы на Фортране и С, строго говоря, не входят в Scilab: они вызываются интерпретатором. Большая их часть доступна независимо, например из сетевой базы данных Netlib: http://www.netlib.org/. Некоторые подпрограммы слегка модифицированы для лучшей совместимости с интерпретатором Scilab. Основным достоинством Scilab, как и системы MATLAB, является возможность легко оперировать с числовыми матрицами: такие операции, как сложение, умножение, транспонирование, слияние, вырезание, выполняются непосредственно. Scilab также дает возможность манипулирования с более сложными объектами (такими как полиномы, полиномиальные или рациональные матрицы передаточных функций) при помощи обработки списков и типизированных списков. Это в свою очередь открывает пути естественного символического представления таких объектов, как линейные системы и графы. При этом синтаксис операций идентичен используемому при работе с постоянными векторами и матрицами и очень похож на синтаксис языка MATLAB.

Scilab предоставляет широкий спектр процедур анализа нелинейных систем: решение уравнений, численное инте-

203

грирование явных или неявных систем, методы оптимизации (в том числе недифференцируемой). В состав системы Scilab входят сотни математических функций; встроенные библиотеки по линейной алгебре (включая работу с разреженными матрицами, кронекеровским и шуровским представлениями); по системам управления (классической, линейноквадратичной и #°°-оптимизации, редукции моделей, идентификации и др.); пакет программ по линейным матричным неравенствам (LMI), по обработке сигналов и фильтрации и пакет анализа и оптимизации сетей Metanet.

Прилагаемый к пакету Scilab тулбокс Scicos предоставляет средства графического описания и моделирования взаимосвязанных дискретных и непрерывных систем (аналогично системе SIMULINK) [128].

Имеются развитые графические возможности (двумерная и трехмерная графика, анимация). Возможности символических вычислений реализуются через интерфейс с Maple. Наконец, имеется система параллельных вычислений Parallel Scilab.

Общая философия Scilab состоит в создании вычислительной среды, предоставляющей пользователю:

-гибкий набор типов данных с естественным и простым синтаксисом (например, поддерживается тип данных "строка символов", а функции считаются объектами и могут определяться в системе Scilab и передаваться как входные или выходные аргументы других функций);

-достаточный набор вычислительных процедур как осно-

ву для выполнения разнообразных вычислений;

-открытую среду программирования, в которую легко добавляются новые примитивы;

-поддержку разработки новых библиотек и тулбоксов, в частности добавление новых модулей на Фортране и С.

Приведем некоторые примеры конструкций и текстов на Scilab.

А.2. Основные конструкции языка

Скалярные объекты

Синтаксические правила определения скалярных объектов и действий с ними на Scilab, как и многие другие конструкции,

204

идентичны соответствующим правилам языка MATLAB. Рассмотрим некоторые примеры определения таких объектов.

а=1	-	вещественная константа;
1==1	-	логическое выражение;
' s t r i n g '	-	строка	символов;
z=poly(0, ' z ' )	- многочлен от z, имеющий один ну-
		левой	корень;
p=l+3*z+4.5*z~2	-	многочлен от z.
Для задания рациональных			функций (например, передаточ-

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

р = 1 + 3z + 4.5z~2 r=z/p

В результате выполнения этих операторов получим:

г =

z

l + 3z + 4.5z 2

Матрицы

Матричные объекты вводятся операторами вида

А= [а+1 2 3; 0	0	a t a n ( l ) ;	5 9 - 1 ] ;
Здесь определена ЗхЗ-матрица вещественных				констант.
Чтобы задать матрицу логических			констант размера 1x2 ис-
пользуем выражение
Ъ= [e/.t, т .
Матрицы строк задаются операторами вида
M c = [ ' t h i s ' , ' i s ' ;		' а ' , ' m a t r i x ' ]
(здесь определена 2х2-матрица строк).
Может быть введена и матрица			многочленов.	Выражение
Мр=[р, 1 - z; 1, z*p]
задает 2х2-матрицу многочленов вида
Мр =
! l + 3 z	+ 4.5z~2		1 - z	!
!	1	z + 3z~2 +4.5z~3		!

205

Аналогично матрица			рациональных			функций
F =
,	1 + 3z + 4J>z 2 ^					- 1		,
	- 2 + 4г - Зг 2 + г 3				2 - 2г + Г 2 ^
I	1			г + 3г^2 + 4.5г~\				,
	- 2 + 4г - Зг 2 + г 3 - 2 + 4г - Зг 2 + г 3
задается	оператором	F=Mp/poly([l+'/.i				1-JJi	l ] , ' z ' ) .
Можно задавать и разреженные					матрицы.		Выполнив опе-
ратор
Sp=sparse( [ 1 , 2 ; 4 , 5 ; 3 , 1 0 ] , [ 1 , 2 , 3 ] )
получим:
	sP =
	( 4, 10) sparse matrix
	(	1, 2)		1.
	( 3,		10) 3.
	( 4, 5)			2.

Разреженные матрицы могут	иметь	и другой тип элементов,
например - логический: S p	( l , 1 0 )	~ S p ( l , l ) .

Списки

Лля образования списка используется оператор list. Выполнив

L = l i s t ( a , - ( 1 : 5 ) , М р , [ ' t h i s ' , ' i s ' ; ' а ' , ' l i s t ' ] ) ,

получим список элементов разного типа:
L =
L(l)
1.
L(2)
! - 1. - 2. - 3. - 4. - 5. !
L(3)
!	1 + Зг + 4.5z 2		1 - г	!
!		1	г + 3z~2 + 4.5z~3	!
L(4)
!	this is	!
!	a list	!

206

Можно также создать типизированный список с помощью оператора tlist:

L t = t l i s t ( [ ' m y l i s t ' , ' c o l o r ' , ' p o s i t i o n ' , ' w e i g h t ' ] , ' b l u e ' , [ 0 , 1 ] , 1 0 ) >

Для извлечения элемента из списка используются операторы вида L t С c o l o r ' ) .

Динамические системы

Линейную динамическую систему ("ЛС-систему", см. 3.3.2) можно задать стандартными уравнениями состояния. Операторы:

	A=diag([2, 3, 4]); В=[1 0; 0						1;	0 0];
	С= [1 -1	0];		D=0*C*B;	х0=[0;		0;	0];
	Sl=syslin('c',A,B,C,D,x0)
	- приводят к
	S1 =
	Sl(l) (state-space system:)
	lss					2.	0.	0.
	Sl(2)		=	A matrix	=
						0.	3.	0.
	Sl(3)		=	В matrix	=	0.	0.	4.
		1 .		0 .
		0.		1.
		0.		0.
	Sl(4) =			С matrix	=
	! 1. - 1 . 0. !
	Sl(5) =			D matrix	=
	!	0.		0. !

Sl(6) = XO (initial state) =

0.!

0.!

0.!

Sl(7) = Time domain =

с

Структура SI с точки зрения языка является типизированным списком, доступ к элементам которого осуществляется операторами вида

207

S l ( ' A ' ) , S l ( ' C ' )

Для преобразования уравнений состояния к передаточной функции используется оператор ss2tf. Выполнив

Slt=ss2tf(S1)

получим:
Sit	=	- 1
,	1
•	- 2 + 5	- 3 + 5

Операторы Sit ('num'), Sit ( ' d e n ' ) позволяют извлечь матрицы числетелей и знаменателей передаточной функции.

Операторы
v=l:5;	W=v'*v	- умножение матрицы-константы;
W(l,:)		-	извлечение	первого столбца;
W(:,\$)		-	извлечение	последней строки;
Mp'*Mp+eye		-	матрица многочленов.

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

ans	=
column 1		^	^	^
!	3 + 6z + 18*		2 + 27г 3 + 20.25*	4 !
!		1 + Зг + 4.5z~2		!
column 2			2
!

!^ 1 + Зг + 4.5г

!2 — 2г + 2 + 6-Гз + 18г~4 + 27г~5 + 20.25г~6

Mpl=Mp(l ,l)+4.5*'/,i					-	комплексная		матрица;
	Fi=C*(z*eye-A)~(-l)*B;				-	вычисление		передаточной
						функции;
	F ( : , l ) * F i				- действия		с	рациональными
						функциями.
	В результате выполнения двух последних операторов по-
лучим:
ans						1	+ Зг + 4.5г~2
	1 + Зг + 4.5г~2
4	- Юг + Юг	2 -	5г	3 + г	4	6 - 14г + 13г		г ^ б г 3 + г 4
		1
4	- Юг + Юг	2 -	5г	3 + г 4		6 — 14г + 13г 2 — 6г 3 + г 4

208

М= [Мр -Мр; Мр'Мр+еуе] ; - сцепление многочленных матриц; [Fi, F i ( : ,1)] - сцепление рациональных функций; F=syslin( ' с ' , F ) ; - определение линейной системы через передаточную функцию;

Num=F( 'num');		Den=F('den') ; - извлечение числителя и зна-
менателя передаточной функции.
Некоторые числовые операции
inv(A)	-	обращение числовой матрицы;
inv(Mp)	-	обращение многочленной	матрицы;
inv(Sl*Sl')	-	соединение двух линейных систем и обраще-
		ние;
w=ss2tf (ans)	-	получение передаточной	функции.
Результатом выполнения последнего оператора будет:
w _	18 - 30418 .55^2 - 5д"з +		О . бА
		6 . 5 - 5 5 + 5 2

wl=inv(ss2tf (Sl)*ss2tf (SI*)) - произведение двух передаточных функций и обращение;

A=rand(3,3); B=rand(3,l); n=contr(A,B) - проверка управляемости системы третьего порядка, заданной матрицами уравнений состояния Л, В со случайными значениями элементов;

K-ppol(A,B,[-l-Xi -1+y.i - 1]);

poly(A-B*K, 'z') - poly([ - l - e / . i -1+УД - 1 3 , ' z ' ) - синтез модального регулятора по состоянию (см. [4, 6]) и проверка правильности вычислений.

Врезультате выполнения этих операторов получим:

К= - 113.28616 111.62667 33.092441

s=sin(0:0.1:5*e /,pi);
s s = f f t ( s ( l : 1 2 8 ) , - 1 ) ; - быстрое	преобразование Фурье
xbascO; plot2d3(H enn" , l , a b s ( s s ) ' ) ; - вывод графика
Определение функции в режиме интерпретации
d e f f ( ' [ x ] = f a c t ( n ) ' , ' if n=0 then	x=l,
else x=n*fact(n - l),end')
10+fact(5)

209

ы ScilabQrapbicl	|.|у|к! и SciUbOrvhicZ			I.
	в» го2ю	1Mb» х>ш
	\	<с	„ д »	\

томи IK

а)

б)

Рис. А.1. Примеры графических изображений пакета Scilab. а - трехмерный график; 6 - моделирование велосипеда.

Оптимизация

d e f f ( ' [ f , g , i n d ] = r o s e n b r o ( x , i n d ) ' , 'a=x(2)-x(l)~2, b=l - x(2),

f=100.*a~2 + b~2		,
g(l)= - 400 . *x(l)*a		, g(2)=200.*a -2.*b ' ) ;
[f ^jgDsoptimCrosenbro, [2;2] , ' q n ' )
Моделирование
a=rand(3,3);	e=expm(a);

d e f f ( ' [ y d o t ] = f ( t , у ) ' , ' y d o t = a * y ' ) ; e ( : , l ) - o d e ( [ l ; 0 ; 0 ] , 0 , l , f )

Описание систем

Введем скалярную переменную 5 : s = p o l y ( 0 , ' s ' ) ;

Построим матрицу h рациональных функций от s : h=[1/s,1/(s+1); ..

1/s/(s+1),1/(s+2)/(s+2)]

Определим систему с именем w, заданную передаточной функцией с матрицей Л, и получим описание этой системы в форме уравнений состояния:

w=tf2ss(h); ss2tf(w)

210