网站域名 空间申请表免费创建个人博客网站

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

目录

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

 二、程序代码

 三、程序运行结果

---

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

方法	递推公式	性能指标
经典最小二乘法（CAR）
递推最小二乘法 （CAR）
带遗忘因子递推最小二乘法 （CAR）
递推增广最小二乘法 （CARMA）

 二、程序代码

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); %输入采用白噪声系列，方差为1
xi_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' %计算参数估计值theta
end%%%%%%%%%%%%%%%%%%%%%%%%%%递推最小二乘法%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

 三、程序运行结果

 

参考文献：系统辨识与自适应控制MATLAB仿真（第3版）庞中华、崔红编著