Модальный метод синтеза в MATLAB

Частым заданием в различных курсах по теории автоматического управления является нахождение матрицы K для модального управления системой вида dx/dt = Ax+Bu y = Cx.

Такой тип задач легко решается в среде MATLAB.

Сперва наперво требуется задать нашу систему. Для примера возьмем типовую модель электродвигателя:

J = 0.01;   % Момент инерции ротора kg.m^2
b = 0.1;    % Коэффициент затухания мех. системы Nms
K = .01;    % ЭДС константа Nm/A
R = 1.0;    % Сопротивление Ohms
L = 0.5;    % Индуктивность Henrys
A = [-b/J K/J; -K/L -R/L];
B = [0; 1/L];
C = [1 1];
D = 0;
W = ss(A,B,C,D); % создаем модель в пространстве состояний

Далее проверим критерий устойчивости и управляемость системы

pole(W) % найдем полюса системы
% ans = 2×1 -9.9975 -2.0025
rank(ctrb(A,B)) % если ранг равен порядку системы то она полностью управляема
% ans = 2

Полюса все оба два отрицательны, а значит система устойчива. Ранг матрицы управления равен порядку системы, а значит система полностью управляема.

Найдем реакцию системы на ступенчатый сигнал:

step(W)
grid on

Устоявшееся конечное значение чуть больше единицы, попытаемся исправить это замкнув систему обратной связью -Kx так, чтобы полюса итоговой системы оказались равны p = 4.45*[-1 -1.1] и конечное значение оказалось близко к единице.

В пакете MATLAB это делается буквально парой строк:

p = 4.45*[-1 -1.1]; % заданные полюса
K = place(A,B,p) % матрица К
% K = 1×2 14.1564 -1.3275

% Проверим полюса замкнутой системы
AClosed = A - B*K;
BClosed = B;
CClosed = C;
DClosed = D;
Wclosed = ss(AClosed,BClosed,CClosed,DClosed);
pole(Wclosed)
% ans = 2×1 -4.8950 -4.4500

Как можно видеть, полюса замкнутой системы точно такие же, как в векторе p.

Найдем реакцию замкнутой системы на ступенчатый сигнал:

step(Wclosed)
grid on

Устоявшееся значение замкнутой системы около единицы, и сама реакция системы на ступенчатое воздействие раза в 2 быстрее.

Теперь, для сравнения, можно найти ту же матрицу K классическим алгоритмом:

% Запишем матрицу замкнутой системы
L = length(B);
k = sym('k',[1,L]);
F = A - B*k % найдем матрицу F
% Найдем характеристический полином
syms lamda
polyF = det(lamda*eye(L)-F) % характеристический полином
CoefF = fliplr(coeffs(polyF,lamda)) % коэффициенты полинома
CoefF(1) = [];
% Назначим желаемый характеристический полином
polyL = (lamda-p(1))*(lamda-p(2)) % желаемый полином
CoefL = fliplr(coeffs(polyL,lamda)) % коэффициенты полинома
CoefL(1) = [];
% Найдем элементы матрицы К
K = solve(CoefF == CoefL);
K = double([K.k1 K.k2]) % матрица К
% K = 1×2 14.1564 -1.3275

Строк кода много больше, но ответ вышел точно такой же.

UPD1

В MATLAB есть еще одна функция для вычисления матрицы K, которая использует формулу Аккермана, но которая, в отличие от функции place, работает только если система имеет один вход.

K = acker(A,B,p) % матрица К
% K = 1×2 14.1564 -1.3275

Еще один способ найти коэффициенты матрицы усиления K это решение матричного уравнения Сильвестра.

Запишем матрицы описания эталонной модели:

G = [0 -CoefL(2); 1 -CoefL(1)]; % матрица nxn определяющая требуемые свойства
H = [0 1]; % матрица mxn выхода эталонной модели

Решим уравнение Сильвестра:

% найдем такую матрицу М, чтобы выполнялось M*G-A*M-B*H == 0
M = lyap(-A,G,-B*H)
% ans M = 2×2
% [ -0.0265 -0.0176;
% -0.2826 0.5657]

Найдем матрицу усиления:

K = -H*M^-1
% K = 1×2 14.1564 -1.3275

Найдем полюса системы:

eig(A - B*K)
% ans = 2×1 -4.8950 -4.4500

Как можно увидеть, результат одинаков для каждого из способов.

Ссылка на исходники