【MPC学习笔记】02：MPC详细简介（Lecture 1_1 Unconstrained MPC）
本笔记来自北航诸兵老师的课程
课程地址：模型预测控制（2022春）lecture 1-1 Unconstrained MPC
接上一篇：【MPC学习笔记】01：MPC简介（Lecture 1_1 Unconstrained MPC）
1 详细介绍
1.1 状态方程
对 LTI 离散系统：
x
(
k
+
1
)
=
A
x
(
k
)
+
B
u
(
k
)
x
∈
R
n
,
u
∈
R
p
x(k+1) = Ax(k) + Bu(k)\quad x\in\R^n, u\in R^p
x(k+1)=Ax(k)+Bu(k)x∈Rn,u∈Rp
对传统控制系统，连续系统是好处理的，离散系统是要额外考虑其他因素的
对MPC，则是反过来，离散系统是好处理的，连续系统是要额外考虑其他因素的
假设：
- A , B A,B A,B 可控（Stablizable）
- 状态和控制输入不存在约束（本节讨论无约束MPC）
Define：
x
(
i
∣
k
)
,
u
(
i
∣
k
)
x(i|k), u(i|k)
x(i∣k),u(i∣k): Prediction of
i
i
isteps ahead from time
k
k
k (比如，在时刻
k
k
k 预测下一时刻的状态，记为
x
(
1
∣
k
)
x(1|k)
x(1∣k)，当前时刻的输入，记为
u
(
0
∣
k
)
u(0|k)
u(0∣k))
预测：
x
(
1
∣
k
)
=
A
x
(
0
∣
k
)
+
B
u
(
0
∣
k
)
x
(
2
∣
k
)
=
A
x
(
1
∣
k
)
+
B
u
(
1
∣
k
)
=
A
[
A
x
(
0
∣
k
)
+
B
u
(
0
∣
k
)
]
+
B
u
(
1
∣
k
)
=
A
2
x
(
0
∣
k
)
+
A
B
u
(
0
∣
k
)
+
B
u
(
1
∣
k
)
⋮
⋮
⋮
x
(
i
∣
k
)
=
A
i
x
(
0
∣
k
)
+
A
i
−
1
B
u
(
0
∣
k
)
+
A
i
−
2
B
2
u
(
1
∣
k
)
+
⋯
+
B
u
(
i
−
1
∣
k
)
\begin{aligned} x(1|k) &= Ax(0|k) + Bu(0|k) \\ x(2|k) &= Ax(1|k) + Bu(1|k)=A[Ax(0|k) + Bu(0|k)] + Bu(1|k) \\ &=A^2x(0|k) + ABu(0|k) + Bu(1|k) \\ \quad&\quad \quad \quad \quad \quad \vdots\quad \quad \quad \vdots\quad \quad \quad \vdots \\ x(i|k) &= A^ix(0|k) + A^{i-1}Bu(0|k) + A^{i-2}B^2u(1|k) + \cdots + Bu(i-1|k) \end{aligned}
x(1∣k)x(2∣k)x(i∣k)=Ax(0∣k)+Bu(0∣k)=Ax(1∣k)+Bu(1∣k)=A[Ax(0∣k)+Bu(0∣k)]+Bu(1∣k)=A2x(0∣k)+ABu(0∣k)+Bu(1∣k)⋮⋮⋮=Aix(0∣k)+Ai−1Bu(0∣k)+Ai−2B2u(1∣k)+⋯+Bu(i−1∣k)
In compact form:
X
(
k
)
=
F
x
(
k
)
+
Φ
U
(
k
)
X(k) = Fx(k) + \Phi U(k)
X(k)=Fx(k)+ΦU(k)
X
(
k
)
≜
[
x
(
1
∣
k
)
x
(
2
∣
k
)
⋮
x
(
N
∣
k
)
]
U
(
k
)
≜
[
u
(
0
∣
k
)
u
(
1
∣
k
)
⋮
u
(
N
−
1
∣
k
)
]
X(k)\triangleq \begin{bmatrix} x(1|k)\\ x(2|k)\\ \vdots\\ x(N|k) \end{bmatrix} \quad\quad U(k)\triangleq \begin{bmatrix} u(0|k)\\ u(1|k)\\ \vdots\\ u(N-1|k) \end{bmatrix}
X(k)≜
x(1∣k)x(2∣k)⋮x(N∣k)
U(k)≜
u(0∣k)u(1∣k)⋮u(N−1∣k)
X ( k ) X(k) X(k) 式中的 x ( k ) x(k) x(k) 也即 x ( 0 ∣ k ) x(0|k) x(0∣k)
≜ \triangleq ≜ : 表示定义为
N N N : Control/Predictive horizon，实际上二者有区别，但这里不做区分
1.2 Cost Function
这里cost function 的控制/预测时域是一个有限的数
J
(
k
)
=
∑
i
=
1
N
∣
∣
x
(
i
∣
k
)
∣
∣
Q
2
+
∣
∣
u
(
i
−
1
∣
k
)
∣
∣
R
2
=
X
T
(
k
)
Q
X
(
k
)
+
U
T
(
k
)
R
U
(
k
)
\begin{aligned} J(k) &= \sum^{N}_{i=1}||x(i|k)||_Q^2 + ||u(i-1|k)||_R^2 \\ &= X^T(k)\mathcal{Q}X(k) + U^T(k)\mathcal{R}U(k) \end{aligned}
J(k)=i=1∑N∣∣x(i∣k)∣∣Q2+∣∣u(i−1∣k)∣∣R2=XT(k)QX(k)+UT(k)RU(k)
假设
Q
Q
Q和
R
R
R是正定的，是权重
Q
=
[
Q
Q
⋱
Q
]
R
=
[
R
R
⋱
R
]
\mathcal{Q} = \begin{bmatrix} Q\\ &Q\\ &&\ddots\\ &&&Q \end{bmatrix} \quad\quad \mathcal{R} = \begin{bmatrix} R\\ &R\\ &&\ddots\\ &&&R \end{bmatrix}
Q=
QQ⋱Q
R=
RR⋱R
将
X
(
k
)
=
F
x
(
k
)
+
Φ
U
(
k
)
X(k) = Fx(k) + \Phi U(k)
X(k)=Fx(k)+ΦU(k) 代入
J
(
k
)
J(k)
J(k) 有
J
(
k
)
=
(
F
x
(
k
)
+
Φ
U
(
k
)
)
T
Q
(
F
x
(
k
)
+
Φ
U
(
k
)
)
+
U
T
(
k
)
R
U
(
k
)
=
(
x
(
k
)
T
F
T
+
U
(
k
)
T
Φ
T
)
(
Q
F
x
(
k
)
+
Q
Φ
U
(
k
)
)
+
U
T
(
k
)
R
U
(
k
)
=
x
(
k
)
T
F
T
Q
F
x
(
k
)
+
U
(
k
)
T
Φ
T
Q
F
x
(
k
)
+
x
(
k
)
T
F
T
Q
Φ
U
(
k
)
+
U
(
k
)
T
Φ
T
Q
Φ
U
(
k
)
+
U
T
(
k
)
R
U
(
k
)
=
x
(
k
)
T
F
T
Q
F
x
(
k
)
+
2
x
(
k
)
T
F
T
Q
Φ
U
(
k
)
+
U
(
k
)
T
(
Φ
T
Q
Φ
+
R
)
U
(
k
)
\begin{equation*} \begin{aligned} J(k) &= (Fx(k) + \Phi U(k))^{T} \mathcal{Q} (Fx(k) + \Phi U(k)) + U^T(k)\mathcal{R}U(k) \\ &= (x(k)^TF^T + U(k)^T\Phi^T)(\mathcal{Q}Fx(k) + \mathcal{Q}\Phi U(k))+ U^T(k)\mathcal{R}U(k) \\ &= \textcolor{green}{x(k)^TF^T\mathcal{Q}Fx(k)}+ \textcolor{red}{U(k)^T\Phi^T\mathcal{Q}Fx(k) + x(k)^TF^T\mathcal{Q}\Phi U(k)} + \textcolor{blue}{U(k)^T\Phi^T\mathcal{Q}\Phi U(k)+ U^T(k)\mathcal{R}U(k)} \\ &=\textcolor{green}{x(k)^TF^T\mathcal{Q}Fx(k)}+ \textcolor{red}{2x(k)^TF^T\mathcal{Q}\Phi U(k)} +\textcolor{blue}{U(k)^T(\Phi^T\mathcal{Q}\Phi+\mathcal{R})U(k)} \end{aligned} \end{equation*}
J(k)=(Fx(k)+ΦU(k))TQ(Fx(k)+ΦU(k))+UT(k)RU(k)=(x(k)TFT+U(k)TΦT)(QFx(k)+QΦU(k))+UT(k)RU(k)=x(k)TFTQFx(k)+U(k)TΦTQFx(k)+x(k)TFTQΦU(k)+U(k)TΦTQΦU(k)+UT(k)RU(k)=x(k)TFTQFx(k)+2x(k)TFTQΦU(k)+U(k)T(ΦTQΦ+R)U(k)
F x ( k ) Fx(k) Fx(k) 和 Φ U ( k ) \Phi U(k) ΦU(k) 维数相同（是系统状态的维数*N），而 Q \mathcal{Q} Q是一个对角方阵，故 U ( k ) T Φ T Q F x ( k ) = x ( k ) T F T Q Φ U ( k ) = 一个标量 {U(k)^T\Phi^T\mathcal{Q}Fx(k) = x(k)^TF^T\mathcal{Q}\Phi U(k)}=一个标量 U(k)TΦTQFx(k)=x(k)TFTQΦU(k)=一个标量，故红色部分相加相当于其中一个乘2
1.3 状态变量 u ( k ) u(k) u(k) 的求解
Minimize the control function by predictive control series:
（可以不严谨地理解为：让
J
J
J最小，相当于求
J
J
J在
J
J
J对
U
U
U导数为0 点的值）
∇
U
∣
U
=
U
∗
=
∂
J
∂
U
∣
U
=
U
∗
=
0
\nabla_U{|}_{U=U^*}=\frac{\partial{J}}{\partial U}{|}_{U=U^*}=0
∇U∣U=U∗=∂U∂J∣U=U∗=0
∂
J
∂
U
=
0
+
2
x
(
k
)
T
F
T
Q
Φ
+
2
U
(
k
)
T
(
Φ
T
Q
Φ
+
R
)
\frac{\partial{J}}{\partial U}=\textcolor{green}{0} +\textcolor{red}{2x(k)^TF^T\mathcal{Q}\Phi}+\textcolor{blue}{2U(k)^T(\Phi^T\mathcal{Q}\Phi+\mathcal{R})}
∂U∂J=0+2x(k)TFTQΦ+2U(k)T(ΦTQΦ+R)
令
∂
J
∂
U
=
0
\frac{\partial{J}}{\partial U}=0
∂U∂J=0，可得：
x
(
k
)
T
F
T
Q
Φ
+
U
(
k
)
T
(
Φ
T
Q
Φ
+
R
)
=
0
(
x
(
k
)
T
F
T
Q
Φ
+
U
(
k
)
T
(
Φ
T
Q
Φ
+
R
)
)
T
=
0
(
x
(
k
)
T
F
T
Q
Φ
)
T
+
(
Φ
T
Q
Φ
+
R
)
T
U
(
k
)
=
0
Φ
T
Q
F
x
(
k
)
+
(
Φ
T
Q
Φ
+
R
)
U
(
k
)
=
0
(
Φ
T
Q
Φ
+
R
)
U
(
k
)
=
−
Φ
T
Q
F
x
(
k
)
U
(
k
)
=
−
(
Φ
T
Q
Φ
+
R
)
−
1
Φ
T
Q
F
x
(
k
)
\begin{equation*} \begin{aligned} \textcolor{red}{x(k)^T F^T \mathcal{Q} \Phi}+\textcolor{blue}{U(k)^T (\Phi^T\mathcal{Q}\Phi+\mathcal{R})}= 0 \\ (\textcolor{red}{x(k)^TF^T\mathcal{Q}\Phi}+\textcolor{blue}{U(k)^T(\Phi^T\mathcal{Q}\Phi+\mathcal{R})})^T=0 \\ \textcolor{red}{(x(k)^TF^T\mathcal{Q}\Phi)^T}+\textcolor{blue}{(\Phi^T\mathcal{Q}\Phi+\mathcal{R})^TU(k)}=0 \\ \textcolor{red}{\Phi^T \mathcal{Q} Fx(k)} + \textcolor{blue} {(\Phi^T\mathcal{Q}\Phi+\mathcal{R})U(k)}=0 \\ \textcolor{blue}{(\Phi^T\mathcal{Q}\Phi+\mathcal{R})U(k)}=-\textcolor{red}{\Phi^T\mathcal{Q}Fx(k)} \\ \textcolor{blue}{U(k)}=-\textcolor{blue}{(\Phi^T\mathcal{Q}\Phi+\mathcal{R})^{-1}}\textcolor{red}{\Phi^T\mathcal{Q}Fx(k)} \end{aligned} \end{equation*}
x(k)TFTQΦ+U(k)T(ΦTQΦ+R)=0(x(k)TFTQΦ+U(k)T(ΦTQΦ+R))T=0(x(k)TFTQΦ)T+(ΦTQΦ+R)TU(k)=0ΦTQFx(k)+(ΦTQΦ+R)U(k)=0(ΦTQΦ+R)U(k)=−ΦTQFx(k)U(k)=−(ΦTQΦ+R)−1ΦTQFx(k)
即：
U
(
k
)
=
−
(
Φ
T
Q
Φ
+
R
)
−
1
Φ
T
Q
F
x
(
k
)
(
R
>
0
,
Q
≥
0
;
o
r
R
≥
0
,
Q
>
0
,
a
n
d
Φ
i
s
f
u
l
l
y
r
a
n
k
e
d
)
\begin{equation*} \begin{aligned} \textcolor{blue}{U(k)}=- \textcolor{blue}{(\Phi^T\mathcal{Q}\Phi+\mathcal{R})^{-1}} \textcolor{red}{\Phi^T\mathcal{Q}Fx(k)}\\ (R>0,Q\ge0;or R\ge0,Q\gt0, and\ \Phi\ is\ fully\ ranked) \end{aligned} \end{equation*}
U(k)=−(ΦTQΦ+R)−1ΦTQFx(k)(R>0,Q≥0;orR≥0,Q>0,and Φ is fully ranked)
满足括号里的条件， ( Φ T Q Φ + R ) (\Phi^T\mathcal{Q}\Phi+\mathcal{R}) (ΦTQΦ+R)才可逆
u
∗
(
k
)
=
−
[
I
p
×
p
0
⋯
0
]
(
Φ
T
Q
Φ
+
R
)
−
1
Φ
T
Q
F
x
(
k
)
=
−
K
m
p
c
x
(
k
)
\begin{equation*} \begin{aligned} \begin{aligned} \textcolor{blue}{u^*(k)}&=- \begin{bmatrix} I_{p\times p}&0&\cdots 0 \end{bmatrix} \textcolor{blue}{(\Phi^T\mathcal{Q}\Phi+\mathcal{R})^{-1}} \textcolor{red}{\Phi^T\mathcal{Q}Fx(k)} \\ &=-K_{mpc}x(k) \end{aligned} \end{aligned} \end{equation*}
u∗(k)=−[Ip×p0⋯0](ΦTQΦ+R)−1ΦTQFx(k)=−Kmpcx(k)
取
u
∗
(
k
)
u^*(k)
u∗(k) ，则在
U
(
k
)
U(k)
U(k) 前乘一个分块矩阵，对角线上的第一个分块是一个单位阵，维度为控制输入的维度
p
p
p ,如果不是多输入而是单输入，则
p
=
1
p=1
p=1 。
u
∗
(
k
)
u^*(k)
u∗(k) 最后化简为一个常数矩阵
K
m
p
c
K_{mpc}
Kmpc （因为
Φ
\Phi
Φ,
Q
\mathcal{Q}
Q,
R
\mathcal{R}
R,
F
F
F这些全部是已知量）乘上
k
k
k 时刻的状态变量，从形式上看是状态反馈。
因此，无约束线性MPC实际上是一个线性反馈控制。
1.4 举例
写一段matlab程序，即可求解
F
F
F，
Φ
\Phi
Φ 和
u
∗
(
k
)
u^*(k)
u∗(k)。但问题是，求解出来的
u
∗
(
k
)
u^*(k)
u∗(k)是否能保证系统是稳定的？
对于有稳定性的判定，有李雅普诺夫直接法和李雅普诺夫间接法，见下：
- > 0 和 < 0 >0和<0 >0和<0 分别指的是正定和负定
- 验证稳定性的前提是 K m p c K_{mpc} Kmpc 存在
- 可优化性 并不决定 可稳定性（两种可能的原因见下a和b），所以这里验证稳定性的操作是必要的。
- 优化是在一段时间上进行的，在这段时间内， x x x不一定由大变小，也可能先变小再变大，从而不收敛
- 对于非最小相位系统，系统响应方向可能相反，N取得不够大时，预测不能反映真实运动趋势，那么优化的不是系统真正的性能，导致不稳定的情况发生
本例采用李雅普诺夫间接法在离散时间系统下的判据
当
N
N
N的取值越来越小，如下面这个例子所示， 特征值超出了单位圆，从而不稳定。
知道如何判断稳定性后，现在的问题变成了：每一次优化后，都要去验证一下系统的稳定性呢？有没有一种机制，保证每一次算出来的
K
K
K都保证系统的稳定性？