非线性模型辨识

非线性模型辨识
张慧勇 S082121 Identification of a nonlinear system represented by the following simple Hammerstein model.
G(s)
1 1 2s v(t ) 2 u (t ) 0.5u 2 (t )
五、模型验证
在 matlab 的 simulink 中构造给定模型和假定模型（如图 6 ） ，对给定模型与假 定模型的输出进行比较，来测试假定模型是否合适。
图 6 模型验证的 simulink 图
图 6 中 零阶保持器（Zero-Order Hold）中的 Sample time 设置为 0.4；离散传 递函数中的 Sample time 设置为 0.2；To Workspace 中的 Sample time 设置为 0.2；PRTS 测试信号要从 workspace 中导入。
p 为输入u为参数的函数。
二、 稳态分析
在不同的阶跃输入下，测设系统的稳态输出，记录如下表；
阶跃输入u -20 -15 -10 -5 -3 -2 -1 0 稳态输出p 182.00 99.50 42.00 9.50 3.50 2.00 1.50 2.00 表 1 不同阶跃输入对应的稳态输出 阶跃输入u 20 15 10 5 3 2 1 稳态输出p 222.00 129.50 62.00 19.50 9.50 6.00 3.50
ˆ (k ) 表示假定模型输出的采样数 用 y (k ) 表示系统给定模型输出的采样数据 , y ˆ (k ) 绘制如图 7 据， y (k ) 与 y
图 7 系统模型与假定模型的输出曲线
ˆ (k ) 表示系统模型与假定模型的输出误差， e(k ) 绘制如图８ e(k ) y(k ) y
图 8 系统模型与假定模型的输出误差曲线
三、动态分析
y(t ) p(1 e qt ) 可以视为阶跃输入为 u(t ) p 的输出函数，u(t ) p 的拉普拉斯变
换 为 U ( s)
G(s)
p p ， 输 出 y(t) 的 拉 普 拉 斯 变 换 为 Y (s) ， s s( s / q 1)
Y ( s) 1 ，该模型的动态部分是一阶系统，因此假设系统动态部分 U ( s) s / q 1
由最小二乘法可得 C ( AT A)1 AT Y ，求的 c1 0.9048, c2 0.0952 离散模型为 y(k ) 0.9048 y(k 1) 0.0952 p(k 1) ， p(k ) 0.5u 2 (k ) u(k ) 2
模型动态部分的 Z 变换为
四、模型测试信号
前文假设的模型的多项式稳态只有 2 阶，因此构造 PRTS 作为模型的测试信 号，来验证假定模型。 设 PRTS 的时钟周期为 Te 0.4s ，信号的周期数 Np 3T95 / Te ； 取 Np 80 ， 80 r nr 1 3nr 1 nr 4 ，需要 4 个移位寄存器； PRTS 信号的时间周期长 Np Te 80*0.4s 32s ； PRTS 信号的差分方程为：
y(t)
两模型输出比较的 Matlab 程序
t=simout(:,1);%从 workspace 中读入时间 u=simout(:,2);% 从 workspace 中读入输入
y1=simout(:,3);% 从 workspace 中读入系统模型输出 y2=simout(:,4);% 从 workspace 中读入假定模型输出 figure(1) plot(t,y1,'r'); hold on plot(t,y2,'*') grid on xlabel('时间') ylabel('输出') legend('系统模型输出','假定模型输出') figure(2) e=y1-y2; plot(t,e,'r') grid on xlabel('时间') ylabel('误差') title('误差曲线')
Gz=zpk([],0.9048,0.0952,0.2) Gs=d2c(Gz)
模型动态部分为： G( s) 系统连续模型为
0.50021 1 s 0.5002 1.9992s 1.0000
u(t)
0.5u (t ) u(t ) 2
2
p(t)
1 1.9992s 1.0000
一、 数据采样分析
在matlab的simulink中构造给定模型(如图1)， 并利用给定模型产生数据 （如图2） , 将采样数据导入Workspace,通过分析该数据来辨识系统模型。
图1 给定Hammerstein model 的simulink图
图2 simulink中数据采样图
采样时间取为0.2s，将图2中To Workspace中的sample time设置为0.2即可。 输入 u0 (t ) 0 时，输出的采样数据记为 y0 (k ) ； 输入 u1 (t ) 1 时，输出的采样数据记为 y1 (k ) ； 输入 u2 (t ) 2u1 (t ) 2 时，输出的采样数据记为 y2 (k ) 。
的离散模型为 y(k ) c1 y(k 1) c2 p(k 1) ，其中 p(k ) 0.5u 2 (k ) u(k ) 2
令 C [c1
c2 ]T
p(1) M p(k 1) M p ( N 1)
y (2) y (1) M M 令 Y y ( k ) ， A y (k 1) M M y( N ) y ( N 1)
where G( s) is the transfer function of the linear dynamic part, v(t ) is the nonlinear static part. 1) Design the experiment: test signal, sample time etc. 2) Use the given model to generate data. 3) Estimate the parameters. 4) Validate the model. 5) Write a report on the identification procedure, results and analysis including the programs.
x(k ) c1 3 x(k 1) 3 c2 3 x(k 2) 3 c3 3 x(k 3) 3 c4 3 x(k 4) ， k 4
其中 c1 0, c2 0, c3 1, c4 1 初值 x(0) 0, x(1) 1, x(2) 2, x(3) 1
根据上表绘出不同阶跃输入对应的稳态输出图，如图 4 左，由图可以看出像一抛 物线，因此进行二次曲线拟合。 设 p c1u 2 c2u c3 ，用最小二乘法求系数 c1 c2 令 C [c1 c2
c3
c3 ]T ， A [u 2
u 1]
1.0000
（u,p 都为列向量） 2.0000]，拟合曲线如图 4 右；由图可以
u=[-20,-15,-10,-5,-3,-2,-1,0,1,2,3,5,10,15,20]'; p=[182.00 99.50 42.00 9.50 3.50 2.00 1.50 2.00 3.50 6.00 9.50 19.50 62.00 129.50 222.00]'; subplot(1,2,1) plot(u,p,'*') grid on xlabel('阶跃输入 ') ylabel('稳态输出 ') A=[u.^2,u,ones(15,1)]; c=inv(A'*A)*A'*p %c=polyfit(u,p,2) u1=-22:1:22; p1=polyval(c,u1); subplot(1,2,2) plot(u,p,'*') hold on plot(u1,p1,'r') grid on axis([-22 22 -10 250]) xlabel('阶跃输入 ') ylabel('稳态输出 ')
(k )
y2 (k ) y0 (k ) 2.667 2 y1 (k ) y0 (k )
因此该系统模型为非线性模型。 采样数据 y0 (k ) ， y1 (k ) ， y2 (k ) 绘制如图3
ቤተ መጻሕፍቲ ባይዱ
图3、阶跃输入的系统输出图
由图3可以看出输出曲线也可用函数 y p(1 e qt ) 来描述，其中 q 0 描述输 出的动态部分，当时间 t增大时，输出y将趋近于 p ， 因此 p 描述输出的稳态部分，
图 5 周期数为 80 的 PRTS 信号
PRTS 信号产生的 Matlab 程序
%PRTS N=160; Ts=0.4;%时钟周期 x=zeros(1,N); c1=0;c2=0;c3=1;c4=1; x(1)=0;x(2)=1;x(3)=2;x(4)=1; for k=5:N x(k)=c1*x(k-1); x(k)=mod(x(k),3)+mod(c2*x(k-2),3); x(k)=mod(x(k),3)+mod(c3*x(k-3),3); x(k)=mod(x(k),3)+mod(c4*x(k-4),3); x(k)=mod(x(k),3); end time=0:Ts:(N-1)*Ts; signal=[time',x'];
G( z )
c z 1 y( z ) 0.0952 z 1 0.0952 ，采样时间为 0.2s 2 1 1 p( z ) 1 c1 z 1 0.9048 z z 0.9048