系统辨识.docx
《系统辨识.docx》由会员分享，可在线阅读，更多相关《系统辨识.docx（25页珍藏版）》请在冰点文库上搜索。
系统辨识
课程实践报告
——《系统辩识与自适应控制》
2.3方法一：
递推最小二乘法
%递推y（k）=（z^-1+0.5*z^-2）/（1-1.3*z^-1+0.3*z^-2）*u（k）+xi（k）
clearall;closeall;
a=[1-1.30.3]',b=[10.5]';d=1;
na=length（a）-1,nb=length（b）-1;
L=400;
uk=zeros（d+nb,1）;
yk=zeros（na,1）;
theta=[a（2:
na+1）;b];
thetagj_1=zeros（na+nb+1,1）;
P=10^6*eye（na+nb+1）;
u=randn（L,1）;%输入采用白噪声序列
xi=randn（L+2,1）;
fork=1:
L
%系统的差分方程
phi=[-yk;uk（d:
d+nb）];
y（k）=phi'*theta+xi（k+2）-1.3*xi（k+1）+0.3*xi（k）;
%递推最小二乘法
K=P*phi/（1+phi'*P*phi）;
thetagj（:
k）=thetagj_1+K*（y（k）-phi'*thetagj_1）;
P=（eye（na+nb+1）-K*phi'）*P;
%数据更新
thetagj_1=thetagj（:
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],thetagj）;
title（'2.3递推最小二乘法'）
xlabel（'time'）;ylabel（'参数估计a、b'）;
legend（'a_1','a_2','b_0','b_1'）;axis（[0L-22]）;
2.3方法二：
具有遗忘因子的递推最小二乘法
%***********遗忘因子****
%y（k）=（z^-1+0.5*z^-2）/（1-1.3*z^-1+0.3*z^-2）*u（k）+xi（k）
clear
a=[1-1.30.3]',b=[10.5]';d=1;
na=length（a）-1,nb=length（b）-1;
L=400;
namuda=0.95
uk=zeros（d+nb,1）;
yk=zeros（na,1）;
theta=[a（2:
na+1）;b]
thetagj_1=zeros（na+nb+1,1）;
P=10^6*eye（na+nb+1）;
u=randn（L,1）;%输入采用白噪声序列
xi=randn（L+2,1）;
fork=1:
L
%系统的差分方程
phi=[-yk;uk（d:
d+nb）];
y（k）=phi'*theta+xi（k+2）-1.3*xi（k+1）+0.3*xi（k）;
%具有遗忘因子的递推最小二乘法
K=P*phi/（namuda+phi'*P*phi）;
thetagj（:
k）=thetagj_1+K*（y（k）-phi'*thetagj_1）;
P=1/namuda*（eye（na+nb+1）-K*phi'）*P;
%数据更新
thetagj_1=thetagj（:
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],thetagj）;
xlabel（k）;ylabel（'参数估计a、b'）;
legend（'a_1','a_2','b_0','b_1'）;axis（[0L-22]）;
2.4增广最小二乘
%增广最小二乘（1-0.8*z^-1+0.15*z^-2）*y（k）=（z^-2+0.5*z^-3）
*u（k）+（1-0.65*z^-1+0.1*z^-2）*xi（k）
clearall;closeall;
a=[1-0.80.15]',b=[10.5]',c=[1-0.650.1]';d=2;
na=length（a）-1,nb=length（b）-1,nc=length（c）-1;
L=1000;
uk=zeros（d+nb,1）;
yk=zeros（na,1）;
xik=zeros（nc,1）;
xik_gj1=zeros（nc,1）;
theta=[a（2:
na+1）;b;c（2:
nc+1）];
theta_gj_1=zeros（na+nb+nc+1,1）;
I=eye（na+nb+nc+1）
P=10^6*I;
u=randn（L,1）;%输入采用白噪声序列
xi=randn（L+2,1）;
fork=1:
L
%系统的差分方程
phi=[-yk;uk（d:
d+nb）;xik];
y（k）=phi'*theta+xi（k）;
phi_gj=[-yk;uk（d:
d+nb）;xik_gj1];
%增广最小二乘法
K=P*phi_gj/（1+phi_gj'*P*phi_gj）;
theta_gj（:
k）=theta_gj_1+K*（y（k）-phi_gj'*theta_gj_1）;
P=（I-K*phi_gj'）*P;
xik_gj=y（k）-phi_gj'*theta_gj（:
k）;
%数据更新
theta_gj_1=theta_gj（:
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）;
fori=nc:
-1:
2
xik（i）=xik（i-1）;
xik_gj1（i）=xik_gj1（i-1）;
end
xik
（1）=xi（k）;
xik_gj1
（1）=xik_gj;
end
figure
（1）
plot（[1:
L],theta_gj（1:
na,:
）,[1:
L],theta_gj（na+1:
na+nb+1,:
）,[1:
L],theta_gj（na+nb+2:
na+nb+nc+1,:
））;
xlabel（'time'）;ylabel（'参数估计a'）;
title（'2.4增广最小二乘'）
legend（'a_1','a_2','b_0','b_1','c_1','c_2'）;axis（[0L-22]）;
3.5极点配置间接自校正控制
%极点配置间接自校正控制
clearall;closeall;
a=[1-1.70.6]';b=[11.2]';d=2;%对象参数
am=[1-0.60.08]';%期望传递函数分母多项式
na=length（a）-1;nb=length（b）-1;nam=length（am）-1;%na、nb为多项式A、B阶次
nF=nb+d-1;%nf为多项式F的阶次
L=400;%控制步数
namuda=0.95
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）;
yk=zeros（na,1）;%输出初值
yrk=zeros（nam,1）;
yr=20*[ones（L/4,1）;-ones（L/4,1）;ones（L/4,1）;-ones（L/4+d,1）];%期望输出
%RELS初值设置
theta=[a（2:
na+1）;b]
thetagj_1=0.001*ones（na+nb+1,1）;%非常小的正数（这里不能为0）
P=10^6*eye（na+nb+1）;
fork=1:
L
time（k）=k;
phi=[-yk;uk（d:
d+nb）]
y（k）=phi'*theta
%遗忘递推最小二乘法
K=P*phi/（namuda+phi'*P*phi）;
thetagj（:
k）=thetagj_1+K*（y（k）-phi'*thetagj_1）;
P=1/namuda*（eye（na+nb+1）-K*phi'）*P;
%提取辨识参数
agj=[1thetagj（1:
na,k）']
bgj=thetagj（na+1:
na+nb+1,k）'
c=0
BN=1
BY=1
rootb=roots（bgj）
fori=1:
nb
ifabs（rootb（i））>=1
BN=conv（BN,[1-rootb（i）]）
else
BY=conv（BY,[1-rootb（i）]）
end
end
%确定F1和G的阶次
nF1=length（BN）-1+d-1
nG=na-1
maxA0=na+nF1-nam
minA0=2*na-1-nam-length（BY）+1
nA0=minA0
%为了简便这里就取nA0=1设A0=1+0.1z^-1
A0=1;
fori=1:
nA0
A0=conv（A0,[1-0.9]）;
end
Am0=conv（A0,am）
[F1,G]=sindiophantine（agj,bgj,d,Am0）%求解单步Diophantine方程
F=conv（F1,BY）
%无静差取Bm’=Am/BY
sum（am）
Bm1=sum（am）/sum（BN）
R=Bm1*A0
u（k）=（-F（2:
nF+1）*uk（1:
nF）+R*yrk（1:
length（R））-G*yk（1:
na））;%求控制量
%更新数据
thetagj_1=thetagj（:
k）;
fori=d+nb:
-1:
2
uk（i）=uk（i-1）;
end
uk
（1）=u（k）;
fori=nam:
-1:
2
yrk（i）=yrk（i-1）;
end
yrk
（1）=yr（k）;
fori=na:
-1:
2
yk（i）=yk（i-1）;
end
yk
（1）=y（k）;
end
figure
（1）;
subplot（2,1,1）;
plot（time,yr（1:
L）,'r:
',time,y）;
title（'3.5极点配置间接自校正控制'）
xlabel（'time'）;ylabel（'y_r（k）、y（k）'）;
legend（'y_r（k）参考输入','y（k）期望输出'）;axis（[0L-120120]）;
subplot（2,1,2）;
plot（time,u）;
xlabel（'time'）;ylabel（'u（k）'）;axis（[0L-4040]）;
figure
（2）
subplot（211）
plot（[1:
L],thetagj（1:
na,:
））;
title（'3.5极点配置间接自校正控制'）
xlabel（'time'）;ylabel（'参数估计a'）;
legend（'a_1','a_2'）;axis（[0L-32]）;
subplot（212）
plot（[1:
L],thetagj（na+1:
na+nb+1,:
））;
xlabel（'time'）;ylabel（'参数估计b'）;
legend（'b_0','b_1'）;axis（[0L01.5]）;
function[F1,G]=sindiophantine（A,B,d,c）
dB=[zeros（1,d）B];
na=length（A）-1;nd=length（dB）-1;
nc=length（c）;
T=[c;zeros（na+nd-nc,1）];
AB=zeros（na+nd）;
fori=1:
na+1
forj=1:
nd
AB（i+j-1,j）=A（i）;
end
end
fori=1:
nd+1
forj=1:
na
AB（i+j-1,j+nd）=dB（i）;
end
end
L=（AB）\T;
F1=L（1:
nd）';
G=L（nd+1:
na+nd）';
end
4.1最小方差控制
%最小方差控制（MVC）
clearall;closeall;
a=[1-1.30.4];b=[1-0.5];c=[1-0.650.1];d=2;%对象参数
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%na、nb、nc为多项式A、B、C阶次
nf=nb+d-1;%nf为多项式F的阶次
L=400;%控制步数
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）;
yk=zeros（na,1）;%输出初值
yrk=zeros（nc,1）;%期望输出初值
xik=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）;%白噪声序列
[e,f,g]=diophantine（a,b,c,d）;%求解单步Diophantine方程
fork=1:
L
time（k）=k;
y（k）=-a（2:
na+1）*yk+b*uk（d:
d+nb）+c*[xi（k）;xik];%采集输出数据
u（k）=（-f（2:
nf+1）*uk（1:
nf）+c*[yr（k+d:
-1:
k+d-min（d,nc））;yrk（1:
nc-d）]-g*[y（k）;yk（1:
na-1）]）/f
（1）;%求控制量
%更新数据
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）;
fori=nc:
-1:
2
yrk（i）=yrk（i-1）;
xik（i）=xik（i-1）;
end
ifnc>0
yrk
（1）=yr（k）;
xik
（1）=xi（k）;
end
end
subplot（2,1,1）;
plot（time,yr（1:
L）,'r:
',time,y）;
title（'4.1最小方差控制'）
xlabel（'time'）;ylabel（'y_r（k）、y（k）'）;
legend（'y_r（k）','y（k）'）;
subplot（2,1,2）;
plot（time,u）;
xlabel（'time'）;ylabel（'u（k）'）;
function[e,f,g]=diophantine（a,b,c,d）
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%A、B、C的阶次
ne=d-1;ng=na-1;%E、G的阶次
ad=[a,zeros（1,ng+ne+1-na）];cd=[c,zeros（1,ng+d-nc）];
e
（1）=1;
fori=2:
ne+1
e（i）=0;
forj=2:
i
e（i）=e（i）+e（i+1-j）*ad（j）;
end
e（i）=cd（i）-e（i）;%计算ei
end
fori=1:
ng+1
g（i）=0;
forj=1:
ne+1
g（i）=g（i）+e（ne+2-j）*ad（i+j）;
end
g（i）=cd（i+d）-g（i）;%计算gi
end
f=conv（b,e）;%计算F
5.5状态变量可测时的模型参考自适应控制
%状态变量可测时的模型参考自适应控制
clearall;closeall;
h=0.01;L=40/h;
Am=[01;-10-5];Bm=[1;2];
Ap=[20;-6-7];Bp=[2;4];
sz=size（Bp）;n=sz
（1）;m=sz
（2）
P=[31;11];
R1=1*eye（m）;R2=1*eye（m）;
F0=zeros（m,n）;G0=eye（m,m）;
yr0=zeros（m,1）;u0=zeros（m,1）;e0=zeros（n,1）;
xp0=zeros（n,1）;xm0=zeros（n,1）;
fork=1:
L
time（k）=k*h;
yr（k）=sin（pi*time（k））;
xp（:
k）=xp0+h*（Ap*xp0+Bp*u0）;
xm（:
k）=xm0+h*（Am*xm0+Bm*yr0）;
e（:
k）=xm（:
k）-xp（:
k）;
F=F0+h*（R1*Bm'*P*e0*xp0'）;
G1=inv（G0）-h*（R2*Bm'*P*e0*（F*xp0+yr0）'*G0'）;
G=inv（G1）;
u（:
k）=G*（yr（k）+F*xp（:
k））;
parae（1,k）=norm（Am-Ap-Bp*G*F）;
parae（2,k）=norm（Bm-Bp*G）;
%数据更新
yr0=yr（:
k）;
xp0=xp（:
k）;
xm0=xm（:
k）;
e0=e（:
k）;
F0=F;
G0=G;
u0=u（:
k）;
end
figure
（1）
subplot（2,1,1）;
plot（time,xm（1,:
）,'r',time,xp（1,:
）,':
'）;
title（'5.5状态变量可测时的模型参考自适应控制'）
xlabel（'t'）;ylabel（'x_m_1（t）、x_p_1（t）'）;
legend（'x_m_1（t）','x_p_1（t）'）;
subplot（2,1,2）;
plot（time,xm（2,:
）,'r',time,xp（2,:
）,':
'）;
xlabel（'t'）;ylabel（'x_m_2（t）、x_p_2（t）'）;
legend（'x_m_2（t）','x_p_2（t）'）;
figure
（2）
plot（time,e（1,:
）,':
',time,e（2,:
）,'r'）;
title（'5.5状态变量可测时的模型参考自适应控制'）
xlabel（'t'）;ylabel（'状态输出误差'）;
legend（'状态变量x1的误差','状态变量x2的误差'）;
5.7输入输出信息的模型参考自适应控制
%n-m=2Narendra自适应控制律
clearall;closeall;
h=0.02;M=40/h;
nump=[2.5];denp=[1125];[Ap,Bp,Cp,Dp]=tf2ss（nump,denp）;n=length（denp）-1;
numm=[16];denm=[16.416];[Am,Bm,Cm,Dm]=tf2ss（numm,denm）;
L=[11];
Df=conv（L,numm）;
Af=[[zeros（n-2,1）,eye（n-2）];-Df（n:
-1:
2）];
Bf=[zeros（n-2,1）;1];
yr0=0;yp0=0;u0=0;e0=0;
v10=zeros（n-1,1）;v20=zeros（n-1,1）;
xp0=zeros（n,1）;xm0=zeros（n,1）;
theta0=zeros（2*n,1）;
zeta0=zeros（2*n,1）;
r=1;yr=r*[ones（1,M/4）-ones（1,M/4）ones（1,M/4）-ones（1,M/4）];
Gamma=550*eye（2*n）;
fork=1:
M
time（k）=k*h;
xp（:
k）=xp0+h*（Ap*xp0+Bp*u0）;
yp（k）=Cp*xp（:
k）+Dp*u0;
xm（:
k）=xm0+h*（Am*xm0+Bm*yr0）;
ym（k）=Cm*xm（:
k）+Dm*yr0;
e（k）=ym（k）-yp（k）;
v1=v10+h*（Af*v10+Bf*u0）;
v2=v20+h*（Af*v20+Bf*yp0）;
phi0=[yr0;v10;yp0;v20];
zeta=zeta0+h*（-L
（2）*zeta0+phi0）;
theta（:
k）=theta0+h*e0*Gamma*zeta0;
phi=[yr（k）;v1;yp（k）;v2];
u（k）=theta（:
k）'*phi+e（k）*zeta'*Gamma*zeta;
%更新数据
yr0=yr（k）;yp0=yp（k）;u0=u（k）;e0=e（k）;
v10=v1;v20=v2;
xp0=xp（:
k）;xm0=xm（:
k）;
phi0=phi;theta0=theta（:
k）;zeta0=zeta;
end
subplot（3,2,1）;
plot（time,ym,'r',time,yp,':
'）;
title（'方波输入时系统输出跟踪输入变化'）
xlabel（'t'）;ylabel（'y_m（t）、y_p（t）'）;
legend（'y_m（t）','y_p（t）'）;
subplot（3,2,2）;
plot（time,e）;
title（'方波输入时输出跟踪输入的误差仿真'）
xlabel（'t'）;ylabel（'e（t）'）;
%输入为阶跃信号的仿真
clearall;
h=0.02;M=40/h;
nump=[2.5];denp=[1125];[Ap,Bp,Cp,Dp]=tf2ss（nump,denp）;n=length（denp）-1;
numm=[16];denm=[16.416];[Am,Bm,Cm,Dm]=tf2ss（numm,denm）;
L=[11];
Df=conv（L,numm）;
Af=[[zeros（n-2,1）,eye（n-2）];-Df（n:
-1:
2）];
Bf=[zeros（n-2,1）;1];
%初值
yr0=0;yp0=0;u0=0;e0=0;
v10=zeros（n-1,1）;v20=zeros（n-1,1）;
xp0=zeros（n,1）;xm0=zeros（n,1）;
theta0=zeros（2*n,1）;
zeta0=zeros（2*n,1）;
r=1;yr=r*ones（M）;
Gamma=550*eye（2*n）;
fork=1:
M
time（k）=k*h;
xp（:
k）=xp0+h*（Ap*xp0+Bp*u0）;
yp（k）=Cp*xp（:
k）+Dp*u
升级会员 