《系统辨识》作业 1题目最小二乘参数估计仿真学院信息科学与工程学院专业控制科学与工程姓名王瑞荣学号 030120694教师顾幸生2012年12月26日一：一般最小二乘法 考虑仿真系统：式中，)(k v 为方差为1的白噪声。

选用幅值为1的逆M 序列作为输入信号)(k u ，利用LS 算法进行参数估计，仿真结果如下表所示。

表1.1 批处理最小二乘法的参数估计结果仿真程序：%批处理最小二乘参数估计（LS ） clear all;a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; %na 、nb 为A 、B 阶次L=500; %数据长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值x1=1; x2=1; x3=1; x4=0; S=1; %移位寄存器初值、方波初值 xi=randn(L,1); %白噪声序列theta=[a(2:na+1);b]; %对象参数真值 for k=1:Lphi(k,:)=[-yk;uk(d:d+nb)]'; %此处phi(k,:)为行向量，便于组成phi 矩阵 y(k)=phi(k,:)*theta+xi(k); %采集输出数据IM=xor(S,x4); %产生逆M 序列 if IM==0 u(k)=-1; elseu(k)=1; endS=not(S); M=xor(x3,x4); %产生M 序列%更新数据)()4(5.0)3()2(7.0)1(5.1)(k v k u k u k y k y k y +-+-=-+--x4=x3; x3=x2; x2=x1; x1=M;for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetae=inv(phi'*phi)*phi'*y' %计算参数估计值thetae二：递推最小二乘算法 算法：已知：n a 、n b 和d 。

Step1 设置初值)0(P 和)0(∧θ，输入初始数据； Step2 采样当前输出)(k y 和输入)(k u ；Step3 利用递推公式，计算K(k)、)(k ∧θ和)(k P ； Step4 1+→k k ,返回Step2,继续循环。

考虑仿真系统：式中，)(k v 为方差为1的白噪声。

取初值I P 610)0(=、0)0(=∧θ。

选择方差为1的白噪声作为输入信号)(k u ，采用RLS 算法进行参数估计，仿真结果如图1.2所示。

5101520253035-2.5-2-1.5-1-0.500.511.52k参数估计a 、b图1.2 递推最小二乘法的参数估计结果当k=35时，参数估计值为0.48780.9276,0.6389,-1.4399,1021====∧∧∧∧b b a a 。

)()4(5.0)3()2(7.0)1(5.1)(k v k u k u k y k y k y +-+-=-+--仿真程序：%递推最小二乘参数估计（RLS）clear all; close all;a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=35; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列theta=[a(2:na+1);b]; %对象参数真值thetae_1=zeros(na+nb+1,1); %thetae初值P=10^6*eye(na+nb+1);for k=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'*thetae_1);P=(eye(na+nb+1)-K*phi')*P;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=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','a_2','b_0','b_1');%axis([0 L -2 2]);三：遗忘因子递推最小二乘算法 算法：已知：n a 、n b 和d 。

Step1 设置初值)0(P 和)0(∧θ及遗忘因子λ，输入初始数据； Step2 采样当前输出)(k y 和输入)(k u ；Step3 利用遗忘因子递推公式，计算K(k)、)(k ∧θ和)(k P ； Step4 1+→k k ,返回Step2,继续循环。

考虑仿真系统：式中，)(k v 为方差为1的白噪声。

对象时变参数[]Tb b a a k 1021)(=θ为[]⎩⎨⎧-≤-=500,]2.05.14.01[500,5.017.05.1)( k k k TT θ 取初值I P 610)0(=、0)0(=∧θ。

选择方差为1的白噪声作为输入信号)(k u 。

遗忘因子最小二乘法、取遗忘因子λ=0.98，采用FFRLS 算法进行参数估计，仿真结果如图1.3(1)当[]；）（时，参数估计值为T k 4835.00198.17103.05156.1500500-==∧θ 当[]；）（时，参数估计值为T k 1026.05147.14609.00669.110001000-==∧θ可以看出，几十对于参数突变的系统，采用FFRLS 算法也能够有效地进行参数估计。

)()4()3()2()1()(1021k v k u b k b k y a k y a k y +-+-=-+-+5001000-2-1.5-1-0.500.511.52k参数估计a500100000.511.522.533.54k参数估计b图1.3(1) 遗忘因子递推最小二乘法的参数估计结果（λ=0.98）当取遗忘因子λ=1时，FFRLS 将退化为普通的RLS 算法，仿真结果如图1.3(2)所示。

5001000-2-1.5-1-0.500.511.52k参数估计a500100000.20.40.60.811.21.41.61.8k参数估计b图1.3(2) 递推最小二乘法的参数估计结果（λ=1）可以看出，RLS对于参数时变系统，即使增加数据长度，也不能有效的跟踪参数的变化（除非增加代数到30000代左右时才可较好跟踪（本例））。

仿真程序：%遗忘因子递推最小二乘参数估计（FFRLS）clear all; close all;a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=1000; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列thetae_1=zeros(na+nb+1,1); %thetae初值P=10^6*eye(na+nb+1);lambda=0.98; %遗忘因子范围[0.9 1]for k=1:Lif k==501a=[1 -1 0.4]';b=[1.5 0.2]'; %对象参数突变endtheta(:,k)=[a(2:na+1);b]; %对象参数真值phi=[-yk;uk(d:d+nb)];y(k)=phi'*theta(:,k)+xi(k); %采集输出数据%遗忘因子递推最小二乘法K=P*phi/(lambda+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P/lambda;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endsubplot(1,2,1)plot([1:L],thetae(1:na,:)); hold on; plot([1:L],theta(1:na,:),'k:');xlabel('k'); ylabel('参数估计a');legend('a_1','a_2'); axis([0 L -2 2]);subplot(1,2,2)plot([1:L],thetae(na+1:na+nb+1,:)); hold on; plot([1:L],theta(na+1:na+nb+1,:),'k:'); xlabel('k'); ylabel('参数估计b');legend('b_0','b_1'); %axis([0 L -0.5 2]);。