系统辨识与建模实验报告.docx

系统辨识与建模实验 电加热炉动态特性辨识实验系统辨识与建模实验报告电加热炉动态特性辨识实验姓名学号： 张春燕 312102332同组同学：沈 剑 312102331序号： 81指导老师：郭毓实验时间：2013年5月目 录一 实验目的32.2 数据获取32.3 离线辨识33.1 数据预处理33.2 结构辨识43.3.1 RLS辨识参数63.3.2 RELS辨识参数73.3.3 RIV辨识参数83.3.4 RML辨识参数93.4 模型验证103.4.1 输入阶跃响应比较103.4.2 比较残差11四实验结果分析12五实验心得12附录1：13附录2：13附录3：14附录4：16171 实验目的通过实验了解辨识方法在工程应用中的一些实际问题；了解数据获取和数据处理的各种方法和手段，掌握各种辨识方法的应用特点。2 实验内容及其步骤2.1 编写M序列的产生程序在实验参数设定时选择加热电压60V，采样周期为3S，所以加入的M序列电压最好为加热电压的10-20%，M序列的采样周期为数据采样周期的整数倍，因为实验时间有限，选择了2组数据，即M序列信号为6V、10V，采样周期60s.2.2 数据获取高温老化试验温箱，以控制电压作为炉温控制系统的输入控制变量，即，设备的输入量是燃料供给量或电压、电流，而输出量是炉膛内腔的温度。在热稳定工况的基础上，在电压稳定值上再附加一个辨识信号，即M序列电压信号。加热炉热惯性大，升温过程较长，所以采样周期较长，M序列的周期也较长。这样施加M序列周期信号之后，记录几个周期的炉温试验数据。2.3 离线辨识利用处理过的数据，选择某种辨识方法；如RLS、RELS、RIV或RML等参数估计方法计算，以及F检验方法或AIC定阶法。离线估计出参数模型参数，并计算相应的模型静态增益，同时比较利用不同方法所得到的辨识结果。最后，模型验证。三数据处理数据处理主要包括输入输出数据，模型结构确定，然后辨识参数。3.1 数据预处理在实验中采集了四组数据，仿真时选择M序列为10的数据组。其中采样，为隔20点采样。在离线辨识时发现，如果数据全部运用则采样时间长度过短，则计算量大并容易产生病态方程，所以数据采用隔点取数据的方法。程序见附录1。图3.1 M序列及其温度变化3.2 结构辨识残差校验估计方差最小定阶。利用最小二乘法求出值估值，计算残差。不考虑相关误差可以假设系统的传递函数如下：即需要估计的参数：已知数据构成的向量：根据最小二乘法进行系统定阶时的误差平方和：在F检验法中代入数据由matlab计算得：表1 阶数与误差对应关系12345J86.28670.23760.29758.58355.195T6.39794.45090.760931.5342AIC27.79919.45214.29616.56616.992图3.2 残差校验定阶有数据及图可以看出，当n=3时，AIC最小，而且此时t3，与工程实践结合，辨识时暂定系统为三阶系统。另外通过辨识，因为系统采样时为隔点采样，所以辨识时滞为0-1，在下面的参数估计中时滞d=1。程序见附录2。3.3 利用参数估计方法计算参数模型参数首先，根据最小二乘法辨识出参数theta_ls = -0.4890 -0.1501 -0.2295 0.0141 0.0543 0.0050 0.0154在参数辨识时，由于采用隔点取样数据量较少，可能产生误差较大。3.3.1 RLS辨识参数递推最小二乘法不考虑相关噪声，辨识模型与最小二乘法相似。根据递推最小二乘法公式：RLS步骤：1） 输入系统阶次n，遗忘因子和数据；2） 置初值和，输入初始数据；3） 采样当前输入和输出并构成；4） 按递推公式计算，和；5） 返回2）直到收敛或满足要求为止。由matlab仿真得到结果：theta_rls = -0.4865 -0.1510 -0.2305 0.0141 0.0538 0.0058 0.0157图3.3 递推最小二乘法估计参数3.3.2 RELS辨识参数增广矩阵法是一种用于实时过程控制中系统参数估计的较好方法，可同时获得系统参数和噪声模型参数。改写LS模型为需要估计的参数：已知数据构成的向量：其中，为噪声项。则有：由matlab仿真得到结果：theta_rels = -0.5341 -0.2569 -0.2151 0.0239 0.0606 0.0036 0.0113 -0.0922 -0.1563 -0.0757图3.4 增广矩阵法估计参数3.3.3 RIV辨识参数辅助变量法是在算法中引入一个辅助变量。与RLS相比，其中由matlab仿真得到结果：theta_riv = -0.3645 -0.3645 -0.3645 0.0588 0.0196 0.0304 0.0183图3.5 递推辅助变量法估计参数3.3.4 RML辨识参数极大似然估计的基本思想是构造一个与观测数据和未知参数有关的似然函数，并极大化这个函数获得的模型的参数估计。其公式：其中：由matlab仿真得到结果：theta_rml = -0.3223 -0.4026 -0.1489 0.0022 0.0010 0.0196 -0.0077 -0.0870 0.1454 -0.0334图3.6 递推极大似然估计参数3.4 模型验证常用的模型辨识方法：比较实际系统的输出与各个模型的输出之间的残差，选出残差较小的模型。比较实际系统于各个模型的阶跃响应，选出误差较小的模型。本节以RLS的模型验证为例。程序见附录43.4.1 输入阶跃响应比较在实验时，温箱初始温度30度左右，系统加热电压为60V，一段时间后温度大约稳定在147。给模拟的系统一个60V的阶跃响应，系统最终温度142.55，误差在可接受范围内，基本符合要求。 图3.7 阶跃响应3.4.2 比较残差给系统一个按实验时的M序列，如图为模拟的系统与真实系统的输出反应。图3.8 真实值与理论值比较四实验结果分析首先在选取数据采样频率时，放弃大部分的采样数据以免选择全部数据造成病态方程，但是也样做的同时也可能导致Ts过大，才而降低系统动态特性的精确度。其次在定阶时，从分析数据可以看出模拟系统为3阶比较合适，其实，系统的残差在7,8阶的时候更小，但是理论结合实际，选择3阶，由于选择隔点取样，采样周期较大，时滞相对来说变小。分析温箱系统是一个慢时变系统，所以在参数估计中加入了遗忘因子，从第三章三种递推算法中可以看出随着迭代次数增加，最终收敛到一定范围内。比较参数估计，由于最小二乘法模型与递推最小二乘法相似，可以看出旗鼓基础的参数也是相似的，后来模型验证得出，RLS误差总体小于LS。RELS算法能辨识出系统参数以及噪声模型。RIV引进一个新的相关噪声，理论上误差应该较小，但是实际模型验证时误差较大，可能是数据较少导致的。对于RML，由于未知噪声的概率密度按函数，辨识出误差较大。五实验心得在离线辨识时发现，如果数据全部运用则采样时间长度过短，则计算量大并容易产生病态方程，所以数据采用隔点取数据的方法。因采样时，M序列设置的为15，采样两个周期即为60组，数据偏少，导致辨识时结果误差较大。造成这些问题的主要原因是，在实验前对实验的了解不充分，验前知识准备不足。从而导致实验时虽然记录的实验数据不足，无法做当场辨识出实验模型从而无法运用实验装置及其PID参数辨识模型。如果验前知识准备充足，采样时采集数据足够可以更精确的辨识出实验模型。这次实验充分暴露出我在实验上的缺点，在真正动手实验时浮躁的问题，努力在毕业前改正。最后，感谢郭老师的认真教导。附录1： clear all;%-导入数据M10=csvread(D:wenxiangM55.csv);%温度154左右%-去直流分量L=60;z=0;u=0;for i=1:L; %求直流分量 z=M10(30*(i-1)+19,3)+z; u=M10(30*(i-1)+2,1)+u;endu0=u/L;z0=z/L;z=zeros(L,1);u=zeros(L,1);for i=1:L; % 隔点取值 z(i)=M10(30*(i-1)+19,3)-z0; u(i)=M10(30*(i-1)+2,1)-60;endKs=z0/u0;u=u;z=z;figure(1);%画第一个图形k=1:1:L;subplot(2,1,1); %画第一个图形的第一个子图stem(u),grid on%画出M序列输入信号title(移位寄存器产生的M序列);ylabel(电压);subplot(2,1,2); %画第一个图形的第一个子图plot(k,z,k,z,rx);ylabel(温度);title(温箱温度变化);附录2：%-定阶n=6;J=zeros(n,1);AIC=zeros(n,1);e=zeros(L,1);n=1; %假设为1阶for i=2:L H(i,:)=-z(i-1) u(i) u(i-1);endtheta_ls=inv(H*H)*H*z;for i=2:L e(i)=z(i)-H(i,:)*theta_ls; J(1)=J(1)+e(i)2;endAIC(n)=L*log(e*e/L)+4*n+2;n=2; %假设为2阶H=zeros(L,5);e=zeros(L,1);for i=3:L H(i,:)=-z(i-1) -z(i-2) u(i) u(i-1) u(i-2);endtheta_ls=inv(H*H)*H*z;for i=3:L e(i)=z(i)-H(i,:)*theta_ls; J(2)=J(2)+e(i)2;endAIC(n)=L*log(e*e/L)+4*n+2;n=3; %假设为3阶H=zeros(L,7);e=zeros(L,1);for i=4:L H(i,:)=-z(i-1) -z(i-2) -z(i-3) u(i) u(i-1) u(i-2) u(i-3);endtheta_ls=inv(H*H)*H*z;for i=4:L e(i)=z(i)-H(i,:)*theta_ls; J(3)=J(3)+e(i)2;endAIC(n)=L*log(e*e/L)+4*n+2;n=4; %假设为4阶H=zeros(L,9);e=zeros(L,1);for i=5:L H(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) u(i) u(i-1) u(i-2) u(i-3) u(i-4);endtheta_ls=inv(H*H)*H*z;for i=5:L e(i)=z(i)-H(i,:)*theta_ls; J(4)=J(4)+e(i)2;endAIC(n)=L*log(e*e/L)+4*n+2;n=5; %假设为5阶H=zeros(L,11);e=zeros(L,1);for i=6:L H(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) -z(i-5) u(i) u(i-1) u(i-2) u(i-3) u(i-4) u(i-5);endtheta_ls=inv(H*H)*H*z;for i=6:L e(i)=z(i)-H(i,:)*theta_ls; J(5)=J(5)+e(i)2;endAIC(n)=L*log(e*e/L)+4*n+2;n=6; %假设为6阶H=zeros(L,13);e=zeros(L,1);for i=7:L H(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) -z(i-5) -z(i-6) u(i) u(i-1) u(i-2) u(i-3) u(i-4) u(i-5) u(i-6);endtheta_ls=inv(H*H)*H*z;for i=6:L e(i)=z(i)-H(i,:)*theta_ls; J(6)=J(6)+e(i)2;endAIC(n)=L*log(e*e/L)+4*n+2;t=zeros(n,1);for i=1:n-1 t(i)=J(i)-J(i+1)/J(i+1)*(L-2*(i+1)/2;endn=6;figure(2);subplot(3,1,1); %画第二个图形的第一个子图plot(1:n),J);xlabel(n);ylabel(J);title(估计误差方差最小定阶);subplot(3,1,2); %画第二个图形的第二个子图plot(1:n),t);%画出taxis(1 5 0 3);xlabel(n);ylabel(t);title(F检验定阶);subplot(3,1,3); %画第二个图形的第三个子图plot(1:n),AIC);%画出txlabel(n);ylabel(AIC);title(F检验定阶);附录3：%LSd=1; %时滞H=zeros(L,7);e_ls=zeros(L,1);for i=(4+d):L H(i,:)=-z(i-1) -z(i-2) -z(i-3) u(i-d) u(i-1-d) u(i-2-d) u(i-3-d);endtheta_ls=inv(H*H)*H*z;for i=(4+d):L e_ls(i)=z(i)-H(i,:)*theta_ls;end%RLSe_rls=zeros(L,1);theta_rls=0.001*ones(7,L);%a1;a2;a3;b0;b1;b2;b3P=106*eye(7,7);delta=0.998;for i=(4+d):L H=-z(i-1) -z(i-2) -z(i-3) u(i-d) u(i-1-d) u(i-2-d) u(i-3-d); K=P*H/(delta+H*P*H); e_rls(i)=z(i)-H*theta_rls(:,i-1); theta_rls(:,i)=theta_rls(:,i-1)+K*e_rls(i); e_rls(i)=z(i)-H*theta_rls(:,i); P=(P-K*H*P)/delta;endfigure(3);i=1:L;plot(i,theta_rls(:,i);xlabel(t)ylabel(辨识参数ai,bi)legend(a1,a2,a3,b0,b1,b2,b3)title(RLS);%RELSP=106*eye(10);theta_rels=0.001*ones(10,L);%a1;a2;a3;b0;b3;b1;b2;c1;c2;c3e_rels=zeros(L,1);for i=(4+d):L; H=-z(i-1) -z(i-2) -z(i-3) u(i-d) u(i-1-d) u(i-2-d) u(i-3-d) e_rels(i-1) e_rels(i-2) e_rels(i-3); K=P*H/(delta+H*P*H); e_rels(i)=z(i)-H*theta_rels(:,i-1); theta_rels(:,i)=theta_rels(:,i-1)+K*e_rels(i); P=(P-K*H*P)/delta;endfigure(4);i=1:L;plot(i,theta_rels(:,i);xlabel(t)ylabel(辨识参数ai,bi,ci)title(RELS);legend(a1,a2,a3,b0,b1,b2,b3,c1,c2,c3)%RIVe_riv=zeros(L,1);P=106*eye(7);theta_riv=0.001*ones(7,L);y=zeros(L,1);for i=(4+d):L; H=-z(i-1) -z(i-2) -z(i-3) u(i-d) u(i-1-d) u(i-2-d) u(i-3-d); w=-y(i-1) -y(i-2) -y(i-3) u(i-d) u(i-1-d) u(i-2-d) u(i-3-d); y(k)=w*theta_riv(:,i-1); K=P*w*inv(H*P*w+1); theta_riv(:,i)=theta_riv(:,i-1)+K*(z(i)-H*theta_riv(:,i-1); P=(eye(7)-K*H)*P; e_riv(i)=z(i)-H*theta_riv(:,i);endfigure(5);i=(4+d):Lplot(i,theta_riv(:,i);title(RIV);xlabel(t)ylabel(辨识参数ai,bi)legend(a1,a2,a3,b0,b1,b2,b3)%RMLe_rml=zeros(L,1);I=eye(10,10);theta_rml=0.001*ones(10,L);P=I;zf(1)=0.1;zf(2)=0.1;zf(3)=0.1;zf(4)=0.1;vf(1)=0.1;vf(2)=0.1;vf(3)=0.1;vf(4)=0.1;uf(1)=0.1;uf(2)=0.1;uf(3)=0.1;uf(4)=0.1;v=randn(L,1);for k=(4+d):L; H_f=-zf(k-1);-zf(k-2);-zf(k-3);uf(k-d);uf(k-1-d);uf(k-2-d);uf(k-3-d);vf(k-1);vf(k-2);vf(k-3); K=P*H_f/(H_f*P*H_f+1); P=I-K*H_f*P; v(k)=z(k)-H_f*theta_rml(:,k-1); theta_rml(:,k)=theta_rml(:,k-1)+K*v(k); zf(k)=z(k)-theta_rml(8,k-1)*zf(k-1)-theta_rml(9,k-1)*zf(k-2)-theta_rml(10,k-1)*zf(k-3); uf(k)=u(k)-theta_rml(8,k-1)*uf(k-1-d)-theta_rml(9,k-1)*uf(k-2-d