系统辨识与自适应控制实验

中南大学

系统辨识及自适应控

制实验

指导老师贺建军

姓名史伟东

专业班级测控1102班0909111814号

实验日期2014年11月

实验一 递推二乘法参数辨识

设被辨识系统的数学模型由下式描述：

234123123

2.0 1.51

()()()1 1.50.70.11 1.50.70.1z z z y k u k k z z z z z z ξ---------++=+-++-++ 式中ξ(k )为方差为0.1的白噪声。要求：

（1） 当输入信号u (k )是方差为1的白噪声序列时，利用系统的输入

输出值在线辨识上述模型的参数；

（2）

当输入信号u (k )是幅值为1的逆M 序列时，利用系统的输入输出值在线辨识上述模型的参数；

分析比较在不同输入信号作用下，对系统模型参数辨识精度的影响。 （1）clear all; close all;

a=[1 -1.5 0.7 0.1]';b=[1 2 1.5]';d=3; %对象参数 na=length(a)-1;nb=length(b)-1; %计算阶次 L=500; %数据长度

uk=zeros(d+nb,1);yk=zeros(na,1); %输入输出初值 u=randn(L,1); %输入采用方差为1的白噪声序列 xi=sqrt(0.1)*randn(L,1); % 方差为0.1的白噪声干扰序列 theta=[a(2:na+1);b]; %对象参数真值 thetae_1=zeros(na+nb+1,1); %参数初值 P=10^6*eye(na+nb+1); for k=1:L

phi=[-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:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=na:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

end

plot([1:L],thetae); %line([1:L],[theta,theta]); xlabel('k');ylabel('参数估计a,b');

legend('a_1','a_2','a_3','b_0','b_1','b_2');

axis([0 L -2 2]);

（2）clear all;

a=[1 -1.5 0.7 0.1]';b=[1 2 1.5]';d=2; %对象参数

na=length(a)-1;nb=length(b)-1; %计算阶次

L=20; %数据长度

uk=zeros(d+nb,1);yk=zeros(na,1); %输入初值

x1=1;x2=1;x3=1;x4=0;S=1;%移位寄存器初值，方波初值

xi=rand(L,1);%白噪声序列

theta=[a(2:na+1);b]; %对象参数真值

for k=1:L

phi(k,:)=[-yk;uk(d:d+nb)]'; % phi(k,:)为行向量，便于组成phi矩阵y(k)=phi(k,:)*theta+xi(k); %采集输出数据

IM=xor(S,x4);

if IM==0

u(k)=-1;

else

u(k)=1;

end

S=not(S);M=xor(x3,x4); %产生M序列

%更新数据

x4=x3;x3=x2;x2=x1;x1=M;

for i=nb+d:-1:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=na:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

End

实验二 最小方差自校正控制实验

设二阶纯滞后被控对象的数学模型参数未知或慢时变，仿真实验时用下列模型：

3411212

2.5 1.510.5()()()1 1.50.71 1.50.7z z z y k u k k z z z z ξ-------++=+-+-+ 式中ξ(k )为方差为0.1的白噪声。要求：

（1） 当设定输入y r (k )为幅值是10的阶跃信号时，设计最小方差直接

自校正控制算法对上述对象进行闭环控制；

（2） 1）当设定输入y r (k )为幅值是10的方波信号时，设计最小方差

直接自校正控制算法对上述对象进行闭环控制； （3） 如果被控对象模型改为：

341

1212

0.5 1.510.5()()()1 1.50.71 1.50.7z z z y k u k k z z z z

ξ-------++=+-+-+

重复上述（1）、（2）实验，控制结果如何？分析原因。

（1）clear all;close all;

a=[1 -1.5 0.7]; b=[2.5 1.5]; c=[1 0.5]; d=4; %对象参数

na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次

nh=nb+d-1; ng=na-1; %nh 为多项式H 的阶次,ng 为多项式G 的阶次L=400;

uk=zeros(d+nh,1);

yk=zeros(d+ng,1);

yek=zeros(nc,1); %最优输出预测估计初值

yrk=zeros(nc,1);

xik=zeros(nc,1);

%xiek=zeros(nc,1); %白噪声估计值

yr=10*[ones(L/4,1);ones(L/4,1);ones(L/4,1);ones(L/4+d,1)];% 期望输出xi=sqrt(0.1)*randn(L,1);%方差为0.1 的白噪声序列

thetaek=ones(na+nb+d+nc,d);

P=10^6*eye(na+nb+d+nc);

for k=1:L

time(k)=k;

y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*[xi(k);xik];%采集输出数据

phie=[yk(d:d+ng);uk(d:d+nh);-yek(1:nc)];

K=P*phie/(1+phie'*P*phie);

thetae(:,k)=thetaek(:,1)+K*(y(k)-phie'*thetaek(:,1));

P=(eye(na+nb+d+nc)-K*phie')*P;

ye=phie'*thetaek(:,d);%预测输出估计值

%提取辨识参数

ge=thetae(1:ng+1,k)';

he=thetae(ng+2:ng+nh+2,k)';

ce=[1 thetae(ng+nh+3:ng+nh+nc+2,k)'];

if abs(ce(2))>0.9

ce(2)=sign(ce(2))*0.9;

end

if he(1)<0.1

he(1)=0.1;%设h0 的下界为0.1

end

u(k)=(-he(2:nh+1)*uk(1:nh)+ce*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-ge*[y(k);yk(1:na-1)])/he(1);%求控制量

%更新数据

for i=d:-1;2

thetaek(:,i)=thetaek(:,i-1);

end

thetaek(:,1)=thetae(:,k);

for i=d+nh:-1:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=d+ng:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

for i=nc:-1:2

yek(i)=yek(i-1);

yrk(i)=yrk(i-1);

xik(i)=xik(i-1);

end

if nc>0

yek(1)=ye;

yrk(1)=yr(k);

xik(1)=xi(k);

end

end

figure(1);

subplot(2,1,1);

plot(time,yr(1:L),'r:',time,y); xlabel('k');ylabel('y_r(k)、y(k)');

legend('y_r(k)','y(k)');axis([0 L -20 20]);

subplot(2,1,2);

plot(time,u);

xlabel('k');ylabel('u(k)');axis([0 L -40 40]);

figure(2);

subplot(2,1,1);

plot([1:L],thetae(1:ng+1,:),[1:L],thetae(ng+nh+3:ng+2+nh+nc,:)); xlabel('k');ylabel('参数估计g,c');

legend('g_0','g_1','c_1');axis([0 L -3 4]);

subplot(2,1,2);

plot([1:L],thetae(ng+2:ng+2+nh,:));

xlabel('k');ylabel('参数估计h');

legend('h_0','h_1','h_2','h_3','h_4');axis([0 L 0 4]);

（2）clear all;close all;

a=[1 -1.5 0.7]; b=[2.5 1.5]; c=[1 0.5]; d=4; %对象参数

na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次

nh=nb+d-1; ng=na-1; %nh 为多项式H 的阶次,ng 为多项式G 的阶次L=400;

uk=zeros(d+nh,1);

yk=zeros(d+ng,1);

yek=zeros(nc,1); %最优输出预测估计初值

yrk=zeros(nc,1);

xik=zeros(nc,1);

%xiek=zeros(nc,1); %白噪声估计值

yr=10*[ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1)];% 期望输出

xi=sqrt(0.1)*randn(L,1);%方差为0.1 的白噪声序列

thetaek=zeros(na+nb+d+nc,d);

P=10^6*eye(na+nb+d+nc);

for k=1:L

time(k)=k;

y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*[xi(k);xik];%采集输出数据

phie=[yk(d:d+ng);uk(d:d+nh);-yek(1:nc)];

K=P*phie/(1+phie'*P*phie);

thetae(:,k)=thetaek(:,1)+K*(y(k)-phie'*thetaek(:,1));

P=(eye(na+nb+d+nc)-K*phie')*P;

ye=phie'*thetaek(:,d);%预测输出估计值

%提取辨识参数

ge=thetae(1:ng+1,k)';

he=thetae(ng+2:ng+nh+2,k)';

ce=[1 thetae(ng+nh+3:ng+nh+nc+2,k)'];

if abs(ce(2))>0.9

ce(2)=sign(ce(2))*0.9;

end

if he(1)<0.1

he(1)=0.1;%设h0 的下界为0.1

end

u(k)=(-he(2:nh+1)*uk(1:nh)+ce*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-ge*[y(k);yk(1:na-1)])/he(1);%求控制量

%更新数据

for i=d:-1;2

thetaek(:,i)=thetaek(:,i-1);

end

thetaek(:,1)=thetae(:,k);

for i=d+nh:-1:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=d+ng:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

for i=nc:-1:2

yek(i)=yek(i-1);

yrk(i)=yrk(i-1);

xik(i)=xik(i-1);

end

if nc>0

yek(1)=ye;

yrk(1)=yr(k);

xik(1)=xi(k);

end

end

figure(1);

subplot(2,1,1);

plot(time,yr(1:L),'r:',time,y);

xlabel('k');ylabel('y_r(k)、y(k)'); legend('y_r(k)','y(k)');axis([0 L -20 20]); subplot(2,1,2);

plot(time,u);

xlabel('k');ylabel('u(k)');axis([0 L -40 40]);

figure(2);

subplot(2,1,1);

plot([1:L],thetae(1:ng+1,:),[1:L],thetae(ng+nh+3:ng+2+nh+nc,:)); xlabel('k');ylabel('参数估计g,c');

legend('g_0','g_1','c_1');axis([0 L -3 4]);

subplot(2,1,2);

plot([1:L],thetae(ng+2:ng+2+nh,:));

xlabel('k');ylabel('参数估计h');

legend('h_0','h_1','h_2','h_3','h_4');axis([0 L 0 4]);

（3-1）clear all;close all;

a=[1 -1.5 0.7]; b=[5 1.5]; c=[1 0.5]; d=4; %对象参数

na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次

nh=nb+d-1; ng=na-1; %nh 为多项式H 的阶次,ng 为多项式G 的阶次L=400;

uk=ones(d+nh,1);

yk=ones(d+ng,1);

yek=ones(nc,1); %最优输出预测估计初值

yrk=ones(nc,1);

xik=ones(nc,1);

%xiek=zeros(nc,1); %白噪声估计值

yr=10*[ones(L/4,1);ones(L/4,1);ones(L/4,1);ones(L/4+d,1)];% 期望输出xi=sqrt(0.1)*randn(L,1);%方差为0.1 的白噪声序列

thetaek=ones(na+nb+d+nc,d);

P=10^6*eye(na+nb+d+nc);

for k=1:L

time(k)=k;

y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*[xi(k);xik];%采集输出数据phie=[yk(d:d+ng);uk(d:d+nh);-yek(1:nc)];

K=P*phie/(1+phie'*P*phie);

thetae(:,k)=thetaek(:,1)+K*(y(k)-phie'*thetaek(:,1));

P=(eye(na+nb+d+nc)-K*phie')*P;

ye=phie'*thetaek(:,d);%预测输出估计值

%提取辨识参数

ge=thetae(1:ng+1,k)';

he=thetae(ng+2:ng+nh+2,k)';

ce=[1 thetae(ng+nh+3:ng+nh+nc+2,k)'];

if abs(ce(2))>0.9

ce(2)=sign(ce(2))*0.9;

end

if he(1)<0.1

he(1)=0.1;%设h0 的下界为0.1

end

u(k)=(-he(2:nh+1)*uk(1:nh)+ce*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-ge*[y(k);yk(1:na-1)])/he(1);%求控制量

%更新数据

for i=d:-1;2

thetaek(:,i)=thetaek(:,i-1); end

thetaek(:,1)=thetae(:,k); for i=d+nh:-1:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=d+ng:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

for i=nc:-1:2

yek(i)=yek(i-1);

yrk(i)=yrk(i-1);

xik(i)=xik(i-1);

end

if nc>0

yek(1)=ye;

yrk(1)=yr(k);

xik(1)=xi(k);

end

end

figure(1);

subplot(2,1,1);

plot(time,yr(1:L),'r:',time,y);

xlabel('k');ylabel('y_r(k)、y(k)');

legend('y_r(k)','y(k)');axis([0 L -20 20]);

subplot(2,1,2);

plot(time,u);

xlabel('k');ylabel('u(k)');axis([0 L -40 40]);

figure(2);

subplot(2,1,1);

plot([1:L],thetae(1:ng+1,:),[1:L],thetae(ng+nh+3:ng+2+nh+nc,:)); xlabel('k');ylabel('参数估计g,c');

legend('g_0','g_1','c_1');axis([0 L -3 4]);

subplot(2,1,2);

plot([1:L],thetae(ng+2:ng+2+nh,:));

xlabel('k');ylabel('参数估计h');

legend('h_0','h_1','h_2','h_3','h_4');axis([0 L 0 4]);

（3-2）clear all;close all;

a=[1 -1.5 0.7]; b=[5 1.5]; c=[1 0.5]; d=4; %对象参数

na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次