系统辨识与自适应作业.docx
《系统辨识与自适应作业.docx》由会员分享，可在线阅读，更多相关《系统辨识与自适应作业.docx（19页珍藏版）》请在冰豆网上搜索。
系统辨识与自适应作业
系统辨识与自适应控制
学院：
电气工程学院
专业：
控制理论与控制工程
学号：
2010020119
姓名：
玉海超
系统辨识与自适应控制作业
一、系统方程为如下方程，e（k）为零均值白噪声
要进行系统辨识并讨论：
定常系统a=0.8,b=0.5参数递推估计
时变系统，λ取不同值时递推辨识结果并讨论
解答：
方法分析
利用批处理法得到系统参数的最小二乘法估计
式中Y=
采用遗忘因子递推最小二乘法，有所学理论知，其参数估计公式如下：
其中：
程序代码如下：
%批处理最小二乘参数估计（LS）
clearall;
a=[10.8]';b=[0.5]';d=1;%对象参数
na=length（a）-1;nb=length（b）-1;%na、nb为A、B阶次
L=300;%数据长度
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）
yk=zeros（na,1）;%输出初值
x1=1;x2=1;x3=1;x4=0;S=1;%移位寄存器初值、方波初值
xi=randn（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）;%产生逆M序列
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=d+nb:
-1:
1
%uk（i）=uk（i-1）;
%end
uk
（1）=u（k）;
%fori=na:
-1:
1
%yk（i）=yk（i-1）;
%end
yk
（1）=y（k）;
end
thetae=inv（phi'*phi）*phi'*y'%计算参数估计值thetae
运行结果如下：
程序代码如下：
%遗忘因子递推最小二乘参数估计（FFRLS）
clearall;closeall;
a=[10.8]';b=[0.5]';d=1;%对象参数
na=length（a）-1;nb=length（b）-1;%na、nb为A、B阶次
L=600;%仿真长度
uk=zeros（d+nb,1）;%输入初值：
uk（i）表示u（k-i）
yk=zeros（na,1）;%输出初值
u=randn（L,1）;%输入采用白噪声序列
xi=sqrt（0.1）*randn（L,1）;%白噪声序列
thetae_1=zeros（na+nb+1,1）;%thetae初值
P=10^6*eye（na+nb+1）;
lambda=0.95;%遗忘因子范围[0.91]
fork=1:
L
ifk==301
a=[10.6]';b=[0.3]';%对象参数突变
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_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
subplot（1,2,1）
plot（[1:
L],thetae（1:
na,:
））;holdon;plot（[1:
L],thetae（1:
na,:
）,'k:
'）;
xlabel（'k'）;ylabel（'参数估计a'）;axis（[0L-0.52]）;
subplot（1,2,2）
plot（[1:
L],thetae（na+1:
na+nb+1,:
））;holdon;plot（[1:
L],thetae（na+1:
na+nb+1,:
）,'k:
'）;
xlabel（'k'）;ylabel（'参数估计b'）;axis（[0L-0.52]）;
仿真结果：
依次为λ=0.95、λ=0.98、λ=1的结果
λ=0.95
当遗忘因子λ=0.95时，从仿真图可以看出，辨识结果很明显。
开始K=0，估计值a的值比较小，而b的值比较大。
随着K的增大估计值a增大，b减小，并且变化率很大，当K=20左右变化比较缓慢了，一直到k=300，a围绕0.8上下波动，b围绕0.6上下波动，波动范围0.2左右；当k>300的时候，由于a和b变化了，所以当变化稳定的时候a就围绕0.6变化，b围绕0.3变化。
总的来说当λ=0.95时，系统达到辨识的要求，但是辨识效果不太理想。
λ=0.98
当遗忘因子λ=0.98时，可以从图中看出和λ=0.95的时候存在一定的差别。
首先是当K=0的时候，a,b的值都很大，随着K的增大，a,b开始减小，当K=300的时候，a,b接近真值，当K大于300的时候，两者围绕变化的数值变为0.6和0.3。
可以看出当λ=0.98的时候辨识结果明显。
λ=1
当λ=1的时候我们称为此时为递推最小二乘法，此时就不带遗忘因子了。
从图中可以从看出当K=0的时候出现的结果与λ=0.98时的一样，但是当K>100后变化缓慢，最终a达到0.8左右，b达到0.5左右。
但是当K>300系统估计值a和b却没有变化到0.6和0.3。
所以可以看出系统辨识结果也不太理想。
二、已知系统方程为
其中e（K）为白噪声，在输入信号为方波时，分析（白噪声为信号幅值的10%，最大为20%）
系统开环响应
在PID控制下系统闭环响应
在最小方差控制下，系统闭环响应（要有分析、比较说明）
解：
增量式PID控制算法，如下式所示
△u（k）=kp（e（k）-e（k-1））+kie（k）+kd（e（k）-2e（k-1）+e（k-2））
u（k）=u
（1）+△u（k）;
式中，kp，ki，kd为PID调节参数，且
△u（k）=u（k）-u（k-1）
e（k）=yr（k）-y（k）
最小方差控制算法：
1）输入初始数据
2）采样当前实际输出y（k）和期望输出Yr（k+d）（或Yr（k））;
3）求解Diophantine方程
，得到多项式E、F、G的系数
4）采用式
计算并实施u（k）;
5）返回第二步（k→k+1）,继续循环
程序代码如下：
c%白噪声及有色噪声序列的产生
clearall;closeall;
L=400;%仿真长度
d=[11.20.35];j=[10.4];c=[11.10.28];a=2;%D、C多项式的系数（可用roots命令求其根）
nd=length（d）-1;nc=length（c）-1;nj=length（j）-1;%nd、nc为D、C的阶次
xik=zeros（nc,1）;%白噪声初值，相当于ξ（k-1）...ξ（k-nc）
ek=zeros（nd,1）;%输出Y初始值
eiu=zeros（a+nj,1）;%输入信号的初始值
xi=randn（L,1）;%randn产生均值为0，方差为1的高斯随机序列（白噪声序列）
err=zeros（L,1）;
u=10*[ones（L/4,1）;-ones（L/4,1）;ones（L/4,1）;-ones（L/4,1）];%产生方波
fork=1:
L
time（k）=k;
e（k）=-d（2:
nd+1）*ek+c*[xi（k）;xik]+j*eiu（a:
a+nj）;%系统的实际输出
err=u（k）-e（k）;%计算偏差
%数据更新
fori=nd:
-1:
2
ek（i）=ek（i-1）;
end
ek
（1）=e（k）;
fori=nc:
-1:
2
xik（i）=xik（i-1）;
end
xik
（1）=xi（k）;
fori=nj+a:
-1:
2
eiu（i）=eiu（i-1）;
end
eiu
（1）=u（k）;
end
subplot（2,1,1）;
plot（time,u（1:
L）,'r:
',time,e）;
xlabel（'t'）;ylabel（'y_r（t）、y（t）'）;
legend（'y_r（t）','y（t）'）;
subplot（2,1,2）;
plot（time,err,'r'）;%输出偏差图形
xlabel（'t'）;ylabel（'e（t）'）;
仿真结果：
程序代码：
%PID控制
clearall;closeall;
a=[11.20.35];b=[10.4];c=[11.10.28];d=2;%对象参数
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%na、nb、nc为多项式A、B、C阶次
kp=0.2;ki=0.5;kd=0.1;%PID控制器参数（试凑法
L=400;%控制步数
uk=zeros（d+nb,1）;%控制初值：
uk（i）表示u（k-i）;
yk=zeros（na,1）;%输出初值
xik=zeros（nc,1）;%白噪声初值
ek=zeros（2,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）;%白噪声序列
fork=1:
L
time（k）=k;
y（k）=-a（2:
na+1）*yk+b*uk（d:
d+nb）+c*[xi（k）;xik];%采集输出数据
e（k）=yr（k）-y（k）;
%增量式PID控制律
du=kp*（e（k）-ek
（1））+ki*e（k）+kd*（e（k）-2*ek
（1）+ek
（2））;
u（k）=uk
（1）+du;
%更新数据
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）;
end
ifnc>0
xik
（1）=xi（k）;
end
fori=d+nb:
-1:
2
ek（i）=ek（i-1）;
end
ek
（1）=e（k）;
end
subplot（2,1,1）;
plot（time,yr（1:
L）,'r:
',time,y）;
xlabel（'t'）;ylabel（'y_r（t）、y（t）'）;
legend（'y_r（t）','y（t）'）;
subplot（2,1,2）;
plot（time,e）;
xlabel（'t'）;ylabel（'e（t）'）;
仿真结果：
程序代码
%最小方差控制（MVC）
clearall;closeall;
a=[11.20.35];b=[10.4];c=[11.10.28];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）;%白噪声序列
ek=zeros（L,1）;
[e,f,g]=sindiophantine（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];%采集输出数据
ek（k）=yr（k）-y（k）;
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）;
xlabel（'k'）;ylabel（'y_r（k）、y（k）'）;
legend（'y_r（k）','y（k）'）;
subplot（2,1,2）;
plot（time,ek）;
xlabel（'k'）;ylabel（'e（k）'）;
function[e,f,g]=sindiophantine（a,b,c,d）
%*********************************************************
%功能：
单步Diophanine方程的求解
%调用格式：
[e,f,g]=sindiophantine（a,b,c,d）
%输入参数：
多项式A、B、C系数（行向量）及纯滞后（共4个）
%输出参数：
Diophanine方程的解e,f,g（共3个）
%*********************************************************
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）];%令a（na+2）=a（na+3）=...=0
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
仿真结果如下：
最小方差控制放大图
分析如下：
从仿真图可以看出，系统的开环响应，实际的输出与期望的误差另很大，超调量也比较高，这主要是开环模型的一个效果。
当系统引进闭环的PID控制时，系统的超调量比较小，上升时间也比较快，最后稳态误差几乎为零，但在信号改变的那一刻，也即波峰与波谷转变的那瞬间，系统的误差也比较大。
PID控制能够达到一定的控制效果，但PID参数的调整比较麻烦，需要反复的试奏才能达到一个满意的控制效果。
从最小方差控制的仿真图来看，刚开始超调量有点大，误差也比较大。
但很快的，随着控制的深入，系统的误差很小，几乎与输入信号的波形重合，在信号发生变化的瞬间，控制效果也很好，但将图形放大，还是有点小误差，但总体上，控制效果还是比较好的。
三、进行基于波波夫稳定性理论的MARS设计及算法仿真
算法的模型图如下
程序算法：
1）选择参考模型，即Gm（s）;
2）选择输入信号yr（t）
3）采样当前参考模型输出ym（t）和系统实际输出yp（t）;
4）利用式
;
计算u（t），其中r为输入信号，v=e*（d0+d1*s）,
为
的初始状态
5）t→t+h,返回第3步，继续循环
考虑如下被控对象模型：
kp未知（仿真时取kp=1）
选择参考模型为：
Km=1
程序代码如下：
%波波夫程序
clearall;closeall;
h=0.01;L=100/h;%数值积分步长和仿真步数（减小h，可以提高积分精度）
num=[21];den=[121];n=length（den）-1;%对象参数（严格正实）\
d0=1;d1=5;
ek=zeros（n,1）;
kp=1;[Ap,Bp,Cp,Dp]=tf2ss（kp*num,den）;%对象参数（传递函数型转换为状态空间型）
km=1;[Am,Bm,Cm,Dm]=tf2ss（km*num,den）;%参考模型参数
yr0=0;u0=0;e0=0;v0=0;%初值
xp0=zeros（n,1）;xm0=zeros（n,1）;%状态向量初值
kc0=0;%可调增益初值
r=2;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）;%计算yp
xm（:
k）=xm0+h*（Am*xm0+Bm*yr0）;
ym（k）=Cm*xm（:
k）;%计算ym
e（k）=ym（k）-yp（k）;%e=ym-yp
v（k）=d0*ek（n）+d1*（（ek（n）-ek（n-1））/h）;
kc=kc0+1*h*v0*yr0+2*v0*yr0;%波波夫算法
u（k）=kc*yr（k）;%控制量
%更新数据
yr0=yr（k）;u0=u（k）;v0=v（k）;
fori=n:
-1:
2
ek（n）=ek（n-1）;
end
ek
（1）=e（k）;
xp0=xp（:
k）;xm0=xm（:
k）;
kc0=kc;
end
subplot（2,1,1）;
plot（time,ym,'r',time,yp,':
'）;
xlabel（'t'）;ylabel（'y_m（t）、y_p（t）'）;
legend（'y_m（t）','y_p（t）'）;
subplot（2,1,2）;
plot（time,u）;
xlabel（'t'）;ylabel（'u（t）'）;
仿真结果：
升级会员 