电加热炉的系统辨识与自适应控制

1、电加热炉的系统辨识与自适应限制一、电加热炉的先验知识1.1.1 电加热炉的工作原理1.1.2 电加热炉温度限制系统的硬件构成2二、电加热炉系统辨识3.2.1 电加热炉温度系统模型32.2 最小二乘估计的递推算法.4.2.3 最小二乘估计的递推算法辨识及仿真5三、电加热炉系统的自适应限制算法及仿真83.1 电加热炉系统限制问题的提出83.2 广义最小方差间接自校正限制算法83.3 广义最小方差间接自校正限制仿真9参考资料1.5电加热炉的系统辨识与自适应限制一、电加热炉的先验知识1.1 电加热炉的工作原理我选择电加热炉作为辨识和自适应限制设计与仿真实验的对象.电加热炉的工作原理为：布置在炉内的加热

2、元件将电能转化为热能,通过辐射或对流的方式将热能传递给加热对象,从而改变对象的温度.通常的工业过程都对炉温的限制提出了一定的要求,这就需要对电加热炉的进行限制,调节它的通电时间或通电强度来改变它输出的热能.传统的限制方法有两种：第一种就是手动调压法,即是依靠人的经验直接改变电加热炉的输入电压,其控温效果依赖于人为的调节,限制精度不高,且浪费人力资源.第二种限制方法在主回路中采取可控硅装置,并结合一些简单的仪表,保温阶段自动调节,升温过程仍依赖于试验者的调节,它属于半自动限制.随着微型计算机、可编程逻辑限制器的出现和迅速更新换代,智能温度限制仪表、工业限制计算机在电加热炉温度限制领域日益得到广泛

3、地应用.借助计算机强大的数据处理和运算水平,引入反应的思想,运用现代限制理论,实现对炉温的全自动化限制1.以常用的恒温箱式电加加热炉为例,采用反应限制.该限制系统的目的是要实现炉内的温度与给定温度值一致,即保持温度恒定,是一个典型的自动限制系统.当然,系统给定的不是具体的期望温度值,而是通过给定电位器给定一个电压UsT.电加热炉内的实际温度由热电偶转换为对应的电压UfT.给定电压信号Ust与实际温度所对应的电压UfT比拟得温度偏差信号4U经放大器放大后,用以驱动执行电动机,并通过传动机构拖动调压器动触头.当温度偏高时,动触头向减小电压的方向运动,反之加大电压,直到温度到达给定值为止,此时,偏差

4、AU=0,电机停止转动.上面只是一个比拟简单的闭环温度限制系统.1.2 电加热炉温度限制系统的硬件构成电加热炉温度限制系统框图如图1.1所示+图1.1电加热炉温度限制系统根据信号的流动,其工作原理大致是：首先将热电偶传来的带有温度信号的毫伏级电压滤波、放大,送至A/D转换器,这样通过采样和A/D转换,就将所检测的炉温对应的电压信号转换成数字量送入了限制装置如微机、智能仪表的处理器等；在限制装置内计算出该电压信号对应的温度值,然后将它与给定的温度值进行比拟,并按一定的限制算法进行运算；运算结果通过限制品闸管在限制周期内的触发角,也就是限制电加热炉的平均功率的大小来到达温度限制的目的.设计温度限制

5、算法时还需要将上述的原理图简化成模型如图1.2所示,以便于系统进行分析.图1.2系统简化模型图模型中的限制器就是广义的加载到计算机或微处理器上的限制算法,品闸管模块、电加热炉、加热对象一起归为限制对象,而A/D转换器、热电偶那么构成反应回路.限制器给定的温度作为系统输入信号r,传感器检测到的温度作为输出信号y,误差e、限制信号u均在限制装置里通过计算得到,最后经过一系列转换实现对热工对象温度的限制,这就形成了一个典型的反应限制系统2o二、电加热炉系统辨识2.1 电加热炉温度系统模型电加热炉的温度限制是典型的过程限制.由于传热问题的复杂性,电加热炉系统具有非线性、时变性、大滞后、不对称等特点.它

6、的滞后主要是容积滞后,炉体的结构、容量、测温元件及其安装的位置都影响着滞后的大小；而在使用过程中,随着温度的升降,加热元件的特性发生变化,保温绝热材料会逐渐老化,环境也在不断变化,因而炉温特性是时变的；又由于绝大多数电加热炉都是在温度上升时强迫加热,而温度下降时那么自然冷却,所以其温度特性是不对称的；另外由于炉温取决于加热元件的发热量、散热量和负荷的情况,发热时间总比传热时间短得多,所以炉温动态特性主要由传热过程决定,传导、对流、辐射三种的传热方式都在起作用,只是在不同温区所占比例不同,三者中只有传导是线性的,辐射是绝对温度的四次方,对流那么更加复杂,故电加热炉是一个本质非线性的系统.由于电加

7、热炉可认为是一个大容积滞后加纯滞后的对象(容积滞后比纯滞后大得多),故在其整个温度工作区域,对象动态参数是随炉温变化的,而每个炉子都有一个设定的工作温区,在工作点附近的小范围内,炉子的动态特性可看成近似线性.在过程限制中,为了方便,通常把电加热炉温控系统看成是一个线性系统,其模型可以定性描述为：dY2T+Y=K0U2(t-t)(2-1)dt式中,Y为加热对象温度,t为加热时间,T为系统时间常数,Ko为放大倍数,U为限制电压,丁为纯滞后时间.如果设定限制器输出为u,而u正比于U2,即KoU2=Ku,对式(2-1)作拉氏变换,可得:TsY(s)Y(s)=KU(s)e-s(2-2)所以Y(s)Ke-

8、sU(s)一Ts1故系统的传递函数为一阶惯性加纯滞后环节,其中G(s)=Ke-sTs1(2-3)K为静态增益3.(2-4)2.2 最小二乘估计的递推算法最小二乘法由于原理简明、收敛较快、易于理解、易于编程实现等特点、在系统参数估计中应用相当广泛.而最小二乘法中的递推算法由于能对对象参数在线实时估计,从而改善了估计精度,在系统辨识中倍受青睐.最小二乘估计递推算法的根本思想可以为：新的估计值0(k)=旧的估计值(k-1)+修正项.最小二乘估计递推算法的公式为4-I-,T.(k)-?(k-1)K(k)y(k)-(k)(k-1)K(k).P(K-1)%)(2-5)1T(k)P(k-1)(k)P(k)=

9、I-K(k)T(k)P(k-1)最小二乘估计递推算法如下5：na、叫和d.1设置初值其0)和P(0),输入初始数据；2采样当前输出y(k)和输入u(k)；3利用式(2-5),计算K(k)、的k)和P(k)；4kTk+1,返回2,继续循环.2.3最小二乘估计的递推算法辨识及仿真电加热炉的参数模型,其传递函数表示为60.44e20sG(s)=520s1(2-6)可以考虑利用最小二乘估计递推算法对上述系统进行在线辨识.为了得到系统的离散模型,可利用MATLAB对系统的传递函数(2-6)进行离散化.采样时间取T=20so在MATLAB中输入命令G=tf(0.44,5201,'inputdela

10、y',20),Gz=c2d(G20,'z').得到系统的离散系统的传递函数为G(z)=z-0.9623z将电加热炉转换为差分方程,有y(k)-0.9623y(k-1)-0.0166u(k-2)(2-7)(2-8)在实际的电加热炉模型中,一定存在噪声干扰.往差分方程里添加噪声干扰项后,得到电加热炉的动态方程为y(k)-0.9623y(k-1)=0.0166u(k-2)(k)(2-9)式中£(k)为白噪声.取初值P(0)=107I,$(0)=0,选才¥M序列为输入信号u(k严.取1500个采样点.采用最小二乘估计递推算法进行参数估计,具体程序如下(程序参

11、考?系统辨识与自适应限制MATLAB仿真?,并进行了修改)：%递推最小二乘参数估计(RLS)clearall;closeall;a=1-0.9623;b=0.0166;d=2;%对象参数na=length(a)-1;nb=length(b)-1;%nanb为A、B阶次L=1500;%仿真长度x1=1;x2=1;x3=1;x4=0;%产生M序列的移位存放器初值s=1;%方波初值uk=zeros(d+nb,1);%俞入初值：uk(i)表示u(k-i)yk=zeros(na,1);%输出初值fork1=1:Ls=not(s);%产生方波m(k1)=xor(x3,x4);%进行异或运算,产生M序列x4

12、=x3;x3=x2;x2=x1;x1=m(k1);%存放器移位endxi=sqrt(0)*randn(L,1);%白噪声序列,修改括号中数值可产生不同方差的白噪声theta=a(2:na+1);b;%对象参数真值thetae_1=zeros(na+nb+1,1);%theta初值P=10A7*eye(na+nb+1);fork=1:Lphi=-yk;uk(d:d+nb);%止匕处phi为列向量y(k)=phi'*theta+xi(k);%采集输出数据%递推最小二乘法K=P*phi/(1+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi&#

13、39;*thetae_1);P=(eye(na+nb+1)-K*phi')*P;%更新数据thetae_1=thetae(:,k);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=m(k);%M序列!！fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot(1:L,thetae);%line(1,L,theta,theta);xlabel('k');ylabel('参数估计a、b');legend('a_1','b_0');axis(0L-1.51.5);分

14、别设置不同方差的白噪声,运行程序,具体辨识结果如图2.1-2.4所示.1.51.5abQ0.5ba计估数-0.5-1-1.5050010001500图2.1白色噪声方差为0时参数辨识结果b-0.5-1.510001500计估数图2.3白色噪声方差为0.1时参数辨识结果0.5ba计估数-0.5-1.5050010001500图2.2白色噪声方差为0.01时参数辨识结果1.510.50-0.5-1-1.5kb、a计估数图2.4白色噪声方差为1时参数辨识结果不同白噪声下的具体数据如表2.1所示:表2.1参数b0真值-0.96230.0166白噪声方差为0.01-0.96280.0174白噪声力差为0

15、.1-0.95750.0199白噪声力差为1-0.96470.0238以上参数数据和曲线说明随着扰动强度的增强,辨识的效果会有所下降.当扰动十分剧烈的时候,系统的辨识效果会十分差.三、电加热炉系统的自适应限制算法及仿真3.1 电加热炉系统限制问题的提出第一章中已经介绍了电加热炉的工作原理,电加热炉通常是通过限制系统中传动机构拖动调压器动触头来改变输入电压,保持炉温的稳定平衡,从而保证电加热炉系统稳定正常的工作.在实际工业生产中,一般要求电加热炉炉温的超调量尽可能小.因此,对于电加热炉的限制问题是：给定电加热炉的输入电压大小,通过某种限制方式来调节电加热炉温度,使电加热炉温度的超调量尽可能小并到

16、达目标温度数值.3.2 广义最小方差间接自校正限制算法自适应限制算法包括有自校正调节器、自校正限制器、自适应极点配置PID限制器、自校正PID限制等.我选用广义最小方差自校正限制器对电加热炉系统的温度进行自适应限制.最小方差限制具有算法简单、易于理解、易于实现等优点,是其他自校正限制算法的根底.其根本思想是：由于一般工业对象存在纯延时d,当前的限制作用要滞后d个采样周期才能影响输出.因此,要使输出方差最小,就必须提前d步对输出量作出预测,然后根据所得的预测值来设计所需的限制律.这样,通过连续不断的预测和限制,就能保证稳态输出方差最小.由此可见,实现最小方差限制的关键在于输出预测.由于最小方差限

17、制不适用非最小相位系统,且输入限制量未受到约束,因此,在此根底上,出现了广义最小方差限制算法.其根本思想为：在求解限制律的性能指标中引入对限制量的加权项,从而限制限制作用过于剧烈变化；另外,只要适中选择性能指标中的各加权多项式,广义最小方差限制可以可以适用于非最小相位系统.广义最小方差间接自校正限制算法如下5:：na、nb、nc及纯延时d.1设置初值政0和P0,输入初始数据,并设置加权多项式Pz"、Rz"和Q(z)；2采样当前实际输出y(k)和期望输出y(k+d);3利用递推增广最小二乘法在线实时估计被控对象参数即A、2和6;-1'A/-lx.-/-1.-d-14求

18、解Diophantine方程Z1Z1Z1ZZ,得到多项式E、F、F(z-)=B(z-)E(z-)和G的系数；5利用式u(k)=C(z-1)R(z-1)y-G(z)P(z-1Mk)计算并实施u(k)；q0C(z-1)Q(z-1)F(z-1)P(z-1)b06返回2(kTk+1),继续循环.3.3 广义最小方差间接自校正限制仿真在上面系统辨识中已经知道电加热炉转的差分方程,故自校正限制中,被控对象电加热炉的动态方程同为y(k)-0.9623y(k-1)=0.0166u(k-2)(k)(2-9)式中,“k)为白噪声.取初值P(0)=1071,&0)=0.001;设置加权多项式P(z,)=1、

19、R(z,)=1和Q(z,)=0.01.期望输出y.(k)为幅值从0.5-10依次递增变化的方波.具体程序如下(程序参考?系统辨识与自适应限制MATLAB仿真?,并进行了修改)：%广义最小方差自校正限制(问接算法)clearall;closeall;clc;a=1-0.9623;b=0.0166;c=1;d=2;%对象参数na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%nanb、nc为多项式A、B、C阶次nf=nb+d-1;ng=na-1;%nf、ng为多项式F、G的阶次Pw=1;R=1;Q=0.01;%加权多项式P、R、Qnp=length(Pw)-

20、1;nr=length(R)-1;nq=length(Q)-1;L=1200;%限制步数uk=zeros(d+nb,1);%俞入初值：uk(i)表示u(k-i);yk=zeros(na,1);%输出初值yrk=zeros(nc,1);%期望输出初值yr=10*0.15*ones(L/8,1);0.05*ones(L/8,1);0.4*ones(L/8,1);0.05*ones(L/8,1);0.7*ones(L/8,1);0.05*ones(L/8,1);1*ones(L/8,1);0.05*ones(L/8+d,1);%期望输出vki=zeros(nc,1);%白噪声初值vkie=zeros

21、(nc,1);%白噪声估计初值vk=sqrt(0.01)*randn(L,1);%白噪声序列,修改括号中数值可产生不同方差的白噪声%RLS初值thetae_1=0.001*ones(na+nb+1+nc,1);%£常小的正数,此处不能为0P=10A7*eye(na+nb+1+nc);fork=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*vk(k);%采集输出数据%c=1,故用递推最小二乘法phie=-yk(1:na);uk(d:d+nb);K=P*phie/(1+phie'*P*phie);thetae(:,k)=thetae

22、_1+K*(y(k)-phie'*thetae_1);P=(eye(na+nb+1+nc)-K*phie')*P;%提取辨识参数ae=1thetae(1:na,k)'be=thetae(na+1:na+nb+1,k)'ce=1;e,f,g=sindiophantine(ae,be,ce,d);%Jt解单步Diophantine方程CQ=conv(ce,Q);FP=conv(f,Pw);CR=conv(ce,R);GP=conv(g,Pw);%CQ=Ce*Qu(k)=(-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/be(1)-FP(2:np+nf

23、+1)*uk(1:np+nf).+CR*yr(k+d:-1:k+d-min(d,nr+nc);yrk(1:nr+nc-d).-GP*y(k);yk(1:np+ng)/(Q(1)*CQ(1)/be(1)+FP(1);%求限制量%更新数据thetae_1=thetae(:,k);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endfigure(1);subplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k

24、');ylabel('y_r(k)、y(k)');legend('y_r(k)','y(k)');axis(0L-212);subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)');axis(0L-2030);figure(2)plot(1:L,thetae);%line(1,L,theta,theta);xlabel('k');ylabel('参数估计a1、b0');legend('a1','b0&

25、#39;);axis(0L-1.50.2);其中sindiophantine函数程序为：functione,f,g=sindiophantine(a,b,c,d)%*%功能：单步Diophanine方程的求解%调用格式：e,f,g=sindiophantine(a,b,c,d)%输入参数：多项式A、B、C系数(行向量)及纯滞后(共4个)%输出参数：Diophanine方程的解e,f,g(共3个)%*na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%A、B、C的阶次ne=d-1;ng=na-1;%E、G的阶次ad=a,zeros(1,ng+ne+1-na);cd=c,zeros(1,ng+d-nc);%令a(na+2)=a(na