自适应控制程序.docx
《自适应控制程序.docx》由会员分享，可在线阅读，更多相关《自适应控制程序.docx（38页珍藏版）》请在冰豆网上搜索。
自适应控制程序
%M序列及其逆序列的产生
设M序列{M（k）}由如下4位移位寄存器产生：
{S（k）}为方波序列，逆M序列{IM（k）={M（k）S（k）}
clearall;closeall;
L=60;%序列长度
x1=1;x2=1;x3=1;x4=0;%移位寄存器初值
S=1;%方波初值
fork=1:
L
IM=xor（S,x4）;%进行异或运算，产生逆M序列
ifIM==0
u（k）=-1;
else
u（k）=1;
end
S=not（S）;%产生方波
M（k）=xor（x3,x4）;%产生M序列
x4=x3;x3=x2;x2=x1;x1=M（k）;%寄存器移位
end
subplot（2,1,1）;
stairs（M）;grid;
axis（[0L/2-0.51.5]）;xlabel（'k'）;ylabel（'M序列幅值'）;title（'M序列'）;
subplot（2,1,2）;
stairs（u）;grid;
axis（[0L-1.51.5]）;xlabel（'k'）;ylabel（'逆M序列幅值'）;title（'逆M序列'）;
%白噪声及有色噪声序列的产生
设（k）为均值为0，方差为1的高斯白噪声序列，e（k）为有色噪声序列：
高斯白噪声序列（k）在Matlab中由rand（）函数产生，程序如下：
clearall;closeall;
L=500;%仿真长度
d=[1-1.50.70.1];c=[10.50.2];%分子分母多项式系数
nd=length（d）-1;nc=length（c）-1;%阶次
xik=zeros（nc,1）;%白噪声初值
ek=zeros（nd,1）;
xi=randn（L,1）;%产生均值为0，方差为1的高斯白噪声序列
fork=1:
L
e（k）=-d（2:
nd+1）*ek+c*[xi（k）;xik];%产生有色噪声
%数据更新
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）;
end
subplot（2,1,1）;
plot（xi）;
xlabel（'k'）;ylabel（'噪声幅值'）;title（'白噪声序列'）;
subplot（2,1,2）;
plot（e）;
xlabel（'k'）;ylabel（'噪声幅值'）;title（'有色噪声序列'）;
%批处理最小二乘参数估计（LS）
考虑如下系统：
式中（k）为方差为1的白噪声。
clearall;
a=[1-1.50.7]';b=[10.5]';d=3;%对象参数
na=length（a）-1;nb=length（b）-1;%计算阶次
L=500;%数据长度
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
thetaevaluation=inv（phi'*phi）*phi'*y'%计算参数估计值
thetaevaluation=
-1.5362
0.6802
1.0068
0.4864
%遗忘因子递推最小二乘参数估计（FFRLS）
考虑如下系统：
式中（k）为均值为0、方差为0.1的白噪声，
对象时变参数
为：
取遗忘因子=0.98，
clearall;closeall;
a=[1-1.50.7]';b=[10.5]';d=3;%对象参数
na=length（a）-1;nb=length（b）-1;%计算阶次
L=1000;%数据长度
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）;
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_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（2,1,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（2,1,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]）;
%增广递推最小二乘参数估计（ERLS）
考虑如下系统：
式中（k）为方差为0.1的白噪声。
选择方差为1的白噪声作为输入信号u（k）.
clearall;closeall;
a=[1-1.50.7]';b=[10.5]';c=[1-10.2]';d=3;%对象参数
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）;%噪声初值
xiek=zeros（nc,1）;%噪声估计初值
u=randn（L,1）;%输入采用方差为1的白噪声序列
xi=sqrt（0.1）*randn（L,1）;%方差为0.1的白噪声干扰序列
theta=[a（2:
na+1）;b;c（2:
nc+1）];%对象参数真值
thetae_1=zeros（na+nb+1+nc,1）;%参数初值,na+nb+1+nc为辨识参数个数
P=10^6*eye（na+nb+1+nc）;
lambda=0.98;%遗忘因子范围[0.91]
fork=1:
L
phi=[-yk;uk（d:
d+nb）;xik];%测量向量
y（k）=phi'*theta+xi（k）;%采样输出数据
phie=[-yk;uk（d:
d+nb）;xiek];%估计的测量向量
%增广递推最小二乘公式
K=P*phie/（1+phie'*P*phie）;
thetae（:
k）=thetae_1+K*（y（k）-phie'*thetae_1）;
P=（eye（na+nb+1+nc）-K*phie'）*P;
xie=y（k）-phie'*thetae（:
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）;
end
xik
（1）=xi（k）;
xiek
（1）=xie;
end
figure
（1）
plot（[1:
L],thetae（1:
na,:
））;
xlabel（'k'）;ylabel（'参数估计a'）;
legend（'a_1','a_2'）;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]）;
figure（3）
plot（[1:
L],thetae（na+nb+2:
na+nb+nc+1,:
））;
xlabel（'k'）;ylabel（'参数估计c'）;
legend（'c_1','c_2'）;axis（[0L-22]）;
递推最小二乘参数估计（RLS）
考虑如下系统：
式中（k）为方差为0.1的白噪声。
clearall;closeall;
a=[1-1.50.7]';b=[10.5]';d=3;%对象参数
na=length（a）-1;nb=length（b）-1;%计算阶次,na=2,nb=1
L=500;%数据长度（仿真长度）
uk=zeros（d+nb,1）;yk=zeros（na,1）;%输入输出初值uk:
4x1,ykx1
u=randn（L,1）;%输入采用方差为1的白噪声序列
xi=sqrt（0.1）*randn（L,1）;%方差为0.1的白噪声干扰序列
theta=[a（2:
na+1）;b];%对象参数真值theta=[-1.5,0.7;1,0.5]
thetae_1=zeros（na+nb+1,1）;%参数初值为4x1的全零矩阵
P=10^6*eye（na+nb+1）;
fork=1:
L
phi=[-yk;uk（d:
d+nb）];%此处phi为列向量4x1
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]）;
%最小方差控制（MVC）
考虑如下系统：
式中（k）为方差为0.1的白噪声。
取期望输出yr（k）为幅值为10的方波信号。
clearall;closeall;
a=[1-1.70.7];b=[10.5];c=[10.2];d=4;%对象参数
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%计算阶次
nh=nb+d-1;%nh为多项式H的阶次
L=400;
uk=zeros（d+nb,1）;
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）;%方差为0.1的白噪声序列
[h,f,g]=singlediophantine（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）=（-h（2:
nh+1）*uk（1:
nh）+c*[yr（k+d:
-1:
k+d-min（d,nc））;yrk（1:
nc-d）]-g*[y（k）;yk（1:
na-1）]）/h
（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,u）;
xlabel（'k'）;ylabel（'u（k）'）;
单步Diophantine方程求解
求解下列系统的单步Diophantine方程：
（1）
（2）
（3）
%单步Diophantine方程的求解
clearall;
a=[1-1.50.7];b=[10.5];c=[1];d=3;%例4.1
（1）
%a=[1-0.95];b=[12];c=[1-0.7];d=2;%例4.1
（2）
%a=[1-1.70.7];b=[0.91];c=[12];d=4;%例4.1（3）
[e,f,g]=sindiophantine（a,b,c,d）%调用函数sindiophantine
function[e,f,g]=singlediophantine（a,b,c,d）%单步Diophantine方程求解
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%计算阶次
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-1）*ad（i+j）;
end
g（i）=cd（i+d）-g（i）;%计算gi
end
f=conv（b,e）;%计算F
e=
1.00001.50001.5500
f=
1.00002.00002.30000.7750
g=
1.2750-1.0850
多步Diophantine方程求解
求解如下系统的多步Diophantine方程，并去预测长度N=3
%多步Diophantine方程的求解
clearall;
a=[1-2.72.4-0.7];b=[0.91];c=[12];
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%A、B、C的阶次
N=3;%预测步数
[E,F,G]=multidiophantine（a,b,c,N）%调用函数multidiophantine
function[E,F,G]=multidiophantine（a,b,c,N）
%********************************************************
%功能：
多步Diophanine方程的求解
%调用格式：
[E,F,G]=sindiophantine（a,b,c,N）（注：
d=1）
%输入参数：
多项式A、B、C系数向量及预测步数（共4个）
%输出参数：
Diophanine方程的解E、F、G（共3个）
%********************************************************
na=length（a）-1;nb=length（b）-1;nc=length（c）-1;%A、B、C的阶次
%E、F、G的初值
E=zeros（N）;E（1,1）=1;F（1,:
）=conv（b,E（1,:
））;
ifna>=nc
G（1,:
）=[c（2:
nc+1）zeros（1,na-nc）]-a（2:
na+1）;%令c（nc+2）=c（nc+3）=...=0
else
G（1,:
）=c（2:
nc+1）-[a（2:
na+1）zeros（1,nc-na）];%令a（na+2）=a（na+3）=...=0
end
%求E、G、F
forj=1:
N-1
fori=1:
j
E（j+1,i）=E（j,i）;
end
E（j+1,j+1）=G（j,1）;
fori=2:
na
G（j+1,i-1）=G（j,i）-G（j,1）*a（i）;
end
G（j+1,na）=-G（j,1）*a（na+1）;
F（j+1,:
）=conv（b,E（j+1,:
））;
end
%最小方差自校正控制
用最小方差自校正控制算法对以下系统进行闭环控制：
式中（k）为方差为0.1的白噪声。
取期望输出yr（k）为幅值为10的方波信号。
%最小方差间接自校正控制
clearall;closeall;
a=[1-1.70.7];b=[10.5];c=[10.2];d=4;%对象参数
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）;%白噪声初值
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）;%白噪声序列
%RELS初值设置
thetae_1=0.001*ones（na+nb+1+nc,1）;%非常小的正数（这里不能为0）
P=10^6*eye（na+nb+1+nc）;
fork=1:
L
time（k）=k;
y（k）=-a（2:
na+1）*yk+b*uk（d:
d+nb）+c*[xi（k）;xik];%采集输出数据
%递推增广最小二乘法
phie=[-yk;uk（d:
d+nb）;xiek];
K=P*phie/（1+phie'*P*phie）;
thetae（:
k）=thetae_1+K*（y（k）-phie'*thetae_1）;
P=（eye（na+nb+1+nc）-K*phie'）*P;
xie=y（k）-phie'*thetae（:
k）;%白噪声的估计值
%提取辨识参数
ae=[1thetae（1:
na,k）'];be=thetae（na+1:
na+nb+1,k）';ce=[1thetae（na+nb+2:
na+nb+1+nc,k）'];
ifabs（be
（2））>0.9
be
（2）=sign（ce
（2））*0.9;%MVC算法要求B稳定
end
ifabs（ce
（2））>0.9
ce
（2）=sign（ce
（2））*0.9;
end
[e,f,g]=sindiophantine（ae,be,ce,d）;%求解单步Diophantine方程
u（k）=（-f（2:
nf+1）*uk（1:
nf）+ce*[yr（k+d:
-1:
k+d-min（d,nc））;yrk（1:
nc-d）]-g*[y（k）;yk（1:
na-1）]）/f
（1）;%求控制量
%更新数据
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
yrk（i）=yrk（i-1）;
xik（i）=xik（i-1）;
xiek（i）=xiek（i-1）;
end
ifnc>0
yrk
（1）=yr（k）;
xik
（1）=xi（k）;
xiek
（1）=xie;
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（
升级会员 