系统辨识大作业

自动化08-2班罗山08051220系统辨识作业学院：信息与控制工程学院专业班级：姓名：学号：2011年4月15日自动化08-2班罗山08051220-2-一、（1）题目要求一钢球以零初速由电视天线上自由下落，实验测得下落时间和距离的数据如下：t(s)123456l(m)8.9420.0550.6572.19129.85171.56假设时间t无量测误差，而距离l有量测误差。用最小二乘法确定重力加速度g。分析：对于初速度为零，从初始位置下落，下落距离等于重力加速度乘以下落时间t的平方的二分之一，因此，此处使用最小二乘法的基本算法，即可求得g的估计值，同时经画图可在坐标上描出对应的点，通过其曲线对结果进行相应的验证。作业程序：clearallcloseallclc%%%%第一题（1）程序%%%%t=[0123456]';L=[08.9420.0550.6572.19129.85171.56]';T=0.5*t.^2;A=polyfit(T,L,1);z=polyval(A,T);figureplot(T,z);figureplot(T,L,'b+')g=A(1)%辨识参数holdonplot(T,z,'r');xlabel('T/s^2');ylabel('L/m');title('最小二乘法求重力加速度g');holdoff运行结果：a.最小二乘法拟合曲线图b．参数辨识：g=9.5793即求得重力加速度值为9.5793。自动化08-2班罗山080512203（2）题目要求设两国军备竞赛模型为x(k)=ax(k-1)+by(k-1)+fy(k)=cx(k-1)+dy(k-1)+g式中x(k)和y(k)为两国军事费用（单位：百万美元），已知数据如下：k伊朗伊拉克北约华约1972289190921647811289319733982112321114611502019748801221021226711716919751123022472105251196121976121782204205717121461197798672303212009123561197891652179215988125498197950802675218561127185198040400225411129000198100233957131595试用最小二乘法确定模型参数a,b,c,d,f和g。分析同上题，仍然利用最小二乘法，只不过此题不是单个数，而是矩阵。作业程序伊朗和伊拉克clearallcloseallclc%%%%第一题（2）1程序%%%%y=[9091123221022472204230321792675];x=[2891398288011123012178986791655080];fork=2:8h(k,:)=[x(k-1)y(k-1)1];z(k,:)=[x(k)y(k)];endest=inv(h'*h)*h'*z;%算出a、b、f、c、d、g估计值，为三行两列的矩阵a=est(1,1),b=est(2,1),f=est(3,1),c=est(1,2),d=est(2,2),g=est(3,2)%分别取到a、b、f、c、d、gplot([1972:1979],x,'r');holdon%画出实际与各自的估计曲线hg=h*est;abf=(hg(:,[1]))';plot([1972:1979],abf,'b');holdonplot([1972:1979],y,'r');holdonhg=h*est;cdg=(hg(:,[2]))';plot([1972:1979],cdg,'g');xlabel('年代');ylabel('军事费用/百万美元');%写出坐标表头legend('伊朗实际值','伊朗估计值','伊拉克实际值','伊拉克估计值');title('伊朗与伊拉克军事费用对比');holdoff自动化08-2班罗山080512204北约和华约clearallcloseallclcx=[216478211146212267210525205717212009215988218561255411233957];y=[112893115020117169119612121461123561125498127185129000131595];fork=2:10h(k,:)=[x(k-1)y(k-1)1];z(k,:)=[x(k)y(k)];endest=inv(h'*h)*h'*z;%用最小二乘法求出估计值a=est(1,1),b=est(2,1),f=est(3,1),c=est(1,2),d=est(2,2),g=est(3,2)plot([1972:1981],x,'r');holdonhg=h*est;abf=(hg(:,[1]))';%取guji矩阵的第一列即a，b,f的值plot([1972:1981],abf,'b');holdonplot([1972:1981],y,'g');holdoncdg=(hg(:,[2]))';plot([1972:1981],cdg,'b');holdonxlabel('年代');ylabel('军事费用/百万美元');%写出坐标表头legend('北约实际值','北约估计值','华约克实际值','华约克估计值');title('北约与华约事费用对比');holdoff自动化08-2班罗山080512205最小二乘法辨识结果最小二乘法辨识结果参数AbfcdG伊朗与伊拉克0.5481-0.20344447.5-0.03060.7443987.1565北约与华约0.15401.7438-707700.01810.96752075.6结果分析由以上对比曲线可以看出，经过最小二乘法估计得到的数据与实际数据之间虽然存在区别，但是基本符合要求,故可取该组辨识参数数据。二、题目要求考虑理想数学模型为选择如下的辨识模型进行增广递推最小二乘参数辨识。()1.5(1)0.7(2)(1)0.5(2)1.2()(1)0.2(2)zkzkzkukukvkvkvk式中，)(kV是服从正态分布的白噪声)1,0(N。输入信号采用4阶M序列，其幅值为1。1212123()(1)(2)(1)(2)()(1)(2)zkazkazkbukbukcvkcvkcvk自动化08-2班罗山080512206给出各参数的辨识曲线和辨识误差曲线。分析应该使用增广递推最小二乘算法。作业程序clearallcloseallclc%%%%第二题程序%%%%N=15;%4位移位寄存器产生的M序列的周期y1=1;y2=1;y3=1;y4=0;fori=1:N;x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5,u(i)=-1;elseu(i)=1;endy1=x1;y2=x2;y3=x3;y4=x4;end%白噪声的产生A=19;x0=12;M=1024;fork=1:Nx=A*x0;x1=mod(x,M);v(k)=x1/512;x0=x1;end%辨识主程序z=zeros(7,N);zs=zeros(7,N);zm=zeros(7,N);zmd=zeros(7,N);z(1)=0;z(2)=0;zs(1)=0;zs(2)=0;zm(1)=0;zm(2)=0;zmd(1)=0;zmd(2)=0;theta=zeros(7,1);p0=10^6*eye(7,7);the=[theta,zeros(7,14)];e=zeros(7,15);fork=3:N;z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);h=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';x=inv(h'*p0*h+1);q=p0*h*x;d1=z(k)-h'*theta;theta1=theta+q*d1;zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2);zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[theta1(1);theta1(2);theta1(3);theta1(4)];zmd(k)=h'*theta1;e(:,k)=theta1-theta;theta=theta1;the(:,k)=theta1;p1=p0-q*q'*[h'*p0*h+1];p0=p1;endfigure(1);i=1:N;plot(i,the(1,:),'r',i,the(2,:),'r:',i,the(3,:),'b',i,the(4,:),'b:',i,the(5,:),'g',i,the(6,:),'g:',i,the(7,:),'g+')title('辨识曲线')figure(2);i=1:N;plot(i,e(1,:),'r',i,e(2,:),'r:',i,e(3,:),'b',i,e(4,:),'b:',i,e(5,:),'g',i,e(6,:),'g:',i,e(7,:),'r+')自动化08-2班罗山080512207title('辨识误差曲线')运行程序，得到以下曲线自动化08-2班罗山080512208得到参数辨识结果增广最小二乘参数辨识结果参数a1a2b1b2d1d2d3真值1.50.71.00.51.2-1.00.2估计值1.50.71.00.51.2-1.00.2三、题目要求完善试验2，给出由相关分析法给出的脉冲响应函数辨识曲线，在试验2的基础上，给出拟合二阶及三阶系统的最小二乘辨识参数曲线，及由脉冲响应函数反求传函()Gz。分析利用带遗忘因子的递推算法。作业程序clearallcloseallclc%%%%第三题程序%%%%n=2;Ts=2;np=15;%n阶次，Ts采样时间U=UY(61:300,1);%输入初值UY=UY(61:300,2);%输入初值Yf=zeros(120,4);f(:,1)=-1*Y(61:180);f(:,2)=-1*Y(60:179);f(:,3)=U(61:180);f(:,4)=U(60:179);yy=Y(62:181);q=inv(f'*f)*f'*yy;a=q(1:2);b=q(3:4);G=tf(b',[1a'],2)%传递函数GLS(z)yc=zeros(n+np,1);%脉冲响应uc=zeros(n+np,1);uc(n+1)=1/Ts;fori=1:np-1yc(n+i)=-1*[yc(i+n-1:-1:i)]'*a+[uc(i+n-1:-1:i)]'*b;endyimp=yc(n+1:n+np);figure(1);stem(0:Ts:(np-1)*Ts,2*yimp,'r')holdonplot(0:Ts:(np-1)*Ts,2*yimp,'r^-')impulse(G)自动化08-2班罗山080512209程序运行结果：（n=2）Transferfunction:0.3933z+0.2639-------------------------------z^2-0.9544z+0.2004Samplingtime:2得到辨识曲线如图所示：（n=3）Transferfunction:0.3954z^2+0.323z+0.05386--------------------------------------------------z^3-0.8075z^2+0.08785z+0.004068自动化08-2班罗山0805122010得到的传递函数为：Transferfunction:0.3882z+0.2645-----------------------z^2-0.9641z+0.2137三阶系统的脉冲响应函数四、题目要求查资料，比较最小二乘法及其改进