В этом примере показано, как использовать встроенный решатель QP активного набора для реализации пользовательского алгоритма MPC, поддерживающего генерацию кода C в MATLAB.
Модель установки представляет собой дискретно-временную систему состояния-пространства и является нестабильной с разомкнутым контуром. Мы предполагаем, что все состояния растений поддаются измерению. Поэтому мы избегаем необходимости проектирования оценщика состояния, что выходит за рамки данного примера.
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 с выходным весовым коэффициентом. Этот контроллер служит базовой линией для сравнения с пользовательским алгоритмом 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-контроллера с использованием веса терминала на последнем шаге прогнозирования.
Предсказанные последовательности состояний 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
Ограничьте выход контроллера 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-контроллера является то, что он явно обрабатывает ограничения ввода и вывода, решая задачу оптимизации на каждом интервале управления.
Используйте встроенный решатель 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 производит отклик по замкнутому контуру с более быстрым временем оседания и меньшим колебанием.
mpcActiveSetSolver может использоваться внутри функционального блока MATLAB для моделирования и генерации кода в среде Simulink.
mdl = 'mpc_activesetqp';
open_system(mdl)
Пользовательский блок контроллера MPC является функциональным блоком MATLAB. Для проверки кода MATLAB дважды щелкните блок. С тех пор Linv, F, 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)
mpcActiveSetSolver поддерживает генерацию кода C с кодером MATLAB. Предположим, что у вас есть функция, 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);
Для создания кода C с помощью кодера MATLAB можно использовать следующую команду:
fun = 'mycode'; Cfg = coder.config('mex'); % or 'lib', 'dll', etc. codegen('-config',Cfg,fun,'-o',fun);
Этот пример вдохновлен записями лекций профессора Марка Кэннона для класса модельного предиктивного контроля в Оксфордском университете. Модель установки аналогична модели, использованной в примере 2.1 в разделе «Прогнозирование и оптимизация».
bdclose(mdl)
mpcqpsolver | mpcqpsolverOptions