预测控制1

分数: ___________华北电力大学研究生作业学年学期：课程名称：预测控制学生姓名：学号：提交时间：目录引言 (1)1.预测控制 (2)1.1.单步预测控制 (2)1.2.多步预测控制 (4)1.2.1.对差分方程的思考 (5)1.2.2.差分方程与传递函数关系 (6)2.内膜控制(IMC) (7)2.1.基本原理 (7)2.1.1.模型获取 (7)2.1.2.控制器计算 (7)2.1.3.性质 (8)2.2.基于IMC的PID设计 (9)2.2.1.IMC变形 (9)2.2.2.IMC控制器与PID控制器的联系 (9)2.3.效果 (10)3.动态模型控制（DMC） (10)附录A (13)附录B (16)附录C (20)附录D (22)引言本学期开了新的一门控制学科——预测控制，预测控制是一种基于模型的先进控制技术，它源于工业实践，最大限度地结合了工业实际地要求，并且在实际中取得了许多成功的应用，是一类新型的计算机控制算法。

因为工业生产的过程相当复杂，通过“黑白箱”建立起的模型并不完善，即使使用理论非常复杂的现代控制理论，其控制效果往往也不尽如人意，甚至在某些方面还不及传统的PID 控制。

随后通过老师课堂上对预测控制学科的介绍和讲解，我们渐渐认识到了这门新的学科，它不同与之前学过的PID控制器的思想，也和前馈系统控制思想有很大的不同，预测控制采用多步预测、滚动优化和反馈校正等控制策略实现控制被控对象，适用于控制不易建立精确数字模型并且比较复杂的工业生产过程。

并且由于当今计算机技术的快速发展，Matlab、C、C++的出现极大的方便了我们的理论分析和数学运算。

本课程一方面通过老师课堂上的讲解，另一方面通过课下练习老师布置的习题，一步一步的深入了解了预测控制的原理和运行机制。

1. 预测控制预测控制大致分为以下几个过程。

1) 模型预测对于线性模型有：y (k )=a 1y (k −1)+a 2y (k −2)+bu(k −1)，此模型中u(k −1)为一步固有延迟。

其传递函数可以写为G (s )=kT 1s 2+T 2s+1e −τs （τ=0，T 1，T 2与a 1，a 2有关，k =b1−a 1−a 2）y ̂(k +p )=a 1y ̂(k +p −1)+a 2y ̂(k +p −2)+bu(k +p −1)2) 滚动优化目标函数：J min =∑q i [y r (k +i )−y ̂(k +i)]2P i=1+∑λj Δu 2(k +j −1)M j=1滚动概念：ΔU (k )=[Δu (k )⋯Δu (k +m −1)]T ΔU (k +1)=[Δu (k +1)⋯Δu (k +m )]T3) 反馈校正ê(k +j )=y (k )−y ̂(k |k −j )其中ê(k +j )——预测误差、y (k )——当前时刻输出、y ̂(k |k −j )——第j 个时刻前对当前时刻做的预测y p ̂(k +j |k )=ŷ(k +j |k )+ℎj ê(k +j) 其中y p ̂(k +j |k )——反馈校正之后的预测输出。

4) 参考轨迹y r (k +i )=∂y r (k +i −1)+(1−ð)r(k +j) 1.1. 单步预测控制对于单步预测控制的学习通过例子进行学习，例子如下：对y (k )=a×y (k −1)+b×u (k −l )作d 步预测，极小化如下指标J d=12(y r(k+d)−y p(k+d))2+12λΔu2(k)传递函数为Y(s)U(s)=bs+(1−a)，其中时间常数T s=11−a，因此当a值越小时，系统响应时间会越慢。

此外，b(1−a)×U是系统的最终输出值。

参考轨迹为Y r(k)=α×Y r(k−1) + (1−α)×r(k−1) ( 0<α<1 )，可以理解为Y rr =(1−α)s+(1−α)，此时调节α意味着调节该传递函数的时间参数T r=11−α。

因此，当α接近0时，预测轨迹跟踪r会快。

同时，由于该传递函数为一阶系统，因此参考轨迹不会出现超调震荡的现象。

最后参数λ是控制量增量的约束参数，当λ越大时，目标函数会相应的增大，当然，课堂上老师通过公式推导得出如下公式：dU =b×Yr−a×Y−b×Uλ+b2由上述公式可知，当系统模型确定时，λ越大，控制量的增量dU会越小，这样最直观的现象是控制量U会缓慢的增加。

由于在实际过程中，执行器的控制量的响应速度有限制，同时控制输出也不是无限大，因此从考虑执行器的响应速度的角度分析得出λ可以反映执行器的响应速度。

通过Matlab编程得出如下控制效果图：通过以上各图可以清晰地发现，当系统的时间常数T s较大时（即系统为慢速系统），执行器的控制器输出较为平缓，并且系统输出能够良好的跟踪参考轨迹；当系统是的时间常数T s较小时（即系统为快速系统），执行器控制量的输出会有较大波动，同时系统的输出出现了一定的超调。

通过本次实验发现单步预测控制通过改变离散控制过程中的目标值，使其逐步逼近已设定的目标值，同时利用目标函数得出有极小目标函数值时的控制量增量进而实现控制。

因为，当将参考轨迹的时间常数设为0时，就相当于直接将目标值施加系统中，然后用目标函数计算出控制量的增量即可。

1.2.多步预测控制多步预测控制的学习通过例子进行学习，例子如下：对y(k)=a×y(k −1)+b×u(k −L)作d 步预测，极小化如下指标J d =12(y r (k +d )−y p (k +d ))2+12λΔu 2(k ) 其中:L =1，d =2；L =2，d =2,3。

1.2.1. 对差分方程的思考当L =1时，差分方程为y(k)=a×y(k −1)+b×u(k −1)，可推出其参考轨迹递推式如下：y p (k +1)=a×y (k )+b×u (k )⋯⋯y p (k +d )=a×y p (k +d −1)+b×u (k )化简上述各式可得：y p (k +d )=a d ×y (k )+b×∑a i d−1i=0×u (k ) 又u (k )=u (k −1)+∆u (k )∴ y p (k +d )=a d ×y (k )+b×∑a i d−1i=0×u (k −1)+b×∑a i d−1i=0×∆u(k) (1)由题目知，最小化J d ，且公式(1)中的自变量是∆u ，所以用求导公式令∂J d ∂∆u =0即可。

此时求得：∆u (k )=b×∑a i d−1i=0×(y r (k +d )−a d ×y (k )−b×∑a i d−1i=0×u(k −1))(b×∑a i i=0)2×λ上述公式是L =1时的通式。

同理，可以推导出当L =2时：∆u (k )=b×∑a i d−2i=0×(y r (k +d )−a d ×y (k )−b×∑a i d−1i=0×u(k −1))(b×∑a i i=0)2×λ1.2.2.差分方程与传递函数关系由数字仿真理论可知，差分方程与传递函数存在如下关系式：Y(s)=kT×s+1×U(s)y(k+1)=a×y(k)+bu(k−1)其中：a=e−∆tT，b=k×(1−a)通过Matlab编程得出如下控制效果图：2. 内膜控制(IMC)2.1. 基本原理2.1.1. 模型获取G m (s)可分解为G m−可逆部分（非最小相位，零点在左半平面）和G m+其余部分（稳态增益为1）.2.1.2. 控制器计算G c (s )=G m−−1(s )∙F (s )F(s)=1(T f s+1)n f其中F(s)是为了使G c(s)物理可实现而引入的低通滤波器，n f与G m−有关，T f 是整个控制器时长的10%左右。

例如对于下述被控对象：G m(s)=25s+1e−2s有G m−(s)=25s+1G m+(s)=e−2s2.1.3.性质将内膜控制器构成的闭环控制器进行化简得G B(s)=y(s)r(s)=G p G c1+G c(G p−G m)理想情况下，G m=G p并且G c=G m−1、y(s)=r(s)。

同时实现了稳态无差，因为控制器的稳态增益等于模型稳态增益的逆。

2.2.基于IMC的PID设计2.2.1.IMC变形其中C(s)=G c1−G c G m =G m−(T f s+1)n f⁄1−G m+G m−∙G m−(T f s+1)n f⁄=G m−−1(T f s+1)n f−G m+2.2.2.IMC控制器与PID控制器的联系对于一阶系统G p(s)=k Ts+1G m(s)=k mT m s+a此时G m+=1,G m−=G m,T f=λT m,n f=1，所以C(s)=T m s+1K m(λT m s+1)−1=T m s+1λT m k m s=1λk m(1+1T m s)对于二阶系统G p (s )=1(T 1s +1)(T 2s +1)G m =k m(T m1s +1)(T m2s +1)此时G m+=1,G m−=G m ,T f =λ(T m1+T m2)2,n f =2，所以C (s )=(T m1s +1)(T m2s +1)k m (λT m s +1)2−1=T m1T m2s 2+(T m1+T m2)s +1k m ((λT m s )2+2λT m s ) (1)当n f =1时C (s )=T m1T m2s 2+(T m1+T m2)s +1λT m k m s =T m1+T m2λT m k m +1λT m k m s +T m1T m2s λT m k m(2)观察(1)、(2)式，发现n f =2时相当于串联了一阶滤波器。

2.3. 效果以一阶控制系统为例，通过Matlab 编程得出如下控制效果图：3. 动态模型控制（DMC ）动态矩阵预测控制算法是一种基于系统的单位阶跃响应序列非参数模型的控制策略。

动态矩阵通过被控对象的阶跃响应来建立系统的非参数数学模型，可用于渐进稳定的线性系统。

动态矩阵预测控制算法通常三部分：预测模型、滚动优化、反馈校正。

动态矩阵控制算法结构图图3-1如所示。

图3-1 DMC控制算法结构图运行中以加热炉为例进行学习，通过Matlab编程得出如下控制效果图：输出结果如下:Percent error in the last step response coefficientof output yi for input uj is∶44%Percent error in the last step response coefficient of output yi for input uj is∶44%Time remaining 500/500Time remaining 400/500Time remaining 300/500Time remaining 200/500Time remaining 100/500Time remaining 0/500Simulation time is 0.032 seconds.附录A% ============================================ % 单步预测控制% by: 2017−03−01% ============================================ % function PreSP =3;% 仿真时间100ST =0.1;% 仿真步长LP =SP/ST;% 仿真循环次数% 模型参数% Y(k)=a∗Y(k−1) +b∗U(k)% Y(s)/U(s)=b/(s+(1−a))a =0.5;b =2;r =1;% 目标值L =1;% 延迟环节系数d =1;% 预测步长% 参考轨迹方程Yr =Alpha∗Yr + (1−Alpha)∗r ( 0<Alpha<1 )% Alpha 接近1 参考轨迹可信度越高% 接近0 原目标值可信度越高Tr =2;% 参考轨迹时间参数Ts =1/(1−a);% 模型时间常数Alpha =exp(−Ts/Tr);% 参考轨迹参数% Lamda 为目标函数对控制量增量约束的参数，是一个给定参数Lambda =0.5;% 控制量增量约束参数% 初始化参数Yo =zeros(1,LP+1);% 模型输出值YrRes =zeros(1,LP+1);% 参考轨迹Uo =zeros(1,LP+1);% 控制量输出T =zeros(1,LP+1);% 仿真时间% 初值初始化Y0 =0;U0 =0;dU =0;Y =Y0;U =U0 +dU;Yr =Y0;for i=1:LPYr =Alpha∗Yr+(1−Alpha)∗r;% 参考轨迹dU =b∗(Yr−a∗Y−b∗U)/(Lambda+b∗b);% 控制量增量U =U +dU;% 控制量Y =a∗Y +b∗U;% 模型输出% 存储历史值T(i+1)=i;Yo(i+1)=Y;Uo(i+1)=U;YrRes(i+1)=Yr;end% 图形显示plot(T,Yo,′ro−′,T,Uo,′k.−′,T,YrRes,′∗−′);legend(′模型输出Yo′,′控制量输出Uo′,′参考轨迹Yr′);% ShowTitle =sprintf(′a =%1.3f,b =%3.3f,\alpha =%1.3f,\Lamda =%4.3f′,...% a,b,Alpha,Lamda);title([′a =′,num2str(a),′,b =′,num2str(b),...′,\alpha =′,num2str(Alpha),′,\lambda =′,num2str(Lambda)]);grid minor附录B% ============================================ % 多步预测控制% by: 2017−03−07% ============================================ function Pre0307ST =3;% 仿真时间100DT =0.1;% 仿真步长LP =ST/DT;% 仿真循环次数R =1;% 目标值% 模型参数% Y(k)=a∗Y(k−1) +b∗U(k)% a =exp(−DT/Ts) −>Ts =−DT/ln(a)% b =K∗(1−a) −>K =b/(1−a)% Y(s)/U(s)=K/(Ts∗s+1)% ============================================ % a =0.5;% Ts =−DT/log(a);% 模型时间常数b =2;% k =b/(1−a);% 模型增益% 为便于结合传递函数模型，将上述表达式进行改写如下% ============================================ Ts =10;% 模型时间常数a =exp(−DT/Ts);% k =200;% b =k∗(1−a);% ============================================ L =1;% 延迟环节系数d =1;% 预测步长% 参考轨迹方程Yr =Alpha∗Yr + (1−Alpha)∗r ( 0<Alpha<1 )% Alpha 接近1 参考轨迹可信度越高% 接近0 原目标值可信度越高Tr =10;% 参考轨迹时间参数Alpha =exp(−Ts/Tr);% 参考轨迹参数% Lamda 为目标函数对控制量增量约束的参数，是一个给定参数Lambda =0.7;% 控制量增量约束参数% 初始化参数Yo =zeros(1,LP+1);% 模型输出值YrRes =zeros(1,LP+1);% 参考轨迹Uo =zeros(1,LP+1);% 控制量输出T =zeros(1,LP+1);% 仿真时间Uc =zeros(1,L);% 初值初始化Y0 =0;Y =Y0;Yr =Y0;% 求取控制量增量所用参数A =0;if L ==1for i=0:d−1A =A+power(a,i);endA =A∗b;B =A;elseif L==2for i=0:d−2A =A+power(a,i);endB =b∗(A+power(a,d−1));A =A∗b;endfor i=1:d−1Yr =Alpha∗Yr+(1−Alpha)∗R;% 参考轨迹初始化endYrRes(1)=Yr;for i=1:LPYrRes(i+1)=Alpha∗YrRes(i)+(1−Alpha)∗R;% 参考轨迹dU =A∗(YrRes(i+1)−power(a,d)∗Y−B∗Uc(1))/...(Lambda+A∗A);% 求取控制量增量通式if L>1Uc(2:L)=Uc(1:L−1);% 延迟环节endUc(1)=Uc(1) +dU;% 控制量Y =a∗Y +b∗Uc(L);% 模型输出% 存储历史值T(i+1)=i∗DT;Yo(i+1)=Y;Uo(i+1)=Uc(L);end% 图形显示plot(T,Yo,′ro−′,T,Uo,′k.−′,T,YrRes,′∗−′);legend(′模型输出Yo′,′控制量输出Uo′,′参考轨迹Yr′); title([′\lambda =′,num2str(Lambda),...′,T_s =′,num2str(Ts),...′,T_r =′,num2str(Tr)]);grid minor附录Cfunction [yimc]=yijieIMC(alfa1,k,t1,r1,dt,st,c1,c2)K=k;T=t1;r=r1;DT=dt;ST=st;LP=ST/DT;alfa=alfa1;x(1:4)=0;e1=0;e2(1:LP)=0;for i=1:LPe1=r−x(4);e2(i)=e1+x(3);if i>1x(1)=e2(i)∗(DT+T)/K/DT−T/DT/K∗e2(i−1);x(2)=x(2)+DT∗(−x(2)+x(1))/T/alfa;x(3)=x(3)+DT∗(−x(3)+c1∗K∗x(2))/T/c2;x(4)=x(4)+DT∗(−x(4)+K∗x(2))/T;yimc(i)=x(4);t(i)=i∗DT;endendfunction [ypid,u,t]=yijiePID(alfa1,k,t1,r1,dt,st,c1,c2)K=k;T=t1;r=r1;DT=dt;ST=st;LP=ST/DT;alfa=alfa1;x1=0;xp=0;xi=0;u0=0;for i=1:LPe=r−x1;xp=e/alfa/K/c1;xi=xi+DT∗e/T/c1/alfa/K/c2;u0=xp+xi;x1=x1+DT∗(−x1+K∗u0)/T;ypid(i)=x1;t(i)=i∗DT;u(i)=u0;end%% 一阶系统IMC−PIDclc;clear all;close all;K=5;T=3;r=1;DT=1;ST=500;alfa=18;c1=1;c2=1;[yimc]=yijieIMC(alfa,K,T,r,DT,ST,c1,c2); [ypid,u,t]=yijiePID(alfa,K,T,r,DT,ST,c1,c2); plot(yimc)hold on;plot(ypid,′r′)附录D% 加热炉DMC设计Km = 2.5;Tp1 =20;Tp2 =120;dp =80;Kp =Km;Tm1 =Tp1;Tm2 =Tp2;dm =dp;% modelGm =poly2tfd(Kp,conv([Tp1,1],[Tp2,1]),0,dp);% Sample Time Ts.Model length N.Output number ny Ts =10;N =200;ny =1;% Stop Response Coefficientsmodel =tfd2step(N,Ts,ny,Gm);plant =tfd2step(N,Ts,ny,Gp);% else Parametorsp =20;M =5;Q =1;R =1;Kmpc =mpccon(model,Q,R,M,p);%Step point and simulationyr =1;tend =500;[y,u,ym]=mpcsim(plant,model,Kmpc,tend,yr); plotall(y,u,Ts);。