【控制】模型预测控制，公式推导，数值仿真，有程序有图

1 从连续时间模型到增广模型

1.1 从连续时间模型到离散模型

连续时间 (continuous time) 模型的状态空间表达式为

x ˙ m ( t ) = A c x m ( t ) + B c u ( t ) y m ( t ) = C c x m ( t ) (1) ˙xm(t)=Acxm(t)+Bcu(t)ym(t)=Ccxm(t)

\tag{1} x˙m​(t)ym​(t)​=Ac​xm​(t)+Bc​u(t)=Cc​xm​(t)​(1)

离散时间 (discrete time) 模型的状态空间表达式为
x m ( k + 1 ) = A d x m ( k ) + B d u ( k ) y m ( k ) = C d x m ( k ) (2) xm(k+1)=Adxm(k)+Bdu(k)ym(k)=Cdxm(k)

\tag{2} xm​(k+1)ym​(k)​=Ad​xm​(k)+Bd​u(k)=Cd​xm​(k)​(2)

在Matlab中可以通过 c2dm() 函数，同时指定一个采样时间间隔 $dT$，将连续时间模型转换成离散模型。

数值例子

给定连续时间状态空间模型如下，求出对应的离散时间状态空间模型，采样时间间隔为1。

x ˙ m ( t ) = [ 0 1 0 3 0 1 0 1 0 ] x m ( t ) + [ 1 1 3 ] u ( t ) y m ( t ) = [ 0 1 0 ] x m ( t ) ˙xm(t)=[010301010]xm(t)+[113]u(t)ym(t)=[010]xm(t)

x˙m​(t)ym​(t)​=⎣ ⎡​030​101​010​⎦ ⎤​xm​(t)+⎣ ⎡​113​⎦ ⎤​u(t)=[0​1​0​]xm​(t)​

Matlab程序为

A_c = [0 1 0; 3 0 1; 0 1 0 ];
B_c = [1; 1; 3];
C_c = [0 1 0];
D_c = zeros(1,1);
dT = 1;
[A_d, B_d, C_d, D_d] = c2dm(A_c, B_c, C_c, D_c, dT)
1
2
3
4
5
6

得到的结果如下

A_d =
    3.0716    1.8134    0.6905
    5.4403    3.7622    1.8134
    2.0716    1.8134    1.6905

B_d =
    2.9107
    5.9567
    4.9107

C_d =
     0     1     0
   
D_d =
     0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

---

1.2 从离散模型到增广模型

令
$\Delta x_m(k+1)=x_m(k+1)-x_m(k)$，
$\Delta x_m(k)=x_m(k)-x_m(k-1)$，
$\Delta u(k)=u(k)-u(k-1)$，
$\Delta y_m(k+1)=y_m(k+1)-y_m(k)$。
于是，有如下增广模型

[ Δ x m ( k + 1 ) y m ( k + 1 ) ] ⏟ x ( k + 1 ) = [ A d 0 C d A d I ] ⏟ A [ Δ x m ( k ) y m ( k ) ] ⏟ x ( k ) + [ B d C d B d ] ⏟ B Δ u ( k ) y m ( k ) = [ 0 I ] ⏟ A [ Δ x m ( k ) y m ( k ) ] ⏟ x ( k ) (3) [Δxm(k+1)ym(k+1)]⏟x(k+1)=[Ad0CdAdI]⏟A[Δxm(k)ym(k)]⏟x(k)+[BdCdBd]⏟BΔu(k)ym(k)=[0I]⏟A[Δxm(k)ym(k)]⏟x(k)

\tag{3} x(k+1) [Δxm​(k+1)ym​(k+1)​]​​ym​(k)​=A [Ad​Cd​Ad​​0I​]​​x(k) [Δxm​(k)ym​(k)​]​​+B [Bd​Cd​Bd​​]​​Δu(k)=A [0​I​]​​x(k) [Δxm​(k)ym​(k)​]​​​(3)

其中 $0,I$ 分别是适当维度的全零矩阵和单位矩阵。

根据式 (3)，我们就有了如下的增广矩阵
x ( k + 1 ) = A x ( k ) + B Δ u ( k ) y ( k ) = C x ( k ) (4) x(k+1)=Ax(k)+BΔu(k)y(k)=Cx(k)

\tag{4} x(k+1)y(k)​=Ax(k)+BΔu(k)=Cx(k)​(4)

其中 $x,y,\Delta u$ 和 $A,B,C$ 的定义参考式 (3)。

数值例子

给定离散时间状态空间模型如下，

x m ( k + 1 ) = [ 1 1 0 1 ] x m ( k ) + [ 0.5 1 ] u ( k ) y m ( k ) = [ 1 0 ] x m ( k ) xm(k+1)=[1101]xm(k)+[0.51]u(k)ym(k)=[10]xm(k)

xm​(k+1)ym​(k)​=[10​11​]xm​(k)+[0.51​]u(k)=[1​0​]xm​(k)​

求对应的增广模型矩阵。

Matlab程序为

A_d = [1 1;0 1];
B_d = [0.5;1];
C_d = [1 0];

[BRow, BCol] = size(B_d);
[CRow, CCol] = size(C_d);

A = eye(BRow+CRow, BCol+CCol);
A(1:BRow, 1:CCol) = A_d;
A(BRow+1:BRow+CRow, 1:CCol) = C_d * A_d

B = zeros(BRow+CRow, BCol);
B(1:BRow, 1:BCol) = B_d;
B(BRow+1:BRow+CRow, 1:BCol) = C_d*B_d

C = zeros(CRow, BRow + CRow);
C(CRow, CCol+1:BRow+CRow) = eye(CRow)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

得到的结果如下

A =
     1     1     0
     0     1     0
     1     1     1

B =
    0.5000
    1.0000
    0.5000

C =
     0     0     1
1
2
3
4
5
6
7
8
9
10
11
12

---

通过上述转换为增广矩阵的过程，我们对增广矩阵式 (4) 中变量 $x(k)$ 的研究，就等价于对离散时间模型式 (2) 中变量 $x_m(k)$ 的下一时刻与当前时刻差值的研究。

2 状态预测与优化

2.1 状态预测

假设当前时刻为 $k_i$，$x(k_i)$ 为通过检测得到的系统状态变量，这里假设状态变量可以检测到。那么 $N_u$ 个未来的控制变量序列表示为

$$\Delta u(k_i),\Delta u(k_i+1),\cdots ,\Delta u(k_i+N_u-1)\text{\hspace{1em}(5)}$$

$N_x$ 个预测的状态变量表示为
$$x(k_i+1|k_i),x(k_i+2|k_i),\cdots ,x(k_i+m|k_i),\cdots ,x(k_i+N_x|x_i)\text{\hspace{1em}(6)}$$

注意，这里仅使用了当前时刻的状态值 $x(k_i)$ 就预测了之后的 $N_p$ 个时刻

代入增广模型中，得到 $N_x$ 个预测状态变量的表达式为
x ( k i + 1 ∣ k i ) = A x ( k i ) + B Δ u ( k i ) x ( k i + 2 ∣ k i ) = A x ( k i + 1 ∣ k i ) + B Δ u ( k i + 1 ) = A ( A x ( k i ) + B Δ u ( k i ) ) + B Δ u ( k i + 1 ) = A 2 x ( k i ) + A B Δ u ( k i ) + B Δ u ( k i + 1 ) ⋮ x ( k i + N x ∣ k i ) = A N x x ( k i ) + A N x − 1 B Δ u ( k i ) + A N x − 2 B Δ u ( k i + 1 ) + ⋯ + A N x − N u B Δ u ( k i + N u − 1 ) (7) x(ki+1|ki)=Ax(ki)+BΔu(ki)x(ki+2|ki)=Ax(ki+1|ki)+BΔu(ki+1)=A(Ax(ki)+BΔu(ki))+BΔu(ki+1)=A2x(ki)+ABΔu(ki)+BΔu(ki+1)⋮x(ki+Nx|ki)=ANxx(ki)+ANx−1BΔu(ki)+ANx−2BΔu(ki+1)+⋯+ANx−NuBΔu(ki+Nu−1)

\tag{7} x(ki​+1∣ki​)x(ki​+2∣ki​)x(ki​+Nx​∣ki​)​=Ax(ki​)+BΔu(ki​)=Ax(ki​+1∣ki​)+BΔu(ki​+1)=A(Ax(ki​)+BΔu(ki​))+BΔu(ki​+1)=A2x(ki​)+ABΔu(ki​)+BΔu(ki​+1)⋮=ANx​x(ki​)+ANx​−1BΔu(ki​)+ANx​−2BΔu(ki​+1)+⋯+ANx​−Nu​BΔu(ki​+Nu​−1)​(7)

得到 $N_x$ 个预测输出变量的表达式为
y ( k i + 1 ∣ k i ) = C A x ( k i ) + C B Δ u ( k i ) y ( k i + 2 ∣ k i ) = C A x ( k i + 1 ∣ k i ) + C B Δ u ( k i + 1 ) = C A ( A x ( k i ) + B Δ u ( k i ) ) + C B Δ u ( k i + 1 ) = C A 2 x ( k i ) + C A B Δ u ( k i ) + C B Δ u ( k i + 1 ) ⋮ y ( k i + N x ∣ k i ) = C A N x x ( k i ) + C A N x − 1 B Δ u ( k i ) + C A N x − 2 B Δ u ( k i + 1 ) + ⋯ + C A N x − N u B Δ u ( k i + N u − 1 ) (8) y(ki+1|ki)=CAx(ki)+CBΔu(ki)y(ki+2|ki)=CAx(ki+1|ki)+CBΔu(ki+1)=CA(Ax(ki)+BΔu(ki))+CBΔu(ki+1)=CA2x(ki)+CABΔu(ki)+CBΔu(ki+1)⋮y(ki+Nx|ki)=CANxx(ki)+CANx−1BΔu(ki)+CANx−2BΔu(ki+1)+⋯+CANx−NuBΔu(ki+Nu−1)

\tag{8} y(ki​+1∣ki​)y(ki​+2∣ki​)y(ki​+Nx​∣ki​)​=CAx(ki​)+CBΔu(ki​)=CAx(ki​+1∣ki​)+CBΔu(ki​+1)=CA(Ax(ki​)+BΔu(ki​))+CBΔu(ki​+1)=CA2x(ki​)+CABΔu(ki​)+CBΔu(ki​+1)⋮=CANx​x(ki​)+CANx​−1BΔu(ki​)+CANx​−2BΔu(ki​+1)+⋯+CANx​−Nu​BΔu(ki​+Nu​−1)​(8)

可以看出，所有预测的变量都只与当前观测到的系统状态 $x(k_i)$ 和未来的变量序列 $\Delta u(k_i+j),j=1,2,\cdots ,N_u-1$ 有关。

将未来的输出变量序列 $y$ 的和未来的控制变量序列 $\Delta u$ 写成向量的形式

Y = [ y ( k i + 1 ∣ k i ) y ( k i + 2 ∣ k i ) ⋮ y ( k i + N x ∣ k i ) ] N x × 1 ,      Δ U = [ Δ u ( k i ) Δ u ( k i + 1 ) ⋮ Δ u ( k i + N u − 1 ) ] N u × 1 (9) Y = \left[y(ki+1|ki)y(ki+2|ki)⋮y(ki+Nx|ki)

\right]_{N_x \times 1}, ~~~~ \Delta U = \left[Δu(ki)Δu(ki+1)⋮Δu(ki+Nu−1)\right]_{N_u \times 1} \tag{9} Y=⎣ ⎡​y(ki​+1∣ki​)y(ki​+2∣ki​)⋮y(ki​+Nx​∣ki​)​⎦ ⎤​Nx​×1​,    ΔU=⎣ ⎡​Δu(ki​)Δu(ki​+1)⋮Δu(ki​+Nu​−1)​⎦ ⎤​Nu​×1​(9)

在单输入单输出 SISO 系统中，$Y$ 的维度是 $N_x\times 1$，$\Delta U$ 的维度是 $N_u\times 1$。

合并式 (8) 和 式 (9)，我们有
[ y ( k i + 1 ∣ k i ) y ( k i + 2 ∣ k i ) ⋮ y ( k i + N x ∣ k i ) ] = [ C A C A 2 ⋮ C A N x ] x ( k i ) + [ C B 0 0 ⋯ 0 C A B C B 0 ⋯ 0 C A 2 B C A B C B ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ C A N x − 1 B C A N x − 2 B C A N x − 3 B ⋯ C A N x − N u B ] [ Δ u ( k i ) Δ u ( k i + 1 ) ⋮ Δ u ( k i + N u − 1 ) ] Y = F x ( k i ) + Φ Δ U (10) [y(ki+1|ki)y(ki+2|ki)⋮y(ki+Nx|ki)]=[CACA2⋮CANx]x(ki)+[CB00⋯0CABCB0⋯0CA2BCABCB⋯0⋮⋮⋮⋱⋮CANx−1BCANx−2BCANx−3B⋯CANx−NuB][Δu(ki)Δu(ki+1)⋮Δu(ki+Nu−1)]Y=Fx(ki)+ΦΔU

\tag{10} ⎣ ⎡​y(ki​+1∣ki​)y(ki​+2∣ki​)⋮y(ki​+Nx​∣ki​)​⎦ ⎤​Y​=⎣ ⎡​CACA2⋮CANx​​⎦ ⎤​x(ki​)+⎣ ⎡​CBCABCA2B⋮CANx​−1B​0CBCAB⋮CANx​−2B​00CB⋮CANx​−3B​⋯⋯⋯⋱⋯​000⋮CANx​−Nu​B​⎦ ⎤​⎣ ⎡​Δu(ki​)Δu(ki​+1)⋮Δu(ki​+Nu​−1)​⎦ ⎤​=Fx(ki​)+ΦΔU​(10)

2.2 优化

在 $k_i$ 时刻的期望输入值为 $r(k_i)$，需要通过MPC控制，将系统输出尽可能的接近系统输入。假设在优化窗口内，输入是恒定值。问题就转换成了需要在优化窗口内，找到一个最佳的控制变量序列，使得系统输出与输入的误差最小。

将输入 $r(k_i)$ 写成向量的形式有
R s = r ( k i ) [ 1 1 ⋮ 1 ] N x × 1 (11) R_s = r(k_i) \left[11⋮1

\right]_{N_x \times 1} \tag{11} Rs​=r(ki​)⎣ ⎡​11⋮1​⎦ ⎤​Nx​×1​(11)

定义代价函数为
$$J=(R_s-Y{)}^{\text{T}}(R_s-Y)+\Delta U^{\text{T}}\bar{R}\Delta U\text{\hspace{1em}(12)}$$

其中 $\bar{R}$ 是权重矩阵。上述代价函数是一个典型的二次型代价函数，第一项反映了输入信号 $R_s$ 与预测的输出信号 $Y$ 之间的误差，第二项则是为了限制控制量 $\Delta U$，避免控制量过大。

关于权重矩阵 $\bar{R}$ 的选择，这里我们使用
$$\bar{R}=r_w{I}_{N_u\times N_u},r_w\ge 0\text{\hspace{1em}(13)}$$

其中 $r_w$ 是影响闭环性能的可调参数。如果 $r_w=0$，则说明忽略 $\Delta U$ 的大小，只考虑控制误差。

根据上述分析过程，代价函数 $J$ 可以写成
J = ( R s − Y ) T ( R s − Y ) + Δ U T R ˉ Δ U = [ ( R s − F x ( k i ) ) − Φ Δ U ] T [ ( R s − F x ( k i ) ) − Φ Δ U ] + Δ U T R ˉ Δ U = [ ( R s − F x ( k i ) ) T − Δ U T Φ T ] [ ( R s − F x ( k i ) ) − Φ Δ U ] + Δ U T R ˉ Δ U = ( R s − F x ( k i ) ) T ( R s − F x ( k i ) ) − ( R s − F x ( k i ) ) T Φ Δ U − Δ U T Φ T ( R s − F x ( k i ) ) + Δ U T Φ T Φ Δ U + Δ U T R ˉ Δ U = ( R s − F x ( k i ) ) T ( R s − F x ( k i ) ) − Δ U T Φ T ( R s − F x ( k i ) ) − Δ U T Φ T ( R s − F x ( k i ) ) + Δ U T Φ T Φ Δ U + Δ U T R ˉ Δ U = ( R s − F x ( k i ) ) T ( R s − F x ( k i ) ) − 2 Δ U T Φ T ( R s − F x ( k i ) ) + Δ U T ( Φ T Φ + R ˉ ) Δ U (14) J=(Rs−Y)T(Rs−Y)+ΔUTˉRΔU=[(Rs−Fx(ki))−ΦΔU]T[(Rs−Fx(ki))−ΦΔU]+ΔUTˉRΔU=[(Rs−Fx(ki))T−ΔUTΦT][(Rs−Fx(ki))−ΦΔU]+ΔUTˉRΔU=(Rs−Fx(ki))T(Rs−Fx(ki))−(Rs−Fx(ki))TΦΔU−ΔUTΦT(Rs−Fx(ki))+ΔUTΦTΦΔU+ΔUTˉRΔU=(Rs−Fx(ki))T(Rs−Fx(ki))−ΔUTΦT(Rs−Fx(ki))−ΔUTΦT(Rs−Fx(ki))+ΔUTΦTΦΔU+ΔUTˉRΔU=(Rs−Fx(ki))T(Rs−Fx(ki))−2ΔUTΦT(Rs−Fx(ki))+ΔUT(ΦTΦ+ˉR)ΔU

\tag{14} J​=(Rs​−Y)T(Rs​−Y)+ΔUTRˉΔU=[(Rs​−Fx(ki​))−ΦΔU]T[(Rs​−Fx(ki​))−ΦΔU]+ΔUTRˉΔU=[(Rs​−Fx(ki​))T−ΔUTΦT][(Rs​−Fx(ki​))−ΦΔU]+ΔUTRˉΔU=(Rs​−Fx(ki​))T(Rs​−Fx(ki​))−(Rs​−Fx(ki​))TΦΔU−ΔUTΦT(Rs​−Fx(ki​))+ΔUTΦTΦΔU+ΔUTRˉΔU=(Rs​−Fx(ki​))T(Rs​−Fx(ki​))−ΔUTΦT(Rs​−Fx(ki​))−ΔUTΦT(Rs​−Fx(ki​))+ΔUTΦTΦΔU+ΔUTRˉΔU=(Rs​−Fx(ki​))T(Rs​−Fx(ki​))−2ΔUTΦT(Rs​−Fx(ki​))+ΔUT(ΦTΦ+Rˉ)ΔU​(14)

对 $J$ 关于控制量 $\Delta U$ 求一阶偏导数
∂ J ∂ Δ U = − 2 Φ T ( R s − F x ( k i ) ) + 2 ( Φ T Φ + R ˉ ) Δ U (15) ∂J∂ΔU=−2ΦT(Rs−Fx(ki))+2(ΦTΦ+ˉR)ΔU

\tag{15} ∂ΔU∂J​​=−2ΦT(Rs​−Fx(ki​))+2(ΦTΦ+Rˉ)ΔU​(15)

当一阶导数 $\frac{\partial J}{\partial \Delta U}=0$ 时，有最优的控制量为
− 2 Φ T ( R s − F x ( k i ) ) + 2 ( Φ T Φ + R ˉ ) Δ U = 0 ( Φ T Φ + R ˉ ) Δ U = Φ T ( R s − F x ( k i ) ) Δ U = ( Φ T Φ + R ˉ ) − 1   Φ T ( R s − F x ( k i ) ) (16) −2ΦT(Rs−Fx(ki))+2(ΦTΦ+ˉR)ΔU=0(ΦTΦ+ˉR)ΔU=ΦT(Rs−Fx(ki))ΔU=(ΦTΦ+ˉR)−1 ΦT(Rs−Fx(ki))

\tag{16} −2ΦT(Rs​−Fx(ki​))+2(ΦTΦ+Rˉ)ΔU(ΦTΦ+Rˉ)ΔUΔU​=0=ΦT(Rs​−Fx(ki​))=(ΦTΦ+Rˉ)−1 ΦT(Rs​−Fx(ki​))​(16)

这里我们用了一个假设，也就是 $({\Phi }^{\text{T}}\Phi +\bar{R})$ 的逆一定存在。

换种形式表示最优控制量 $\Delta U^{\ast }$ 还可以为
Δ U ∗ = ( Φ T Φ + R ˉ ) − 1   Φ T ( R s − F x ( k i ) ) = ( Φ T Φ + R ˉ ) − 1   Φ T ( r ( k i ) [ 1 1 ⋮ 1 ] − F x ( k i ) ) (17) ΔU∗=(ΦTΦ+ˉR)−1 ΦT(Rs−Fx(ki))=(ΦTΦ+ˉR)−1 ΦT(r(ki)[11⋮1]−Fx(ki))

\tag{17} ΔU∗​=(ΦTΦ+Rˉ)−1 ΦT(Rs​−Fx(ki​))=(ΦTΦ+Rˉ)−1 ΦT⎝ ⎛​r(ki​)⎣ ⎡​11⋮1​⎦ ⎤​−Fx(ki​)⎠ ⎞​​(17)

数值例子

给定离散时间状态空间模型如下，

x m ( k + 1 ) = a x m ( k ) + b u ( k ) y m ( k ) = c x m ( k ) xm(k+1)=axm(k)+bu(k)ym(k)=cxm(k)

xm​(k+1)ym​(k)​=axm​(k)+bu(k)=cxm​(k)​

其中 $a=0.8,b=0.1,c=1$，求
(1) 求对应的增广模型矩阵，
(2) 计算矩阵 $F,\Phi$ 和 ${\Phi }^{\text{T}}\Phi ,{\Phi }^{\text{T}}F$，
(3) 假设在 $k_i=10$ 时刻有 $r(k_i)=1$， x ( k i ) = [ 0.1 0.2 ] T x(k_i) = \left[0.10.2

\right]^\text{T} x(ki​)=[0.1​0.2​]T，计算当权重矩阵参数分别 $r_w=0$ 和 $r_w=10$ 的最优控制解 $\Delta U^{\ast }$，并比较效果。

(1) 增广状态空间方程为
[ Δ x m ( k + 1 ) y m ( k + 1 ) ] = [ a 0 c a 1 ] [ Δ x m ( k ) y m ( k ) ] + [ b c b ] Δ u ( k ) y m ( k ) = [ 0 1 ] [ Δ x m ( k ) y m ( k ) ] [Δxm(k+1)ym(k+1)]=[a0ca1][Δxm(k)ym(k)]+[bcb]Δu(k)ym(k)=[01][Δxm(k)ym(k)]

[Δxm​(k+1)ym​(k+1)​]ym​(k)​=[aca​01​][Δxm​(k)ym​(k)​]+[bcb​]Δu(k)=[0​1​][Δxm​(k)ym​(k)​]​

其中增广矩阵分别为 A = [ a 0 c a 1 ] = [ 0.8 0 0.8 1 ] A=\left[a0ca1

\right]=\left[0.800.81\right] A=[aca​01​]=[0.80.8​01​]，$B=\left[\begin{array}{c}b\\ cb\end{array}\right]=\left[\begin{array}{c}0.1\\ 0.1\end{array}\right]$，$C=\left[\begin{array}{cc}0& 1\end{array}\right]$。

(2) 矩阵 $F$ 和 $\Phi$ 和下，
F = [ C A C A 2 C A 3 C A 4 C A 5 C A 6 C A 7 C A 8 C A 9 C A 10 ] ,      Φ = [ C B 0 0 0 C A B C B 0 0 C A 2 B C A B C B 0 C A 3 B C A 2 B C A B C B C A 4 B C A 3 B C A 2 B C A B C A 5 B C A 4 B C A 3 B C A 2 B C A 6 B C A 5 B C A 4 B C A 3 B C A 7 B C A 6 B C A 5 B C A 4 B C A 8 B C A 7 B C A 6 B C A 5 B C A 9 B C A 8 B C A 7 B C A 6 B ] F=\left[CACA2CA3CA4CA5CA6CA7CA8CA9CA10

\right], ~~~~ \Phi=\left[CB000CABCB00CA2BCABCB0CA3BCA2BCABCBCA4BCA3BCA2BCABCA5BCA4BCA3BCA2BCA6BCA5BCA4BCA3BCA7BCA6BCA5BCA4BCA8BCA7BCA6BCA5BCA9BCA8BCA7BCA6B\right] F=⎣ ⎡​CACA2CA3CA4CA5CA6CA7CA8CA9CA10​⎦ ⎤​,    Φ=⎣ ⎡​CBCABCA2BCA3BCA4BCA5BCA6BCA7BCA8BCA9B​0CBCABCA2BCA3BCA4BCA5BCA6BCA7BCA8B​00CBCABCA2BCA3BCA4BCA5BCA6BCA7B​000CBCABCA2BCA3BCA4BCA5BCA6B​⎦ ⎤​

这里我们预测了 $k_i$ 时刻之后10个状态变量和未来4个控制变量。矩阵 $F$ 中的系数可以按如下方式计算
C A = [ s 1 1 ] = [ a 1 ] C A 2 = [ s 2 1 ] = [ a 2 + s 1 1 ] C A 3 = [ s 3 1 ] = [ a 3 + s 2 1 ] ⋮ C A k = [ s k 1 ] = [ a k + s k − 1 1 ] CA=[s11]=[a1]CA2=[s21]=[a2+s11]CA3=[s31]=[a3+s21]⋮CAk=[sk1]=[ak+sk−11]

CACA2CA3CAk​=[s1​​1​]=[a​1​]=[s2​​1​]=[a2+s1​​1​]=[s3​​1​]=[a3+s2​​1​]⋮=[sk​​1​]=[ak+sk−1​​1​]​

矩阵 $\Phi$ 中的系数可以按如下方式计算
C B = g 0 = b C A B = g 1 = a b + g 0 C A 2 B = g 2 = a 2 b + g 1 C A 3 B = g 3 = a 3 b + g 2 ⋮ C A k B = g k = a k b + g k − 1 CB=g0=bCAB=g1=ab+g0CA2B=g2=a2b+g1CA3B=g3=a3b+g2⋮CAkB=gk=akb+gk−1

CBCABCA2BCA3BCAkB​=g0​=b=g1​=ab+g0​=g2​=a2b+g1​=g3​=a3b+g2​⋮=gk​=akb+gk−1​​

将 $a=0.8,b=0.1$ 代入后有
Φ T Φ = [ 1.1541 1.0407 0.9116 0.7726 1.0407 0.9549 0.8475 0.7259 0.9116 0.8475 0.7675 0.6674 0.7726 0.7259 0.6674 0.5943 ] ,      Φ T F = [ 9.2325 3.2147 8.3259 2.7684 7.2927 2.3355 6.1811 1.9194 ] ΦTΦ=[1.15411.04070.91160.77261.04070.95490.84750.72590.91160.84750.76750.66740.77260.72590.66740.5943],    ΦTF=[9.23253.21478.32592.76847.29272.33556.18111.9194]

ΦTΦ​=⎣ ⎡​1.15411.04070.91160.7726​1.04070.95490.84750.7259​0.91160.84750.76750.6674​0.77260.72590.66740.5943​⎦ ⎤​,    ΦTF​=⎣ ⎡​9.23258.32597.29276.1811​3.21472.76842.33551.9194​⎦ ⎤​​

(3) 根据之前的运算，我们知道最优控制为
Δ U ∗ = ( Φ T Φ + R ˉ ) − 1   Φ T ( R s − F x ( k i ) ) ΔU∗=(ΦTΦ+ˉR)−1 ΦT(Rs−Fx(ki))

ΔU∗​=(ΦTΦ+Rˉ)−1 ΦT(Rs​−Fx(ki​))​

当 $r_w=0$ 时，有 $\bar{R}=0$，计算得到的 $N_u=4$ 个控制量序列为
Δ U ∗ = [ 7.2 − 6.4 0 0 ] T \Delta U^* = \left[7.2−6.400

\right]^\text{T} ΔU∗=[7.2​−6.4​0​0​]T

当 $r_w=10$ 时，有 $\bar{R}=10I$，计算得到的 $N_u=4$ 个控制量序列为
Δ U ∗ = [ 0.1269 0.1034 0.0829 0.065 ] T \Delta U^* = \left[0.12690.10340.08290.065

\right]^\text{T} ΔU∗=[0.1269​0.1034​0.0829​0.065​]T

对比两次的控制量，可以看出来第一种情况控制量的变化较大，而第二种则平稳的多。

clear
clc

a = 0.8;
b = 0.1;
c = 1;
N_u = 4;
N_x = 10;

% discrete model
A_d = a;
B_d = b;
C_d = c;

% augmented model
[BRow, BCol] = size(B_d);
[CRow, CCol] = size(C_d);

A = eye(BRow+CRow, BCol+CCol);
A(1:BRow, 1:CCol) = A_d;
A(BRow+1:BRow+CRow, 1:CCol) = C_d * A_d;

B = zeros(BRow+CRow, BCol);
B(1:BRow, 1:BCol) = B_d;
B(BRow+1:BRow+CRow, 1:BCol) = C_d*B_d;

C = zeros(CRow, BRow + CRow);
C(CRow, CCol+1:BRow+CRow) = eye(CRow);

% matrix F
F = zeros(N_x, 2);
for i=1:N_x
    F(i,:) = C * A^i;
end

% matrix Phi
Phi = zeros(N_x, N_u);
for i=1:N_x
    for j=1:N_u
        if i < N_u && i