DMC算法MATLAB编程及仿真

1一、前言工业生产的过程是复杂的，建立起来的模型也是不完善的。即使是理论非常复杂的现代控制理论，其效果也往往不尽人意，甚至在一些方面还不及传统的PID控制。20世纪70年代，人们除了加强对生产过程的建模、系统辨识、自适应控制等方面的研究外，开始打破传统的控制思想，试图面向工业开发出一种对各种模型要求低、在线计算方便、控制综合效果好的新型算法。在这样的背景下，预测控制的一种，也就是动态矩阵控制（DMC）首先在法国的工业控制中得到应用。因此预测控制不是某种统一理论的产物，而是在工业实践中逐渐发展起来的。预测控制中比较常见的三种算法是模型算法控制（MAC），动态矩阵控制（DMC）以及广义预测控制。本篇所采用的是动态矩阵控制，其采用增量算法，因此在消除稳态余差方面非常有效。二、控制系统设计方案2.1控制系统方案设计图动态矩阵控制是基于系统阶跃响应模型的算法，隶属于预测控制的范畴。它的原理结构图如下图(1)所示：图(1)预测控制原理结构图上图就是预测控制原理结构图，从图中我们可以看到，预测控制的主要特点。即建立预测模型；采用滚动优化策略，采用模型误差反馈矫正。这也是预测控制的本质所在。下面将对这三个特点一一说明。2.2预测控制基本原理1、预测模型：预测模型的功能是根据对象历史信息和未来输入对对象输出进行预测，它是被控对象的准确模型。这里只强调模型的功能而不强调其结构形式。因此，预测模型可以是状态方程、传递函数等传统的参数模型，对于线性稳2定对象，阶跃响应、脉冲响应这类非参数模型也可以作为预测模型使用。预测模型具有展示系统未来动态行为的功能，这样，就可以利用预测模型来预测未来时刻被控对象的输出变化及被控变量与其给定值的偏差，作为控制作用的依据，使之适应动态系统所具有的因果性的特点，得到比常规控制更好的控制效果。2、滚动优化：预测控制是一种优化控制算法，它是通过某一性能指标的最优来确定未来的控制作用。这一性能指标涉及到系统未来的行为。如，通常可取被控对象输出在未来的采样点上跟踪某一期望轨迹的方差最小，性能指标中涉及到的系统未来的行为，是根据预测模型由未来的控制策略决定的。但是，预测控制中的优化与传统意义下的离散最优控制有很大的区别，这主要表现在预测控制中的优化是一种有限时段的滚动优化。每一采样时刻，优化性能指标只涉及到从该时刻起未来的有限时段，到下一采样时刻，这一优化时段同时向前推移。因此，预测控制不是用一个对全局相同的优化性能指标，而是在每一采样时刻有一个相对于该时刻的优化性能标。不同时刻优化性能指标的相对形式是相同的，但其绝对形式，即，所包含的时间区域是不同的。在预测控制中，优化不是一次离线进行，而是反复在线进行，这就是滚动优化的含义。虽然各个不同时刻的优化指标的相对形式是相同的，每一步实现的是静态参数的优化，但因为其包含的时间区域不同，从控制的全过程看，实现的是动态的优化。因而预测控制被称为是“滚动优化”的算法。与一般最优控制中的全局优化相比，预测控制中的滚动优化只能得到全局的次优解，但由于它的优化始终建立在实际过程的基础上，使控制结果达到实际意义上的最优控制，能够有效地克服工业过程控制中的模型不精确、非线性、时变等不确定性的影响。3、反馈校正：预测控制是一种闭环控制算法，预测算法在进行滚动优化时，优化的基点应与系统实际一致。但作为基础的预测模型，只是对象动态特性的粗描述，由于实际系统中存在的非线性、模型失配、干扰等因素，基于不变模型的预测不可能和实际情况完全相符，这就需要用附加的预测手段补充模型预测的不足，或者对基础模型进行在线修正。滚动优化只有建立在反馈校正的基础上，才能体现出其优越性。因此，预测控制算法在通过优化确定了一系列未来的控制作用后，为了防止模型失配或环境干扰引起控制对理想状态的偏离，并不是把这些控制作用逐一全部实施，而只是实现本时刻的控制作用。到下一采样时刻，则首先检测对象的实际输出，并利用这一实时信息、对基于模型的预测进行修正，然后再进行新的优化。反馈校正的形式是多样的，可以在保持预测模型不变的基础上，对未来的误差作出预测并加以补偿，也可以根据在线辨识的原理直接修改预测模型。不论取何种校正形式，预测控制都把优化建立在系统实际的基础上，并力图在优化时对系统未来的动态行为作出较准确的预测。因此，预测控制中的优3化不仅基于模型，而且利用了反馈信息，因而构成了闭环优化。2.3DMC控制策略DMC算法的基本原理DMC（DynamicMatrixControl）算法是基于被控对象单位阶跃响应非参数模型的预测控制算法。通过反馈校正和滚动优化当前和未来时刻的控制量，使得输出响应符合预先设定的轨迹。DMC动态矩阵控制是一种用被控对象的阶跃响应特性来描述系统动态模型的预测控制算法。三、模型建立3.1模型和辨识算法的选择已知加热炉温度被控系统的模型可以表达为如下形式，且已知参数范围：1K5,100T2500,50T1200,8.5D11.5现从生产现场采集到的输入输出运行数据存放在文件u和yi中（采样间隔是5秒）。DsesTsTKsG1)(122辨识算法分别用一次最小二乘算法和递推最小二乘算法对被控对象进行辨识。最小二乘法是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。3.2辨识算法介绍1、一次最小二乘算法数学原理：观测矩阵H和LL=Tnlll)](),.....,4(),3([H=)2(...)2()1()1(...)3()2()2(...)2()1()1(...)3()2(nuuunuuunlllnlll被辨识参数矩阵C=Tbbaa]2,1,2,1[公式1：LHHHCTT1)(用MATLAB程序编程如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%loady22;loadu;%导入输入输出序列u=u;y=y22;%给u，y赋值D=1;%设置纯迟延为14N=length(u);%设置步长fori=D+3:NH1(i-D-1,:)=[-y(i-D),-y(i-D-1),u(i-D-1),u(i-D-2)];end%此循环为求观测矩阵HH3=H1'*H1;H4=pinv(H3);%矩阵求逆C1=H4*H1'*y(D+2:N)';%由公式1求得被辨识参数矩阵Ca11=C1(1),a21=C1(2),b11=C1(3),b21=C1(4),%分离参数z=tf('z',5);%设置算子Gz1=(b11*z+b21)/(z^2+a11*z+a21)%系统的离散时间域模型Gc1=d2c(Gz1)%系统的连续时间域模型%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%其程序的流程框图如图(2)所示：辨识结果如下：a11=-1.1505a21=0.1779b11=0.0365b21=0.0208Transferfunction:0.03651z+0.02077----------------------z^2-1.151z+0.1779Samplingtime:5Transferfunction:1.375e-013s+0.004834-------------------------s^2+0.3453s+0.002314图(2)一次最小二乘辨识程序流程框图2、递推最小二乘算法数学原理：公式(2)：1])(1)()1()()[()1()(iihiPihihiPiKT5公式(3)：)]1()()()[()1()(iCihkliKiCiCT公式(4)：])(1)()1()()[()()1()(iihiPihiKiKiPiPTT注：其中的H阵和L阵形式见一次辨识算法部分。利用MATLAB程序编程如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%loady22;loadu;%导入输入输出数据u=u;y=y22;D=1;N=length(u);%数据赋值，设置步长C1=[0.001,0.001,0.001,0.001];%给被辨识参数赋初值P1=10^6*eye(4,4);%给初始状态赋值e=10^(-6);%给收敛判定校对的误差赋值Cx=[C1;zeros(N-4,4)];%设置辨识参数矩阵并赋初值E=zeros(N-3,4);%设置相对误差矩阵并赋初值fori=D+3:N%由公式2,3,4求K，C，P阵H1=[-y(i-D-1),-y(i-D-2),u(i-D-1),u(i-D-2)];x1=H1*P1*H1'+1;x2=pinv(x1);%加权矩阵是单位阵。K1=P1*H1'*x2;%由公式2求出K值t1=y(i)-H1*C1';C2=C1+(K1*t1)';%由公式3求出被辨识参数E1=C2-C1;%求参数当前值与上一次值的误差E2=E1./C1;%求参数的相对变化E3=0;%给收敛判定判定误差赋初值E3=abs(E2(1,1))+abs(E2(1,2))+abs(E(1,3))+abs(E(1,4));%求收敛判定误差E(i-D-1,:)=E2;%把当前相对变化的行向量加入误差矩阵的最后一行C1=C2;%把新获得的参数作为下一次递推的旧参数Cx(i-D-1,:)=C2;%把辨识参数C的行向量加入辨识参数矩阵的最后一行P2=P1-K1*K1'*(H1*P1*H1'+1);%由公式4求出P(i)的值P1=P2;%给下次用ifE3=ebreak;%如果参数收敛满足要求，终止计算endend6a1x=Cx(:,1)';a2x=Cx(:,2)';b1x=Cx(:,3)';b2x=Cx(:,4)';%分离辨识参数矩的列向量fori=4:Nifa1x(i)==0&&b1x(i)==0&&a2x(i)==0&&b2x(i)==0a12=a1x(i-1);a22=a2x(i-1);b12=b1x(i-1);b22=b2x(i-1);break;endend%此处循环选择语句是从参数列向量中分离出最终所需要的辨识参数a12,a22,b12,b22%显示辨识参数z=tf('z',5);%设置算子Gz2=(b12*z+b22)/(z^2+a12*z+a22)%系统的离散时间域模型Gc2=d2c(Gz2)%系统的连续时间域模型%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序设计的流程图如图(3)所示：辨识结果如下所示：a12=-1.5294a22=0.5759b12=0.0371b22=0.0600Transferfunction:0.03715z+0.06001----------------------z^2-1.529z+0.5759Samplingtime:5Transferfunction:7-0.004176s+0.005082-------------------------s^2+0.1104s+0.002436图(3)递推最小二乘算法程序流程框图四、动态矩阵控制算法（DMC）预测模型：)1()()1(0kUAkUAkYmTmmmmkpkykkykkykY)]/(,.......,/)2(),/)1([)1(TMkukukukU)]1(,),........1(),([)(TkuNkuNkukU)]1(),.......,2(),1([)1(其中0A是由阶跃响应系数组成的矩阵。并且以递推辨识参数作为预测模型参数，而且递推辨识参数作为被控对象参数。基于MATLAB程序M文件编程如下所示：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%D=1;%纯迟延CD=1;%采样时间8g=poly2tfd([b12,b22],[1,a12,a22],CD,D);%将预测模型转换为阶跃响应模型p=poly2tfd(