系统辨识仿真作业.docx
- 文档编号:23422615
- 上传时间:2023-05-16
- 格式:DOCX
- 页数:11
- 大小:70.77KB
系统辨识仿真作业.docx
《系统辨识仿真作业.docx》由会员分享,可在线阅读,更多相关《系统辨识仿真作业.docx(11页珍藏版)》请在冰豆网上搜索。
系统辨识仿真作业
考虑图1所示仿真对象,其中
为服从
正态分布的不相关随机噪声;输入信号采用周期大于500的逆M序列,幅度为1;λ可用于控制数据的信噪比。
图1仿真对象框图
选择如下的模型结构
数据长度取
;初始条件取
,要求:
1)根据仿真模型编程生成辨识所需数据,绘出各数据随时间的变化曲线(
可取定值,如
等);
2)利用广义最小二乘递推算法对模型参数进行辨识,绘出各参数随时间的变化曲线;
3)试分析不同信噪比的数据对辨识结果的影响。
第一问解答:
取
=1,实际上c1、c2的真值应该为零,所以可取c1=0,c2=0,可得e(k)=v(k),采用最小二乘递推辨识方法。
程序如下:
function[seq]=mseries(con,reg,len);
con=[11000011];%多项式系数方程
m=length(con);
reg=[10100011];%寄存器初始值
len=810;%一个周期
fan=0;
fori=1:
len
seq(i)=reg(m);
forj=1:
m
fan=fan+con(j)*reg(j);
fan=(mod(fan,2));
end
fort=m:
-1:
2%寄存器移位
reg(t)=reg(t-1);
end
reg
(1)=fan;
fan=0;
end
s=[];
fori0=1:
(len/2)
s=[s,[1,0]];
end
u=xor(seq,s);%生成一个周期的m逆序列
v=randn(1,800);%产生均值为1方差为0的高斯白噪声
e=v;
z=zeros(802,1);
z
(1)=0;
z
(2)=0;
fork=3:
802
z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+e(k-2);
end
p=10^6*eye(4);%给出初始状态p
pstore=zeros(4,802);
pstore(:
1)=[p(1,1),p(2,2),p(3,3),p(4,4)];
theta=zeros(4,802);
theta(:
1)=[0.0010.0010.0010.001]';%辨识参数初始值
K=zeros(4,800);
fori1=3:
802
h=[-z(i1-1);-z(i1-2);u(i1-1);u(i1-2)];
K=p*h*inv(h'*p*h+1);
theta(:
i1-1)=theta(:
i1-2)+K*(z(i1)-h'*theta(:
i1-2));
p=(eye(4)-K*h')*p;
pstore(:
1)=[p(1,1),p(2,2),p(3,3),p(4,4)];%存储上一次数据
end
k=1:
1:
802;
figure
(1);
plot(k,theta(1,:
),'b');
holdon;
gridon;
k=1:
1:
802;
plot(k,theta(2,:
),'r');
holdon;
gridon;
k=1:
1:
802;
plot(k,theta(3,:
),'g');
holdon;
gridon;
k=1:
1:
802;
plot(k,theta(4,:
),'y');
holdon;
gridon;
title('最小二乘法递推辨识结果');
text(780,1.0,'b1');
text(780,0.7,'a2');
text(780,0.5,'b2');
text(780,-1.5,'a1');
xlabel('k');
ylabel('theta');
仿真结果如下:
第二问解答:
程序如下:
function[seq]=mseries(con,reg,len);
con=[11000011];%多项式系数方程
m=length(con);
reg=[10100011];%寄存器初始值
len=810;%一个周期
fan=0;
fori=1:
len
seq(i)=reg(m);
forj=1:
m
fan=fan+con(j)*reg(j);
fan=(mod(fan,2));
end
fort=m:
-1:
2%寄存器移位
reg(t)=reg(t-1);
end
reg
(1)=fan;
fan=0;
end
s=[];
fori0=1:
(len/2)
s=[s,[1,0]];
end
uf=xor(seq,s);%生成一个周期的m逆序列
v=randn(1,802);%产生均值为1方差为0的高斯白噪声
zf=zeros(802,1);
zf
(1)=0;
zf
(2)=0;
fork=3:
802
zf(k)=1.5*zf(k-1)-0.7*zf(k-2)+uf(k-1)+0.5*uf(k-2)+v(k);%在此处修改信噪比可观
%测出信噪比对辨识结果的影响
end
e=zeros(802,1);
p=10^6*eye(4);%给出初始状态p
pstore=zeros(4,802);
pstore(:
1)=[p(1,1),p(2,2),p(3,3),p(4,4)];
thetak=zeros(4,802);
thetak(:
1)=[0.0010.0010.0010.001]';%辨识参数初始值
Kf=zeros(4,800);
pe=eye
(2);
pstoree=zeros(2,802);
pstoree(:
1)=[pe(1,1),pe(2,2)];
thetae=zeros(2,802);
thetae(:
1)=[0.0010.001]';
Ke=zeros(2,800);
fori1=3:
802
hf=[-zf(i1-1);-zf(i1-2);uf(i1-1);uf(i1-2)];
Kf=p*hf*inv(hf'*p*hf+1);
thetak(:
i1-1)=thetak(:
i1-2)+Kf*(zf(i1)-hf'*thetak(:
i1-2));
p=(eye(4)-Kf*hf')*p;
pstore(:
1)=[p(1,1),p(2,2),p(3,3),p(4,4)];%存储上一次数据
e(k)=zf(k)-hf'*thetak(:
i1-1);
he=[-e(i1-1);-e(i1-2)];
Ke=pe*he*inv(he'*pe*he+1);
thetae(:
i1-1)=thetae(:
i1-2)+Ke*(e(i1)-he'*thetae(:
i1-2));
pe=(eye
(2)-Ke*he')*pe;
pstoree(:
1)=[pe(1,1),pe(2,2)];
end
k=1:
1:
802;
figure
(2);
plot(k,thetak(1,:
),'b');
holdon;
gridon;
k=1:
1:
802;
plot(k,thetak(2,:
),'r');
holdon;
gridon;
k=1:
1:
802;
plot(k,thetak(3,:
),'g');
holdon;
gridon;
k=1:
1:
802;
plot(k,thetak(4,:
),'y');
holdon;
gridon;
k=1:
1:
802;
plot(k,thetae(1,:
),'b');
holdon;
gridon;
k=1:
1:
802;
plot(k,thetae(2,:
),'r');
holdon;
gridon;
title('广义最小二乘法辨识结果');
text(780,1.0,'b1');
text(780,0.7,'a2');
text(780,0.5,'b2');
text(780,-1.5,'a1');
text(760,0,'c1');
text(700,0,'c2');
xlabel('k');
ylabel('参数');
仿真结果如下:
第三问解答:
程序如下:
可在第二问的基础上修改信噪比观测信噪比对辨识结果的影响。
由如下仿真结果可知信噪比越小,参数收敛速度越快,且与真实值越接近,信噪比越大,收敛速度越慢,与真实值的接近状况不如信噪比小的时候。
仿真结果如下:
为0时
为0.5时
为1时
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 系统 辨识 仿真 作业