模型参考自适应控制—MIT法

北京化工大学信息学院[e-mail]1一原理及方法模型参考自适应系统，是用理想模型代表过程期望的动态特征，可使被控系统的特征与理想模型相一致。一般模型参考自适应控制系统的结构如图1所示。图1一般的模型参考自适应控制系统其工作原理为，当外界条件发生变化或出现干扰时，被控对象的特征也会产生相应的变化，通过检测出实际系统与理想模型之间的误差，由自适应机构对可调系统的参数进行调整，补偿外界环境或其他干扰对系统的影响，逐步使性能指标达到最小值。基于这种结构的模型参考自适应控制有很多种方案，其中由麻省理工学院科研人员首先利用局部参数最优化方法设计出世界上第一个真正意义上的自适应控制律，简称为MIT自适应控制，其结构如图2所示。图2MIT控制结构图系统中，理想模型Km为常数，由期望动态特性所得，被控系统中的增益Kp在外界环境发生变化或有其他干扰出现时可能会受到影响而产生变化，从而使其动态特征发生偏离。而Kp的变化是不可测量的，但这种特性的变化会体现在广义误差e上，为了消除或降低由于Kp的变化造成的影响，在系统中增加一个可调增益Kc，来补偿Kp的变化，自适应机构的任务即是依据误差最小指标及时调整Kc，使得Kc与Kp的乘积始终与理想的Km一致，这里使用的优化方法为最优梯度法，自适应律为：tmdyeBKctKc0)0()(YpYme+__+R参考模型调节器被控对象适应机构可调系统———kmq(s)p(s)KcKpq(s)-----p(s)适应律Rymype+-北京化工大学信息学院[e-mail]2MIT方法的优点在于理论简单，实施方便，动态过程总偏差小，偏差消除的速率快，而且用模拟元件就可以实现；缺点是不能保证过程的稳定性，换言之，被控对象可能会发散。二对象及参考模型该实验中我们使用的对象为：122)()()(2ssspsqKsGpp参考模型为：121)()()(2ssspsqKsGmm用局部参数最优化方法设计一个模型参考自适应系统，设可调增益的初值Kc(0)=0.2，给定值r(t)为单位阶跃信号，即r(t)=A×1(t)。A取1。三自适应过程将对象及参考模型离散化，采样时间取0.1s，进而可得对象及参考模型的差分方程分别为：)2(0044.0)1(0047.0)2(8187.0)1(8079.1)(krkrkykykym)2(0088.0)1(0094.0)2(8187.0)1(8097.1)(kukukykykyppp其中u为经过可调增益控制器后的信号。编程进行仿真，经大量实验发现，取修正常数B为0.3，可得较好的动态过度过程，如下图3所示：北京化工大学信息学院[e-mail]3图3仿真结果由图3中第一个图形可以看出，在阶跃扰动后，经过一段时间对象的输出完全跟踪上了理想模型的值，系统最终趋于稳定；由第二个图可以看出，当系统稳定后，Kp*Kc等于Km，说明补偿环节达到了期望的补偿效果，这与系统设计的目标一致；由第三个图可以看出，在控制的动态过程中，偏差的总和是比较小的，而且偏差的消除是很快的，这是由于所选用的优化方法为最有梯度法的结果。在1中我们已经得到一个能使对象得到较好控制的参数B=0.3，在此情况下，我们将Kp取为1，对应于实际中即指对象增益发生漂移，再做仿真，结果如图4所示。图4对象增益变化后的仿真图由图4我们可以看出，在一个适当的修真参数B下，当对象的特性参数Kp发生漂移后，控制器依然能很好的控制对象，这也证明了MIT方法的自适应特性。而且我们发现，当Kp由2变为1后，控制器的控制效果更好了，具体表现为振荡减弱，过渡过程有所加快，关于导致这一现象发生的原因，我们会在第四部分中做详细的分析与说明。四研究分析1对于一个被控过程，系统能稳定运行是设计与控制的首要指标，然而如前所述，依据最优控制的原则设计出来的MIT自适应控制器却可能会使得系统不稳定，输出发散，以下我们对此做一研究，以期找出其中的相关信息。我们设某连续二阶对象为：北京化工大学信息学院[e-mail]41)()()(122sbsbKpspsqKsGp则有：pmmmmmcpppppyyeRKyybybRKKuKyybyb1212控制律为：mBeycKR为一阶跃信号，即R(t)=A×1(t)，则偏差的动态方程为：0212eAKBKeebebmp根据劳斯稳定判据，列出劳斯行列式：10101221121223sbAKBKbbsAKBKbsbsmpmp得知，对于该连续系统，当212/bbAKBKmp时会不稳定。试验中Kp=2，Km=1，A=1，1,221bb，因而对于连续系统，可求得当B=1时，系统将会等幅振荡。现取B=1，得仿真曲线如图5所示。图5显示，当B为1的时候系统发散，另取原使系统稳定的B=0.3，计算出此时可使系统振荡的阶跃幅值A=sqrt(1/0.3)做仿真，结果如图6所示：图6显示结果与图5一样，系统也发散。图5、6过程中所取参数均为由劳斯判据所得临界值，然而系统并未做等幅振荡，而是发散，这似乎使得理论计算与仿真结果不符。但稍作分析我们就会发现，问题在于我们仿真时用的是离散化的模型，而所用参数为由连续系统计算所得。我们知道用连续系统分析的结论是不完全适用于离散系统的，这是因为随着采样时间取不同的值，同一对象的连续特性和离散特性会不同。因而对于离散系统，我们对其做稳定性分析时还需考虑采样时间的影响。正确的做法应是：将连续开环对象做Z变换，进而得到闭环的Z域特征方程，对此方程做双线性变化，然后对所得w域方程列出其劳斯阵列，应用劳斯稳定判据即可得到使离散系统做等幅振荡的相关参数。本实验中广义偏差方程为三阶系统，在应用采样系统的劳斯稳定判据时需要求解含有参变量的三解方程的解析解，运算量较大，因而这里未做相应的求解。只是对其做一些定性的分析，指出对于同一对象，使得连续系统和离散系统做等幅振荡的参数B是不一样的，因而仿真的结果并没有问题。北京化工大学信息学院[e-mail]5在对离散系统进行大量的仿真实验后发现当B取0.8367左右的时候，离散系统会发生等幅振荡。如图7所示：图5B=1时的仿真结果图6A=sqrt(1/0.3)时的仿真北京化工大学信息学院[e-mail]6图7离散系统等幅振荡2为了更进一步的了解该实验的相关特征，我们设计以下实验，来分别研究该仿真中Kc的初始值、阶跃信号的幅值对实验的影响。(1)我们得知当B=0.3时原系统是稳定的，这里我们逐步改变阶跃信号的幅值，使A分别取1、1.3、1.6、1.9、2.3来观察其结果，如图8所示：北京化工大学信息学院[e-mail]7图8A取不同值的过渡过程由图8可以看出，当A由小逐渐增大时，系统将由稳定转向发散，因为在设计实验或真实过程中，该扰动的幅值不可太大，否则将使得系统发散。其原因已在1中做过说明。(2)在原系统稳定的情况下，我们改变修改常数Kc的初始值，分别取Kc=-5、-0.5、-0.2、0、0.2、0.5、5来进行观测，结果如图9所示：图9取不同Kc初始值的仿真由实验我们得到Kc的最终稳定值为0.49921，有图可以看出，当Kc的初始值取得离此稳态值越远的话，过程的初始超调越大，但最终过程都能趋于稳定。因而在设计实验的时候，Kc的初始值应根据先验知识或粗略计算去一个与其稳定值较为相近的值为宜。五结论：由以上的推导及仿真结果可以看出，依据最优控制的方法设计出的MIT控制律并不能保证控制器在任何情况下都能很好的工作，换言之，对于连续系统当212/bbAKBKmp时系统会不稳定。对于离散系统随未给出准确的解析表达式，但从定性的角度来说，各参数的影响是相似的。因而在实验设计中为避免系统发散，阶跃信号的幅值不可选择太大。同时由此式可以看出，当Kp减小后，对应使系统振荡的B将增大，这就说明了第三部分中我们在没有改变B的情况下将Kp由2变为1后系统的性能为什么会得到提高。北京化工大学信息学院[e-mail]8参考文献[1]韩曾晋.自适应控制[M].北京:清华大学出版社,1995:148-151[2]厉玉鸣,等.自动控制原理[M].北京：化学工业出版社，2005:279-280附程序clcclearts=0.1;B=0.8367;ei=0;%临界值B=0.836Kp=2;Km=1;Kc(1)=0.2;num1=[Km];den1=[121];%参考模型sys1=tf(num1,den1);dsys1=c2d(sys1,ts,'z')[num11,den11]=tfdata(dsys1,'v');num2=[Kp];den2=[121];%对象模型sys2=tf(num2,den2);dsys2=c2d(sys2,ts,'z')[num22,den22]=tfdata(dsys2,'v');ym_1=0;ym_2=0;r_1=0;r_2=0;yp_1=0;yp_2=0;u_1=0;u_2=0;fori=1:1:500time(i)=i*ts;rin(i)=1;r(i)=rin(i);u(i)=Kc(i)*rin(i);ym(i)=-den11(2)*ym_1-den11(3)*ym_2+num11(2)*r_1+num11(3)*r_2;yp(i)=-den22(2)*yp_1-den22(3)*yp_2+num22(2)*u_1+num22(3)*u_2;error(i)=ym(i)-yp(i);Err(i)=error(i)^2;北京化工大学信息学院[e-mail]9gain(i)=Kc(i)*Kp;ei=ei+error(i)*ym(i)*ts;Kc(i+1)=Kc(1)+B*ei;ym_2=ym_1;ym_1=ym(i);r_2=r_1;r_1=r(i);yp_2=yp_1;yp_1=yp(i);u_2=u_1;u_1=u(i);endsubplot(1,3,1)plot(time,rin,'r',time,ym,'g',time,yp,'b')legend('R','ym','yp',4)title('inputandoutput')xlabel('time/second')subplot(1,3,2)plot(time,Km,'r',time,gain,'b')legend('Km','Kp*Kc')title('changeofKc')xlabel('time/second')subplot(1,3,3)plot(time,Err)title('changeoferror^2')xlabel('time/second')