Задайте модель объекта управления

Модель объекта управления является системой в пространстве состояний дискретного времени, и это - нестабильный разомкнутый контур. Мы принимаем, что все состояния объекта измеримы. Поэтому мы избегаем потребности в разработке средства оценки состояния, которое выходит за рамки этого примера.

A = [1.1 2; 0 0.95];
B = [0; 0.0787];
C = [-1 1];
D = 0;
Ts = 1;
sys = ss(A,B,C,D,Ts);
x0 = [0.5;-0.5]; % initial states at [0.5 -0.5]

Спроектируйте неограниченный Линейный квадратичный регулятор (LQR)

Спроектируйте неограниченный LQR с выходным взвешиванием. Этот контроллер служит базовой линией, чтобы соответствовать пользовательскому алгоритму MPC. Законом о надзоре LQ является u(k) = -K_lqr*x(k).

Qy = 1;
R = 0.01;
K_lqr = lqry(sys,Qy,R);

Запустите симуляцию с начальными состояниями в [0.5 - 0.5]. Ответ с обратной связью устойчив.

t_unconstrained = 0:1:10;
u_unconstrained = zeros(size(t_unconstrained));
Unconstrained_LQR = tf([-1 1])*feedback(ss(A,B,eye(2),0,Ts),K_lqr);
lsim(Unconstrained_LQR,'-',u_unconstrained,t_unconstrained,x0);
hold on

Спроектируйте пользовательский контроллер MPC с терминальным весом

Спроектируйте пользовательский контроллер MPC с терминальным весом, примененным на последнем шаге предсказания.

Предсказанные последовательности состояния, X (k), сгенерированный линейной моделью и входной последовательностью, U (k), могут быть сформулированы как: X(k) = M*x(k) + CONV*U(k). В этом примере используйте четыре шага предсказания (N = 4).

M = [A;A^2;A^3;A^4];
CONV = [B zeros(2,1) zeros(2,1) zeros(2,1);...
        A*B B zeros(2,1) zeros(2,1);...
        A^2*B A*B B zeros(2,1);...
        A^3*B A^2*B A*B B];

Целевой функцией MPC является J(k) = sum(x(k)'*Q*x(k) + u(k)'*R*u(k) + x(k+N)'*Q_bar*x(k+N)). Гарантировать, что целевая функция MPC имеет ту же квадратичную стоимость как бесконечный горизонт квадратичная стоимость, используемая LQR, терминальный вес Q_bar получен путем решения следующего уравнения Ляпунова:

Q = C'*C;
Q_bar = dlyap((A-B*K_lqr)', Q+K_lqr'*R*K_lqr);

Преобразуйте проблему MPC в стандартную проблему QP, которая имеет целевую функцию J(k) = U(k)'*H*U(k) + 2*x(k)'*F'*U(k).

Q_hat = blkdiag(Q,Q,Q,Q_bar);
R_hat = blkdiag(R,R,R,R);
H = CONV'*Q_hat*CONV + R_hat;
F = CONV'*Q_hat*M;

Когда нет никаких ограничений, оптимальной предсказанной входной последовательностью U (k) сгенерированный контроллером MPC является -K*x, где K = inv(H)*F.

K = H\F;

На практике, только первое перемещение управления u(k) = -K_mpc*x(k) применяется к объекту (отступающий управление горизонтом).

K_mpc = K(1,:);

Запустите симуляцию с начальными состояниями в [0.5 - 0.5]. Ответ с обратной связью устойчив.

Unconstrained_MPC = tf([-1 1])*feedback(ss(A,B,eye(2),0,Ts),K_mpc);
lsim(Unconstrained_MPC,'*',u_unconstrained,t_unconstrained,x0)
legend show

LQR и контроллеры MPC приводят к тому же результату, потому что законы о надзоре являются тем же самым.

K_lqr
K_mpc

K_lqr =

    4.3608   18.7401


K_mpc =

    4.3608   18.7401

Эффективность управления LQR ухудшается когда применение ограничений

Ограничьте контроллер выход, u (k), чтобы быть между-1 и 1. Контроллер LQR генерирует медленный и колебательный ответ с обратной связью из-за насыщения.

x = x0;
t_constrained = 0:40;
for ct = t_constrained
    uLQR(ct+1) = -K_lqr*x;
    uLQR(ct+1) = max(-1,min(1,uLQR(ct+1)));
    x = A*x+B*uLQR(ct+1);
    yLQR(ct+1) = C*x;
end
figure
subplot(2,1,1)
plot(t_constrained,uLQR)
xlabel('time')
ylabel('u')
subplot(2,1,2)
plot(t_constrained,yLQR)
xlabel('time')
ylabel('y')
legend('Constrained LQR')

Контроллер MPC решает задачу QP онлайн когда применение ограничений

Одно из главных преимуществ использования контроллера MPC - то, что оно обрабатывает ограничения ввода и вывода явным образом путем решения задачи оптимизации в каждом контрольном интервале.

Используйте встроенный решатель KWIK QP, mpcActiveSetSolver, реализовывать пользовательский контроллер MPC, спроектированный выше. Ограничительные матрицы заданы как Ac*x> =b0.

Ac = [1 0 0 0;...
      -1 0 0 0;...
       0 1 0 0;...
       0 -1 0 0;...
       0 0 1 0;...
       0 0 -1 0;...
       0 0 0 1;...
       0 0 0 -1];
b0 = [1;1;1;1;1;1;1;1];

Поскольку в этом случае матрица H Гессиана является постоянной, можно предварительно вычислить инверсию ее нижнего треугольного разложения Холесского, и затем передать ее mpcActiveSetSolver функция, вместо того, чтобы передать матрицу Гессиана непосредственно. В результате mpcActiveSetSolver может постараться не выполнять этот расчет на каждом временном шаге.

L = chol(H,'lower');
Linv = L\eye(size(H,1));

Запустите симуляцию путем вызова mpcActiveSetSolver в каждом шаге симуляции. Первоначально все неравенства неактивны (холодный запуск).

x = x0;
iA = false(size(b0));

% create options for the solver, and specify non-hessian first input
opt = mpcActiveSetOptions;
opt.IntegrityChecks = false;
opt.UseHessianAsInput = false;

for ct = t_constrained
    [u,status,iA] = mpcActiveSetSolver(Linv,F*x,Ac,b0,[],zeros(0,1),iA,opt);
    uMPC(ct+1) = u(1);
    x = A*x+B*uMPC(ct+1);
    yMPC(ct+1) = C*x;
end
figure
subplot(2,1,1)
plot(t_constrained,uMPC)
xlabel('time')
ylabel('u')
subplot(2,1,2)
plot(t_constrained,yMPC)
xlabel('time')
ylabel('y')
legend('Constrained MPC')

Контроллер MPC производит ответ с обратной связью с более быстрым временем урегулирования и меньше колебания.

Симулируйте пользовательский MPC Используя блок MATLAB function в Simulink

mpcActiveSetSolver может использоваться в блоке MATLAB function, чтобы обеспечить симуляцию и генерацию кода в окружении Simulink.

mdl = 'mpc_activesetqp';
open_system(mdl)

Блоком Custom MPC Controller является блок MATLAB function. Чтобы исследовать код MATLAB, дважды кликните блок. Начиная с LinvF, Ac, b0 матрицы и opt структура является постоянной, они передаются в блок MATLAB function как параметры.

Запустите симуляцию в Simulink. Закрытые ответы LQR и контроллеров MPC идентичны своим дубликатам в симуляции MATLAB.

open_system([mdl '/u_lqr'])
open_system([mdl '/y_lqr'])
open_system([mdl '/u_mpc'])
open_system([mdl '/y_mpc'])
sim(mdl)

Генерация кода в MATLAB

mpcActiveSetSolver генерация кода C поддержек с MATLAB Coder. Примите, что у вас есть функция, mycode, это совместимо со стандартом генерации кода.

function [x,iter,iA1,lam] = mycode()
%#codegen
n = 5;
m = 10;
q = 2;
H = diag(10*rand(n,1));
f = randn(n,1);
A = randn(m,n);
b = randn(m,1);
Aeq = randn(q,n);
beq = randn(q,1);
Linv = chol(H,'lower')\eye(n);
iA = false(m,1);
Opt = mpcActiveSetOptions();
[x,iter,iA1,lam] = mpcActiveSetSolver(Linv,f,A,b,Aeq,beq,iA,Opt);

Можно использовать следующую команду, чтобы сгенерировать код С с MATLAB Coder:

fun = 'mycode';
Cfg = coder.config('mex'); % or 'lib', 'dll', etc.
codegen('-config',Cfg,fun,'-o',fun);

Подтверждение

Этот пример вдохновлен примечаниями лекции профессора Марка Кэннона для Прогнозирующего класса Управления Модели в Оксфордском университете. Модель объекта управления - тот же используемый в Примере 2.1 в разделе "Prediction and optimization".

bdclose(mdl)