0、前言
动态矩阵控制(Dynamic Matrix Control,DMC)是一种典型的模型预测控制方法,其不需要被控对象的数学模型,只需要获取被控对象的阶跃响应序列即可实现控制效果,但其需要被控对象是渐近稳定的。
1、稳定 SISO 系统的阶跃响应模型
考虑单输入单输出(Signal Input Signal Output,SISO)系统,其传递函数为
G
(
s
)
=
y
(
s
)
u
(
s
)
(1)
G(s)=\frac{y(s)}{u(s)}\tag{1}
G(s)=u(s)y(s)(1)
首先考虑当输入不变时,系统的非零初始状态响应
如图所示(系统过渡过程时间为
N
N
N 个采样间隔)
系统在
k
+
N
−
2
k+N-2
k+N−2 以后进入稳态,输出保持不变,定义
k
−
1
k-1
k−1 时刻的未来
N
N
N 步输出为
Y
(
k
−
1
)
=
{
[
y
(
k
−
1
)
y
(
k
)
⋯
y
(
k
+
N
−
3
)
y
(
k
+
N
−
2
)
]
T
Δ
u
(
k
+
i
)
=
0
,
i
=
−
1
,
0
,
⋯
,
N
−
2
}
(2)
Y(k-1)=\Bigg\lbrace \begin{matrix} \left[ \begin{matrix} y(k-1) & y(k) & \cdots & y(k+N-3) & y(k+N-2)\\ \end{matrix} \right]^T\\[2ex] \Delta u(k+i)=0,\quad i=-1,0,\cdots,N-2 \end{matrix} \Bigg\rbrace\tag{2}
Y(k−1)={[y(k−1)y(k)⋯y(k+N−3)y(k+N−2)]TΔu(k+i)=0,i=−1,0,⋯,N−2}(2)
我们称 Y ( k − 1 ) Y(k-1) Y(k−1) 为系统在 k − 1 k-1 k−1 时刻的 “ 状态 ”。它的物理意义为:在输入不变条件下系统的非零初始状态响应的 N N N 步输出,即系统无外部输入的 “ 自由响应 ” 的 N N N 步输出为状态变量的 N N N 个分量。
因此定义
k
k
k 时刻的状态变量为
Y
(
k
)
=
{
[
y
(
k
)
y
(
k
+
1
)
⋯
y
(
k
+
N
−
2
)
y
(
k
+
N
−
1
)
]
T
Δ
u
(
k
+
i
)
=
0
,
i
=
0
,
1
,
⋯
,
N
−
1
}
(3)
Y(k)=\Bigg\lbrace \begin{matrix} \left[ \begin{matrix} y(k) & y(k+1) & \cdots & y(k+N-2) & y(k+N-1)\\ \end{matrix} \right]^T\\[2ex] \Delta u(k+i)=0,\quad i=0,1,\cdots,N-1 \end{matrix} \Bigg\rbrace\tag{3}
Y(k)={[y(k)y(k+1)⋯y(k+N−2)y(k+N−1)]TΔu(k+i)=0,i=0,1,⋯,N−1}(3)
当
Δ
u
(
k
−
1
)
=
0
\Delta u(k-1)=0
Δu(k−1)=0 时
Y
(
k
)
Y(k)
Y(k) 与
Y
(
k
−
1
)
Y(k-1)
Y(k−1) 之间的关系为
Y
(
k
)
=
M
s
s
Y
(
k
−
1
)
(4)
Y(k)=M_{ss}Y(k-1)\tag{4}
Y(k)=MssY(k−1)(4)
其中
M
s
s
=
[
0
1
0
⋯
0
0
0
0
1
⋯
0
0
⋮
⋮
⋮
⋮
⋮
0
0
0
⋯
0
1
0
0
0
⋯
0
1
]
(5)
M_{ss}=\left[ \begin{matrix} 0 & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 0 & 1\\ 0 & 0 & 0 & \cdots & 0 & 1\\ \end{matrix} \right]\tag{5}
Mss=
00⋮0010⋮0001⋮00⋯⋯⋯⋯00⋮0000⋮11
(5)
这就是稳定 SISO 系统的不变输入非零初始状态响应。
接着考虑当初始状态为零时,系统的输入响应
在系统平衡状态下做一个单位阶跃响应实验,得到系统的输出响应如图所示。
采样得到系统零初始条件下单位阶跃响应序列为
{
0
,
s
1
,
s
2
,
⋯
,
s
N
,
s
N
,
⋯
}
(6)
\lbrace 0,s_1,s_2,\cdots,s_N,s_N,\cdots\rbrace\tag{6}
{0,s1,s2,⋯,sN,sN,⋯}(6)
由线性系统的齐次性,对任意的输入变化
Δ
u
(
k
−
1
)
\Delta u(k-1)
Δu(k−1),系统的响应序列为
{
0
,
s
1
,
s
2
,
⋯
,
s
N
,
s
N
,
⋯
}
⋅
Δ
u
(
k
−
1
)
(6)
\lbrace 0,s_1,s_2,\cdots,s_N,s_N,\cdots\rbrace\cdot\Delta u(k-1)\tag{6}
{0,s1,s2,⋯,sN,sN,⋯}⋅Δu(k−1)(6)
记单位阶跃响应系数矩阵为
S
=
[
s
1
s
2
⋮
s
N
]
N
×
1
(7)
S=\left[ \begin{matrix} s_1 \\[1ex] s_2 \\[1ex] \vdots \\[1ex] s_N \\ \end{matrix} \right]_{N\times 1}\tag{7}
S=
s1s2⋮sN
N×1(7)
则零初始状态下,系统对任意输入的响应可以描述为
Y
(
k
)
=
S
Δ
u
(
k
−
1
)
(8)
Y(k)=S\Delta u(k-1)\tag{8}
Y(k)=SΔu(k−1)(8)
由于线性系统满足叠加性,因此,由(4)和(8)得到,在非零初始状态下系统对任意输入变化的响应为
Y
(
k
)
=
M
s
s
Y
(
k
−
1
)
+
S
Δ
u
(
k
−
1
)
(9)
\boxed{Y(k)=M_{ss}Y(k-1)+S\Delta u(k-1)}\tag{9}
Y(k)=MssY(k−1)+SΔu(k−1)(9)
其中初始条件为
Y
(
0
)
=
[
y
(
0
)
⋮
y
(
0
)
]
N
×
1
(10)
Y(0)=\left[ \begin{matrix} y(0) \\[1ex] \vdots \\[1ex] y(0) \\ \end{matrix} \right]_{N\times 1}\tag{10}
Y(0)=
y(0)⋮y(0)
N×1(10)
而系统在
k
k
k 时刻的输出为
y
(
k
)
=
C
Y
(
k
)
(11)
\boxed{y(k)=CY(k)}\tag{11}
y(k)=CY(k)(11)
其中
C
=
[
1
0
⋯
0
]
1
×
N
(12)
C=\left[ \begin{matrix} 1 & 0 & \cdots & 0 \\ \end{matrix} \right]_{1\times N}\tag{12}
C=[10⋯0]1×N(12)
综上,(9)和(11)就是系统的单位阶跃响应模型。
2、SISO 系统的动态矩阵控制(DMC)
2.1、被控系统描述
如图所示为被控系统,
u
u
u 为控制输入,
y
y
y 为输出,
d
d
d 为可以测量的外部干扰,
w
w
w 为不能测量的外部干扰,
P
u
P_u
Pu 为输入
u
u
u 到输出
y
y
y 的传递函数,
P
d
P_d
Pd 为可测量干扰
d
d
d 到输出
y
y
y 的传递函数。设控制输入
u
u
u 和可测量干扰
d
d
d 对输出
y
y
y 的单位阶跃响应系数矩阵分别为
S
u
=
[
s
1
u
s
2
u
⋮
s
N
u
]
,
S
d
=
[
s
1
d
s
2
d
⋮
s
N
d
]
,
(13)
S_u=\left[ \begin{matrix} s_1^u \\[1ex] s_2^u \\[1ex] \vdots \\[1ex] s_N^u \\ \end{matrix} \right],\quad S_d=\left[ \begin{matrix} s_1^d \\[1ex] s_2^d \\[1ex] \vdots \\[1ex] s_N^d \\ \end{matrix} \right], \tag{13}
Su=
s1us2u⋮sNu
,Sd=
s1ds2d⋮sNd
,(13)
由于线性系统满足齐次性和叠加性,因此,由(9)可以得出带可测干扰的单位阶跃响应模型为
Y
(
k
)
=
M
s
s
Y
(
k
−
1
)
+
S
u
Δ
u
(
k
−
1
)
+
S
d
Δ
d
(
k
−
1
)
y
(
k
)
=
C
Y
(
k
)
(14)
\begin{aligned} Y(k)&=M_{ss}Y(k-1)+S_u\Delta u(k-1)+S_d\Delta d(k-1)\\[1ex] y(k)&=CY(k)\tag{14} \end{aligned}
Y(k)y(k)=MssY(k−1)+SuΔu(k−1)+SdΔd(k−1)=CY(k)(14)
2.2、状态估计
由于单位阶跃响应模型(14)的状态不是全部可以测量的(只有第一个分量是可以测量的)。因此,需要对状态进行估计,用估计的状态作为初始条件预测系统未来的动态。
在
k
−
1
k-1
k−1 时刻,由(14)计算
k
k
k 时刻的状态,记为
Y
(
k
∣
k
−
1
)
Y(k|k-1)
Y(k∣k−1),即
Y
(
k
∣
k
−
1
)
=
M
s
s
Y
^
(
k
−
1
)
+
S
u
Δ
u
(
k
−
1
)
+
S
d
Δ
d
(
k
−
1
)
(15)
Y(k|k-1)=M_{ss}\hat Y(k-1)+S_u\Delta u(k-1)+S_d\Delta d(k-1)\tag{15}
Y(k∣k−1)=MssY^(k−1)+SuΔu(k−1)+SdΔd(k−1)(15)
其中,
Y
^
(
k
−
1
)
\hat Y(k-1)
Y^(k−1) 是
k
−
1
k-1
k−1 时刻对状态的估计。计算
k
k
k 时刻的输出为
y
(
k
∣
k
−
1
)
=
C
Y
(
k
∣
k
−
1
)
(16)
y(k|k-1) = CY(k|k-1)\tag{16}
y(k∣k−1)=CY(k∣k−1)(16)
在
k
k
k 时刻的测量值为
y
(
k
)
y(k)
y(k),与计算值之差为
y
(
k
)
−
y
(
k
∣
k
−
1
)
y(k)-y(k|k-1)
y(k)−y(k∣k−1)。以这个误差作为校正量,得到校正后的状态分量如下:
y
^
(
k
∣
k
)
=
y
(
k
∣
k
−
1
)
+
[
y
(
k
)
−
y
(
k
∣
k
−
1
)
]
,
y
^
(
k
+
1
∣
k
)
=
y
(
k
+
1
∣
k
−
1
)
+
[
y
(
k
)
−
y
(
k
∣
k
−
1
)
]
,
⋮
,
y
^
(
k
+
N
−
1
∣
k
)
=
y
(
k
+
N
−
1
∣
k
−
1
)
+
[
y
(
k
)
−
y
(
k
∣
k
−
1
)
]
.
(17)
\begin{aligned} \hat y(k|k) & = y(k|k-1) + [y(k)-y(k|k-1)],\\[1ex] \hat y(k+1|k) & = y(k+1|k-1) + [y(k)-y(k|k-1)],\\[1ex] \vdots,\\[1ex] \hat y(k+N-1|k) & = y(k+N-1|k-1) + [y(k)-y(k|k-1)]. \end{aligned}\tag{17}
y^(k∣k)y^(k+1∣k)⋮,y^(k+N−1∣k)=y(k∣k−1)+[y(k)−y(k∣k−1)],=y(k+1∣k−1)+[y(k)−y(k∣k−1)],=y(k+N−1∣k−1)+[y(k)−y(k∣k−1)].(17)
记
Y
^
(
k
)
=
[
y
^
(
k
∣
k
)
y
^
(
k
+
1
∣
k
)
⋮
y
^
(
k
+
N
−
1
∣
k
)
]
N
×
1
(18)
\hat Y(k)=\left[ \begin{matrix} \hat y(k|k)\\[1ex] \hat y(k+1|k)\\[1ex] \vdots\\[1ex] \hat y(k+N-1|k) \end{matrix} \right]_{N\times 1}\tag{18}
Y^(k)=
y^(k∣k)y^(k+1∣k)⋮y^(k+N−1∣k)
N×1(18)
则上式变为
Y
^
(
k
)
=
Y
(
k
∣
k
−
1
)
+
K
I
(
y
(
k
)
−
y
(
k
∣
k
−
1
)
)
,
其中
K
I
=
[
1
⋮
1
]
N
×
1
(19)
\hat Y(k) = Y(k|k-1) + K_I\big(y(k)-y(k|k-1)\big),\quad 其中\quad K_I=\left[ \begin{matrix} 1\\ \vdots\\[1ex] 1\\ \end{matrix} \right]_{N\times 1}\tag{19}
Y^(k)=Y(k∣k−1)+KI(y(k)−y(k∣k−1)),其中KI=
1⋮1
N×1(19)
将(15)和(16)代入(19),得
Y
^
(
k
)
=
Y
(
k
∣
k
−
1
)
+
K
I
(
y
(
k
)
−
y
(
k
∣
k
−
1
)
)
=
(
I
−
K
I
C
)
Y
(
k
∣
k
−
1
)
+
K
I
y
(
k
)
=
(
I
−
K
I
C
)
M
s
s
Y
^
(
k
−
1
)
+
K
I
y
(
k
)
+
(
I
−
K
I
C
)
S
u
Δ
u
(
k
−
1
)
+
(
I
−
K
I
C
)
S
d
Δ
d
(
k
−
1
)
(20)
\begin{aligned} \hat Y(k) &= Y(k|k-1) + K_I\big(y(k)-y(k|k-1)\big)\\[1ex] &=(I-K_IC)Y(k|k-1) + K_Iy(k)\\[1ex] &=(I-K_IC)M_{ss}\hat Y(k-1) + K_Iy(k)+(I-K_IC)S_u\Delta u(k-1)+(I-K_IC)S_d\Delta d(k-1)\\[1ex] \end{aligned}\tag{20}
Y^(k)=Y(k∣k−1)+KI(y(k)−y(k∣k−1))=(I−KIC)Y(k∣k−1)+KIy(k)=(I−KIC)MssY^(k−1)+KIy(k)+(I−KIC)SuΔu(k−1)+(I−KIC)SdΔd(k−1)(20)
这是一个典型的状态估计器方程,其中
(
I
−
K
I
C
)
M
s
s
=
[
0
0
0
0
⋯
0
0
−
1
1
0
⋯
0
0
−
1
0
1
⋯
0
⋮
⋮
⋮
⋮
⋮
0
−
1
0
0
⋯
1
0
−
1
0
0
⋯
1
]
N
×
N
(21)
(I-K_IC)M_{ss}=\left[ \begin{matrix} 0 & 0 & 0 & 0 & \cdots & 0\\[1ex] 0 & -1 & 1 & 0 & \cdots & 0\\[1ex] 0 & -1 & 0 & 1 & \cdots & 0\\[1ex] \vdots & \vdots & \vdots & \vdots & & \vdots\\[1ex] 0 & -1 & 0 & 0 & \cdots & 1\\[1ex] 0 & -1 & 0 & 0 & \cdots & 1\\ \end{matrix} \right]_{N\times N}\tag{21}
(I−KIC)Mss=
000⋮000−1−1⋮−1−1010⋮00001⋮00⋯⋯⋯⋯⋯000⋮11
N×N(21)
可以证明, ( I − K I C ) M s s (I-K_IC)M_{ss} (I−KIC)Mss 的所有特征值均位于单位圆内,因此,估计器(20)是名义渐近稳定的。
3、预测方程
采用 基于状态空间模型的无约束预测控制 相同的方法推导预测方程。
对系统未来
p
p
p 步输出的预测可以由下面的预测方程计算:
Y
p
(
k
+
1
∣
k
)
=
M
Y
^
(
k
)
+
S
u
Δ
U
(
k
)
+
S
d
Δ
d
(
k
)
Y_p(k+1|k)={\cal M}\hat Y(k)+{\cal S_u}\Delta U(k) + {\cal S_d}\Delta d(k)
Yp(k+1∣k)=MY^(k)+SuΔU(k)+SdΔd(k)
其中,
M
=
[
0
C
c
0
0
0
⋯
0
⋯
0
0
0
C
c
0
0
⋯
0
⋯
0
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
0
0
⋯
0
C
c
⋯
0
⋯
0
]
p
×
N
S
d
=
[
C
c
S
1
d
C
c
S
2
d
⋮
C
c
S
p
d
]
p
×
1
S
u
=
[
C
c
S
1
u
0
0
⋯
0
C
c
S
2
u
C
c
S
1
u
0
⋯
0
⋮
⋮
⋮
⋮
C
c
S
m
u
C
c
S
m
−
1
u
⋯
⋯
C
c
S
1
u
⋮
⋮
⋮
⋮
C
c
S
p
u
C
c
S
p
−
1
u
⋯
⋯
C
c
S
p
−
m
+
1
u
]
p
×
m
\begin{aligned} &{\cal M}=\left[ \begin{matrix} \pmb 0 & \pmb C_c & \pmb 0 & \pmb 0 & \pmb 0 &\cdots & \pmb 0 &\cdots & \pmb 0\\[1ex] \pmb 0 & \pmb 0 & \pmb C_c & \pmb 0 & \pmb 0 &\cdots & \pmb 0 &\cdots & \pmb 0\\[1ex] \vdots & \vdots & \vdots &\vdots &\vdots &\vdots & \vdots & & \vdots\\[1ex] \pmb 0 & \pmb 0 & \cdots & \pmb 0 & \pmb C_c &\cdots & \pmb 0 &\cdots & \pmb 0\\ \end{matrix} \right]_{p\times N}\\ &{\cal S_d}=\left[ \begin{matrix} \pmb C_c\pmb S_1^d\\[1ex] \pmb C_c\pmb S_2^d\\[1ex] \vdots\\[1ex] \pmb C_c\pmb S_p^d\\ \end{matrix} \right]_{p\times1}\\ &{\cal S_u}=\left[ \begin{matrix} \pmb C_c\pmb S_1^u & \pmb 0 & \pmb 0 & \cdots & \pmb 0 \\[1ex] \pmb C_c\pmb S_2^u & \pmb C_c\pmb S_1^u & \pmb 0 & \cdots & \pmb 0 \\[1ex] \vdots & \vdots & \vdots & & \vdots \\[1ex] \pmb C_c\pmb S_m^u & \pmb C_c\pmb S_{m-1}^u & \cdots & \cdots & \pmb C_c\pmb S_{1}^u \\[1ex] \vdots & \vdots & \vdots & & \vdots \\[1ex] \pmb C_c\pmb S_p^u & \pmb C_c\pmb S_{p-1}^u & \cdots & \cdots & \pmb C_c\pmb S_{p-m+1}^u \\[1ex] \end{matrix} \right]_{p\times m} \end{aligned}
M=
00⋮0Cc0⋮00Cc⋮⋯00⋮000⋮Cc⋯⋯⋮⋯00⋮0⋯⋯⋯00⋮0
p×NSd=
CcS1dCcS2d⋮CcSpd
p×1Su=
CcS1uCcS2u⋮CcSmu⋮CcSpu0CcS1u⋮CcSm−1u⋮CcSp−1u00⋮⋯⋮⋯⋯⋯⋯⋯00⋮CcS1u⋮CcSp−m+1u
p×m