系统辨识作业2doc.docx
《系统辨识作业2doc.docx》由会员分享，可在线阅读，更多相关《系统辨识作业2doc.docx（14页珍藏版）》请在冰豆网上搜索。
系统辨识作业2doc
系统辨识作业
学院：
专业：
姓名：
学号：
日期：
系统辨识作业：
以下图为仿真对象
图中，v（k）为服从N（0,1）正态分布的不相关随即噪声，输入信号采用循环周期Np>500的逆M序列，幅值为1，选择辨识模型为：
加权因子
，数据长度L=500，初始条件取
，
要求：
（1）采用一次完成最小二乘法对系统进行辨识，给出数据u（k）和z（k），及
，
和
和
的值。
（2）采用递推最小二乘法进行辨识，要给出参数收敛曲线以及新息
，残差
，准则函数
随着递推次数K的变化曲线。
（3）对仿真对象和辨识出的模型进行阶跃响应对比分析以检验辨识结果的实效。
1、一次完成法对系统进行辨识：
估计值
，其中
一次完成算法对系统辨识的Matlab程序见附录:
部分输入、输出数据如下，全部的输入输出数据用图1.1所示
输入数据u（k）=Columns1through16
00111100001
00110
输出数据z（k）=Columns1through9
001.23723.93316.49877.99097.71326.59475.4523
Columns10through18
3.22120.84190.62431.71100.71261.07122.80373.80474.6372
图1.1输入数据u（k）和输出数据z（k）
的值为（部分）：
HL=
-5.4523-6.594700
-3.2212-5.452300
-0.8419-3.22121.00000
…………
-3.5706-5.194400
-0.6787-3.570600
2.3020-0.678700
的值为（部分）：
ZL=
3.2212
0.8419
0.6243
1.7110
0.7126
…
0.6787
-2.3020
-3.8270
一次完成辨识后
的值为：
theta=-1.49180.69321.05410.4691
的值为：
J（theta）=493.1837
2、递推最小二乘法辨识系统：
初始条件：
，
递推最小二乘法辨识系统的Matlab程序见附录:
其参数收敛曲线如图2.1
图2.1参数收敛曲线
新息
，残差
，准则函数
随着递推次数K的变化曲线如图2.2中依次所示:
图2.2新息、残差、准则函数变化曲线
3、仿真对象和辨识出的模型进行阶跃响应对比分析
仿真对象和辨识模型阶跃响应对比Matlab程序见附录:
图3.1分别给出了一次完成算法辨识出来系统和辨识对象的阶跃响应对比图，递推算法辨识出来系统和辨识对象的阶跃响应对比图。
附录：
一次完成和递推法系统辨识Matlab程序
%%最小二乘法辨识系统;
%叠加噪声为1/（1-1.5z^（-1）+0.7z^（-2））;
%化为差分方程形式为;
%e（k）=v（k）+1.5e（k-1）-0.7e（k-2）;
%仿真对象为（1z^-1+0.5z^-2）/（1-1.5z^-1+0.7z^-2）;
%化为差分方程形式为;
%y（k）-1.5y（k-1）+0.7y（k-2）=u（k-1）+0.5u（k-2）;
%辨识模型为;
%z（k）=-a1z（k-1）-a2z（k-2）+b1u（k-1）+b2u（k-2）+v（k）;
%================================
%主函数;
functionmain
closeall;
clc;
clear;
%================
%产生逆M序列;
X=[0,1,1,0,1,0,0,1];%寄存器初始值;
F=[0,1,1,1,0,0,0,1];%特征多项式;
N=1000;%生成长度;
M=[];%存放产生的M序列;
%产生逆M序列函数调用;
out=IMxulie（X,F,N）;
%阶梯图输出逆M序列;
figure
（1）;
M=out;
subplot（2,1,1）;
stairs（M）;
xlabel（'k'）
ylabel（'逆M序列'）
title（'移位寄存器产生的逆M序列'）;
gridon;
%=================
%一次完成最小二乘法辨识系统;
%e（k）=v（k）+1.5e（k-1）-0.7e（k-2）;
%y（k）=1.5y（k-1）-0.7y（k-2）+u（k-1）+0.5u（k-2）;
%z（k）=y（k）+e（k）;
%即z（k）=1.5y（k-1）-0.7y（k-2）+u（k-1）+0.5u（k-2）+v（k）+1.5e（k-1）-0.7e（k-2）;
%产生均值为0，方差为1的白噪声;
v=[];
v=randn（1,600）;
%产生系统辨识所需要数据;
y
（1）=0;
y
（2）=0;
e
（1）=0;
e
（2）=0;
fori=3:
600
y（i）=1.5*y（i-1）-0.7*y（i-2）+M（i-1）+0.5*M（i-2）;
e（i）=v（i）+1.5*e（i-1）-0.7*e（i-2）;
z（i）=y（i）+e（i）;
end
subplot（2,1,2）;
stem（z）;
xlabel（'k'）;
ylabel（'幅度值'）;
title（'输出数据'）;
gridon;
%辨识模型为;
%z（k）=-a1z（k-1）-a2z（k-2）+b1u（k-1）+b2u（k-2）+v（k）;
fori=10:
509
H（i-9,:
）=[-z（i-1）,-z（i-2）,M（i-1）,M（i-2）];
end
ZL=（z（10:
509））';
ES=inv（H'*H）*H'*ZL;%inv用来求矩阵的逆矩阵;
%输出辨识的参数;
disp（'输入数据u（k）='）;
disp（M）;
disp（'输出数据z（k）='）;
disp（z）;
disp（'HL='）;
disp（H）;
disp（'ZL='）;
disp（ZL）;
disp（'theta='）;
disp（ES）;
%由估计值计算准则函数的值;
J=（ZL-H*ES）'*（ZL-H*ES）;
disp（'J（theta）='）;
disp（J）;
%================================
%递推最小二乘法辨识系统;
%一式s（k）=s（k-1）+k（k）（z（k）-h'（k）s（k-1））;
%二式k（k）=p（k-1）h（k）[h'（k）p（k-1）h（k）+1/w（k）]';
%三式p（k）=[I-k（k）h'（k）]p（k-1）;
%s（0）=[0.01,0.01,...0.01]';
%p（0）=10^6I;
%新息ZX=z（k）-H'（k）S（k-1）;
%残差E=z（k）-H'（k）S（k）;
%准则函数J（k）
%递推求解;
z1=z（10:
509）;
M1=M（10:
509）;
P=10^6*eye（4）;
S（:
1）=[0.01;0.01;0.01;0.01]';
S（:
2）=[0.01;0.01;0.01;0.01]';
J1
（1）=0;
J1
（2）=0;
fori=3:
500
h=[-z1（i-1）;-z1（i-2）;M1（i-1）;M1（i-2）];
K=P*h*inv（h'*P*h+1）;
S（:
i）=S（:
i-1）+K*（z1（i）-h'*S（:
i-1））;%待估计参数变化;
zx（i）=z1（i）-h'*S（:
i-1）;%新息变化;
E（i）=z1（i）-h'*S（:
i）;%残差变化;
J1（i）=J1（i-1）+（zx（i）*zx（i））/（h'*P*h+1）;%准则函数变化;
P=（eye（4,4）-K*h'）*P;
end
%待估计参数过度曲线;
figure
（2）;
i=1:
500;
plot（i,S（1,:
）,i,S（2,:
）,i,S（3,:
）,i,S（4,:
））;
title（'待估计参数过度过程'）;
legend（'参数a1的过渡过程','参数a2的过渡过程','参数b1的过渡过程','参数b2的过渡过程'）;
gridon;
%disp（S（:
500））;
%新息变化曲线;
figure（3）;
subplot（3,1,1）;
i=1:
500;
plot（i,zx（i））;
title（'新息变化曲线'）;
gridon;
%残差变化曲线;
subplot（3,1,2）;
i=1:
500;
plot（i,E（i））;
title（'残差变化曲线'）;
gridon;
%准则函数变化曲线;
subplot（3,1,3）;
i=1:
500;
plot（i,J1（i））;
title（'准则函数变化曲线'）;
gridon;
%disp（J1（500））;
%================================
%分别对一次完成，递推和辨识对象阶跃响应对比;
%%一次完成法辨识阶跃响应对比;
%一次辨识后对仿真对象和辨识模型进行阶跃响应对比;
figure（4）;
subplot（2,1,1）;
%辨识模型阶跃响应;
%z（k）+a1z（k-1）+a2z（k-2）=b1u（k-1）+b2u（k-2）;
a=[1,ES
（1）,ES
（2）];
b=[0,ES（3）,ES（4）];
%求离散系统阶跃响应;
gn=dstep（b,a,30）;
stem（gn,'ro'）;
gridon;
holdon;
%加阶跃响应后辨识对象的系统函数为：
%（1+1z^-1+0.5z^-2）/（1-1.5z^-1+0.7z^-2）;
b1=[0,1.0,0.5];
a1=[1,-1.5,0.7];
gn1=dstep（b1,a1,30）;
stem（gn1,'bx'）;
holdoff;
xlabel（'k'）;
ylabel（'幅值h（k）'）;
title（'辨识对象和一次完成法辨识系统阶跃响应对比'）;
legend（'估计值','理论值'）;
%%递推法辨识阶跃响应对比;
subplot（2,1,2）;
%辨识模型阶跃响应;
%z（k）+a1z（k-1）+a2z（k-2）=b1u（k-1）+b2u（k-2）;
a=[1,S（1,500）,S（2,500）];
b=[0,S（3,500）,S（4,500）];
%求离散系统阶跃响应;
gn=dstep（b,a,30）;
stem（gn,'ro'）;
gridon;
holdon;
%加阶跃响应后辨识对象的系统函数为：
%（1+1z^-1+0.5z^-2）/（1-1.5z^-1+0.7z^-2）;
b1=[0,1.0,0.5];
a1=[1,-1.5,0.7];
gn1=dstep（b1,a1,30）;
stem（gn1,'bx'）;
holdoff;
xlabel（'k'）;
ylabel（'幅值h（k）'）;
title（'辨识对象和递推法辨识系统阶跃响应对比'）;
legend（'估计值','理论值'）;
%================================
end
%================================
%================================
%子函数，产生周期大于500的逆M序列函数;
functionout=IMxulie（X,F,N）
n=N;
X1=X;
F1=F;
s=1;
fori=1:
n
Y=X1;
t=F1（8）*Y（8）;
forj=7:
-1:
1
t=xor（t,F1（j）*Y（j））;
end
X1
（1）=t;
fork=2:
8
X1（k）=Y（k-1）;
end
U（i）=Y（8）;
V（i）=xor（U（i）,s）;
s=not（s）;
end
out=V;
end
%================================
升级会员 