系统辨识及自适应控制实验报告.docx
《系统辨识及自适应控制实验报告.docx》由会员分享，可在线阅读，更多相关《系统辨识及自适应控制实验报告.docx（25页珍藏版）》请在冰点文库上搜索。
系统辨识及自适应控制实验报告
中南大学
系统辨识及自适应控制实验
指导老师贺建军
姓名
专业班级测控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]）;
实验三模型参考自适应控制实验
设被控对象模型参数未知或慢时变，但其状态变量完全可观测，仿真时取状态方程为：
选择参考模型：
状态完全可观测的模型参考自适应控制系统如下图所示：
控制器自适应规律为：
，
式中：
，为mm矩阵（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];%参考模型参数
Sz=size（Bs）;n=Sz
（1）;m=Sz
（2）;
升级会员 