最小二乘参数辨识方法及应用程序清单.doc

第3章 最小二乘参数辨识方法及应用程序清单一、3.2.2节利用最小二乘法求取模型参数的程序程序1： 热敏电阻例子clear allclose allclcT=20.5 26 32.7 40 51 61 73 80 88 95.7;R=765 790 826 850 873 910 942 980 1010 1032;m,n=size(T);figure plot(T,R,b+)t=0;z=0;tz=0;tt=0;for i=1:n t=t+T(i); tt=tt+T(i)*T(i); z=z+R(i); tz=tz+T(i)*R(i);enda=(tt*z-t*tz)/(n*tt-t*t)b=(n*tz-t*z)/(n*tt-t*t)%A=polyfit(T,R,1);z=polyval(A,T);figureplot(T,z);figureplot(T,R,b+)hold onplot(T,z,r);hold off程序2t=0.1:0.1:1;N=10;y=a+b*t;v=rands(1,10);z=y+v;sum_t1=0;sum_t2=0;sum_tz=0;sum_N=0;sum_z=0;for i=1:1:N sum_t1=sum_t1+t(i); sum_t2=sum_t2+t(i)2; sum_tz=sum_tz+t(i)*z(i); sum_N=sum_N+z(i); sum_z=sum_z+z(i);end ap_num=sum_t2*sum_N-sum_t1*sum_tz; ap_den=N*sum_t2-sum_t12; ap=ap_num/ap_den bp_num=N*sum_tz-sum_t1*sum_z; bp_den=N*sum_t2-sum_t12; bp=bp_num/bp_den二、3.3 加权最小二乘算法程序3：最小二乘参数辨识程序clear all%清理工作间变量close all%关闭所有图形clc%清屏z(1)=440,z(2)=430,z(3)=420,z(4)=380,z(5)=370,z(6)=360,z(7)=320,z(8)=310,z(9)=300,z(10)=260,z(11)=250,z(12)=240,z(13)=220,z(14)=210,z(15)=170,z(16)=160;u(1)=3,u(2)=2.7,u(3)=2.4,u(4)=2.1,u(5)=2.0,u(6)=1.9,u(7)=1.6,u(8)=1.54,u(9)=1.48,u(10)=1.2,u(11)=1.14,u(12)=1.08,u(13)=0.95,u(14)=0.9,u(15)=0.7,u(16)=0.6;HL=-z(1) u(1);-z(2) u(2);-z(3) u(3);-z(4) u(4);-z(5) u(5);-z(6) u(6); -z(7) u(7); -z(8) u(8);-z(9) u(9); -z(10) u(10);-z(11) u(11);-z(12) u(12);-z(13) u(13);-z(14) u(14) %给样本矩阵HL赋值ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)% 给样本矩阵zL赋值%calculating parameters%计算参数c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3 %计算并显示 %DISPLAY PARAMETERSa2=c(1), b2=c(2)三、3.4 递推最小二乘算法程序4：递推最小二乘参数辨识程序clear all%清理工作间变量close all%关闭所有图形clc%清屏% 产生白噪声的程序%A=6;x0=1;M=255;for k=1:10000 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*2; x0=x1; v0=v1;endnum=v;k1=k;% num=zeros(1,1000);% randn(seed,100);% num=randn(1,1000);% 白噪声程序结束 % M序列产生程序 %L=15;% M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值for i=1:L;%开始循环，长度为L x1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或” x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出 x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出 x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出 y(i)=y4;%取出第四个移位积存器幅值为0和1的输出信号， if y(i)0.5,u(i)=-1;%如果M序列的值为1时,辨识的输入信号取“-0.03” else u(i)=1;%当M序列的值为0时,辨识的输入信号取“0.03” end%小循环结束 y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号u figure(1);%第1个图形stem(u),grid on%以径的形式显示出输入信号并给图形加上网格% M序列产生程序程序结束 % 递推最小二乘辨识程序%lamt=1;z(2)=0;z(1)=0;%取z的前两个初始值为零for k=3:15;%循环变量从3到15 z(k)=-1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+0*num(1,k);%给出理想的辨识输出采样信号 end%RLS递推最小二乘辨识c0=0.001 0.001 0.001 0.001;%直接给出被辨识参数的初始值,即一个充分小的实向量p0=103*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=c0,zeros(4,14);%被辨识参数矩阵的初始值及大小e=zeros(4,15);%相对误差的初始值及大小for k=3:15; %开始求K h1=-z(k-1),-z(k-2),u(k-1),u(k-2); x=h1*p0*h1+1*lamt; x1=inv(x); %开始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1*c0; c1=c0+k1*d1;%求被辨识参数c % p1=p0-k1*k1*(h1*p0*h1+1);%求出 p(k)的值 p1=1/lamt*(eye(4)-k1*h1)*p0; e1=c1-c0;%求参数当前值与上一次的值的差值 e2=e1./c0;%求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 c0=c1;%新获得的参数作为下一次递推的旧参数 p0=p1;%给下次用 if e20.5,u(i)=-1;%M序列的值为1时,辨识的输入信号取“-1” else u(i)=1;%M序列的值为0时,辨识的输入信号取“1” end y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号作准备endfigure(1);%画第一个图形subplot(2,1,1); %画第一个图形的第一个子图stem(u),grid on%画出M序列输入信号v=randn(1,60); %产生一组60个正态分布的随机噪声subplot(2,1,2); %画第一个图形的第二个子图plot(v),grid on;%画出随机噪声信号R=corrcoef(u,v);%计算输入信号与随机噪声信号的相关系数r=R(1,2);%取出互相关系数u%显示输入型号v%显示噪声型号% M序列、噪声信号产生及其显示程序结束 % 增广递推最小二乘辨识程序 %z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=zeros(7,60);%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出矩阵的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%给输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出赋初值%增广递推最小二乘辨识c0=0.001 0.001 0.001 0.001 0.001 0.001 0.001;%直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(7,7);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=5.0e-15;%取相对误差Ec=c0,zeros(7,59);%被辨识参数矩阵的初始值及大小e=zeros(7,60);%相对误差的初始值及大小for k=3:60; %开始求K z(k)=-1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+1.2*v(k)-v(k-1)+0.2*v(k-2);%系统在M序列输入下的输出采样信号 h1=-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2);%为求K(k)作准备 x=h1*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K d1=z(k)-h1*c0; c1=c0+k1*d1;%辨识参数c zs(k)=-1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%系统在M序列的输入下不考虑扰动时的输出响应 zm(k)=-z(k-1),-z(k-2),u(k-1),u(k-2)*c1(1);c1(2);c1(3);c1(4);%模型在M序列的输入下不考虑扰动时的的输出响应 zmd(k)=h1*c1;%模型在M序列的输入下的的输出响应 e1=c1-c0; e2=e1./c0; %求参数的相对变化 e(:,k)=e2; c0=c1;%给下一次用 c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵 p1=p0-k1*k1*h1*p0*h1+1;%find p(k) p0=p1;%给下次用 if e2=E break;%若收敛情况满足要求，终止计算 end%判断结束end%循环结束c, e, %显示被辨识参数及参数收敛情况z,zmd %显示输出采样值、系统实际输出值、模型输出值%分离变量 a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);%分离出a1、 a2、 b1、 b2 d1=c(5,:); d2=c(6,:); d3=c(7,:); %分离出d1、 d2、 d3 ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); %分离出a1、 a2、 b1、 b2的收敛情况 ed1=e(5,:); ed2=e(6,:); ed3=e(7,:)