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

стей дифференциальных уравнений 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 (£',*»]). На каждом участке движение определяется сво-

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

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

рует вектор индексов элементов массивов по указанному в качестве аргумента логическому выражению. Этим логическим выражением здесь служат операции сравнения 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'. Появление такой заметной ошибки связано с тем, что рассматриваемая система "неудобна" с точки зрения численного интегрирования, так как правые части

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

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

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 с и получается в окрестности

48

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

0.09с.

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

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

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

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

Объединим эти матрицы в модели некоторой динамической системы с именем 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

для так называемых жестких систем.1 Сведения о жестких ОДУ и методах их интегрирования приведены в [74, 87, 97].

2.2.4.Применение системы SIMULINK

Рассмотрим решение задачи о постановке машины в гараж (п. 1.3) в среде системы визуального моделирования SIMULINK. Программа моделирования SIMULINK, хоть и имеет текстовое представление, создается и используется исследователем в виде блок-схемы (S-моделъ), в которой содержится описание рассматриваемой системы и задаются требуемые действия

Рис. 2.4. 5-модель системы (1.4).

нала. Указаны следующие значения параметров настройки:

Step time (длительность ступеньки) - t l ; Initial value (начальное значение) равно b; Final value (конечное значение) равно

-а.

Блок Gain - усилитель. Параметр Gain (коэффициент передачи) установлен равным 1/т.

Блок Integrator - интегрирующий блок. Приняты заданные по умолчанию значения параметров, в частности нулевые начальные условия (параметр Initial condition).

Блок Scope - индикатор. Служит для вывода результатов в процессе интегрирования в "электронное окно". Параметры блока могут быть настроены после моделирования с помощью выводимых на панели окна кнопок [29].

Диалоговое окно 5-модели имеет меню Simulation (моделирование), которая включает опцию Parameters (параметры).

1 Следует заметить, что для рассматриваемого примера точность в значительной степени зависит от того, насколько близко лежат точки i', К к дискретной последовательности t. При попытке использовать lsim с неравномерным шагом по времени выводится сообщение об ошибке.

51