系统辨识举例.docx
《系统辨识举例.docx》由会员分享，可在线阅读，更多相关《系统辨识举例.docx（29页珍藏版）》请在冰点文库上搜索。
系统辨识举例
例3.3考虑仿真对象
（3.41）
其中，
是服从正态分布的白噪声N
。
输入信号采用4阶M序列，幅度为1。
选择如下形式的辨识模型
（3.42）
设输入信号的取值是从k=1到k=16的M序列，则待辨识参数
为
=
。
其中，被辨识参数
、观测矩阵zL、HL的表达式为
，
（3.43）
程序框图如图3.2所示。
Matlab6.0仿真程序如下：
%二阶系统的最小二乘一次完成算法辨识程序，在光盘中的文件名：
FLch3LSeg1.m
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1];%系统辨识的输入信号为一个周期的M序列
z=zeros（1,16）;%定义输出观测值的长度
fork=3:
16
z（k）=1.5*z（k-1）-0.7*z（k-2）+u（k-1）+0.5*u（k-2）;%用理想输出值作为观测值
end
subplot（3,1,1）%画三行一列图形窗口中的第一个图形
stem（u）%画输入信号u的径线图形
subplot（3,1,2）%画三行一列图形窗口中的第二个图形
i=1:
1:
16;%横坐标范围是1到16，步长为1
plot（i,z）%图形的横坐标是采样时刻i,纵坐标是输出观测值z,图形格式为连续曲线
subplot（3,1,3）%画三行一列图形窗口中的第三个图形
stem（z）,gridon%画出输出观测值z的径线图形，并显示坐标网格
u,z%显示输入信号和输出观测信号
%L=14%数据长度
HL=[-z
（2）-z
（1）u
（2）u
（1）;-z（3）-z
（2）u（3）u
（2）;-z（4）-z（3）u（4）u（3）;-z（5）-z（4）u（5）u（4）;-z（6）-z（5）u（6）u（5）;-z（7）-z（6）u（7）u（6）;-z（8）-z（7）u（8）u（7）;-z（9）-z（8）u（9）u（8）;-z（10）-z（9）u（10）u（9）;-z（11）-z（10）u（11）u（10）;-z（12）-z（11）u（12）u（11）;-z（13）-z（12）u（13）u（12）;-z（14）-z（13）u（14）u（13）;-z（15）-z（14）u（15）u（14）]%给样本矩阵HL赋值
ZL=[z（3）;z（4）;z（5）;z（6）;z（7）;z（8）;z（9）;z（10）;z（11）;z（12）;z（13）;z（14）;z（15）;z（16）]%给样本矩阵zL赋值
%CalculatingParameters
c1=HL'*HL;c2=inv（c1）;c3=HL'*ZL;c=c2*c3%计算并显示
%DisplayParameters
a1=c
（1）,a2=c
（2）,b1=c（3）,b2=c（4）%从
中分离出并显示a1、a2、b1、b2
%End
程序运行结果：
>>
u=[-1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1]
z=[0，0，0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949，-4.6352，6.2165，-5.5800，-2.5185]
HL=
001.0000-1.0000
-0.50000-1.00001.0000
-0.2500-0.50001.0000-1.0000
-0.5250-0.25001.00001.0000
-2.1125-0.52501.00001.0000
-4.3012-2.11251.00001.0000
-6.4731-4.3012-1.00001.0000
-6.1988-6.4731-1.0000-1.0000
-3.2670-6.1988-1.0000-1.0000
0.9386-3.26701.0000-1.0000
3.19490.9386-1.00001.0000
4.63523.1949-1.0000-1.0000
6.21654.63521.0000-1.0000
5.58006.21651.00001.0000
ZL=[0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949，-4.6352，-6.2165，-5.5800，-2.5185]T
c=[-1.5000，0.7000，1.0000，0.5000]T
a1=-1.5000
a2=0.7000
b1=1.0000
b2=0.5000
>>
从仿真结果表3.1可以看出，由于所用的输出观测值没有任何噪声成分，所以辨识结果也无任何误差。
例3.5考虑图3.6所示的仿真对象，图中，
是服从N
分布的不相关随机噪声。
且
，
，
选择图3.6所示的辨识模型。
仿真对象选择如下的模型结构
（3.68）
其中，
是服从正态分布的白噪声N
。
输入信号采用4位移位寄存器产生的M序
列，幅度为0.03。
按式
（3.69）
构造h（k）；加权阵取单位阵
；利用式（3.61）计算K（k）、
和P（k），计算各次参数辨识的相对误差，精度满足要求式（3.67）后停机。
最小二乘递推算法辨识的Malab6.0程序流程如图3.7所示。
下面给出具体程序。
%最小二乘递推算法辨识程序,在光盘中的文件名：
FLch3RLSeg3.m
%FLch3RLSeg3
clear%清理工作间变量
L=15;%M序列的周期
y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值
fori=1:
L;%开始循环，长度为L
x1=xor（y3,y4）;%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出
y（i）=y4;%取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列
ify（i）>0.5,u（i）=-0.03;%如果M序列的值为"1",辨识的输入信号取“-0.03”
elseu（i）=0.03;%如果M序列的值为"0",辨识的输入信号取“0.03”
end%小循环结束
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备
end%大循环结束，产生输入信号u
figure
（1）;%第一个图形
stem（u）,gridon%显示出输入信号径线图并给图形加上网格
z
（2）=0;z
（1）=0;%设z的前两个初始值为零
fork=3:
15;%循环变量从3到15
z（k）=1.5*z（k-1）-0.7*z（k-2）+u（k-1）+0.5*u（k-2）;%输出采样信号
end
%RLS递推最小二乘辨识
c0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye（4,4）;%直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.000000005;%取相对误差E=0.000000005
c=[c0,zeros（4,14）];%被辨识参数矩阵的初始值及大小
e=zeros（4,15）;%相对误差的初始值及大小
fork=3:
15;%开始求K
h1=[-z（k-1）,-z（k-2）,u（k-1）,u（k-2）]';x=h1'*p0*h1+1;x1=inv（x）;%开始求K（k）
k1=p0*h1*x1;%求出K的值
d1=z（k）-h1'*c0;c1=c0+k1*d1;%求被辨识参数c
e1=c1-c0;%求参数当前值与上一次的值的差值
e2=e1./c0;%求参数的相对变化
e（:
k）=e2;%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1;%新获得的参数作为下一次递推的旧参数
c（:
k）=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p（k）的值
p0=p1;%给下次用
ife2<=Ebreak;%如果参数收敛情况满足要求，终止计算
end%小循环结束
end%大循环结束
c，e%显示被辨识参数及其误差（收敛）情况
%分离参数
a1=c（1,:
）;a2=c（2,:
）;b1=c（3,:
）;b2=c（4,:
）;ea1=e（1,:
）;ea2=e（2,:
）;eb1=e（3,:
）;eb2=e（4,:
）;
figure
（2）;%第二个图形
i=1:
15;%横坐标从1到15
plot（i,a1,'r',i,a2,':
',i,b1,'g',i,b2,':
'）%画出a1，a2，b1，b2的各次辨识结果
title（'ParameterIdentificationwithRecursiveLeastSquaresMethod'）%图形标题
figure（3）;%第三个图形
i=1:
15;%横坐标从1到15
plot（i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
'）%画出a1，a2，b1，b2的各次辨识结果的收敛情况
title（'IdentificationPrecision'）%图形标题
程序运行结果：
>>
c=
0.001000.0010-0.4984-1.2328-1.4951-1.4962-1.4991-1.4998-1.4999
0.001000.00100.0010-0.23500.69130.69410.69900.69980.6999
0.001000.25091.24971.06651.00171.00201.00020.99990.9998
0.00100-0.24890.75000.56680.50200.50160.50080.50020.5002
图3.8最小二乘递推算法的参数辨识仿真
-1.5000-1.5000-1.5000-1.4999-1.4999
0.70000.70000.70000.7000-0.7000
0.99990.99990.99990.99990.9999
0.50000.50000.50000.50000.5000
e=
000-499.42001.47340.21280.00070.00200.00040.0000
0000-235.9916-3.94160.00420.00700.00120.0001
00249.86123.9816-0.1466-0.06070.0003-0.0018-0.0003-0.0001
00-249.8612-4.0136-0.2443-0.1143-0.0007-0.0016-0.0012-0.0001
0.00010.0000-0.0000-0.00000.0000
0.0001-0.00000.00000.00000.0000
0.00010.00000.00000.0000-0.0000
-0.00040.0000-0.00000.0000-0.0000
表3.2最小二乘递推算法的辨识结果
参数a1a2b1b2
真值-1.50.71.00.5
估计值-1.49990.70.99990.5000
仿真结果表明，大约递推到第十步时，参数辨识的结果基本达到稳定状态，即a1=-1.4999,a2=0.7000,b1=0.9999,b2=0.5000。
此时，参数的相对变化量E
0.000000005。
从整个辨识过程来看，精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小，但开始阶段的相对误差会很大，从图3.8（3）图形中可看出，参数的最大相对误差会达到三为数
第4章
第4章例：
考虑图4.2表示的线性时不变SISO系统。
根据卷积定理，系统的输出y（k）与输入序列
的关系可表示成
（4.40）
其中，
组成系统的脉冲响应。
系统在伪随机码输入作用下的输出响应如表4.1所示。
利用这些数据辨识系统的脉冲响应，当N取3时，有
（4.41）
令
（4.42）
若系统的真实脉冲响应记作g0,则有
（4.43）
根据式（4.16），脉冲响应估计值
的递推算式为
（4.44）
其中，脉冲响应估计值的初始值取
，权矩阵
取如下形式
（4.45a）
（4.45b）
确定性问题梯度校正递推算法辨识的流程如图4.3所示。
下面给出Malab6.0程序。
%梯度校正参数辨识程序，在光盘中的文件名：
FLch4GAeg1.m
clear
u=[-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1];
y=[0,-2,-6,-7,-7,-3,5,7,3,-1,5,3,-5,-3,1,-1,1,-5,-7,-7];
%画出u和y图形
figure
（1）,subplot（2,1,1）,stem（u）,subplot（2,1,2）,stem（y）,holdon
k=1:
20；plot（k,y）
%给出初始值
h1=[-1,0,0]';h2=[-1,-1,0]';g=[0,0,0]';I=[1,0,0;0,1/2,0;0,0,1/4];
h=[h1,h2,zeros（3,16）];
%计算样本数据h（k）
fork=3:
18
h（:
k）=[u（k）,u（k-1）,u（k-2）]';
end
%计算权矩阵R（k）和g的估计值
fork=1:
18
a=h（1,k）^2+（h（2,k）^2）/2+（h（3,k）^2）/4;a1=1/a;R=a1*I;%按照式（4.45a）和（4.45b）计算权矩阵
g（:
k+1）=g（:
k）+R*h（:
k）*（y（k+1）-h（:
k）'*g（:
k））;%按照式（4.44）计算脉冲响应的估计值
end
%绘图
g1=g（1,:
）;g2=g（2,:
）;g3=g（3,:
）;
figure
（2）;k=1:
19;subplot（121）;plot（k,g1,'r',k,g2,'g',k,g3,'b'）,gridon
%计算模型输出ym及系统输出与模型输出之间的误差Ey
fork=1:
18
ym（k）=h（:
k）'*g（:
k）;Ey（k）=y（k+1）-ym（k）;
end
k=1:
18;subplot（122）;plot（k,Ey）,gridon
g,ym,Ey%显示脉冲响应估计值、模型输出及系统输出与模型输出之间的误差
figure（3）;x=0:
1:
3;y=[0,g（1,18）,g（2,18）,g（3,18）];xi=linspace（0,3）;yi=interp1（x,y,xi,'cubic'）;
plot（x,y,'o',xi,yi,'m'）,gridon%画出脉冲响应估计值及其三次插值曲线
程序执行结果：
g=
02.00004.66675.23815.23811.53742.14382.23931.96581.9559
001.33331.61901.61903.46943.77263.82043.95713.9621
0000.14290.14291.06800.91640.94031.00871.0062
2.00631.99181.99801.99531.99951.99951.99952.00012.0032
3.98733.99463.99773.99903.99693.99693.99693.99723.9988
0.99360.99720.99570.99630.99740.99740.99740.99720.9980
ym=[0,-2.0000,-6.0000,-7.0000,3.4762,3.9388,6.8328,2.5213,-0.9826,4.9118,2.9746,-4.9891,-2.9953,1.0073,-1.0000,1.0000,-4.9990,-6.9945]
Ey=
[-2.0000,-4.0000,-1.0000,0.0000,-6.4762,1.0612,0.1672,0.4787,-0.0174,0.0882,0.0254,0.0109,-0.0047,-0.0073,
-0.0000,0,-0.0010,-0.0055]
图4.4确定性问题梯度校正参数辨识仿真（被辨识参数的个数为3）
图4.5确定性问题梯度校正参数辨识仿真（被辨识参数的个数为5）
图4.4表明，递推到第10步时，被辨识参数基本上达到了稳定状态，即
;此时系统的输出与模型的输出误差也基本达到稳定状态，即
。
由于被辨识参数的个数较少，递推校正算法的收敛性比较好，也就是说，输入输出的观测数据量已足够。
但从图4.4Figure3中可看出，由于被辨识参数的个数较少，它还不足以充分显示全部的系统脉冲响应。
为了充分地显示出系统的脉冲响应，可以把被辨识参数的个数增加到5个。
图4.5给出了被辨识参数的个数为5时的辨识结果。
图4.5基本显示出了系统的全部脉冲响应过程。
但在图4.5的辨识中，还是利用上面给出的20对输入输出数据。
通过图4.4与图4.5的比较可以看出，在观测数据量一定的情况下，随着被辨识参数的个数的增加，梯度校正辨识算法的收敛性变差，这是由于观测数据不足，递推步数有限造成的。
所以，随着被辨识参数的个数的增加，若要保持梯度校正辨识算法的收敛效果，需要同步增加观测数据量。
第七章的用改进的神经网络MBP算法辨识
例7.1对具有随机噪声的二阶系统的模型辨识
对具有随机噪声的二阶系统的模型辨识，进行标幺化以后系统的参考模型差分方程为
（7.90）
式中，
为随机噪声。
由于神经网络的输出最大为1，所以，被辨识的系统应先标幺化，这里标幺化系数为5。
利用图7.5正向建模（并联辨识）结构，神经网络选用3-9-9-1型，即输入层i，隐层j包括2级，输出层k的节点个数分别为3、9、9、1个；
由于神经网络的最大输出为1，因此在辨识前应对原系统参考模型标么化处理，辨识结束后再乘以标么化系数才是被辨识系统的辨识结果。
1）编程如下:
%w10ij表示第一隐层权值
，w11ij表示
；w120ij表示第二隐层权值%
，w121ij表示
；w20j表示输出层权值
，w21j表示%
；q表示隐层阈值；p表示输出层阈值；置标幺化系数f1=5等
w10ij=[.01.01.02;.1.11.02;.010.1;.11.01.02;.1.1.02;.11.1.1;.1.1.1;0.1.1;.10.1];
w11ij=[.1.2.11;.02.13.04;.09.08.08;.09.1.06;.1.11.02;.060.1;.1.1.1;0.10;.1.1.1];
w20j=[.01;.02;.1;.2;.1;.1;.1;.1;.1];
w21j=[0;0.1;.1;.02;0;.1;.1;.1;.1];
q0j=[.9.8.7.6.1.2.1.1.1];
q120j=q0j;
q11j=[.5.2.3.4.1.2.1.1.1];
q12j=q11j;
w121ij=w20j*q0j;
w120ij=w20j*q11j;
f1=5;
q2j=0;%thresholdvalue
p0=.2;
k1=1;
p1=.3;
w=0;
xj=[111];%inputs
error=0.0001;
a1=[1111];
n=1;
e1=0;
e0=0;
e2=0;
e3=0;
e4=0;
yo=0;
ya=0;
yb=0;
y0=0;
y1=0;
y2=0;
y3=0;
u=0;
u1=0;
u2=0.68;
u3=.780;
u4=u3-u2;
k1=1;
kn=28;
e3=.055;z1=0;
z12=0;q123j=0;t2j=0;o12j=0;r=0;r1=0;s=0.1;d2j=0;
%+++++++++++++++++++++++++++++++++%calculatingoutputofthehiddenlayer
v1=randn（40,1）;
form=1:
40
s1=0.1*v1（m）
yn=.3366*y2+.6634*u1+s*s1;
y1=y2;
y2=yn;
yp=yn;
u0=u1;
u1=u2;
ya（m）=yn;
fork=1:
100
%calculatingoutputofthehiddenlayer
（1）
fori=1:
9
x1=[w11ij（i,1）*xj（:
1）]+[w11ij（i,2）*xj（:
2）]+[w11ij（i,3）*xj（:
3）];
x=x1+q11j（:
i）;
o=1/[1+exp（-x）];
o11j（i）=o;
end
%calculatingoutputofthehiddenlayer
（2）
fori=1:
9
forj=1:
9
z1=z1+w121ij（i,j）*o11j（:
j）;
end
z=z1+q12j（:
i）;
o=
升级会员 