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

а вызов latex (tl) приводит к строке символов, которая после обработки текстовым редактором IATjtX[57] дает выраже-

Теперь перейдем к решению уравнения (1.4) при заданных начальных условиях и найденном через (1.13) управлении. Начнем с символьных вычислений. Для этого продолжим программу следующими операторами:

vy=dsolve('m*Dv=b','Dy=v','у(0)=0','у(0)=0'); v=simple(vy.v)

y=simple(vy.у)

v _ t l = s u b s ( v , ' t ' , t _ l ( l ) ) y _ t l = s u b s ( y , ' t ' , t _ l ( l ) )

В этом фрагменте первая строчка содержит обращение к процедуре символьного (аналитического) решения уравнения (1.4), приведенного предварительно к системе уравнений первого порядка относительно переменных v(t) и s(<), имеющей вид

mv{t) =и{ 0,

t/(0 = v{t).

Символом D в программе обозначается операция дифференцирования. Переменные v, у содержат найденные при = 6 решения. Оператором subs выполняется подстановка в них значения t = t', найденного выше. В формате IATj?X[57] получаются выражения

v

Продолжим вычисления. Рассмотрим следующий интервал времени t 6 [t'}T]. На нем u(t) = —а. Проинтегрируем нашу систему на этом участке. Для этого выполним операторы:

vy2=dsolve(,m*Dv=-a', 'Dy=v', 'v(tl)=vl' , 'y(tl)=yl') v2=simple(vy2.v)

y2=simple(vy2.y)

vT=simple(subs(v2,{'vl','t1','t'>,{v _ tl,t . l(1),T(1)>))

Оператором dsolve находится решение (1.4), для которого в момент t = выполнено v(ti) = vX) y(t 1) = Значения tX)Vi, yx подставляются оператором subs по результатам вычислений на предыдущем участке. Получим

41

v2 =vl - a*(t - tl)/m

y2 =yl+vl*t-vl*tl-l/2*a*(t~2+tl~2-2*t*tl)/m vT =0

Равенство нулю значения скорости на правом конце интервала соответствует условию поставленной задачи.

2.2.2.Применение процедуры численного интегрирования дифференциальных уравнений

Обратимся теперь к численному решению уравнения (1.4). Сначала приведем программу моделирования для более общей ситуации, а именно не будем учитывать, что рассматривается линейное уравнение с постоянными параметрами. Воспользуемся одной из процедур интегрирования обыкновенных дифференциальных уравнений, входящих в систему MATLAB, а именно процедурой ode45, в которой используется метод Рунге-Кутта 4-5-го порядков [74].

Лля интегрирования следует написать головную (вызывающую) программу, содержащую обращение к процедуре ode45 и подпрограмму-функцию, в которой вычисляются правые части системы дифференциальных уравнений. В головной программе должны быть заданы начальные условия для всех интегрируемых переменных, начальное и конечное значения аргумента интегрирования (в рассматриваемых примерах - это время t) и содержаться операторы отображения результатов моделирования. 1 Приведем тексты головной программы и подпрограммы вычисления правых частей.

Головная программа

clear close a l l

global a b m t l a=2e3;

b=5e3;

1 В пакете MATLAB имеется несколько процедур численного интегрирования дифференциальных уравнений. Для каждой из них может быть установлено значительное число параметров (опций). Более подробную информацию можно получить из документации, или выполнив команду help ode45. В нашем примере используются значения параметров, принятые по умолчанию.

42

m=le3;

у.Т=10;

tl=l/(b+a)*2~(l/2)#(m*(b+a)*y_T#a*b)л (1/2)/Ь T=y_T*m*(b+a)*2~(1/2)/(m*(b+a)*y_T*a*b)~(1/2)

x0=zeros(2,1);
[t,x]=ode45('fpl _ 2',[0,			T],x0);
plot(t,x),grid
set(gca,'FontSize',14)
xlabelC't,		с ' )
t i t l e O	v, м/с; у, м		')
l e g e n d ( 1 v ' , ' y ' , 0 )
Подпрограмма		вычисления	правых частей уравнений
function	dx=fpl_2(t,х);

global a b m t l dx=zeros(size(x)); if t<tl

u=b; else

u=-a; end

dx(l)=u/m;

dx(2)=x(l);

Рассмотрим эти программы подробнее. Прежде всего заметим, что если имя rn-файла головной программы безразлично (более того, ее операторы могут быть введены непо-

средственно в Окно управления	MATLAB), то подпрограм-
ма вычисления правых частей	оформляется как т-функция

[35, 74, 84] и помещается в рабочем директории (папке) под определенным именем. В рассматриваемом примере текст этой программы содержится в файле fpl_2.m. Имя программы (совпадающее с именем файла), записано в ее заголовке

function dx=fpl . 2(t,х);

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

Ct,x]=ode45( , fpl . 2 , ,[0, Т],х0);

43

Головная программа начинается операторами clear и close. Эти операторы необязательны, но с них рекомендуется начинать головную программу. Первый из них освобождает рабочую область от всех переменных, которые могут там оказаться после предыдущих вычислений.1 Оператор close закрывает графические окна, если они были открыты перед выполнением программы. Без этого каждый "запуск" программы будет сопровождаться увеличением их числа. В последующих примерах операторы clear и close указываться не будут. В операторе global перечисляются глобальные переменные, значения которых доступны как в головной, так и в вызываемой программах без передачи через параметры процедуры. Такой же оператор имеется и в подпрограмме fpl_2. В следующих операторах (а=2еЗ; Ъ=5еЗ; т=1еЗ; у_Т=10;)

вводятся исходные данные для решаемой задачи.		Далее вы-
числяются значения переменных	t l , Т. 2 Затем	указыва-
ются нулевые начальные условия	на переменные	состояния

(в этом примере вектор переменных состояния х = [v,y] ). Далее производится обращение к процедуре ode45, в котором указаны входные параметры: имя процедуры-функции вычисления правых частей уравнений; вектор [0,Т], содержащий начальный и конечный моменты времени, для которых нужно проинтегрировать уравнение; начальное значение состояния системы. Выходными величинами являются массив моментов времени t, для которых получено решение, и массив результатов интегрирования х. Массив t является вектор-столбцом, число элементов nt которого определяется автоматически процедурой интегрирования. Значение nt можно получить, выполняя по окончании вычислений опе-

раторы l e n g t h ( t ) или s i z e ( t ) . В данном примере nt = 73. Массив результатов имеет размер ntxn, где п - порядок ин-

тегрируемой системы (в нашем случае п = 2). Каждая строка массива х относится к соответствующему моменту вре-

мени.	Чтобы получить решение в заданные моменты		вре-
мени	Jo» • • • )	следует при обращении к программе за-
дать	вектор	•••»*fin] ( к а к э т о можно сделать, будет	рас-

1 Естественно, если программист хочет в дальнейшем использовать эти переменные, то оператор clear выполнять нельзя.

2 Заметим, что в данном примере выражения для этих переменных получены копированием результатов символьных вычислений, описанных выше.

44

смотрено ниже). Далее в головной программе выполняется вывод результатов на график оператором p l o t ( t , x ) . В той же строке через запятую размещен оператор grid, задающий

нанесение	сетки. Этот	и следующие операторы носят не-
обязательный характер.		Оператор		setCgca,'FontSize',14)
выполняет установку заданных свойств объекта axes,							опре-
деляющего	вывод осей	(axis)	на график.		1	Операторы
x l a b e l O t ,	с ' ) и t i t l e d	v,	м/с;	у, и ' )	служат для		выво-

да заголовка графика и надписи по оси аргументов (х-ось). Оператор legend( ' v ' , 'у' ,0) служит для вывода легенды - связи имен переменных с линиями на графике. Цифра "О" в последнем аргументе оператора legend означает, что для меток должно автоматически выбираться "наиболее удобное" (свободное от графика) место. 2

Результат работы программы показан на рис. 2.3, а. Время интегрирования на компьютере с процессором Pentium 166 ММХ составляет 0.15-г 0.17 с. Рассмотрим теперь программу

v. м/с; у. и

Рис. 2.3. Результаты интегрирования и ошибки вычислений.

fpl_2 (стр. 43). Состав формальных параметров этой подпрограммы фиксирован и определяется требованиями процедуры ode45: это входные значения момента времени t и переменных состояния х и выходной массив вектора правых ча-

1Перечень свойств объекта и их значений можно получить оператором set (в нашем случае - вызовом set(gca)) [35, 74, 84]. Отметим, что системе MATLAB 5.3 имеется удобный интерактивный редактор графиков, использование которого позволяет упростить текст программы.

2При выводе на терминал графики имеют разные цвета. Это отличие пропадает при черно-белой печати.

45

стей дифференциальных уравнений dx. Процедуре ode45 требуется, чтобы выходной массив был вектор-столбцом (нужной размерности). Для того чтобы выполнить это требование "автоматически", рекомендуется в начале программы выполнять оператор dx=zeros(size(x)) , который строит массив dx нулевых элементов, размер которого совпадает с размером формального входного массива х (в нашем примере size(x)=(2,l)). Процедура содержит также оператор global для передачи значений между головной программой и подпрограммой, минуя аппарат формальных параметров. В системе MATLAB 5 требуется, чтобы глобальные переменные были указаны до присваивания им значений (по тексту программы). В нашем примере это присваивание происходит в головной программе ниже строки, содержащей global. Далее в подпрограмме вычисляется входное (управляющее) воздействие u(t) как программная функция времени. Для этого используется конструкция if - else - end. Последними операторами программы являются операторы вычисления правых частей интегрируемой системы уравнений и записи результата в выходной массив dx.

Как видно из рис. 2.3, а, полученный в результате интегрирования процесс содержит ошибку (значение у(Т) не равно заданному значению ут = 10). Так как для данной задачи поучены аналитические выражения, можно исследовать эту ошибку более подробно. Получим значения y(t) по точной формуле и найдем относительную ошибку в долях ут :

_ Ш M l } где у - результат интегрирования. Расчеты выполним программой, которая запускается после процесса интегрирования, когда соответствующие переменные имеются в рабочей области MATLAB:

у1 =l/2*t.~2*b/m;

y . t l =1/(b+a)*a*y_T

v_tl=l/(b+a)*2~(l/2)*(m*(b+a)*a*y_T*b)~(l/2)/m y 2 = y _ t l + v _ t l * ( t - t l ) - l / 2 * a * ( t - t l ) . ~ 2 / m ; y=[yl(find(t<=tl)); y2(find(t>tl))]; dy=abs(x(:,2)-y)/y_T;

figure

subplot(211),plot(t,y - x(:,2)),grid title('\Delta y=y(t) - y~~(t)') subplot(212),plot(t,dy*100,'.k').grid

46

x l a b e l O t , с ' )

t i t l e ( ' \ d e l t a у, У.')

Первые четыре оператора программы служат для аналитического вычисления y(t) для моментов времени t) полученных выше при численном интегрировании. Весь процесс разбит на два участка: u(t) = Ь (при t G [0,£']) и u(t) = —а (при t 6 (£',*»]). На каждом участке движение определяется сво-

	bt2
ей формулой: для первого участка y(t) =		для второго -
y(i) = y(t') + v{t'){t - t') -	(t - t')2. Этим	формулам соот-

ветствуют первая и четвертая строки программы. Заметим, что, так как аргумент "время" t является массивом, для поэлементного возведения в степень используется операция .

а не ~ . Значения y(t') и v{t') обозначены соответственно через y . t l и v_tl. Они вычисляются во второй и третьей строках программы. Далее формируется массив у, содержащий искомые значения y(t). Для этого можно использовать оператор цикла for - end в сочетании с условным оператором if

- else -	end. Эта же цель достигается конструкцией
y=[yl(find(t<=tl));		y 2 ( f i n d ( t > t l ) ) ] ;
Смысл	ее состоит в том,	что массив у получается	объеди-
нением	(конкатенацией) частей массивов yl, у2. Эти части
(подмассивы) выделяются оператором And, который			форми-

рует вектор индексов элементов массивов по указанному в качестве аргумента логическому выражению. Этим логическим выражением здесь служат операции сравнения t<=tl для массива yl и t > t l - для массива у2.

Далее в программе вычисляется относительная ошибка 6у, которая помещается в массиве dy. Конструкция х( : ,2) означает "вырезку" массива результатов х по второму измерению, т.е. выделение значений перемещения y(t), полученных интегрированием. Далее следуют операторы вывода на график. С помощью операторов subplot вывод производится в соответствующую часть графического экрана [74]. Оператор title служит для вывода заголовка графика. Результат выполнения программы показан на рис. 2.3, 5, е. Абсолютная ошибка &У = 2/(0 - 2/(0 показана на рис. 2.3, б. Рассмотрим график относительной ошибки (рис. 2.3, в) более подробно. Видно, что она превышает на конце интервала 2.5%, нарастая с

47

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

дифференциальных	уравнений изменяются скачком (в точке
t = t' процесс u(t)	имеет разрыв). Автоматически изменяя

шаг интегрирования, процедура ode45 уменьшает возникающую при этом ошибку, но она все равно остается значительной и имеет тенденцию к накоплению из-за свойств решаемых уравнений. Значения на графике показаны точками именно в те моменты времени, которые получены (автоматически) процедурой интегрирования. Как видно из графика, интервалы между значениями t не одинаковы, наиболее "густо" расположены точки в начале и конце интервала интегрирования, а также в окрестности скачкообразного изменения входно-

го воздействия. Наименьшее значение шага		интегрирования
оказывается равным 10"5 с и имеет место в начале			процесса,
а наибольшее - 0.094 с получается для t = 0.34			с. 1 Мо-
менту времени	= 1.07 с соответствует шаг	интегрирования
0.034 с. Повысить	точность интегрирования	можно заданием
соответствующей	опции программы ode45.	Для этого заме-

ним строку, содержащую обращение к процедуре интегрирования в головной программе, на

opt=odesetСRelTol' ,1е-5); [t,х]=ode45('fр1 . 2',[0; Т],хО,opt);

Здесь с помощью оператора odeset установлено значение параметра относительной ошибки RelTol, содержащееся в ODE-структуре с именем opt, равное 10~5 вместо принятого по умолчанию значения 10~3. Структура opt передается в качестве дополнительного параметра при вызове процедуры ode45. Относительная ошибка интегрирования для рассматриваемой задачи оказывается существенно меньше и достигает 0.04 %. Машинное время интегрирования для того же компьютера составляет около 0.27 с. Задавая RelTol равным 10~7, получим наибольшую относительную ошибку около 4-10~4 % при времени счета 0.4 с. Число шагов интегрирования в последнем случае оказывается равным 125; наименьший шаг имеет длину 3 • 10~6 с и получается в окрестности

1	Для определения этих	значений можно использовать опера-
торы	[mi,ki]«min(diff(t)),	t(ki),	[ma,ka]«max(diff(t)), t(ka) и
kl«find(abs(t-tl)«min(abs(t-		t l ) ) ,	t(kl) - t(kl - 1) .

48

скачка входного воздействия. Наибольший шаг составляет

0.09с.

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

Так как рассматриваемое уравнение (1.4) является линейным ОДУ с постоянными коэффициентами, для его решения можно использовать специальные вычислительные методы, развитые для анализа линейных динамических систем. Эти методы реализованы, в первую очередь, в программах тулбокса CONTROL SYSTEMS, предназначенного для решения задач автоматического управления [6, 74, 61].

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

	х = Ах + Ви,		у = Сх + Du,
где вектор	состояния	x(t) и входное воздействие u(t)		были
определены	выше, а	вектор	выходов y(t) в данном примере
можно считать тождественно			равным я(£), а матрицы	Л, 5,

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

"0	0"	, В =0	,с =	'1	0"
— 1 0				0 1 —

Объединим эти матрицы в модели некоторой динамической системы с именем sys и выполним моделирование этой системы средствами тулбокса CONTROL SYSTEMS. Эти действия реализованы в следующей программе: 1

А= [0 0; 1 0 ] ;

В=[1/т; 0];

С = е у е ( 2 );

D = z e r o s ( 2 , l ) ; sys=ss(A,B,C,D); t=(0:Т/8еЗ:Т)1 ;

1 Моменты переключения t' и окончания процесса tm рассчитываются

так же, как и раньше, поэтому операторы вычисления t l и Т в данной программе не приводятся.

49

u l s b * o n e s ( s i z e ( t ) ); u2=-a*ones(size(t));

u=[ul(find(t<=tl)); u 2 ( f i n d ( t > t l ) ) ] ; x=lsim(sys,u,t);

Первые четыре строчки программы задают значения матриц А} В) С, D. Вызовом функции ss из этих матриц строится структура с именем sys, которая рассматривается дальше как модель линейной стационарной системы (ЛС-системы) в форме уравнений состояния. Оператором t=(0:T/8e3:T) ' ; интервал [O,t0] разбивается равномерно на 8000 точек, которые составляют массив t . 1 В следующих трех строчках формируется входное воздействие u(t): сначала генерируются векторы ul и и2 с одинаковыми значениями элементов, равными соответственно 6 и —а для всех точек t. Затем с помощью вырезки и объединения массивов получается последовательность требуемого вида. Возможно, этот же результат проще достичь оператором цикла

nt=length(t); f o r k=l:nt

if t(k)<=tl u(k)=b;

else

u(k)=-a; end

end

После получения массива u выполняется обращение к процедуре lsim, в которой происходит дискретизация модели системы (см. ниже п. 3.3.3) и рекуррентное вычисление реакции дискретной модели на заданное входное воздействие [4, б, 40, 74, 87]. В результате вычислений по этой программе при времени счета 0.16 с наибольшая относительная ошибка составляет около 0.025 %. Как отмечено в литературе [4, 40, 74, 87], эффективность этого метода интегрирования линейных ОЛУ по сравнению с универсальными методами растет с ростом порядка решаемой задачи и особенно велика

1 Данный оператор заканчивается операцией транспонирования ("штрих") только для удобства сочетания данной программы с ранее написанными фрагментами. Без этого массив t окажется вектор-строкой, а не столбцом, как принято выше.

50