系统辨识与自适应控制实验.docx
《系统辨识与自适应控制实验.docx》由会员分享,可在线阅读,更多相关《系统辨识与自适应控制实验.docx(22页珍藏版)》请在冰豆网上搜索。
系统辨识与自适应控制实验
中南大学
系统辨识及自适应控制实验
指导老师贺建军
姓名史伟东
专业班级测控1102班0909111814号
实验日期2014年11月
实验一递推二乘法参数辨识
设被辨识系统的数学模型由下式描述:
式中ξ(k)为方差为0.1的白噪声。
要求:
(1)当输入信号u(k)是方差为1的白噪声序列时,利用系统的输入输出值在线辨识上述模型的参数;
(2)当输入信号u(k)是幅值为1的逆M序列时,利用系统的输入输出值在线辨识上述模型的参数;
分析比较在不同输入信号作用下,对系统模型参数辨识精度的影响。
(1)clearall;closeall;
a=[1-1.50.70.1]';b=[121.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);
fork=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);
fori=d+nb:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=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([0L-22]);
(2)clearall;
a=[1-1.50.70.1]';b=[121.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];%对象参数真值
fork=1:
L
phi(k,:
)=[-yk;uk(d:
d+nb)]';%phi(k,:
)为行向量,便于组成phi矩阵
y(k)=phi(k,:
)*theta+xi(k);%采集输出数据
IM=xor(S,x4);
ifIM==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;
fori=nb+d:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=na:
-1:
2
yk(i)=yk(i-1);
end
yk
(1)=y(k);
End
实验二最小方差自校正控制实验
设二阶纯滞后被控对象的数学模型参数未知或慢时变,仿真实验时用下列模型:
式中ξ(k)为方差为0.1的白噪声。
要求:
(1)当设定输入yr(k)为幅值是10的阶跃信号时,设计最小方差直接自校正控制算法对上述对象进行闭环控制;
(2)1)当设定输入yr(k)为幅值是10的方波信号时,设计最小方差直接自校正控制算法对上述对象进行闭环控制;
(3)如果被控对象模型改为:
重复上述
(1)、
(2)实验,控制结果如何?
分析原因。
(1)clearall;closeall;
a=[1-1.50.7];b=[2.51.5];c=[10.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);
fork=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=[1thetae(ng+nh+3:
ng+nh+nc+2,k)'];
ifabs(ce
(2))>0.9
ce
(2)=sign(ce
(2))*0.9;
end
ifhe
(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);%求控制量
%更新数据
fori=d:
-1;2
thetaek(:
i)=thetaek(:
i-1);
end
thetaek(:
1)=thetae(:
k);
fori=d+nh:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=d+ng:
-1:
2
yk(i)=yk(i-1);
end
yk
(1)=y(k);
fori=nc:
-1:
2
yek(i)=yek(i-1);
yrk(i)=yrk(i-1);
xik(i)=xik(i-1);
end
ifnc>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([0L-2020]);
subplot(2,1,2);
plot(time,u);
xlabel('k');ylabel('u(k)');axis([0L-4040]);
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([0L-34]);
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([0L04]);
(2)clearall;closeall;
a=[1-1.50.7];b=[2.51.5];c=[10.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);
fork=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=[1thetae(ng+nh+3:
ng+nh+nc+2,k)'];
ifabs(ce
(2))>0.9
ce
(2)=sign(ce
(2))*0.9;
end
ifhe
(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);%求控制量
%更新数据
fori=d:
-1;2
thetaek(:
i)=thetaek(:
i-1);
end
thetaek(:
1)=thetae(:
k);
fori=d+nh:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=d+ng:
-1:
2
yk(i)=yk(i-1);
end
yk
(1)=y(k);
fori=nc:
-1:
2
yek(i)=yek(i-1);
yrk(i)=yrk(i-1);
xik(i)=xik(i-1);
end
ifnc>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([0L-2020]);
subplot(2,1,2);
plot(time,u);
xlabel('k');ylabel('u(k)');axis([0L-4040]);
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([0L-34]);
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([0L04]);
(3-1)clearall;closeall;
a=[1-1.50.7];b=[51.5];c=[10.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);
fork=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=[1thetae(ng+nh+3:
ng+nh+nc+2,k)'];
ifabs(ce
(2))>0.9
ce
(2)=sign(ce
(2))*0.9;
end
ifhe
(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);%求控制量
%更新数据
fori=d:
-1;2
thetaek(:
i)=thetaek(:
i-1);
end
thetaek(:
1)=thetae(:
k);
fori=d+nh:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=d+ng:
-1:
2
yk(i)=yk(i-1);
end
yk
(1)=y(k);
fori=nc:
-1:
2
yek(i)=yek(i-1);
yrk(i)=yrk(i-1);
xik(i)=xik(i-1);
end
ifnc>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([0L-2020]);
subplot(2,1,2);
plot(time,u);
xlabel('k');ylabel('u(k)');axis([0L-4040]);
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([0L-34]);
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([0L04]);
(3-2)clearall;closeall;
a=[1-1.50.7];b=[51.5];c=[10.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);
fork=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=[1thetae(ng+nh+3:
ng+nh+nc+2,k)'];
ifabs(ce
(2))>0.9
ce
(2)=sign(ce
(2))*0.9;
end
ifhe
(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);%求控制量
%更新数据
fori=d:
-1;2
thetaek(:
i)=thetaek(:
i-1);
end
thetaek(:
1)=thetae(:
k);
fori=d+nh:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=d+ng:
-1:
2
yk(i)=yk(i-1);
end
yk
(1)=y(k);
fori=nc:
-1:
2
yek(i)=yek(i-1);
yrk(i)=yrk(i-1);
xik(i)=xik(i-1);
end
ifnc>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([0L-2020]);
subplot(2,1,2);
plot(time,u);
xlabel('k');ylabel('u(k)');axis([0L-4040]);
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([0L-34]);
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([0L04]);
实验三模型参考自适应控制实验
设被控对象模型参数未知或慢时变,但其状态变量完全可观测,仿真时取状态方程为:
选择参考模型:
状态完全可观测的模型参考自适应控制系统如下图所示:
.
.
⨯
⨯
(t)
e(t)
(t)
控制器自适应规律为:
,
式中:
,为m⨯m矩阵(m为输入个数)。
当参考输入为
时,要求选择三组合适的P、R1和R2,实现对被控对象的控制,使被控对象的2个状态变量分别跟踪参考模型的2个状态变量,并分析P、R1和R2对控制系统性能的影响。
clearall;closeall;
h=0.01;L=100/h;
As=[01;-5-3];Bs=[0;6];%对象参数
Am=[01;-10-5];Bm=[0;2