Частым заданием в различных курсах по теории автоматического управления является нахождение матрицы 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
![](https://habrastorage.org/getpro/habr/upload_files/e3e/960/833/e3e960833183e897dd3e05d4c27c1982.jpg)
Устоявшееся конечное значение чуть больше единицы, попытаемся исправить это замкнув систему обратной связью -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
![](https://habrastorage.org/getpro/habr/upload_files/feb/b6d/729/febb6d729c432da2fa3e8c13b883f1c1.jpg)
Устоявшееся значение замкнутой системы около единицы, и сама реакция системы на ступенчатое воздействие раза в 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
Как можно увидеть, результат одинаков для каждого из способов.