系统辨识举例.docx
- 文档编号:30675013
- 上传时间:2023-08-19
- 格式:DOCX
- 页数:29
- 大小:364.88KB
系统辨识举例.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=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 系统 辨识 举例