Главная / Математика / Запись математических моделей систем

Запись математических моделей систем

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

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

При математическом описании систем методами пространства состояний все переменные, характеризующие систему, можно разделить на три группы: 1) входные переменные или входные воздействия Ui(t), I=1,…,m управляемые и wI(t), I=1,2,…,l неуправляемые, генерируемыми устройствами или системами, внешними по отношению к данной; 2) выходные переменные Yj(t), J = 1,…,r, характеризующие реакцию системы и воздействие измерительных шумов
Yj(t), J = 1,…,r и 3) переменные состояния математической модели или промежуточные переменные Xk(t), K=1,n, характеризующие динамическое поведение системы. Для удобства оперирования с многомерными величинами эти величины записываются в виде векторов; — вектор входа, — вектор выхода, — вектор состояния, – вектор возмущения, — вектор шумов измерений, ¢ – значок транспонирования.

Для многосвязного объекта управления могут возникнуть сложности в записи уравнения состояния, даже если есть уравнения его отдельных частей. Для решения задачи целесообразно составить схему аналогового моделирования отдельных частей и всего объекта в целом. В качестве переменных X(t) состояния взять выходные сигналы интеграторов. Тогда правым частям дифференциальных уравнений будут соответствовать сигналы на входах интеграторов, формирование которых легко прослеживается по схеме моделирования.

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

Где – A(n´N) матрица состояния, B(n´M) – матрица входа, C(n´R)— матрица выхода, G(n´L) – Матрица возмущающих воздействий.

Для удобства описания модели вектора U(t) и w(t) Могут быть объединены в один, тогда матрица G Выйдет в B. Объединенная модель может быть составлена из отдельных частей средствами Matlab. Этой цели непосредственно служит функция Blkbuilt, а также функции Append, Cloop, Connect, Feedback, parallel, series.

Уравнение в переменных состояния может быть записано не единственным образом. Это зависит от выбора переменных состояния, что может определяться, в частности, и доступностью их измерения.

Пример. Записать уравнения состояния для двигателя постоянного тока, учитывающие : напряжения U якоря, ток I Якоря, частоту J вращения вала и момент МН нагрузки на валу. Управляющей величиной является напряжение U. Измеряются частота вращения с помощью тахогенератора и ток I Якоря в виде падения напряжения на шунте. Коэффициент передачи КТГ тахогенератора и сопротивление шунта заданы так же, как и сопротивление R и индуктивность L якоря цепи, суммарный момент инерции J двигателя и нагрузки, момент инерции J двигателя, коэффициенты и двигателя.

По закону Ома для якорной цепи можно записать уравнение

По закону Ньютона для вращательного движения уравнение будет иметь вид .

Изобразим оба уравнения в виде схемы моделирования.

Рис.1

Для записи уравнения состояния обозначим

Тогда уравнения состояния будут иметь вид

Следовательно, матрицы A, B, C, D, G можно записать как

.

Пример Для структурной схемы рис. 1 предыдущего примера составить и исследовать модель двигателя МИГ-90А. Результаты представить в форме MATLAB-программы и коментариев к ней..

Решение:

% Расчет и определение коэффициентов модели двигателя МИГ-90А

% по паспортным данным.

Un=27; In=4.5; Fi=315; Rjk=1.1; Tjk=0.005; L=Tjk*Rjk; J=0.2e-4;

Ce=(Un-In*Rjk)/Fi; Rsh=0.01*Rjk; Mn=0.286; Cm=Mn/In; Ktg=5;

% Запись передаточных функций отдельных блоков; ni-используется для

% определения числителя, di-используется для определения знаменателя

% i-го блока.

N1=1; d1=1; n2=1/L; d2=1; n3=1; d3=[1 0]; n4=Ce; d4=1; n5=Rsh; d5=1; n6=Cm; d6=1; n7=1/J; d7=1; n8=1; d8=[1 0]; n9=Ktg; d9=1; n10=Rjk; d10=1; n11=1; d11=1;

Nblocks=11; % Задание общего количества блоков в системе.

% Объединение блоков в систему без их соединения.

Blkbuild; % ***Замечание: Числовые параметры системы будут записаны в матрицы с зарезервированными именами a, b,c, d

% Запись матрицы связей q между блоками.

q=[ 1 0 0 0 % Первый символ обозначает номер описываемого блока,

2 1 -8 -10 % последующие указывают номера блоков, выходы которых

3 2 0 0 % должны быть связаны с входом описываемого блока.

4 8 0 0 % Минус указывает на знак обратной связи.

5 3 0 0

6 3 0 0

7 6 -11 0

8 7 0 0

9 8 0 0

10 3 0 0

11 0 0 0];

Iu=[1 11]; % Определение номеров входных блоков.

Iy=[5 9]; % Определение номеров выходных блоков.

% Соединение блоков согласно матрице связей q в единую систему.

[ac, bc, cc, dc]=connect(a, b,c, d,q, iu, iy); % ac, bc, cc, dc — матрицы системы.

Format short e % Изменение формата на экспоненциальный, короткий.

ac; pause % Вывод матрицы ac c паузой;

% ac = -2.0000e+002 -1.8182e+002

% 3.1778e+003 0.

bc; pause % Вывод матрицы bc c паузой;

% bc = 1.8182e+002 0

% 0 -5.0000e+004.

cc; pause % Вывод матрицы cc c паузой;

% cc = 1.1000e-002 0

% 0 5.0000e+000.

dc; pause % Вывод матрицы dc c паузой;

% dc = 0 0

% 0 0.

% Матрицы A, B, C, D модели двигателя можно записать в общем виде

A = [-Rjk/L — Ce/L; Cm/J 0 ] % здесь ; разделяет строки матрицы,

B = [ 1/L 0; 0 -1/J ] % а пробелы разделяют элементы строк.

C = [ Rsh 0; 0 Ktg ]

D = [ 0 0 ; 0 0 ]

Printsys(A, B, C, D) % Демонстрирует модель системы.

Eig(ac) % Вычислим собственные значения матрицы системы ас.

% Результаты Matlab запомнит в служебной переменной ans;

% ans = -1.0000e+002+ 7.5351e+002i

% -1.0000e+002- 7.5351e+002i.

[V, D]=eig(ac) % Вычислим собственные значения и собственные векторы;

% D-диагональная матрица с собственными значениями;

% V-матрица, столбцы которой — собственные вектора;

% V = -2.3061e-001 — 3.0605e-002i; -2.3061e-001 + 3.0605e-002i;

% 0 + 9.7256e-001i; 0 — 9.7256e-001i;

% D = -1.0000e+002 + 7.5351e+002i 0

% 0 -1.0000e+002 — 7.5351e+002i.

Damp(ac) % Вычислим собственные частоты и коэффициенты

% относительного демпфирования;

% Eigenvalue Damping Freq. (rad/sec)

% -1.0000e+002+ 7.5351e+002i 1.3156e-001 7.6012e+002

% -1.0000e+002 — 7.5351e+002i 1.3156e-001 7.6012e+002

% Переход от модели в пространстве состояний к передаточным

% функциям:

[num1,den1]=ss2tf(ac, bc, cc, dc,1) % от первого входа к обоим выходам.

% В строках матрицы num1 записаны коэффициенты числителя

% передаточной функции по соответствующему выходу в порядке

% убывания степени s ;

% num1 = 0 2.0000e+000 2.3283e-010

% 0 0 2.8889e+006 .

% Коэффициенты знаменателя в векторе den1 также расположены в

% порядке убывания степени s;

% den1 = 1.0000e+000 2.0000e+002 5.7778e+005.

[num2,den2]=ss2tf(ac, bc, cc, dc,2) % от второго входа к обоим выходам

% num2 = 0 0 1.0000e+005

% 0 -2.5000e+005 -5.0000e+007

% den2 = 1.0000e+000 2.0000e+002 5.7778e+005.

% Перейдем от непрерывной модели к дискретной:

% x[k+1] = Ad x[k] + Bd u[k].

[Ad, Bd] = c2d(ac, bc,0.01) % c шагом дискретности 0.01 c.

% Ad = 6.8970e-002 -8.4292e-002

% 1.4732e+000 1.6169e-001;

% Bd = 8.4292e-002 1.3190e+001

% 8.3831e-001 -3.7690e+001 .

% Построим графики АЧХ.

% Определим частотный диапазон от 0.1 до 10000 рад/с в логарифми-

% ческом масштабе с количеством точек равным 100;

W=logspace(-1,4,100);

% Вычислим АЧХ ( mag ) и ФЧХ ( phase ) от первого входа.

[mag, phase]=bode(ac, bc, cc, dc,1,w);

% Построим график АЧХ от первого входа к первому выходу.

Semilogx(w, mag(:,1));pause

% ***Замечание. Запись mag(:,1) означает выбор элементов первого

% столбца, а mag(1,:) первой строки двумерного массива mag.

% Построим график АЧХ от первого входа к второму выходу.

Semilogx(w, mag(:,2));pause

% Построим графики АЧХ от первого входа к обоим выходам.

Semilogx(w, mag);pause % ***Замечание. Графики строятся в одном масштабе

% Построим график АЧХ от первого входа к второму выходу c

% логарифмическим масштабом по обеим осям.

Loglog(w, mag(:,2));pause

% Вычислим и построим ЛАХ и ФЧХ без записи в переменные.

Bode(ac, bc, cc, dc,1) % ЛАХ и ФЧХ по первому входу для обоих выходов;

Margin(ac, bc, cc, dc,1) % к ЛАХ и ФЧХ выведет запасы по фазе и амплитуде.

% Построим графики реакции системы на ступенчатое воздействие

T=0:0.001:0.1; % определим временной диапазон

% Вычислим значения переходных функций, для переменных

% состояния Х и для выходных переменных Y. Ступенчатое

% воздействие подаем на первый вход.

[Y, X] = step(ac, bc, cc, dc,1,t);

Plot(t, X(:,1));pause % график первого компонента вектора состояния;

Plot(t, X(:,2)); pause % график второго компонента вектора состояния;

Plot(t, X); pause % графики обоих компонентов вектора состояния;

Plot(t, Y(:,1)); pause % график первого компонента вектора выхода;

Plot(t, Y(:,2)); pause % график второго компонента вектора выхода;

Plot(t, Y); pause % графики обоих компонентов вектора выхода.

% Вычислим и построим переходную характеристику по первому

% входу без записи в переменные.

Step(ac, bc, cc, dc,1).

% Построим графики на произвольное, гармоническое воздействие

T=0:0.01:3; % во временном диапазоне от 0 до 3.

U=zeros(length(t),2); % Определим нулевую матрицу U для воздействия.

% Вычислим гармоническое воздействие (f=5 гц, А=1), которое подадим

% только на первый вход:

For i=1:length(t)

U(i,1)=sin(t(i)*2*pi*5);

End

Plot(t, U); pause % график возмущений

% Вычислим реакцию системы на гармоническое входное воздействие:

[Y, X] = lsim(ac, bc, cc, dc, U,t);

Plot(t, X(:,1)); pause % график первого компонента вектора состояния;

Plot(t, X(:,2)); pause % график второго компонента вектора состояния;

Plot(t, X); pause % графики обоих компонентов вектора состояния;

Plot(t, Y) % графики обоих компонентов вектора выхода.

Xlabel(‘ time, s’) % Подпишем ось t.

Ylabel(‘output’) % Подпишем ось Y.

% Подпишем название графика:

Title(‘График реакции на гарм. воздействие с 1-го входа (f=5 гц, А=1)’)

Оставить комментарий