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

при формальной подстановке выражения для Р = еАН, полученного из аппроксимаций Тейлора или Паде в (3.31), произвести "сокращение" матрицы Л. Тогда в выражение для Q матрица Л- 1 входить не будет. Например, аппроксимация по методу Эйлера Р = I + Ah приводит к формуле Q = h • В.

Другой способ состоит в расширении уравнений состояния исходной системы (3.21) [6, 74]. Входной процесс u(t) при tk < t < рассматривается как решение некоторого однородного дифференциального уравнения. Тогда расширенная система тоже является однородной и в вычислении по (3.31) нет необходимости. Искомые матрицы Р и Q получаются в качестве подматриц "расширенной" матричной экспоненты. 1 Продемонстрируем этот подход для входного процесса

Соответствующая дискретная модель (аналогично (3.23)) принимает вид

где Р = еА1х. Учитывая блочную структуру матрицы А и формулу (3.32) для Р, непосредственно убеждаемся, что матрица Р имеет следующую блочную структуру:

Р' Q'

Р= 0 1га

1Отметим, что именно этот способ реализован в программе c2d тулбокса CONTROL SYSTEMS пакета MATLAB, специально предназначенной для перехода от непрерывного описания к дискретному.

81

С учетом этого из (3.40) находим, что

Сравнивая (3.41) с (3.23), видим, что матрицы Р} Q в (3.23) совпадают с Р\ Q'. Поэтому они могут быть получены в качестве соответствующих подматриц матрицы Р = eAh .

Если непрерывная система нелинейна, то для перехода к ее дискретному описанию также можно использовать методы численного интегрирования. Например, метод Эйлера дает для системы (3.16) дискретное описание:

Континуализация - это переход от дискретной математической модели системы к непрерывной. Если дискретная модель системы имеет вид (3.23), то перейти к непрерывной модели (3.21) можно по формулам

ция, обратная к экспоненциальной и также определяемая через ряд

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.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 , т.е.

Подставим в полученное выражение начальное я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 .

Итак, получено решение линеаризованного уравнения вбли-

нить графики решений исходного уравнения и линеаризованной модели. Числовые значения можно получить непосредственной подстановкой в найденные выражения массива мо-

Рис. 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 в поле

х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. Сравним полученные численно и аналитически результаты. Для этого выполним операторы

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