基于MATLAB控制系统辨识系列2-最小二乘法

最小二乘算法及优化算法，纵向比较得出结论：只需要在程序中更改响应的迭代公式即可，简单明了。

目录

一.辨识算法之最小二乘法

二、程序代码

三、程序运行结果

一.辨识算法之最小二乘法

二、程序代码

CAR模型：假设控制系统为

CARMA模型：假设控制系统为：

具体代码集合如下，需要运行那一算法，则更改logical(0)为logical(1)，同时保证其余logical(0)

clear all；close all;a=[1 -1.5 0.7]';b=[1 0.5]';d=3;c=[1 -1 0.2]';%对象参数na=length(a)-1;nb=length(b)-1;nc=length(a)-1;%na,nb.nc分别为A,B,C阶次mode=1;%1:递推最小二乘 0：标准最小二乘L=1000;uk=zeros(d+nb,1);%输入初始值：uk(i)表示u(k-i)yk=zeros(na,1);%输出初始值xik=zeros(nc,1);%噪声初值xiek=zeros(nc,1);%噪声估计值lambda=0.977;%遗忘因子0.95-1之间x1=1;x2=1;x3=1;x4=0;S=1;%移位寄存器，方波初始值u=randn(L,1);%输入采用白噪声序列xi=sqrt(1)*randn(L,1); %输入采用白噪声系列，方差为1xi_LSQ=sqrt(0.1)*randn(L,1);%白噪声序列theta=[a(2:na+1);b];%对象参数真值theae_1=zeros(na+nb+1,1);%theae初始值theta_zg=[a(2:na+1);b;c(2:nc+1)];%增广最小二乘对象参数theae_zg_1=zeros(na+nb+1+nc,1);%na+nb+1+nc为辨识参数P=10^6*eye(na+nb+1);%CAR模型P_zg=10^6*eye(na+nb+1+nc);%增广模型%%%%%%%%%%%%%%%%%%%%%%%%%%最小二乘法%%%%%%%%%%%%%%%%%%%%%%%%%%%if logical(0)for k=1:Lphi(k,:)=[-yk;uk(d:d+nb)]';%此处phi(k,:)为行向量，便于组成phi矩阵y(k)=phi(k,:)*theta + xi(k);%采集输出数据M=xor(x3,x4); %产生M序列IM=xor(M,S);%产生逆M序列if IM==0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波%更新数据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' %计算参数估计值thetaend%%%%%%%%%%%%%%%%%%%%%%%%%%递推最小二乘法%%%%%%%%%%%%%%%%%%%%%%%%%%%if logical(0)for k=1:Lphi=[-yk;uk(d:d+nb)];%此处phi(k,:)为列向量，便于组成phi矩阵y(k)=phi'*theta + xi_LSQ(k);%采集输出数据%递推最小二乘公式K=P*phi/(1+phi'*P*phi);thetae(:,k)=theae_1+K*(y(k)-phi'*theae_1);P=(eye(na+nb+1)-K*phi')*P;%更新数据theae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(i);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot([1:L],thetae);xlabel('k');ylabel('参数估计a、b');legend('a_1','a_2','b_0','b_1');axis([0 L -2 1.5]);end%%%%%%%%%%%%%%%%%%%%%%%%%%遗忘因子递推最小二乘法%%%%%%%%%%%%%%%%%%%%%%%%%%%if logical(0)for k=1:Lif k==501a=[1 -1 0.4]';b=[2 0.8]';%对象参数突变endtheta(:,k)=[a(2:na+1);b];%对象参数真值phi=[-yk;uk(d:d+nb)];%此处phi为列向量 y(k)=phi'*theta(:,k) + xi_LSQ(k);%采集输出数据M=xor(x3,x4); %产生M序列IM=xor(M,S);%产生逆M序列if IM==0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波%更新数据x4=x3;x3=x2;x2=x1;x1=M; %输入数据为M序列%遗忘因子递推最小二乘法K=P*phi/(lambda+phi'*P*phi);thetae(:,k)=theae_1+K*(y(k)-phi'*theae_1);P=(eye(na+nb+1)-K*phi')*P/lambda;%更新数据theae_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,:));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,:));xlabel('k');ylabel('参数估计a');legend('b_0','b_1');axis([0 L -1 4]);end%%%%%%%%%%%%%%%%%%%%%%%%%%递推增广最小二乘法%%%%%%%%%%%%%%%%%%%%%%%%%%%if logical(1)for k=1:Lphi=[-yk;uk(d:d+nb);xik];%此处phi(k,:)为列向量，便于组成phi矩阵y(k)=phi'*theta_zg + xi_LSQ(k);%采集输出数据phie=[-yk;uk(d:d+nb);xiek];%组建\hat\phi%递推最小二乘公式K=P_zg*phie/(1+phie'*P_zg*phie);thetae(:,k)=theae_zg_1+K*(y(k)-phie'*theae_zg_1);P_zg=(eye(na+nb+1+nc)-K*phie')*P_zg;xie=y(k)-phie'*thetae(:,k);%白噪声的估计值 %更新数据 theae_zg_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);for i=nc:-1:2xik(i)=xik(i-1);xiek(i)=xiek(i-1);endxik(1)=xi(k);xiek(1)=xie;endfigure(1)plot([1:L],thetae(1:na,:));xlabel('k');ylabel('参数估计a');legend('a_1','a_2');axis([0 L -2 2]); figure(2)plot([1:L],thetae(na+1:na+nb+1,:));xlabel('k');ylabel('参数估计b');legend('b_0','b_1');axis([0 L 0 1.5]);figure(3)plot([1:L],thetae(na+nb+2:na+nb+1+nc,:));xlabel('k');ylabel('参数估计c');legend('c_1','c_2');axis([0 L -2 2]); end

三、程序运行结果