席裕庚《预测控制》学习记录 （一）－基于阶跃响应的动态矩阵控制（DMC）包含仿真

这个博客已经过去四年了，相关内容我已经忘得差不多了，这个文章仅供参考，如果有什么问题，建议问问神奇的chatgpt，不要问我，因为大概率我也不知道（更新于2023/10/8）

• 本文除了黑体部分均为《预测控制》一书中的知识点，感觉像之前的ADRC只放个仿真不太好，就把原理也写上了

简介

模型预测控制（Model Predictive Control,MPC）,或简称预测控制。
　　最早产生于工业过程的预测控制算法，有建立在对象脉冲响应基础上的模型预测启发控制（Model Predictive Heuristic Control,MPHC）,或称模型算法控制（Model Algorithmic Control,MAC）,以及建立在对象阶跃响应基础上的动态矩阵控制(Dynamic Matrix Control,DMC)。它们基于对象的非参数模型，用在线不断进行的有限时域优化，即滚动优化，取代了传统最优控制中的一次性全局优化，在滚动的每一步以实时信息进行反馈校正。
　　除了基于对象非参数模型的预测控制算法外，还有基于辨识模型并带有自校正的预测控制算法，如扩展时域预测自适应控制（EPSAC），广义预测控制（GPC）。
　　为解决预测控制系统的设计参数与系统闭环性能之间的关系问题，80年代出现了内模控制（Internal Model Control,IMC）的概念，但得到的结论有限。
　　从90年代人们对预测控制的理论分析从定量转向定性。记住70年代的滚动时域控制的思路，人们先后提出了终端约束预测控制，双模预测控制，在这一思路下，人们提出了众多具有稳定性保证的预测控制算法，进而研究了不确定系统的鲁棒预测控制。。。
　　预测控制的基本方法原理包括，预测控制，滚动优化，反馈校正。

一. 几种典型的预测控制算法

1.基于阶跃响应的动态矩阵控制（ＤＭＣ）

ＤＭＣ算法包括预测模型，滚动优化和反馈较正３个部分。
　　Ａ．预测模型
　　在ＤＭＣ中首先需要测定对象阶跃响应的采样值$a_i=a(iT),i=1,2,\mathrm{3....}$,其中T为采样周期。对于渐进稳定的对象，在某一时刻t=NT后将趋于平稳，以致$a_i(i>N)$与$a_N$的误差已减小到与量化误差及测量误差有相同的数量级，因而可以认为$a_N$已近似等于阶跃响应在$t\to \infty$时的稳态值$a_{\infty }$。这样对象的动态信息就可近似用有限集合 { a 1 , a 2 , . . . , a n } \lbrace a_1,a_2,...,a_n \rbrace {a1​,a2​,...,an​}加以描述。假定我们的控制对象为二阶欠阻尼系统$$G(s)=\frac{1}{0.002s^2+0.04s+1}$$ ,采样周期T=0.02sec,得其阶跃输出输出为
　　
　 取误差小于0.02稳定，得系统调节时间0.4sec。所以可取N=20。得有限集合 { 0.0867 , 0.2917 , 0.5413 , 0.7851 , 0.9818 , 1.1166 , 1.1882 , 1.2059 , 1.1846 , 1.1407 , 1.0886 , 1.0394 , 1.001 , 0.9738 , 0.9603 , 0.9577 , 0.9627 , 0.9721 , 0.9829 , 0.9929 } \lbrace 0.0867,0.2917,0.5413,0.7851,0.9818,1.1166,1.1882,1.2059,1.1846,1.1407,1.0886,1.0394,1.001,0.9738,0.9603,0.9577,0.9627,0.9721,0.9829,0.9929\rbrace {0.0867,0.2917,0.5413,0.7851,0.9818,1.1166,1.1882,1.2059,1.1846,1.1407,1.0886,1.0394,1.001,0.9738,0.9603,0.9577,0.9627,0.9721,0.9829,0.9929}这个向量称为模型向量，N称为建模时域。
　　阶跃模型是一种非参数模型，参数模型和非参数模型分别是指

   用代数方程、微分方程、微分方程组以及传递函数等描述的模型都是参数模型。建立参数模型
就在于确定已知模型结构中的各个参数。通过理论分析总是得出参数模型。非参数模型是直接或间接
地从实际系统的实验分析中得到的响应，例如通过实验记录到的系统脉冲响应或阶跃响应就是非参数
模型。运用各种系统辨识的方法，可由非参数模型得到参数模型。如果实验前可以决定系统的结构，
则通过实验辨识可以直接得到参数模型。

由于线性系统的比例和叠加性质，故可以利用这组参数预测对象未来时刻的输出值。之所以可以预测是因为当k时刻把控制u(k)作用于对象是，相当于在对象输入端加上了一个幅值为u(k)的阶跃。每个时刻的输入都相当一个阶跃输入，我们可以计算在这些输入作用下系统输出的叠加来预测系统的输出。
　　在ｋ时刻，假定控制作用保持不变，是对未来Ｎ个时刻的输出有初始预测值$\widetilde{y_0}(k+i|k),i=1,....N$,例如稳态启动时便可取$\widetilde{y_0}(k+i|k)=y(k)$,在我们之前假设的被控对象的基础上进一步假设，初始条件为0，即$y(0)=0$,输入信号为单位阶跃信号。 在从当前时刻起 M个连续的控制增量$\Delta u(k),...\Delta u(k+M-1)$的作用下，未来时刻的输出值为：
　　
　　当$i<M$时，即只有$i$个控制量作用在系统上，其中的$\widetilde{y_0}(k+i|k)$在运算中来自上一时刻的$\widetilde{y_M}$再加上误差矫正。在之后的反馈校正中再具体讲它的计算。
　　综上，我们知道了什么是阶跃模型，阶跃模型怎么求，知道了预测输出等于一个从反馈校正中得到的值加上M个控制量。M又称为控制时域，在实际使用时自己整定得出的。
　　
　　B滚动优化
　　既然是优化，就要确定出优化的性能指标，或者说是目标函数.
　　(1)DMC是一种通过优化确定控制输入的算法。在每一个时刻k,要确定从该时刻起的M个控制增量作用下的被控对象未来P个时刻的输出预测值$\widetilde{y_M}(k+i|k)$尽可能接近给定的期望值$w(k+i),i=1,...P$。P称为优化时域，通常规定$M\le P\le N$。
　　(2)在控制过程中还希望控制增量$\Delta u(k)$的变化不要太剧烈。
　　所以k时刻的优化性能指标可取为，
　　
　　其中的$q_i$和$r_i$,为加权系数。分别表示对跟踪误差及控制量变化的抑制程度。
　　在不考虑控制量和输出量的约束的情况下，我们将预测模型的关系式（即第一个式子）写为向量形式：
　　　　　　　　　.
　　其中，
　　
　　
　　
　　A是由单位阶跃响应$a_i$组成的$P\times M$阵，称为动态矩阵。向量$\tilde{y}$的前一个下标P即所预测输出个数，后一个下标M为控制量变化的次数。
　　再将性能指标也写为向量形式，即化为范数形式，有关范数可自行搜索。
　　　　其中，
　　　　　　　　　　　　
　　即P个期望输出组成的列向量。
　　　　　　　　　　　　　　　
　　即加权系数组成的对角矩阵。
　　将式（1）代入式（2）得，
　　
　　即
　　
　　 令
　　　　　　　　　　　　
　　即可以得出使$J（k）$取得极值的$\Delta u(k)$.根据矩阵求导公式$\frac{\partial X^{T}AX}{\partial X}=(A+A^T)X$对（4）式求导化简，得
　　　
　　从这个式子中可以得出M个控制增量的最优值，我们只取其即时控制增量即所得列向量的首元素，它可以表示为
　　　　
其中
　　　　　　　　　　
　　　　　　
求出即时控制增量后，实际控制量为：$u(k)=u(k-1)+\Delta u(k)$

根据前述公式我们暂且定义P=20，M=20。加权系数均相等，其和为1，即Q=[0.025,…0.025],R=[0.025,…,0.025]。 输入为单位阶跃输入，阶跃时间为0，所以$w_p(k)=[1,....,1{]}^T$ Matlab下的计算程序如下所示：

 a=[0.0867,0.2917,0.5413,0.7851,0.9818,1.1166,1.1882,1.2059,1.1846,1.1407,
    1.0886,1.0394,1.001,0.9738,0.9603,0.9577,0.9627,0.9721,0.9829,0.9929;];
 A=zeros(20,20);
  for i=1:20
    for j=1:i
    A(i,j)=a(i-j+1);
    end
 end
 q=0.025*ones(1,20);
 r=0.025*ones(1,20);
 Q=diag(q);
 R=diag(r);
 c=zeros(20,1);
 c(1)=1;
 dt=c'*inv(A'*Q*A+R)*A'*Q;

最后得到$d^T=[0.0849，0.0716，0.0468，0.0234，0.0095，0.0033，0.0011，0.0003，0.00010，0，0，0，0，0，0，0，000]$（保留小数后六位）。
　　综上我们知道了如何得出最优控制量，最优控制量等于$d^T$乘上期望输出减初始预测值。

C.反馈校正
　　 当k时刻把控制u(k)作用于对象是，相当于在对象输入端加上了一个幅值为u(k)的阶跃，利用预测模型可以计算出其作用下未来时刻的输出预测值。式（7）为预测模型的向量形式。
　　　　　　　　　
　　其中$\widetilde{y_{N1}}(k)$为在k时刻对k+1时刻的预测，且是没有加上u(k+1)的预测值。所以这个值可以作为k+1时刻的初始预测值进行新的优化计算。即$\widetilde{y_{N0}}(k)$是上一时刻的预测值即$\widetilde{y_{N1}}(k-1)$。然而由于实际中存在模型失配，环境干扰等未知因素，所以由式（7）给出的值会偏离实际值，所以需要利用实时信息进行反馈校正，以防止偏差的不断累积。
　　在DMC在k+1时刻计算优化量前，需先检测对象的实际输出$y(k+1)$与$\widetilde{y_1}(k+1|k)$相比较，构成输出误差：
　　 　　　　　　　　$e(k+1)=y(k+1)-\widetilde{y_1}(k+1|k)－－（８）$
　　 这一误差反映了模型中未包括的不确定因素对输出的影响，可用来预测未来的输出误差，用以补充基于模型的预测。由于对误差的产生缺乏因果性的描述，故误差预测只能采用启发式的方法。

启发式方法指人在解决问题时所采取的一种根据经验规则进行发现的方法。其特点是在解决问题
时,利用过去的经验,选择已经行之有效的方法，而不是系统地、以确定的步骤去寻求答案

　其中，
　　　　　　　　　
　　ｈ为由权系数组成的Ｎ维向量$h=[h_1,.....,h_N{]}^T$称为校正向量。
　　在k+1时刻，预测的未来时间点变为k+2,…,k+1+N, $\widetilde{y_{cor}}(k)$的元素向前移位一个，作为k+1时刻的预测初始值，原本的第N位变为N-1位，此时的第N位用第N-1位近似。这一运算用向量形式表示为：
　　　　　　　　　　　　
　　其中，
　　 　　　　　　　　　　
　　 有了ｋ＋１时刻的初始预测值 $\widetilde{y_{N0}}(k+1)$，就可以进行k+1时刻的优化运算。得到$\Delta u(k+1)$。
　　 我们暂且假定$h=[0.05，....，0.05]$。
　　 综合三个部分的公式原理，我们可以得到如下的算法结构图。

综合仿真

参数初始化程序如下：

 a=[0.0867,0.2917,0.5413,0.7851,0.9818,1.1166,1.1882,1.2059,1.1846,1.1407,1.0886,1.0394,1.001,0.9738,0.9603,0.9577,0.9627,0.9721,0.9829,0.9929;];
 A=zeros(20,20);
 a1=a';
 for i=1:20
    for j=1:i
    A(i,j)=a(i-j+1);
    end
 end
 %q=0.025*ones(1,20);
 %r=0.025*ones(1,20)*40;
 q=[0.5 0.4 0.3 0.3 0.1 0.1 0.08 0.06 0.04 0.02 0.01 0.008 0.006 0.004 0.002 0.001 0.001 0.001 0.001 0.001];
 r=[0.5 0.4 0.3 0.2 0.2 0.1 0.08 0.06 0.04 0.02 0.01 0.008 0.006 0.004 0.002 0.001 0.001 0.001 0.001 0.001]*10.5;
 Q=diag(q);
 R=diag(r);
 c=zeros(20,1);
 c(1)=1;
 c1=c';
 dt=c'/(A'*Q*A+R)*A'*Q;
 h=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]';
 v=ones(1,19);
 S=diag(v,1);
 S(20,20)=1;

simulink下的仿真框图如下：

经过我一下午参数的整定（第一次调DMC，没啥经验），得到的参数就是上面的初始化程序里的q和r。r要比q大很多，系统才不会在开始的时候震荡的很厉害。得到系统的输出如下图，蓝线是没加预测控制，黄线是加了预测控制，可以看到调节时间是差不多的，但是DMC控制没有超调量。我觉得再好好整定一下参数，DMC的输出曲线应该会更好。

matlab2019A下的仿真文件下载链接，
matlab2019A+matlab2012下的仿真文件下载链接, 注意低版本matlab运行可能会出现问题，仿真文件中的w函数可能需要重新定义，w函数如下：

function y = fcn(u)
persistent x;
persistent N;
if isempty(x)
   x=ones(20,1)
end
if isempty(N)
   N=1;
end
x(N)=u;
N=N+1;
if(N==21)
    N=1;
end
y = x;