Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков)
.pdfпри формальной подстановке выражения для Р = еАН, полученного из аппроксимаций Тейлора или Паде в (3.31), произвести "сокращение" матрицы Л. Тогда в выражение для Q матрица Л- 1 входить не будет. Например, аппроксимация по методу Эйлера Р = I + Ah приводит к формуле Q = h • В.
Другой способ состоит в расширении уравнений состояния исходной системы (3.21) [6, 74]. Входной процесс u(t) при tk < t < рассматривается как решение некоторого однородного дифференциального уравнения. Тогда расширенная система тоже является однородной и в вычислении по (3.31) нет необходимости. Искомые матрицы Р и Q получаются в качестве подматриц "расширенной" матричной экспоненты. 1 Продемонстрируем этот подход для входного процесса
u(t) = ик |
при tk < t |
< tk+\. |
Для указанного промежутка |
вре- |
|||
мени уравнение (3.21) запишем в виде |
|
|
|||||
x(t) |
= |
Ax{t) |
+ Bu{t), |
|
x{tk) = xkl |
tk<t<tk+u |
|
u{t) |
= |
0, |
u(tk)=uk. |
|
|
(3.38) |
|
Введем |
расширенный (n + т)-мерный |
вектор состояния |
|||||
х = со!(х, и) и (п + тп) х (п + |
т)-матрицу |
|
|
||||
|
|
|
А = |
А |
В |
|
|
|
|
|
0 |
Im |
|
|
|
|
|
|
|
|
|
||
Уравнение (3.38) представим в виде |
|
|
|||||
|
|
*(*) = |
Лх(0, |
х(1к) = со!(хьи*). |
(3.39) |
||
Соответствующая дискретная модель (аналогично (3.23)) принимает вид
Zjb+l = Рхк, |
(3.40) |
где Р = еА1х. Учитывая блочную структуру матрицы А и формулу (3.32) для Р, непосредственно убеждаемся, что матрица Р имеет следующую блочную структуру:
Р' Q'
Р= 0 1га
1Отметим, что именно этот способ реализован в программе c2d тулбокса CONTROL SYSTEMS пакета MATLAB, специально предназначенной для перехода от непрерывного описания к дискретному.
81
С учетом этого из (3.40) находим, что
zjt+i = |
+ Q V |
(3.41) |
Сравнивая (3.41) с (3.23), видим, что матрицы Р} Q в (3.23) совпадают с Р\ Q'. Поэтому они могут быть получены в качестве соответствующих подматриц матрицы Р = eAh .
Если непрерывная система нелинейна, то для перехода к ее дискретному описанию также можно использовать методы численного интегрирования. Например, метод Эйлера дает для системы (3.16) дискретное описание:
= |
+ hF(xkl ик) tk)y yk = G(xkl ик) tk). |
(3.42) |
Континуализация - это переход от дискретной математической модели системы к непрерывной. Если дискретная модель системы имеет вид (3.23), то перейти к непрерывной модели (3.21) можно по формулам
Л = |
B = \\nP(P-I)-1Q, |
(3.43) |
вытекающим из (3.31), |
где In Р - логарифм матрицы, |
функ- |
ция, обратная к экспоненциальной и также определяемая через ряд
1п(1 + х ) = х + } £ + - . - + Н Г * 1 ^ , сходящийся при ||Х|| < 1 (здесь X = Р - I). С точностью до величин порядка Л2 можно ограничиться формулами
соответствующими методу Эйлера. Однако удобнее всего переходить от дискретной передаточной функции к непрерывной по формулам (3.33) и (3.36). Например, по методу Эйлера (3.33) достаточно заменить в передаточной функции WR(z~l) переменную г"1 на 1 — hp.
При исследовании линейных систем получили распространение также методы упрощения описаний систем путем редукции (понижение порядка) [64, 115]. Взаимосвязь различных описаний динамических систем представлена на рис. 3.10.
82
Непрерывные |
Лискретизация |
Лискретные |
||||
нелинейные |
нелинейные |
|||||
(3.42) |
||||||
УС |
(3.16) |
УС |
(3.22) |
|||
(3.20) |
Линеаризация |
Линеаризация |
(3.20) |
|||
|
|
Континуализация |
|
|||
Непрерывные |
(3.43) |
Дискретные |
||||
линейные |
линейные |
|||||
Дискретизация |
||||||
УС |
(3.21) |
УС (3.23) |
||||
(3.24) |
|
(3.31) |
|
(3.29) |
||
(3.29) |
|
(3.24) |
||||
|
|
|
||||
|
|
Континуализация |
|
|||
Непрерывная |
(3.33), |
(3.36) |
Дискретная |
|||
ПФ (3.24) |
Дискретизация |
ПФ |
||||
|
|
|
|
|||
|
|
(3.33), |
(3.36) |
|
|
|
Редукция |
|
|
|
|
Редукция |
|
Рис. 3.10. Взаимосвязь различных описаний динамических систем.
83
3.4, Примеры преобразования моделей в среде MATLAB
3.4.1.Линеаризация
Рассмотрим сначала процедуру линеаризации для примера 3.3.1 (с. 76). Эта задача может быть решена в общем виде с помощью тулбокса SYMBOLIC.
Рассмотрим динамическую систему х = — х2и с начальным состоянием х0. Решение этого дифференциального нелинейного уравнения найдем оператором
x=dsolve(, Dx=-x~2*u','х(0)=х0')
Получен ответ: х =l/(u*t+l/xO) , или в формате Mj?X - х = (ut + Xfl"1)"1. Продифференцируем теперь правую часть исходной модели по переменным х}и :
d f x = d i f f ( ' - x ~ 2 * u ' , ' х ' ) dfu=diff 0-x~2*u' , ' u ' )
Получим dfx =-2#x*u , dfu =-x~2 . Далее выполним подстановку: вместо символа х подставим найденное выше решение, а вместо символа и - значение 1.
A=subs(dfxj-t'x', 'u'},-Cx, 1})
I^subsCdfu.-C'x', ' u ' J ^ x , 1})
В результате получаем матрицы (в данном примере - скалярные функции) линеаризованной модели х = A(t)x + B(t)u.:
A=-2/(t+l/x0)
В= - l/(t+l/x0)~2
Вформате IATgX это будут выражения: Л(^) = -2(^+х0 "1 )"1 , B(t) = — (t + x0~l)~2. Видим, что результат совпадает с приведенными в примере 3.3.1 формулами. Продолжим пример.
Пусть задано х0 = l}u(t) = 1. Получим решение исходного уравнения подстановкой этих данных в формулу для общего решения:
x=subs(x,{'xO','u'},{l,l})
Находим х = l / ( t + l ) . Теперь перейдем к решению линеаризованного уравнения. Для этого скопируем строки полученных выше выражений для Л, В в оператор интегрирования дифференциальных уравнений:
84
x _l=s imple( d sol v e ('Dx=-2/ (t+l/ x O)*x-l /(t +l /xO )~2*u', . ..
>x(0)=x0'))
Получим ответ: x_l=-x0*(u*t*x0-l)/(t*x0+l)~2 , т.е.
m_ |
xojutxo- 1) |
X,[t)~ |
(tx0 + 1У • |
Подставим в полученное выражение начальное я0 = 1 и и = 1: x _ l = s i m p l e ( s u b s ( x _ l ' х О ' , ' и ' } , { 1 , 1 » )
Получаем следующий ответ: х_1= - ( t - l ) / ( t + l ) ~ 2 .
Итак, получено решение линеаризованного уравнения вбли-
зи заданной траектории: xi(t) = — ~~ |
Интересно срав- |
(t + 1) |
|
нить графики решений исходного уравнения и линеаризованной модели. Числовые значения можно получить непосредственной подстановкой в найденные выражения массива мо-
ментов времени |
t: |
|
|
|
|
|
|
|
||
|
t=0:0 . 025:0. 5; |
|
|
|
|
|
|
|
||
|
x = s u b s ( x , ' t ' , t ) ; |
|
|
|
|
|
|
|
||
|
xl=subs(x_l, Ч ' |
, t ) ; |
|
|
|
|
|
|
||
|
h = p l o t ( t , x , , t , x _ l , ' - o ' ) ; |
|
|
|
|
|
||||
|
legend(h, ' x ' , ' x . l O |
|
|
|
|
|
|
|||
|
xlabel( 4 ' ) |
, уlabel ( 'x, |
х . Г ) |
|
|
|
|
|
||
Результат вычислений показан на рис. 3.11. |
|
|
|
|||||||
а |
|
|
|
|
б |
|
|
|
|
|
1 |
J |
1 |
1 |
J |
|
• |
1f 1 — |
- 1 |
" • • |
Т •• |
|
|
|
|
|
|
Л |
|
|
|
|
0.8 |
|
|
|
|
-0.5 |
• |
в. |
|
|
|
К |
1 |
|
|
|
А |
|
|
|
||
0.6 |
|
|
— |
|
В |
|
|
|
||
—— X |
|
|
-1 |
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|||
0.4 |
|
|
|
|
|
|
|
|
||
— *m |
|
|
|
|
|
|
|
|
|
|
0.2 |
|
|
|
•1.5 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
i |
= |
0.3 |
0.4 |
t - |
|
0.1 0.2 |
0.3 |
0.4 |
t |
|
0.1 |
0.2 |
|
|||||||
Рис. 3.11. Графики решений исходного (а) и линеаризованного (6) уравнений и параметры линейной модели.
Обратимся теперь к численному решению задачи. Составим S-модель исходной системы (рис. 3.12). Блок и2 явля-
85
Рис. 3.12. 5-модель исходной системы.
ется блоком Math function из раздела Nonlinear. Модель составлена таким образом, что на нее можно подавать внешнее входное воздействие. Вначале используем эту программу для расчета опорной траектории x(t) при х(0) = l}u(t) = 1. Для того чтобы передать эти величины из "рабочей области" (Workspace), укажем хО в строке Initial condition окна параметров интегратора. Кроме того, в меню Simulation (опция
Parameters, окно Workspace I/O, поле Load from workspace) установим параметр Input - и, а в поле Save to workspace - параметры: Time - tout, Output - yout. Под этими именами значения результатов моделирования сохраняются в рабочей области и используются MATLAB-программой. В окне Solver в поле
Solver |
options параметр Max step size (максимальная |
величина |
шага) |
установлен равным 0.005. Выполним теперь |
операторы |
х0=1; u = ' 1 ' ;
sim('S3 . 1')
Результат вычислений обозначен на рис. 3.11 через хш. Соответствующий график полностью сливается с графиком аналитического решения х. Заметим, что входной процесс имеет не числовое, а строковое значение. В строке задается способ нахождения u(t). Это связано с тем, что значение и должно рассчитываться на каждом шаге интегрирования. Можно также задать и массивом соответствующей размерности.
Перейдем теперь к численной линеаризации системы. Для этого используем ту же модель (рис. 3.12), но изменим содержимое поля Load from workspace, исключив из него параметр Input. Это вызвано тем, что процедура линеаризации изменяет входы и состояния системы автоматически. Выполним операторы:
uk=l;
f o r k = l : l e n g t h ( t ) ;
86
хО=х(к);
[А1(к),В1(к),С,D]=linmod(' s3_l',xO,uk); end
В этом фрагменте программы циклически изменяются значения состояния системы, в окрестности которого выполняется линеаризация. Эти значения, обозначенные хО, получаются из рассчитанного выше массива х точек опорной траектории. Устанавливается и значение входного воздействия и = 1, обозначенного uk. Далее в цикле выполняется обращение к процедуре linmod, входящей в состав тулбокса SIMULINK. В результате вычислений получаются коэффициенты (в общем случае - матричные) уравнений состояния линеаризованной системы х = Ах + Ви, у = Сх + Du. Сравним полученные численно и аналитически результаты. Для этого выполним операторы
A=subs(A, Х О ' , 1 ) ; |
A=subs(A, Ч ' ,t) ; |
B=subs(B,'х_0',1); |
B=subs(B,Ч',t); |
p l o t ( t , А 1 , ' о ' , t , B l , , t ,A,t,B);
t i t l e ( ' Параметры линеаризованной модели ') x l a b e l ( Ч ' )
Получим график, на котором сплошными линиями показаны графики изменения параметров модели, полученные аналитически, а символами о и * - графики, рассчитанные процедурой linmod, рис. 3.11,5. Относительная ошибка определения
A(t) найдена оператором max(abs(A-Al))/max(abs(A))*100 и
составляет 5 • 10~4 %. Аналогично для B(t) получим относительную ошибку 5 • Ю"10 %.
3.4.2.Дискретизация моделей
Процедуру дискретизации модели продемонстрируем на примере синтеза цифрового фильтра (ПФ) по аналоговому фильт- ру-прототипу (3.8), принципиальная схема которого приведена ранее на рис. 3.2. Пусть требуется найти передаточную Функцию (другими словами, алгоритм вычислений) цифрового фильтра, который обеспечивает подавление сигнала в заданном диапазоне частот. Перейдем от непрерывной системы, заданной передаточной функцией (ЗЛО), к дискретной Модели. При построении дискретной модели следует задать
87
шаг дискретности Т0 и сделать определенное допущение о виде процесса на входе непрерывной системы (более подробно см., например, [3, 6, 75]). Здесь будем предполагать, что u(t) - кусочно-постоянная внутри интервалов квантования. Запишем фрагмент программы.
Программа синтеза ЦФ по аналоговому прототипу
fl=tf(num,den) Т0=5е-3; dfl=c2d(fl,T0)
Здесь предполагается, что массивы num, den уже заданы так, как это сделано в программах на с. 64, 65. Функция tf относится к тулбоксу CONTROL SYSTEMS и служит для построения линейного стационарного объекта (ЛС-обзект) по массивам числителя и знаменателя передаточной функции [6, 61]. В программе задается значение интервала дискретности Т0 = 5 • 10~3 с и выполняется обращение к функции c2d перехода от непрерывной модели к дискретной. Лалее выполним анализ полученной системы и ее сравнение с исходной. Продолжим программу:
omega=logspace(l,3,300);
[mag,phase]=bode(fl,omega);
[dmag,dphase]=bode(dfl,omega);
Lc=20*logl0(squeeze(mag));
Ld=20*logl0(squeeze(dmag));
cph=squeeze(phase);
dph=squeeze(dphase);
subplot(211)
semilogx(omega,Lc, 'k' ,omega,Ld,' — .k' ) ,grid
t i t l e C Диаграммы |
Боде ') |
|
ylabel('L(\omega), |
дБ') |
|
subplot(212) |
|
|
semilogx(omega, |
cph,'k*,omega, d p h , ' - - . k ' ) , g r i d |
|
xlabel('\omega, |
р а д . / с ' ) |
|
Здесь строятся диаграммы Боде (JIAX, ЛФЧХ) непрерывной и дискретной систем. Вычисления выполняются функцией bode. Для непрерывных систем осуществляется описанная выше подстановка $ = ju> аргумента передаточной функции. Аргумент z дискретной передаточной функции заменяется выражением е*шТ° [6, 75]. Функция squeeze используется
88
для удаления единичных размерностей массивов [84]. В данном примере это массивы mag , dmag , cph, dph, имеющие по три измерения. Без этого действия не удается выполнить операторы semilogx и plot. Частотные характеристики приведены на рис. 3.13, а.
Перейдем к расчету переходных характеристик непрерывной и дискретной систем. Воспользуемся функцией step тулбокса CONTROL SYSTEMS [61]. Лля получения более подробных сведений о виде переходной характеристики непрерывной системы выберем для нее меньший шаг вычислений. Шаг расчета переходной характеристики дискретной системы должен совпадать с ее интервалом квантования Т0. Результаты вычислений выводятся операторами plot - для непрерывной системы и stem - для дискретной. Чтобы они разместились на одном графике, используется оператор hold on (см. рис. 3.13, б). Видно, что в узлах квантования выход непрерывной системы совпадает с выходом ее дискретной модели, что вполне соответствует принятому методу дискретизации (в теории ИФ этот метод так и называется - метод инвариантности переходной характеристики (ИПХ) [19, 24]).
tc=0:5е-4:0.25; td=0:T0:0.25; y c = s t e p ( f l , t c ) ;
yd=step(dfl,td);
p l o t ( t c , у с , ' k ' ) , h o l d on stem(td,yd),grid
З а м е ч а н и е . Можно значительно упростить данную программу, если не задаваться точками диапазона частот, моментов времени и видом графиков. Тогда рассмотренные здесь действия можно выполнить программой:
fl=tf(num,den) T0=5e-3; dfl=c2d(f1,Т0) bode(fl,df1); f i g u r e step(f1,df1);
89
Рис. 3.13. Характеристики аналогового прототипа и цифрового фильтра .
Выбор требуемых параметров и вывод результатов на экран производятся функциями bode, step автоматически.
Синтез по аналоговому прототипу - далеко не единствен-
ный и не самый распространенный в настоящее время |
метод |
||
синтеза ЦФ. Тулбокс SIGNAL PROCESSING содержит доста- |
|||
точно много программ, в которых реализованы |
различные |
||
методы синтеза цифровых фильтров, как рекурсивных |
(филь- |
||
тры с импульсной характеристикой бесконечной |
длительно- |
||
сти, или БИХ-фильтры), так и нерекурсивных |
(фильтры с им- |
||
пульсной характеристикой конечной длительности, |
или |
КИХ- |
|
фильтры) [19, 35, 84]. Продемонстрируем применение |
функ- |
||
ции butter синтеза фильтров Баттерворта |
(Butterworth) для |
||
получения цифрового режекторного фильтра с указанной выше полосой задерживания и интервалом дискретности Т0 = 0.005 с.
При использовании функции butter следует указать границы полосы пропускания (или задерживания - в зависимо-
сти |
от типа фильтра) в единицах частоты |
Найквиста ЦФ |
u N |
= ттг с"1. Например, если полоса пропускания равна 0.2и>уу |
|
|
* о |
|
с~\ |
то соответствующий параметр устанавливается равным |
|
0.2 и т.д. Учитывая это, запишем следующие |
операторы: |
|
|
ошр= [70, 100]/pi*T0; |
|
|
[а,b,с,d]=butter(2,omp,'stop') |
|
|
fb=ss(a,b,c,d,TO) |
|
Параметр omp является вектором, в котором указывается нормированная полоса задерживания. Во второй строке происходит обращение к функции butter, где указаны порядок
90