Модальный метод синтеза в 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

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

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