【预测控制4】基于状态空间模型的预测控制

---

相关文章:
【预测控制1】模型预测控制概论
【预测控制2】预测控制基本原理
【预测控制3】动态矩阵控制（DMC）及Matlab仿真

---

一、基于状态空间模型的预测控制

1、预测模型

  首先考虑SISO线性系统$$\begin{gathered} x(k+1)=Ax(k)+bu(k) \\ y(k)=c^{T}x(k)\text{\hspace{1em}(1-1)} \end{gathered}$$
状态变量$x(k)\in R^n$实时可测，$y(k)、u(k)$分别为系统的输入输出。
  设从$k$时刻起系统输入发生$M$步变化，而后保持不变，则由模型(1-1)，可以预测输出在$u(k),u(k+1),...,u(k+M-1)$作用下未来P个时刻的系统状态$$\begin{array}{rl}x(k+1|k)& =Ax(k|k)+bu(k)\\ x(k+2|k)& =Ax(k+1|k)+bu(k)=A^{2}x(k|k)+Abu(k)+bu(k+1)\\ ⋮& \\ x(k+M|k)& =A^{M}x(k|k)+A^{M-1}bu(k)+A^{M-2}bu(k+1)+\cdots +bu(k+M-1)\\ x(k+M+1|k)& =A^{M+1}x(k|k)+A^{M}bu(k)+A^{M-1}bu(k+1)+\cdots +(Ab+b)u(k+M-1)\\ ⋮& \\ x(k+P|k)& =A^{P}x(k|k)+A^{P-1}bu(k)+A^{P-2}bu(k+1)+\\ & \cdots +(A^{P-M}b+A^{P-M-1}b+\cdots +b)u(k+M-1)\end{array}$$
将上面的式子整合为向量形式，有$$X(k)=F_{x}x(k|k)+G_{x}U(k)\text{\hspace{1em}(1-2)}$$其中$$\begin{gathered} F_x={\left[\begin{array}{c}A\\ A^2\\ ⋮\\ A^P\end{array}\right]}_{nP\times 1},Gx={\left[\begin{array}{ccccc}b& 0& \dots & 0& 0\\ Ab& b& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ A^{P-1}b& A^{P-2}b& \dots & A^{P-M+1}b& {\sum }_{i=0}^{P-M}{A}^{i}b\end{array}\right]}_{nP\times M} \\ X(k)={\left[\begin{array}{c}x(k+1|k)\\ x(k+2|k)\\ ⋮\\ x(k+P|k)\end{array}\right]}_{nP\times 1},U(k)={\left[\begin{array}{c}u(k)\\ u(k+1)\\ ⋮\\ u(k+M-1)\end{array}\right]}_{M\times 1} \end{gathered}$$
类似的，根据输出方程$$y(k)=c^{T}x(k)$$可以得到$$Y(k)=F_{y}x(k|k)+G_{y}U(k)\text{\hspace{1em}(1-3)}$$其中$$\begin{gathered} F_y={\left[\begin{array}{c}{c}^{T}A\\ c^T{A}^2\\ ⋮\\ c^T{A}^P\end{array}\right]}_{P\times n},Gy={\left[\begin{array}{ccccc}{c}^{T}b& 0& \dots & 0& 0\\ c^{T}Ab& c^{T}b& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {c^{T}A}^{P-1}b& {c^{T}A}^{P-2}b& \dots & {c^{T}A}^{P-M+1}b& {\sum }_{i=0}^{P-M}{c^{T}A}^{i}b\end{array}\right]}_{P\times M} \\ Y(k)={\left[\begin{array}{c}y(k+1|k)\\ y(k+2|k)\\ ⋮\\ y(k+P|k)\end{array}\right]}_{P\times 1} \end{gathered}$$   在此总结一下，式（1-2）是状态预测模型，根据状态预测模型，我们可以通过当前时刻的状态变量$x(k|k)$和未来时刻控制量$u(k),u(k+1),...,u(k+M-1)$来预测未来P个时刻的状态。式（1-3）是输出预测模型，根据输出预测模型，我们可以通过当前时刻的状态变量$x(k|k)$和未来时刻控制量$u(k),u(k+1),...,u(k+M-1)$来预测未来P个时刻的输出。

2、滚动优化

  根据预测模型中，当前状态变量$x(k|k)$是可以测得的，而未来M个时刻的控制量$u(k),u(k+1),...,u(k+M-1)$还是未知的，接下来我们通过优化性能指标来求出从当前时刻起的M个控制量$u(k),u(k+1),...,u(k+M-1)$。
性能指标有很多种形式，在这里我们不考虑约束，举例两种简单的性能指标表达形式。

2.1 状态优化问题

  在k时刻的状态优化问题可表述为:确定从该时刻起的M个控制量$u(k),u(k+1),...,u(k+M-1)$，使被控对象在其作用下未来P个时刻的状态得到镇定，及趋近$x=0$，同时一种控制作用的剧烈变化，优化性能指标可以表达为$$minJ(k)={\Vert \begin{array}{c}X(k)\end{array}\Vert }_Q^2+{\Vert \begin{array}{c}U(k)\end{array}\Vert }_R^2\text{\hspace{1em}(2-1)}$$
其中，$Q,R$为状态加权矩阵和控制加权矩阵。接下来把状态预测模型(1-2)带入上式替换掉$X(k)$，然后令$\frac{\partial J}{\partial U}=0$,如下$$\begin{array}{rl}\frac{\partial J(k)}{\partial U(k)}& =2G_x^{T}Q[F_{x}x(k|k)+G_{x}U(k)]+2RU(k)\\ & =2G_x^{T}Q{F}_{x}x(k|k)+2(G_x^{T}Q{G}_x+R)U(k)=0\end{array}$$可以解得$$U(k)=-(G_x^{T}Q{G}_x+R{)}^{-1}{G}_x^{T}Q{F}_{x}x(k|k)$$最后我们只取$U(k)$的首项$u(k)$将其实施到系统上,即$$u(k)=\left[\begin{array}{cccc}1& 0& \dots & 0\end{array}\right]U(k)=-k^{T}x(k|k),\text{\hspace{1em}(2-2)}$$其中反馈增益$$k^T=\left[\begin{array}{cccc}1& 0& \dots & 0\end{array}\right](G_x^{T}Q{G}_x+R{)}^{-1}{G}_x^{T}Q{F}_x$$

2.2 输出优化问题

  现在考虑的输出优化问题，即要确定从 k 时刻起的M个控制量$u(k),u(k+1),...,u(k+M-1)$，使被控对象在其作用下未来P个时刻的输出预测值$y(k+i)$尽可能解近给定的期望值$w(k+i),i=1,\dots ,P$,同时抑制控制作用的剧烈变化，则可类似地写出性能指标的向量形式$$minJ(k)={\Vert \begin{array}{c}W(k)-Y(k)\end{array}\Vert }_Q^2+{\Vert \begin{array}{c}U(k)\end{array}\Vert }_R^2\text{\hspace{1em}(2-3)}$$其中，$Q,R$为输出加权矩阵和控制加权矩阵。接下来把输出预测模型(1-3)带入性能指标中替换$Y(k)$,然后令$\frac{\partial J}{\partial U}=0$,可得$$u(k)=d^T[W(k)-F_{y}x(k|k)],\text{\hspace{1em}(2-4)}$$其中$$d^T=\left[\begin{array}{cccc}1& 0& \dots & 0\end{array}\right](G_y^{T}Q{G}_y+R{)}^{-1}{G}_y^{T}Q$$

3、反馈校正

  由于$x(k)$可测，在每一时刻实测到的$x(k)$可直接用来对该时刻的预测和优化作初始定位，这意味着预测和优化都基于系统的实时反馈信息，从而自然实现了反馈校正，不必再引人额外的校正措施。

二、Matlab仿真示例

例： 有SISO线性时不变系统$$\begin{gathered} x(k+1)=Ax(k)+Bu(k) \\ y(k)=Cx(k) \end{gathered}$$其中$$\begin{gathered} A=\left[\begin{array}{cc}1& 0.1\\ 0& 2\end{array}\right],B=\left[\begin{array}{c}0\\ 0.5\end{array}\right] \\ C=\left[\begin{array}{cc}2& 1\end{array}\right] \end{gathered}$$设计控制器使系统输出趋近$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");

运行结果：