mathworks社群中的這個資料還是值得一說的。

1 openExample('mpc/mpccustomqp')

我們從幾個角度來解析兩者關係，簡單的說就是MPC是帶了約束的LQR.

下面我們從程式碼的角度解析這個問題：
1， 定義被控系統：

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

2，設計無約束LQR：

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

3, 執行模擬閉環結果：

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

4，設計MPC控制器：

 1 %%
 2 % The MPC objective function is |J(k) = sum(x(k)'
 
*Q*x(k) + u(k)'*R*u(k) +
 3 % x(k+N)'*Q_bar*x(k+N))|. To ensure that the MPC objective function has the
 4 % same quadratic cost as the infinite horizon quadratic cost used by LQR,
 5 % terminal weight |Q_bar| is obtained by solving the following Lyapunov
 6 % equation:
 7 Q = C'*C;
 8 Q_bar = dlyap((A-B*K_lqr)'
 
, Q+K_lqr'*R*K_lqr); 
 9 
10 %%
11 % Convert the MPC problem into a standard QP problem, which has the
12 % objective function |J(k) = U(k)'*H*U(k) + 2*x(k)'*F'*U(k)|.
13 Q_hat = blkdiag(Q,Q,Q,Q_bar);
14 R_hat = blkdiag(R,R,R,R);
15 H = CONV'*Q_hat*CONV + R_hat;
16 F = CONV'*Q_hat*M;
17 
18 %%
19 % When there are no constraints, the optimal predicted input sequence U(k)
20 % generated by MPC controller is |-K*x|, where |K = inv(H)*F|.
21 K = H\F;
22 
23 %% 
24 % In practice, only the first control move |u(k) = -K_mpc*x(k)| is applied
25 % to the plant (receding horizon control).
26 K_mpc = K(1,:);
27 
28 %% 
29 % Run a simulation with initial states at [0.5 -0.5]. The closed-loop
30 % response is stable.  
31 Unconstrained_MPC = tf([-1 1])*feedback(ss(A,B,eye(2),0,Ts),K_mpc);
32 lsim(Unconstrained_MPC,'*',u_unconstrained,t_unconstrained,x0)
33 legend show

到這裡，完全可以說明，在無約束前提下，兩種方法是一致的：

1 K_lqr =
2 
3     4.3608   18.7401
4 
5 
6 K_mpc =
7 
8     4.3608   18.7401

5，對LQR施加約束：

 1 x = x0;
 2 t_constrained = 0:40;
 3 for ct = t_constrained
 4 uLQR(ct+1) = -K_lqrx;
 5 uLQR(ct+1) = max(-1,min(1,uLQR(ct+1)));
 6 x = Ax+BuLQR(ct+1);
 7 yLQR(ct+1) = Cx;
 8 end
 9 figure
10 subplot(2,1,1)
11 plot(t_constrained,uLQR)
12 xlabel(‘time’)
13 ylabel(‘u’)
14 subplot(2,1,2)
15 plot(t_constrained,yLQR)
16 xlabel(‘time’)
17 ylabel(‘y’)
18 legend(‘Constrained LQR’)

6，對MPC施加約束：

 1 %% MPC Controller Solves QP Problem Online When Applying Constraints 
 2 % One of the major benefits of using MPC controller is that it handles
 3 % input and output constraints explicitly by solving an optimization
 4 % problem at each control interval.
 5 %
 6 % Use the built-in KWIK QP solver, |mpcqpsolver|, to implement the custom
 7 % MPC controller designed above. The constraint matrices are defined as
 8 % Ac*x>=b0.
 9 Ac = -[1 0 0 0;...
10       -1 0 0 0;...
11        0 1 0 0;...
12        0 -1 0 0;...
13        0 0 1 0;...
14        0 0 -1 0;...
15        0 0 0 1;...
16        0 0 0 -1];
17 b0 = -[1;1;1;1;1;1;1;1];
18 
19 %% 
20 % The |mpcqpsolver| function requires the first input to be the inverse of
21 % the lower-triangular Cholesky decomposition of the Hessian matrix H.
22 L = chol(H,'lower');
23 Linv = L\eye(size(H,1));
24 
25 %%
26 % Run a simulation by calling |mpcqpsolver| at each simulation step.
27 % Initially all the inequalities are inactive (cold start). 
28 x = x0;
29 iA = false(size(b0));
30 opt = mpcqpsolverOptions;
31 opt.IntegrityChecks = false;
32 for ct = t_constrained
33     [u, status, iA] = mpcqpsolver(Linv,F*x,Ac,b0,[],zeros(0,1),iA,opt);
34     uMPC(ct+1) = u(1);
35     x = A*x+B*uMPC(ct+1);
36     yMPC(ct+1) = C*x;
37 end
38 figure
39 subplot(2,1,1)
40 plot(t_constrained,uMPC)
41 xlabel('time')
42 ylabel('u')
43 subplot(2,1,2)
44 plot(t_constrained,yMPC)
45 xlabel('time')
46 ylabel('y')
47 legend('Constrained MPC')

轉載：https://blog.csdn.net/gophae/article/details/104546805/