精华内容
下载资源
问答
  • 2021-04-26 15:05:46

    最早接触到最小二乘算法是在大四上学期的系统建模与辨识的课上,此课也是大四为数不多的认真做作业的课程之一。最小二乘算法原理简单,实用性好,编程很容易实现。

    最小二乘算法有递推最小二乘、辅助变量最小二乘、增广矩阵法等,多阶段最小二乘算法等,大四上学期和研一上学期的系统辨识中对递推、增广的算法曾编程过,这里我们讨论一种很简单的多项式最小二乘的实现。其实这个题目也是袁牛要做的数值分析中的一个。

    题目:给定数表:x=-0.75 -0.5 -0.25 0 0.25 0.5

    0.75,y=0.33 0.88 1.44 2.00 2.56 3.13

    3.71,试分别用一次、二次、三次多项式根据最小二乘原则拟合这些数据,并比较优劣。

    原理:已知X1

    X2...Xm和Y1

    Y2...Ym,求一个简单易算的近似函数y(x)去代替数据反映的关系f(x),对于近似函数,当然可用插值方法进行处理,但是当数据量太大的时候,差值方法有着明显的不足。(1)大量的实验数据难以保证每个数据的精确性,当其中数据存在一定误差的时候,由于插值条件的要求,其误差将完全被插值函数进一步继承。(2)即使所有数据都精确,但是多项式次数过高或者过多的分段处理,都会带来相应的问题。其实,我们没有必要要求近似函数和原数据在采样点值相同,这时需要引入拟合的思想,使得在采样点的误差组合最小即可。

    确定多项式y(x)=C0+C1*X+C2*X^2+...Cn*X^n,使得在采样点处(y(Xi)-Yi)^2的总和最小。对于具体的公式,随便找本书都有的,在这里需要说明的是,所需要的点数即等式的个数要多于需要估计的参数的个数。我一向觉得最重要的是把公式Estimate=inv(H'H)H'y中的H找出来,待估计量Estimate是一个列向量,Yi=[1

    Xi Xi^2...Xi^n]*[C0 C1 C2...Cn]',等式右半部份的前面即是H,从而我们得到如下的程序:

    clc;

    close all;

    %原始数据

    Y=[0.33 0.88 1.44 2.00 2.56 3.13 3.71];%给定Y数值

    X=[-0.75 -0.5 -0.25 0 0.25 0.5 0.75]; %给定X数值

    L=length(Y); %给定的数据个数

    for i=1:1:L

    row1(i)=X(i);

    row2(i)=X(i)^2;

    row3(i)=X(i)^3;

    end

    X1=[ones(L,1),row1']; %构造矩阵L×2

    X2=[ones(L,1),row1',row2']; %构造矩阵L×3

    X3=[ones(L,1),row1',row2',row3']; %构造矩阵L×4

    %一阶估计

    Estimate1=inv(X1'*X1)*X1'*Y'; for i=1:1:L

    EZ1(i)=Estimate1(2,1)*X(i)+Estimate1(1,1); end

    %二阶估计

    Estimate2=inv(X2'*X2)*X2'*Y';

    for i=1:1:L

    EZ2(i)=Estimate2(3,1)*X(i)^2+Estimate2(2,1)*X(i)+Estimate2(1,1);

    end

    %三阶估计

    Estimate3=inv(X3'*X3)*X3'*Y';

    for i=1:1:L

    EZ3(i)=Estimate3(4,1)*X(i)^3+Estimate3(3,1)*X(i)^2+Estimate3(2,1)*X(i)+Estimate3(1,1);

    end

    figure(1);

    plot(X,Y,'+',X,EZ1,'-',X,EZ2,'.',X,EZ3,'^');

    xlabel('时间');

    ylabel('数据曲线');

    title('最小二乘法估计曲线');

    legend('原系统','一阶估计','二阶估计','三阶估计','Location','NorthWest');

    figure(2);

    plot(X,Y,'-',X,EZ1,'^');

    xlabel('时间');

    ylabel('数据曲线');

    title('一阶最小二乘法估计曲线');

    legend('原系统','一阶估计','Location','NorthWest');

    figure(3);

    plot(X,Y,'-',X,EZ2,'+');

    xlabel('时间');

    ylabel('数据曲线');

    title('二阶最小二乘法估计曲线');

    legend('原系统','二阶估计','Location','NorthWest');

    figure(4);

    plot(X,Y,'-',X,EZ3,'+');

    xlabel('时间');

    ylabel('数据曲线');

    title('三阶最小二乘法估计曲线');

    legend('原系统','三阶估计','Location','NorthWest');

    %平方误差信息

    error1=0;

    error2=0;

    error3=0;

    for i=1:1:L

    error1=error1+(Y(i)-EZ1(i))^2;

    error2=error2+(Y(i)-EZ2(i))^2;

    error3=error3+(Y(i)-EZ3(i))^2;

    end

    e1=sqrt(error1)

    e2=sqrt(error2)

    e3=sqrt(error3)

    %输出系数信息

    %digits(n) 控制总位数 vpa(s,n)控制有效位数

    Estimate1

    Estimate2

    Estimate3

    在主窗口得到的输出信息如下:

    e1 =

    0.022677868380554

    e2 =

    0.006172133998484

    e3 =

    0.004629100498863

    Estimate1 =

    2.007142857142857

    2.251428571428571

    Estimate2 =

    1.997619047619047

    2.251428571428571

    0.038095238095238

    Estimate3 =

    1.997619047619047

    2.243650793650794

    0.038095238095238

    0.017777777777781

    分别为平方误差和、系数,可见随着阶数的增加,拟合误差减小。

    a4c26d1e5885305701be709a3d33442f.png

    图1 数据拟合曲线图

    当然了,matlab中也有自带的最小二乘拟合函数polyfit(X,Y,n),拟合所给数据,三次拟合的时候,主窗口如下输出:0.017777777777778

    0.038095238095239 2.243650793650793 1.997619047619048

    系数是从高到低输出,1.997619047619047 2.243650793650794 0.038095238095238 0.017777777777781 系数是从低到高输出,两者几乎完全相同。

    email:dingqian12345@126.com

    2010年07月10日上午于njust automation 101 房间

    参考文献

    [1] 李言俊,张科.系统辨识理论及应用[M].国防工业出版社,2006.

    [2] 李建良,蒋勇,汪光先.计算机数值方法[M].东南大学出版社,2000.

    CopyRight:版权所有若需转载或使用请联系作者

    Email:dingqian12345@126.com

    更多相关内容
  • 自动化专业综合设计报告 自动化专业综合设计报告 设计题目最小二乘递推算法辨识系统参 数 所在实验室 指导教师 学生姓名 班级 学号 自动化系统仿真实验室 计 082-2 班 撰写时间 2012-3-1 成绩评定 自动化专业综合...
  • 1.领域:matlab,最小二乘递推算法的参数估计算法 2.内容:基于最小二乘递推算法的参数估计matlab仿真+操作视频 3.用处:用于最小二乘递推算法的参数估计算法编程学习 4.指向人群:本硕博等教研学习使用 5.运行...
  • 最小二乘递推算法实例应用

    千次阅读 2021-12-10 13:04:05
    利用最小二乘递推算法进行辨识

    1.1 问题描述

    在这里插入图片描述

    1.2 方法思路

    模型类选CAR模型
    在这里插入图片描述

    递推算法

    在这里插入图片描述 1.3 实验结果

    在这里插入图片描述

    1.4 实验代码(Matlab)

    main.m

    %%%%采用递推最小二乘法对四个真值[a1 a2 b1 b2]=[-1.5 0.7 1 0.5]进行估计
    %%%第一步:选用的模型类为CAR模型,由已给的条件可知模型阶数为2
    %产生M序列做为输入,周期为1023,幅值为1
    u=F2(1023,1);%%利用作业2中代码封装成的函数产生周期为1023的M序列
    %产生均值为0,方差为1的高斯白噪声
    V=F1(0.1,0,550);%%利用作业1中代码封装成的函数产生白噪声
    %四个参数估计值
    Y=zeros(4,501);
    %任意给定Y初值
    Y(:,1:2)=[1 2 -1 3;0.8 -2 1 3]';
    %给定P(0)满足课本式5-2-19,使P(0)尽量大
    P=5000*eye(4); 
    %给定增益矩阵初值
    K=[1 1 1 1]';
    %任意给定Z初值
    Z(1:2)=[0 1];
    %根据课本式5-1-1计算z(k)
    for k=3:502
    Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+U(k-1)+0.5*U(k-2)+V(k-2);
    end
    %%%第二步:由课本式5-2-16,另加权因子为1,得到普通最小二乘递推算法RLS,其递推公式如下:
    for k=3:502   %递推500%%由课本式5-1-1可以得到h(k)的表达式,模型阶数为2,na=nb=2
    h=[-Z(k-1) -Z(k-2) U(k-1) U(k-2)]';%此时的h为第k次的值
    %%由课本式5-2-16可得RLS表达式如下:
    K=P*h*inv(h'*P*h+1);%此时的P为第k-1次的值,h为第k次的值
    Y(:,k)=Y(:,k-1)+K*(Z(k)-h'*Y(:,k-1));%此时的K为第k次的值,h为第k次的值
    P=(eye(4)-K*h')*P;%此时的P为第k次的值
    end
    %%%第三步:画图,并比较递推500次后的估计值与真值
    figure
    t=[1:1:502];
    %%画出a1的估计值
    plot(t,Y(1,:),'r')
    e1=-1.5-Y(1,502);%a1与真值的误差
    text(502,Y(1,502),'a1');
    hold on
    %%画出a2的估计值
    plot(t,Y(2,:),'g')
    e2=0.7-Y(2,502);%a2与真值的误差
    text(502,Y(2,502),'a2');
    hold on
    %%画出b1的估计值
    plot(t,Y(3,:),'b')
    e3=1-Y(3,502);%b1与真值的误差
    text(502,Y(3,502),'b1');
    hold on
    %%画出b2的估计值
    plot(t,Y(4,:),'m')
    e4=0.5-Y(4,502);%b1与真值的误差
    text(502,Y(4,502),'b2');
    title('递推500次后的参数估计值的跟踪曲线')
    legend('a1','a2','b1','b2')

    F1.m

    function Y3=F1(v,m,n)
    %利用乘同余法和变换抽样法产生均值为m,方差为v的正态分布的白噪声
    %公式:Xi=A*Xi*mod(M)
    M1=2^11;%2的方幂
    M2=2^17;%2的方幂
    A1=119;%%不能太小
    A2=279;%%不能太小
    x1=1;%伪随机序列1的初值
    x2=11;%伪随机序列2的初值
    X1=x1;
    X2=x2;
    Y1=[];%伪随机序列1
    Y2=[];%伪随机序列1
    % n=700; %伪随机序列长度
    %%%%%采用乘同余法产生伪随机序列%%%%%
    for i=1:n
        %产生伪随机序列1
        X1=mod(A1*X1,M1);
       Y1=[Y1 X1/M1];
        %产生伪随机序列2
        X2=mod(A2*X2,M2);
        Y2=[Y2 X2/M2];
    end
    %正态分布的白噪声
    Y3=m+v*sqrt(-2*log(Y1)).*cos(2*pi*Y2);%正态分布的白噪声
    figure
    plot(Y3,'b');
    title('正态分布的白噪声');
    end
    

    F2.m

    function X=F2(T,a)
    %%%采用伪随机信号发生器生成周期为1023的M序列
    y=[1 1 0 1 1 0 1 0 1 0];%%%初始化10级移位寄存器,保证其初值不全为0
    for i = 1:1:T
        Y(i)=y(10);%%%第n级的输出即是M序列的两种状态
        M(i)=mod(y(10)+y(7),2);%第n级与第N1级模二加结果,此时分别为107
        %%%%移位寄存器工作原理%%%
        y(10)=y(9);
        y(9)=y(8);
        y(8)=y(7);
        y(7)=y(6);
        y(6)=y(5);
        y(5)=y(4);
        y(4)=y(3);
        y(3)=y(2);
        y(2)=y(1);
        y(1)= M(i);
    end
    for i=1:1:T
        if Y(i)==1
            X(i)=a;
        else
            X(i)=-a;
        end
    end
    stairs(X,'b'); %画出伪随机信号
    title('周期为1023、幅值为1的伪随机信号');
    ylim([-1.2 1.2])
    end
    
    展开全文
  • 广义最小二乘递推 matlab实现A. The selected modelB. Identification algorithm usedC.Matlab codeD.The simulation imageE. Make a table to compare the estimation of parameters with the truth valueF. ...

    A. The selected model

    在这里插入图片描述
    在这里插入图片描述

    In this experiment, the specific model is

    在这里插入图片描述

    Where, the truth value of each parameter is:
    a1=-1.8;a2=0.9; b1=1.3;b2=0.7;c1=1.3;c2=0.8;d1=1.1;d2=0.3。
    In the experiment, the parameters are updated once to verify the stability and reliability of the model,after update(when time >=update_time, and we can set update_time):
    a1=-0.9;a2=0.6; b1=1.6;b2=0.3;c1=1.7;c2=1.0;d1=1.5;d2=0.5。

    B. Identification algorithm used

    Recursive generalized extended least squares algorithm,the derivation is as follows:
    在这里插入图片描述

    在这里插入图片描述
    在这里插入图片描述

    在这里插入图片描述

    在这里插入图片描述

    在这里插入图片描述

    在这里插入图片描述

    C.Matlab code

    %广义最小二乘的递推算法
    %z(k)-1.8*z(k-1)+0.9*z(k-2)=1.3*M(k-1)+0.7*M(k-2)+w(k)
    %w(k)+1.3*w(k-1)+0.8*w(k-2)=v(k)+1.1*v(k-1)+0.3*v(k-2)
    %参数更新后:
    %z(k)-0.9*z(k-1)+0.6*z(k-2)=1.6*M(k-1)+0.3*M(k-2)+w(k)
    %w(k)+1.7*w(k-1)+1.0*w(k-2)=v(k)+1.5*v(k-1)+0.5*v(k-2)
    %========================================
    clear all
    close all
    clc
    %===========系统参数真值==============
    a=[-1.8 0.9]';
    b=[1.3 0.7]';
    c=[1.3 0.8]';
    d=[1.1 0.3]';
    na=length(a);nb=length(b);nc=length(c);nd=length(d);  %%na,nb,nc,nd 为输出输入系数矩阵 A,B,C,D 的阶数
    a1=[-0.9 0.6]'; %系统参数更新后的真值
    b1=[1.6 0.3]'; 
    c1=[1.7 1.0]';
    d1=[1.5 0.5]';
    %========================================
    time=1200; %仿真总时间长度
    update_time=500;%更新参数的时间节点
    uk=zeros(nb,1);yk=zeros(na,1);  %输入输出初始化
    u=randn(time,1);  %输入信号,采用方差为1的白噪声序列
    xi=sqrt(0.1)*randn(time,1);  %干扰信号,采用方差为0.1的白噪声序列
    vk=sqrt(0.1)*randn(nd,1);wk=zeros(nc,1);  %w输入输出初始化
    vk_estimate=zeros(nd,1);wk_estimate=zeros(nc,1);  %w输入输出估计初始化
    %========================================
    thetae_1=zeros(na+nb+nc+nd,1);  %参数初值
    P=10^6*eye(na+nb+nc+nd);  %修正系数初值
    
    lambda=1;  %遗忘因子范围[0.9 1],遗忘因子=1时,即为标准递推最小二乘
    
    for k=1:time
        if k>=update_time          %更新参数
            theta_1(:,k)=[a1;b1];  %系统参数真值
            theta_2(:,k)=[c1;d1];  %系统参数真值
        else
            theta_1(:,k)=[a;b];  %系统参数真值
            theta_2(:,k)=[c;d];  %系统参数真值
        end
        phi=[-yk;uk(1:nb)];  %输出输入矩阵
        phi_w=[-wk;vk(1:nd)];  %w输出输入矩阵
        phi_w_estimate=[-wk_estimate;vk_estimate(1:nd)];
        phiall=[phi;phi_w]
        y(k)=phiall'*[theta_1(:,k);theta_2(:,k)]+xi(k);  %系统实际输出
        
        phiall_estimate=[phi;phi_w_estimate]
        %递推公式
        L=P*phiall_estimate/(lambda+phiall_estimate'*P*phiall_estimate);
        thetae(:,k)=thetae_1+L*(y(k)-phiall_estimate'*thetae_1);
        P=(eye(na+nb+nc+nd)-L*phiall_estimate')*P/lambda;
        
        %更新数据
        thetae_1=thetae(:,k);
        for i=nb:-1:2
            uk(i)=uk(i-1);
        end
        uk(1)=u(k);
        for j=nc:-1:2
            wk_estimate(j)=wk_estimate(j-1);
        end
        wk_estimate(1)=y(k)-phi'*thetae_1(1:na+nb);    
        for s=nd:-1:2
            vk_estimate(j)=vk_estimate(j-1);
        end
        vk_estimate(1)=y(k)-phiall_estimate'*thetae_1(:);
        for j=nc:-1:2
            wk(j)=wk(j-1);
        end
        wk(1)=phi_w'*theta_2(:,k)+xi(k);%w系统实际输出    
        for s=nd:-1:2
            vk(j)=vk(j-1);
        end
        vk(1)=sqrt(0.1)*randn(1);
        for i=na:-1:2
            yk(i)=yk(i-1);
        end
        yk(1)=y(k);
    end 
    %=====================输出结果及作图========================= 
    disp('参数 a1 a2 b1 b2 的估计结果(更新参数前):') 
    ab=thetae(1:na+nb,update_time-1)
    disp('噪声传递系数 c1 c2 d1 d2 的估计结果(更新参数前):') 
    cd=thetae(na+nb+1:na+nb+nc+nd,update_time-1)
    disp('参数 a1 a2 b1 b2 的估计结果(更新参数后):') 
    ab_updata=thetae(1:na+nb,time)
    disp('噪声传递系数 c1 c2 d1 d2 的估计结果(更新参数后):') 
    cd_updata=thetae(na+nb+1:na+nb+nc+nd,time)
    disp('参数估计的平方误差总和:') 
    abcd=[a;b;c;d;a1;b1;c1;d1]-[ab;cd;ab_updata;cd_updata];
    erro=abcd'*abcd
    i=1:time;
    figure('name','参数a,b')
    plot(i,thetae(1,:),i,thetae(2,:),i,thetae(3,:),i,thetae(4,:));hold on;plot(i,theta_1(:,:),'k:');
    xlabel('k');ylabel(' 参数估计a1 a2 b1 b2 ');
    legend('a1','a2','b1','b2');
    title(' 参数a1 a2 b1 b2的估计过程 ');
    figure('name',' 参数c,d ')
    plot(i,thetae(5,:),i,thetae(6,:),i,thetae(7,:),i,thetae(8,:));hold on;hold on;plot(i,theta_2(:,:),'k:');
    xlabel('k');ylabel('参数估计c1 c2 d1 d2 ');
    legend('c1','c2','d1','d2');
    title('噪声传递系数c1 c2 d1 d2的估计过程'); 
    

    D.The simulation image

    在这里插入图片描述
    请添加图片描述
    请添加图片描述
    在这里插入图片描述
    请添加图片描述
    请添加图片描述
    在这里插入图片描述

    在这里插入图片描述
    在这里插入图片描述

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    E. Make a table to compare the estimation of parameters with the truth value

    在这里插入图片描述

    F. Conclusion

    As can be seen from the comparison in the table,the experimental results are close to the real value.Compare the experimental results produced by different forgetting factor,when we take the forgetting factor=0.996,sum of squared errors is minimum,about 1.4298,that means parameter estimation works best.

    展开全文
  • 用matlab实现最小二乘递推算法辨识系统参数[参照].pdf
  • 递推最小二乘算法

    2018-05-30 15:44:54
    递推最小二乘算法递推最小二乘算法递推最小二乘算法递推最小二乘算法递推最小二乘算法递推最小二乘算法
  • 此Matlab文件为最小二乘递推算法,可根据需要,修改其中部分代码
  • 增广最小二乘递推算法对应的噪声模型为滑动平均噪声,扩充了参数向量和数据向量H(k)的维数,把噪声模型的辨识同时考虑进去。最小二乘法只能获得过程模型的参数估计,而增广最小二乘法同时又能获得噪声模型的参数...
  • 资源包括递推最小二乘算法程序原代码及word版详细实验报告,程序基于Python开发实现,代码结构清晰、完整。
  • 遗忘因子最小二乘递推算法,参考文献:过程辨识,方崇智,清华大学出版社
  • 3. 限定记忆最小二乘递推算法 9 4. 偏差补偿最小二乘法 11 5. 增广最小二乘法 13 6. 广义最小二乘法 15 7. 辅助变量法 17 8. 二步法 19 9. 多级最小二乘法 21 10. Yule-Walker辨识算法 23 Matlab程序附录 24 附录1、...
  • 递推最小二乘(RLS)算法在迭代过程的每一步都能达到最佳迭代算法。RLS算法可以使输出信号在最小二乘意义上尽可能接近期望信号,因为它可以选择自适应滤波器的权重系数。RLS算法在最小化过程中必须使用所有可用的...
  • 针对有色噪声干扰的输出误差滑动平均系统,将辅助模型与递推增广最小二乘算法相结合:用辅助模型的 输出代替辨识模型信息向量中的未知真实输出项,用估计残差代替信息向量中的不可测噪声项,从而提出了基于辅助模型的...
  • 批次过程时变动态模型的2维递推最小二乘辨识算法
  • 递推最小二乘实现参数估计,对不确定的系统有良好的参数估计效果。
  • 广义最小二乘递推算法,清华大学出版社,方崇智,过程辨识
  • 给出多重线性回归:yi=β0+β|xi|+…+βpxip+εi(i=1,2,…,n)的最小二乘估计的递推算法:β(n)=β(n-1)+Pnxn(yn-xTn β(n-1)) Pn=Pn-1-pn-1xnxTnPn-1/1+xTnPn-1xn β(0)=0,P0=αl(α>>1)
  • 递推最小二乘的matlab算法,真实有效,很好用的 打包好了的
  • 针对探测系统跟踪误差对目标跟踪的影响,提出一种实时递推最小二乘预测跟踪算法。该算法采用平方预测器的估算方式,获得目标运动轨迹的最佳逼近,通过不断更新的历史数据实时递推轨迹参数,预测下一帧目标位置。该...
  • 遗忘因子递推最小二乘参数估计,用于识别系统,MATLAB程序 遗忘因子递推最小二乘参数估计,用于识别系统,MATLAB程序
  • 在已有的偏最小二乘相关算法基础上,提出一种简化的递推最小二乘算法,即直接采用自变量主元的2个回归系数矩阵来取代残差矩阵进行递推计算,进一步简化了递推计算过程,在保证建模精度的同时,使计算速度提高了近一倍。...
  • 基于经典最小二乘法和BP神经网络算法的系统辨识,内含MATLAB算法实例
  • 本人在学习柴天佑老师的自适应控制一书时,在2.3.2参数估计小节遇到了最小二乘法的参数估计,虽然书中给出了详细的推导过程,但是随后的仿真实验中并没有给出详细的代码,亦没有给出很好的解释,只是给出了仿真结果...

    本人在学习柴天佑老师的自适应控制一书时,在2.3.2参数估计小节遇到了最小二乘法的参数估计,虽然书中给出了详细的推导过程,但是随后的仿真实验中并没有给出详细的代码,亦没有给出很好的解释,只是给出了仿真结果,让人摸不着头脑。在参考了大量资料后,我算是有一点点的理解,趁热打铁,赶忙写下来,也算是自己复习了。
    测试使用的时间序列函数为:
    y(k)+1.5y(k-1)+0.6y(k-2) = 2u(k-3)-1.4u(k-4)+V(k)
    1.构造出输入输出序列y(k) 、u(k)
    这种非病态算法 最后肯定会趋向于某个值 所以 k可以取大一些 书中取的是150 这里取400
    输出序列u(k) 可以任意给 输出序列y(k)由输入序列决定
    2.定义 P 数据向量序列seita 、未知参数序列H
    根据下图的算法公式 我们需要先初始化一些向量和矩阵
    在这里插入图片描述
    未知参数H(0)=0
    P(0)=10^6I 其中I是单位阵
    3.根据公式开始计算
    先得到 数据向量序列seita 这里包括 {-y(k-1),-y(k-2),u(k-3),u(k-4)}
    因为k-i大于等于0,所以k应从5开始取值 构造出这种结果
    过程基本就是这样 matlab代码我会附在结尾处

    下面看结果
    测试使用的时间序列函数为:
    y(k)+1.5y(k-1)+0.6y(k-2) = 2u(k-3)-1.4u(k-4)+V(k)
    输出的估计值序列变化过程如下图:

    在这里插入图片描述
    以a0为例 图中红线 基本保持在了 1.5 调出最后一组估计值

    a0 = 1.5037
    a1 = 0.6052
    b0 = 2.0709
    b1 = -1.4479
    对比实际值
    a0 = 1.5
    a1 = 0.6
    b0 = 2.0
    b1 = -1.4
    看起来还不错_
    改变a0的值为1.0 看是否可以继续估计成功
    在这里插入图片描述
    红线是a0的估计值 看起来也还行 调出最后一组数据
    a0 = 1.0286
    a1 = 0.6305
    b0 = 1.9105
    b1 = -1.3208
    对比实际值
    a0 = 1.0
    a1 = 0.6
    b0 = 2.0
    b1 = -1.4
    也算是基本保持跟随了

    好了 基本情况就是这样了 如有错误 请大家批评我!
    附代码:(我不知道如何插入matlab的代码 就用C的格式了)

    ```%递推最小二乘辨识算法仿真   
    %测试使用的时间序列函数为: y(k)+1.5y(k-1)+0.6y(k-2) = 2u(k-3)-1.4u(k-4)+V(k)
    %待估计量 1.5 0.6 2 -1.4
    %y(k) = -1.5y(k-1)-0.6y(k-2)+2u(k-3)-1.4u(k-4)+V(k)
    %阶数 Na = 2 Ma = 4
    %先产生出 y()  u()   u未知  是一个行向量
    clear
    clc
    x = [0 1 0 1 1 0 1 1 1];
    k = 400;    
    num = 4;
    u = [];                  
    for i=1:k
        temp = xor(x(4),x(9));
        u(i) = x(9);
        for j = 9:-1:2
            x(j) = x(j-1);
        end
        x(1) = temp;
    end                      %此处的目的是产生一个无规则的向量  可以任意赋值  此处参考了网上的一种做法
    v = randn(1,400);
    y = zeros(400,1);
    y(1) = -1;
    y(2) = 0;
    y(3) = -1;
    y(4) = 0;
    y(5) = 0;
    for i=5:k
    y(i) = -1.0*y(i-1)-0.6*y(i-2)+2*u(i-3)-1.4*u(i-4)+v(i);
    end
    %y(k) = -1.5y(k-1)-0.6y(k-2)+2u(k-3)-1.4u(k-4)+V(k)
    %y序列 和 u序列 完成后  开始写
    %还需要3个向量
    %P  Seita Data
    P = (10^6)*eye(num);
    P_temp = zeros(num,k);
    P_temp(:,4) = 10^6;%P_temp = { P(1) P(2) P(3).....}
    Data = zeros(num,k);%估计值序列 = {Data(1) Data(2) Data(3) ......}
    seita = zeros(num,1);
    for i = 5:k
        for j = 1:2
            seita(j) = -y(i-j);
        end
        for l = 3:4
            seita(l) = u(i-l);
        end  %h完毕
        K = (P*seita)/(seita'*P*seita+1);
        Data(:,i) = Data(:,i-1) + K*[y(i)-seita'*Data(:,i-1)];
        P = (eye(4)-K*seita')*P; %新的P
       % for kk = 1:4
      %      P_temp(kk,i-1) = P(kk,kk);
     %   end
    end
    d1 = Data(1,:);
    d2 = Data(2,:);
    d3 = Data(3,:);
    d4 = Data(4,:);
    
    i = 1:k;
    plot(i,d1,'r');
    title('辨识结果');
    hold on;
    plot(i,d2,'b');
    hold on;
    plot(i,d3,'g');
    hold on;
    plot(i,d4,'y');
    hold on;
    legend('a0','a1','b0','b1');
    a0 = Data(1,k);
    a1 = Data(2,k);
    b0 = Data(3,k);
    b1 = Data(4,k);
    a0
    a1
    b0
    b1
    
    
    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 898
精华内容 359
关键字:

最小二乘递推算法