Частым заданием в различных курсах по теории автоматического управления является нахождение матрицы 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
Как можно увидеть, результат одинаков для каждого из способов.
