Первый этап составления математической модели объекта управления состоит в выделении объекта из окружающей среды. Что отнести к объекту, а что к окружающей среде, зависит не только от свойств объекта, но и от целей моделирования. Для более точного и тонкого моделирования модель расширяют за счет окружающей среды, учитывают большое число связей и внешних факторов.
Для придания модели математической формы реальные процессы представляются однонаправленными, функциональными преобразователями входных переменных в выходные. Двунаправленные процессы моделируются за счет усложнения структуры.
При математическом описании систем методами пространства состояний все переменные, характеризующие систему, можно разделить на три группы: 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 Якоря в виде падения напряжения на UШ шунте. Коэффициент передачи КТГ тахогенератора и сопротивление RШ шунта заданы так же, как и сопротивление 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)’)