【预测控制4】基于状态空间模型的预测控制
文章目录
- 一、基于状态空间模型的预测控制
- 1、预测模型
- 2、滚动优化
- 2.1 状态优化问题
- 2.2 输出优化问题
- 3、反馈校正
- 二、Matlab仿真示例
相关文章:
【预测控制1】模型预测控制概论
【预测控制2】预测控制基本原理
【预测控制3】动态矩阵控制（DMC）及Matlab仿真
一、基于状态空间模型的预测控制
1、预测模型
? 首先考虑SISO线性系统
x
(
k
+
1
)
=
A
x
(
k
)
+
b
u
(
k
)
y
(
k
)
=
c
T
x
(
k
)
(1-1)
\pmb x(k+1)=\pmb{Ax}(k)+\pmb bu(k)\\y(k)=\pmb{c^Tx}(k) \tag{1-1}
xxx(k+1)=AxAxAx(k)+bbbu(k)y(k)=cTxcTxcTx(k)(1-1)
状态变量
x
(
k
)
∈
R
n
\pmb{x}(k)\in R^n
xxx(k)∈Rn实时可测，
y
(
k
)
、
u
(
k
)
y(k)、u(k)
y(k)、u(k)分别为系统的输入输出。
? 设从
k
k
k时刻起系统输入发生
M
M
M步变化，而后保持不变，则由模型(1-1)，可以预测输出在
u
(
k
)
,
u
(
k
+
1
)
,
.
.
.
,
u
(
k
+
M
?
1
)
u(k),u(k+1),...,u(k+M-1)
u(k),u(k+1),...,u(k+M?1)作用下未来P个时刻的系统状态
x
(
k
+
1
∣
k
)
=
A
x
(
k
∣
k
)
+
b
u
(
k
)
x
(
k
+
2
∣
k
)
=
A
x
(
k
+
1
∣
k
)
+
b
u
(
k
)
=
A
2
x
(
k
∣
k
)
+
A
b
u
(
k
)
+
b
u
(
k
+
1
)
?
x
(
k
+
M
∣
k
)
=
A
M
x
(
k
∣
k
)
+
A
M
?
1
b
u
(
k
)
+
A
M
?
2
b
u
(
k
+
1
)
+
?
+
b
u
(
k
+
M
?
1
)
x
(
k
+
M
+
1
∣
k
)
=
A
M
+
1
x
(
k
∣
k
)
+
A
M
b
u
(
k
)
+
A
M
?
1
b
u
(
k
+
1
)
+
?
+
(
A
b
+
b
)
u
(
k
+
M
?
1
)
?
x
(
k
+
P
∣
k
)
=
A
P
x
(
k
∣
k
)
+
A
P
?
1
b
u
(
k
)
+
A
P
?
2
b
u
(
k
+
1
)
+
?
+
(
A
P
?
M
b
+
A
P
?
M
?
1
b
+
?
+
b
)
u
(
k
+
M
?
1
)
\begin{aligned}\pmb x(k+1|k)&=\pmb{Ax}(k|k)+\pmb bu(k)\\ \pmb x(k+2|k)&=\pmb{Ax}(k+1|k)+\pmb bu(k)=\pmb{A^2x}(k|k)+\pmb{Ab}u(k)+\pmb bu(k+1)\\ \vdots \\\pmb x(k+M|k)&=\pmb{A}^M\pmb{x}(k|k)+\pmb A^{M-1}\pmb bu(k)+\pmb A^{M-2}\pmb bu(k+1)+\dots+\pmb bu(k+M-1) \\\pmb x(k+M+1|k)&=\pmb{A}^{M+1}\pmb{x}(k|k)+\pmb A^{M}\pmb bu(k)+\pmb A^{M-1}\pmb bu(k+1)+\dots+(\pmb{Ab}+\pmb b)u(k+M-1) \\\vdots \\\pmb x(k+P|k)&=\pmb{A}^P\pmb{x}(k|k)+\pmb A^{P-1}\pmb bu(k)+\pmb A^{P-2}\pmb bu(k+1)+\\&\dots+(\pmb A^{P-M}\pmb b+\pmb A^{P-M-1}\pmb b+\dots+\pmb b)u(k+M-1) \end{aligned}
xxx(k+1∣k)xxx(k+2∣k)?xxx(k+M∣k)xxx(k+M+1∣k)?xxx(k+P∣k)?=AxAxAx(k∣k)+bbbu(k)=AxAxAx(k+1∣k)+bbbu(k)=A2xA2xA2x(k∣k)+AbAbAbu(k)+bbbu(k+1)=AAAMxxx(k∣k)+AAAM?1bbbu(k)+AAAM?2bbbu(k+1)+?+bbbu(k+M?1)=AAAM+1xxx(k∣k)+AAAMbbbu(k)+AAAM?1bbbu(k+1)+?+(AbAbAb+bbb)u(k+M?1)=AAAPxxx(k∣k)+AAAP?1bbbu(k)+AAAP?2bbbu(k+1)+?+(AAAP?Mbbb+AAAP?M?1bbb+?+bbb)u(k+M?1)?
将上面的式子整合为向量形式，有
X
(
k
)
=
F
x
x
(
k
∣
k
)
+
G
x
U
(
k
)
(1-2)
\pmb X(k)=\pmb{F_xx}(k|k)+\pmb{G_xU}(k)\tag{1-2}
XXX(k)=Fx?x?Fx?x??Fx?x(k∣k)+Gx?U?Gx?U??Gx?U(k)(1-2)其中
F
x
=
[
A
A
2
?
A
P
]
n
P
×
1
,
G
x
=
[
b
0
…
0
0
A
b
b
…
0
0
?
?
?
?
?
A
P
?
1
b
A
P
?
2
b
…
A
P
?
M
+
1
b
∑
i
=
0
P
?
M
A
i
b
]
n
P
×
M
X
(
k
)
=
[
x
(
k
+
1
∣
k
)
x
(
k
+
2
∣
k
)
?
x
(
k
+
P
∣
k
)
]
n
P
×
1
,
U
(
k
)
=
[
u
(
k
)
u
(
k
+
1
)
?
u
(
k
+
M
?
1
)
]
M
×
1
\pmb{F_x}=\begin{bmatrix}\pmb A \\\pmb A^2 \\\vdots \\\pmb A^{P} \end{bmatrix}_{nP\times 1}, \pmb{Gx}=\begin{bmatrix}\pmb b&0&\dots&0&0 \\\pmb{Ab}&\pmb b&\dots&0&0 \\\vdots&\vdots&\ddots &\vdots&\vdots \\\pmb{A}^{P-1}\pmb b&\pmb{A}^{P-2}\pmb b&\dots&\pmb{A}^{P-M+1}\pmb b&\sum_{i=0}^{P-M}\pmb A^i\pmb b \end{bmatrix}_{nP\times M} \\\pmb X(k)=\begin{bmatrix}\pmb x(k+1|k)\\\pmb x(k+2|k)\\\vdots\\\pmb x(k+P|k)\end{bmatrix}_{nP\times1}, \pmb U(k)=\begin{bmatrix}\pmb u(k)\\\pmb u(k+1)\\\vdots\\\pmb u(k+M-1)\end{bmatrix}_{M\times1}
Fx??Fx???Fx?=??????AAAAAA2?AAAP???????nP×1?,GxGxGx=??????bbbAbAbAb?AAAP?1bbb?0bbb?AAAP?2bbb?……?…?00?AAAP?M+1bbb?00?∑i=0P?M?AAAibbb???????nP×M?XXX(k)=??????xxx(k+1∣k)xxx(k+2∣k)?xxx(k+P∣k)???????nP×1?,UUU(k)=??????uuu(k)uuu(k+1)?uuu(k+M?1)???????M×1?
类似的，根据输出方程
y
(
k
)
=
c
T
x
(
k
)
y(k)=\pmb{c^Tx}(k)
y(k)=cTxcTxcTx(k)可以得到
Y
(
k
)
=
F
y
x
(
k
∣
k
)
+
G
y
U
(
k
)
(1-3)
\pmb Y(k)=\pmb{F_yx}(k|k)+\pmb{G_yU}(k)\tag{1-3}
YYY(k)=Fy?x?Fy?x??Fy?x(k∣k)+Gy?U?Gy?U??Gy?U(k)(1-3)其中
F
y
=
[
c
T
A
c
T
A
2
?
c
T
A
P
]
P
×
n
,
G
y
=
[
c
T
b
0
…
0
0
c
T
A
b
c
T
b
…
0
0
?
?
?
?
?
c
T
A
P
?
1
b
c
T
A
P
?
2
b
…
c
T
A
P
?
M
+
1
b
∑
i
=
0
P
?
M
c
T
A
i
b
]
P
×
M
Y
(
k
)
=
[
y
(
k
+
1
∣
k
)
y
(
k
+
2
∣
k
)
?
y
(
k
+
P
∣
k
)
]
P
×
1
\pmb{F_y}=\begin{bmatrix}\pmb{c^TA} \\\pmb{c^TA^2} \\\vdots \\\pmb{c^TA^{P}} \end{bmatrix}_{P\times n}, \pmb{Gy}=\begin{bmatrix}\pmb{c^Tb}&0&\dots&0&0 \\\pmb{c^TAb}&\pmb{c^Tb}&\dots&0&0 \\\vdots&\vdots&\ddots &\vdots&\vdots \\\pmb{c^TA}^{P-1}\pmb b&\pmb{c^TA}^{P-2}\pmb b&\dots&\pmb{c^TA}^{P-M+1}\pmb b&\sum_{i=0}^{P-M}\pmb{c^TA}^i\pmb b \end{bmatrix}_{P\times M} \\\pmb Y(k)=\begin{bmatrix}\pmb y(k+1|k)\\\pmb y(k+2|k)\\\vdots\\\pmb y(k+P|k)\end{bmatrix}_{P\times1}
Fy??Fy???Fy?=??????cTAcTAcTAcTA2cTA2cTA2?cTAPcTAPcTAP???????P×n?,Gy?Gy??Gy=??????cTbcTbcTbcTAbcTAbcTAb?cTAcTAcTAP?1bbb?0cTbcTbcTb?cTAcTAcTAP?2bbb?……?…?00?cTAcTAcTAP?M+1bbb?00?∑i=0P?M?cTAcTAcTAibbb???????P×M?YYY(k)=??????y?y??y(k+1∣k)y?y??y(k+2∣k)?y?y??y(k+P∣k)???????P×1? ? 在此总结一下，式（1-2）是状态预测模型，根据状态预测模型，我们可以通过当前时刻的状态变量
x
(
k
∣
k
)
\pmb x(k|k)
xxx(k∣k)和未来时刻控制量
u
(
k
)
,
u
(
k
+
1
)
,
.
.
.
,
u
(
k
+
M
?
1
)
u(k),u(k+1),...,u(k+M-1)
u(k),u(k+1),...,u(k+M?1)来预测未来P个时刻的状态。式（1-3）是输出预测模型，根据输出预测模型，我们可以通过当前时刻的状态变量
x
(
k
∣
k
)
\pmb x(k|k)
xxx(k∣k)和未来时刻控制量
u
(
k
)
,
u
(
k
+
1
)
,
.
.
.
,
u
(
k
+
M
?
1
)
u(k),u(k+1),...,u(k+M-1)
u(k),u(k+1),...,u(k+M?1)来预测未来P个时刻的输出。
2、滚动优化
? 根据预测模型中，当前状态变量
x
(
k
∣
k
)
\pmb x(k|k)
xxx(k∣k)是可以测得的，而未来M个时刻的控制量
u
(
k
)
,
u
(
k
+
1
)
,
.
.
.
,
u
(
k
+
M
?
1
)
u(k),u(k+1),...,u(k+M-1)
u(k),u(k+1),...,u(k+M?1)还是未知的，接下来我们通过优化性能指标来求出从当前时刻起的M个控制量
u
(
k
)
,
u
(
k
+
1
)
,
.
.
.
,
u
(
k
+
M
?
1
)
u(k),u(k+1),...,u(k+M-1)
u(k),u(k+1),...,u(k+M?1)。
性能指标有很多种形式，在这里我们不考虑约束，举例两种简单的性能指标表达形式。
2.1 状态优化问题
? 在k时刻的状态优化问题可表述为:确定从该时刻起的M个控制量
u
(
k
)
,
u
(
k
+
1
)
,
.
.
.
,
u
(
k
+
M
?
1
)
u(k),u(k+1),...,u(k+M-1)
u(k),u(k+1),...,u(k+M?1)，使被控对象在其作用下未来P个时刻的状态得到镇定，及趋近
x
=
0
\pmb x=\pmb 0
xxx=000，同时一种控制作用的剧烈变化，优化性能指标可以表达为
m
i
n
J
(
k
)
=
∥
X
(
k
)
∥
Q
2
+
∥
U
(
k
)
∥
R
2
(2-1)
min\pmb J(k)=\begin{Vmatrix} \pmb X(k)\end{Vmatrix}^2_{\pmb Q}+\begin{Vmatrix}\pmb U(k)\end{Vmatrix}^2_{\pmb R}\tag{2-1}
minJJJ(k)=∥∥?XXX(k)?∥∥?Q?Q??Q2?+∥∥?UUU(k)?∥∥?RRR2?(2-1)
其中，
Q
,
R
Q,R
Q,R为状态加权矩阵和控制加权矩阵。接下来把状态预测模型(1-2)带入上式替换掉
X
(
k
)
\pmb X(k)
XXX(k)，然后令
?
J
?
U
=
0
\frac{\partial J }{\partial U}=0
?U?J?=0,如下
?
J
(
k
)
?
U
(
k
)
=
2
G
x
T
Q
[
F
x
x
(
k
∣
k
)
+
G
x
U
(
k
)
]
+
2
R
U
(
k
)
=
2
G
x
T
Q
F
x
x
(
k
∣
k
)
+
2
(
G
x
T
Q
G
x
+
R
)
U
(
k
)
=
0
\begin{aligned}\frac{\partial J(k) }{\partial U(k)}&=2\pmb{G_x^TQ}[\pmb{F_xx}(k|k)+\pmb{G_xU}(k)]+2\pmb{RU}(k) \\&=2\pmb{G_x^TQ}\pmb{F_xx}(k|k)+2(\pmb{G_x^TQ}\pmb{G_x}+\pmb{R})\pmb{U}(k)=0 \end{aligned}
?U(k)?J(k)??=2GxT?Q?GxT?Q??GxT?Q[Fx?x?Fx?x??Fx?x(k∣k)+Gx?U?Gx?U??Gx?U(k)]+2RURURU(k)=2GxT?Q?GxT?Q??GxT?QFx?x?Fx?x??Fx?x(k∣k)+2(GxT?Q?GxT?Q??GxT?QGx??Gx???Gx?+RRR)UUU(k)=0?可以解得
U
(
k
)
=
?
(
G
x
T
Q
G
x
+
R
)
?
1
G
x
T
Q
F
x
x
(
k
∣
k
)
\pmb U(k)=-(\pmb{G_x^TQ}\pmb{G_x}+\pmb{R})^{-1}\pmb{G_x^TQ}\pmb{F_xx}(k|k)
UUU(k)=?(GxT?Q?GxT?Q??GxT?QGx??Gx???Gx?+RRR)?1GxT?Q?GxT?Q??GxT?QFx?x?Fx?x??Fx?x(k∣k)最后我们只取
U
(
k
)
\pmb U(k)
UUU(k)的首项
u
(
k
)
u(k)
u(k)将其实施到系统上,即
u
(
k
)
=
[
1
0
…
0
]
U
(
k
)
=
?
k
T
x
(
k
∣
k
)
,
(2-2)
u(k)=\begin{bmatrix}1&0&\dots&0\end{bmatrix}U(k)=-\pmb{k^Tx}(k|k),\tag{2-2}
u(k)=[1?0?…?0?]U(k)=?kTxkTxkTx(k∣k),(2-2)其中反馈增益
k
T
=
[
1
0
…
0
]
(
G
x
T
Q
G
x
+
R
)
?
1
G
x
T
Q
F
x
\pmb k^T=\begin{bmatrix}1&0&\dots&0\end{bmatrix}(\pmb{G_x^TQ}\pmb{G_x}+\pmb{R})^{-1}\pmb{G_x^TQ}\pmb{F_x}
kkkT=[1?0?…?0?](GxT?Q?GxT?Q??GxT?QGx??Gx???Gx?+RRR)?1GxT?Q?GxT?Q??GxT?QFx??Fx???Fx?
2.2 输出优化问题
? 现在考虑的输出优化问题，即要确定从 k 时刻起的M个控制量 u ( k ) , u ( k + 1 ) , . . . , u ( k + M ? 1 ) u(k),u(k+1),...,u(k+M-1) u(k),u(k+1),...,u(k+M?1)，使被控对象在其作用下未来P个时刻的输出预测值 y ( k + i ) y(k + i) y(k+i)尽可能解近给定的期望值 w ( k + i ) , i = 1 , … , P w(k+ i), i =1,…, P w(k+i),i=1,…,P,同时抑制控制作用的剧烈变化，则可类似地写出性能指标的向量形式 m i n J ( k ) = ∥ W ( k ) ? Y ( k ) ∥ Q 2 + ∥ U ( k ) ∥ R 2 (2-3) min\pmb J(k)=\begin{Vmatrix} \pmb W(k)-\pmb Y(k)\end{Vmatrix}^2_{\pmb Q}+\begin{Vmatrix}\pmb U(k)\end{Vmatrix}^2_{\pmb R}\tag{2-3} minJJJ(k)=∥∥?WWW(k)?YYY(k)?∥∥?Q?Q??Q2?+∥∥?UUU(k)?∥∥?RRR2?(2-3)其中， Q , R Q,R Q,R为输出加权矩阵和控制加权矩阵。接下来把输出预测模型(1-3)带入性能指标中替换 Y ( k ) \pmb Y(k) YYY(k),然后令 ? J ? U = 0 \frac{\partial J }{\partial U}=0 ?U?J?=0,可得 u ( k ) = d T [ W ( k ) ? F y x ( k ∣ k ) ] , (2-4) u(k)=\pmb{d^T}[\pmb W(k)-\pmb{F_yx}(k|k)],\tag{2-4} u(k)=dTdTdT[WWW(k)?Fy?x?Fy?x??Fy?x(k∣k)],(2-4)其中 d T = [ 1 0 … 0 ] ( G y T Q G y + R ) ? 1 G y T Q \pmb{d^T}=\begin{bmatrix}1&0&\dots&0\end{bmatrix}(\pmb{G_y^TQ}\pmb{G_y}+\pmb{R})^{-1}\pmb{G_y^TQ} dTdTdT=[1?0?…?0?](GyT?Q?GyT?Q??GyT?QGy??Gy???Gy?+RRR)?1GyT?Q?GyT?Q??GyT?Q
3、反馈校正
? 由于 x ( k ) \pmb x(k) xxx(k)可测，在每一时刻实测到的 x ( k ) \pmb x(k) xxx(k)可直接用来对该时刻的预测和优化作初始定位，这意味着预测和优化都基于系统的实时反馈信息，从而自然实现了反馈校正，不必再引人额外的校正措施。
二、Matlab仿真示例
例： 有SISO线性时不变系统
x
(
k
+
1
)
=
A
x
(
k
)
+
B
u
(
k
)
y
(
k
)
=
C
x
(
k
)
\pmb x(k+1)=\pmb{Ax}(k)+\pmb Bu(k)\\y(k)=\pmb{Cx}(k)
xxx(k+1)=AxAxAx(k)+BBBu(k)y(k)=CxCxCx(k)其中
A
=
[
1
0.1
0
2
]
,
B
=
[
0
0.5
]
C
=
[
2
1
]
\pmb A=\begin{bmatrix}1&0.1\\0&2\end{bmatrix},\pmb B=\begin{bmatrix}0\\0.5\end{bmatrix}\\\pmb C=\begin{bmatrix}2&1\end{bmatrix}
AAA=[10?0.12?],BBB=[00.5?]CCC=[2?1?]设计控制器使系统输出趋近
y
(
k
)
=
3
y(k)=3
y(k)=3。
注： 该问题是输出优化问题，因此我们通过式（2-4）来计算及时控制量。
参考代码：
close all;
clear all;
k_steps=200;
M=3;P=3;Ts=0.2;
A=[1 0.1;0 2];
B=[0;0.5];
C=[2 1];
x_num=length(A);%状态变量维度
Q=eye(P);%输出权矩阵
R=1;%控制权矩阵
Wk=[3;3;3];%期望值
xk=[0;0];%初始状态
Fy=zeros(P,x_num);
Gy=zeros(P,M);
for i=1:P
Fy(i,1:end)=C*A^i;%计算矩阵Fy
Gy(i,:)=[C*A^(i-1)*B,Gy(max(i-1,1),1:M-1)];%计算矩阵Gy
end
dy=[1 0 0]*(Gy'*Q*Gy+R)^(-1)*Gy'*Q;%计算反馈增益dy，式(2-4)
u_history=zeros(1,k_steps);
y_history=zeros(1,k_steps);
times=zeros(1,k_steps);
%滚动优化
for k=1:k_steps
times(k)=k*Ts;
u_history(k)=dy*(Wk-Fy*xk);%计算及时控制量
y_history(k)=C*xk;%得到真实输出
xk=A*xk+B*u_history(k);
end
figure(1);
plot(times,y_history);
title("输出曲线");
xlabel("t/s")
figure(2);
stairs(times,u_history);
title("控制量变化曲线");
xlabel("t/s");
运行结果：
本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符，请联系我的编程经验分享网邮箱：chenni525@qq.com进行投诉反馈，一经查实，立即删除！
- Python教程
- 深入理解 MySQL 中的 HAVING 关键字和聚合函数
- Qt之QChar编码（1）
- MyBatis入门基础篇
- 用Python脚本实现FFmpeg批量转换
- 【MediaFoundation】OpenCV VideoCapture 读取音频源码
- 电子雨html代码
- Docker：容器化技术的革命者
- HarmonyOS 应用开发学习笔记 ets自定义组件及其引用 @Component自定义组件
- PCA；MDS；TSNE
- 第十六章 ： Spring Boot JWT 集成redis实现分布式token
- 信息安全领域的行业新宠——CISP-DSG
- 数字IC后端设计如何从零基础快速入门？(内附数字IC后端学习视频）
- 将jupyter转换为python文件
- 低代码平台开发 - 物料开发和拓展（二）