系统辨识练习题.docx
《系统辨识练习题.docx》由会员分享，可在线阅读，更多相关《系统辨识练习题.docx（13页珍藏版）》请在冰豆网上搜索。
系统辨识练习题
系统辨识练习题
方法一：
%递推最小二乘参数估计（RLS）
clearall;closeall;
a=[1-1.50.7]';b=[10.5]';d=3;%对象参数
na=length（a）-1;nb=length（b）-1;%na、nb为A、B阶次
L=480;%仿真长度
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）yk=zeros（na,1）;%输出初值
u=randn（L,1）;%输入采用白噪声序列xi=sqrt（O.1）*randn（L,1）;%白噪声序列
theta=[a（2:
na+1）;b];%对象参数真值
thetae_仁zeros（na+nb+1,1）;%thetae初值
P=10A6*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','b_0','b_1'）;axis（[0L-22]）;
FFRLS）
方法三：
%遗忘因子递推最小二乘参数估计（
clearall;closeall;
a=[1-1.50.7]';b=[10.5]';d=3;%对象参数
na=length（a）-1;nb=length（b）-1;%na、nb为A、B阶次
L=1000;%仿真长度
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）yk=zeros（na,1）;%输出初值
u=randn（L,1）;%输入采用白噪声序列xi=sqrt（O.1）*randn（L,1）;%白噪声序列
thetae_仁zeros（na+nb+1,1）;%thetae初值
P=10A6*eye（na+nb+1）;
lambda=0.98;%遗忘因子范围[0.91]
fork=1:
L
ifk==501
a=[1-10.4]';b=[1.50.2]';%对象参数突变
end
theta（:
k）=[a（2:
na+1）;b];%对象参数真值
phi=[-yk;uk（d:
d+nb）];
y（k）=phi'*theta（:
k）+xi（k）;%采集输出数据
%遗忘因子递推最小二乘法
K=P*phi/（lambda+phi'*P*phi）;
thetae（:
k）=thetae_1+K*（y（k）-phi'*thetae_1）;
P=（eye（na+nb+1）-K*phi'）*P/lambda;
%更新数据
thetae_仁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
subplot（1,2,1）
plot（[1:
L],thetae（1:
na,:
））;holdon;plot（[1:
L],theta（1:
na,:
）,'k:
'）;xlabel（'k'）;ylabel（'参数估计a'）;
legend（'a_1','a_2'）;axis（[0L-22]）;
subplot（1,2,2）
plot（[1:
L],thetae（na+1:
na+nb+1,:
））;holdon;plot（[1:
L],theta（na+1:
na+nb+1,:
）,'k:
'）;xlabel（'k'）;ylabel（'参数估计b'）;
legend（'b_0','b_1'）;axis（[0L-0.52]）;
方法四：
%递推极大似然参数估计〔RML〕
clearall;closeall;
a=[1-1.50.7]';b=[10.5]';c=[1-0.5]';d=1;%对象参数
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%na、nb、nc为A、B、C阶次
nn=max（na,nc）;%用于yf（k-i）、uf（k-i）更新
L=480;%仿真长度
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）yk=zeros（na,1）;%输出初值xik=zeros（nc,1）;%白噪声初值xiek=zeros（nc,1）;%白噪声估计初值yfk=zeros（nn,1）;%yf（k-i）ufk=zeros（nn,1）;%uf（k-i）xiefk=zeros（nc,1）;%Ef（k-i）u=randn（L,1）;%输入采用白噪声序列xi=randn（L,1）;%白噪声序列
thetae_仁zeros（na+nb+1+nc,1）;%参数估计初值
P=eye（na+nb+1+nc）;
fork=1:
L
y（k）=-a（2:
na+1）'*yk+b'*uk（d:
d+nb）+c'*[xi（k）;xik];%采集输出数据
%构造向量
phi=[-yk;uk（d:
d+nb）;xiek];xie=y（k）-phi'*thetae_1;
phif=[-yfk（1:
na）;ufk（d:
d+nb）;xiefk];
%递推极大似然参数估计算法
K=P*phif/（1+phif*P*phif）;thetae（:
k）=thetae_1+K*xie;
P=（eye（na+nb+1+nc）-K*phif）*P;
yf=y（k）-thetae（na+nb+2:
na+nb+1+nc,k）'*yfk（1:
nc）;%yf（k）uf=u（k）-thetae（na+nb+2:
na+nb+1+nc,k）'*ufk（1:
nc）;%uf（k）xief=xie-thetae（na+nb+2:
na+nb+1+nc,k）'*xiefk（1:
nc）;%xief（k）
%更新数据
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）;
fori=nc:
-1:
2
xik（i）=xik（i-1）;xiek（i）=xiek（i-1）;
xiefk（i）=xiefk（i-1）;
end
xik
（1）=xi（k）;
xiek
（1）=xie;
xiefk
（1）=xief;
fori=nn:
-1:
2
yfk（i）=yfk（i-1）;
ufk（i）=ufk（i-1）;
end
yfk
（1）=yf;
ufk
（1）=uf;
end
figure
（1）
plot（[1:
L],thetae（1:
na,:
）,[1:
L],thetae（na+nb+2:
na+nb+1+nc,:
））;xlabel（'k'）;ylabel（'参数估计a、c'）;
legend（'a_1','a_2','c_1'）;axis（[0L-22]）;
figure
（2）
plot（[1:
L],thetae（na+1:
na+nb+1,:
））;xlabel（'k'）;ylabel（'参数估计b'）;
legend（'b_0','b_1'）;axis（[0L01.5]）
自适应控制习题
（1）%可调增益MIT-MRAC
clearall;closeall;
h=0.1;L=1OO/h;%数值积分步长、仿真步数
num=[1];den=[111];n=length（den）-1;%对象参数
kp=1;[Ap,Bp,Cp,Dp]=tf2ss（kp*num,den）;%传递函数型转换为状
态空间型
km=1;[Am,Bm,Cm,Dm]=tf2ss（km*num,den）;%参考模型参数
gamma=0.1;%自适应增益
yrO=O;u0=0;eO=O;ymO=O;%初值
xpO=zeros（n,1）;xmO=zeros（n,1）;%状态向量初值
kcO=O;%可调增益初值
r=0.1;yr=r*[ones（1,L/4）-ones（1,L/4）ones（1,L/4）-ones（1,L/4）];%
输入信号
fork=1:
L
time（k）=k*h;
xp（:
k）=xpO+h*（Ap*xpO+Bp*uO）;
yp（k）=Cp*xp（:
k）+Dp*u0;%计算yp
xm（:
k）=xm0+h*（Am*xm0+Bm*yr0）;ym（k）=Cm*xm（:
k）+Dm*yr0;%计算ym
e（k）=ym（k）-yp（k）;%e=ym-ypkc=kcO+h*gamma*eO*ymO;%MIT自适应律u（k）=kc*yr（k）;%控制量
%更新数据
yr0=yr（k）;u0=u（k）;eO=e（k）;ym0=ym（k）;
xp0=xp（:
k）;xm0=xm（:
k）;
kc0=kc;
end
plot（time,ym,'r',time,yp,':
'）;
xlabel（'t'）;ylabel（'y_m（t）、y_p（t）'）;
%axis（[OL*h-1010]）;
legend（'y_m（t）','y_p（t）'）;
⑵%可调增益MIT-MRACclearall;closeall;
h=0.1;L=100/h;%数值积分步长、仿真步数
num=[1];den=[111];n=length（den）-1;%对象参数
kp=1;[Ap,Bp,Cp,Dp]=tf2ss（kp*num,den）;%传递函数型转换为状态空间型km=1;[Am,Bm,Cm,Dm]=tf2ss（km*num,den）;%参考模型参数
gamma=0.1;%自适应增益
yr0=0;u0=0;e0=0;ym0=0;%初值
xp0=zeros（n,1）;xm0=zeros（n,1）;%状态向量初值
kc0=0;%可调增益初值
r=1;yr=r*[ones（1,L/4）-ones（1,L/4）ones（1,L/4）-ones（1,L/4）];%输入信号
fork=1:
L
time（k）=k*h;
xp（:
k）=xp0+h*（Ap*xp0+Bp*u0）;
yp（k）=Cp*xp（:
k）+Dp*u0;%计算yp
xm（:
k）=xm0+h*（Am*xm0+Bm*yr0）;
ym（k）=Cm*xm（:
k）+Dm*yrO;%计算ym
e（k）=ym（k）-yp（k）;%e=ym-ypkc=kcO+h*gamma*eO*ymO;%MIT自适应律u（k）=kc*yr（k）;%控制量
%更新数据
yr0=yr（k）;u0=u（k）;e0=e（k）;ym0=ym（k）;
xp0=xp（:
k）;xm0=xm（:
k）;
kc0=kc;
end
plot（time,ym,'r',time,yp,':
'）;
xlabel（'t'）;ylabel（'y_m（t）、y_p（t）'）;
%axis（[0L*h-1010]）;
legend（'y_m（t）','y_p（t）'）;
⑶
（1）%可调增益MIT-MRAC归一化算法
clearall;closeall;
h=0.1;L=100/h;%数值积分步长和仿真步数
num=[1];den=[111];n=length（den）-1;%对象参数
kp=1;[Ap,Bp,Cp,Dp]=tf2ss（kp*num,den）;%传递函数型转换为状态空间型
km=1;[Am,Bm,Cm,Dm]=tf2ss（km*num,den）;%参考模型参数
gamma=0.1;%自适应增益
alpha=0.01;beta=2;
yr0=0;u0=0;e0=0;ym0=0;%初值
xpO=zeros（n,1）;xm0=zeros（n,1）;%状态向量初值
kc0=0;%可调增益初值
r=0.1;yr=r*[ones（1,L/4）-ones（1,L/4）ones（1,L/4）-ones（1,L/4）];%输入信号
fork=1:
L
time（k）=k*h;
xp（:
k）=xpO+h*（Ap*xpO+Bp*uO）;
yp（k）=Cp*xp（:
k）+Dp*uO;%计算yp
xm（:
k）=xmO+h*（Am*xmO+Bm*yrO）;ym（k）=Cm*xm（:
k）+Dm*yr0;%计算ym
e（k）=ym（k）-yp（k）;%e=ym-yp
DD=e0*ym0/km/（alpha+（ym0/km）A2）;
ifDD<-beta
DD=-beta;
end
ifDD>beta
DD=beta;
end
kc=kcO+h*gamma*DD;%MIT自适应律u（k）=kc*yr（k）;%控制量
%更新数据
yr0=yr（k）;u0=u（k）;e0=e（k）;ym0=ym（k）;xp0=xp（:
k）;xm0=xm（:
k）;
kc0=kc;
end
plot（time,ym,'r',time,yp,':
'）;
xlabel（'t'）;ylabel（'y_m（t）、y_p（t）'）;
（2）%可调增益MIT-MRAC归一化算法
clearall;closeall;
h=0.1;L=1OO/h;%数值积分步长和仿真步数
num=[1];den=[111];n=length（den）-1;%对象参数
kp=1;[Ap,Bp,Cp,Dp]=tf2ss（kp*num,den）;%传递函数型转换为状态空间型km=1;[Am,Bm,Cm,Dm]=tf2ss（km*num,den）;%参考模型参数
gamma=0.1;%自适应增益
alpha=0.01;beta=2;
yrO=O;u0=0;e0=0;ym0=0;%初值
xpO=zeros（n,1）;xmO=zeros（n,1）;%状态向量初值
kc0=0;%可调增益初值
r=1;yr=r*[ones（1,L/4）-ones（1,L/4）ones（1,L/4）-ones（1,L/4）];%输入信号
fork=1:
L
time（k）=k*h;
xp（:
k）=xpO+h*（Ap*xpO+Bp*uO）;
yp（k）=Cp*xp（:
k）+Dp*uO;%计算yp
xm（:
k）=xmO+h*（Am*xmO+Bm*yrO）;
ym（k）=Cm*xm（:
k）+Dm*yr0;%计算yme（k）=ym（k）-yp（k）;%e=ym-yp
DD=e0*ym0/km/（alpha+（ym0/km）A2）;
ifDD<-beta
DD=-beta;
end
ifDD>beta
DD=beta;
end
kc=kc0+h*gamma*DD;%MIT自适应律
u（k）=kc*yr（k）;%控制量
%更新数据
yr0=yr（k）;u0=u（k）;e0=e（k）;ym0=ym（k）;
xp0=xp（:
k）;xm0=xm（:
k）;
kc0=kc;
end
plot（time,ym,'r',time,yp,':
'）;
xlabel（'t'）;ylabel（'y_m（t）、y_p（t）'）;
legend（'y_m（t）','y_p（t）'）;
升级会员 