Модальный метод синтеза в 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Строк кода много больше, но ответ вышел точно такой же.
Ссылка на исходники
