四、神经网络系统辨识分析（25分）
用BP 神经网络进行系统在线逼近的原理框图如图3所示
u (k y n (k )
y (k
图3 图4
假设某控制对象的模型为y (k ) =u (k ) +
3
y (k -1)
，采样时间取t=1ms，输入信号
1+y (k -1) 2
u(k ) =0. 5sin(6πt) 。采用的BP 神经网络结构如图4所示，权值w ij 和w j 2的初值取 [-1，+1]
之间的随机值，权值采用δ学习算法，学习速率η取0.50，动量因子α取0.05。试分析神经网络在线逼近的运行过程，并作Matlab 仿真。
题目四、需要阐述清楚BP 网络逼近控制对象的工作原理和学习过程
BP 算法的基本思想是：对于一个输入样本, 经过权值、阈值和激励函数运算后，得到一个输出y n (k)，然后让它与期望的样本y(k)进行比较，若有偏差，则从输出开始反向传播该偏差，进行权值、阈值调整，使网络输出逐渐与希望输出一致。
BP 算法由四个过程组成：输入模式由输入层经过中间层向输出层的“模式顺传播”过程，网络的希望输出与网络的实际输出之间的误差信号由输出层经过中间层向输入层逐层修正连接权的“误差逆传播”过程，由“模式顺传播”与“误差逆传播”的反复交替进行的网络“记忆训练”过程，网络趋向于收敛即网络的全局误差趋向极小值的 “学习收敛”过程。
BP 网络（Back Propagation），该网络是一种单向传播的多层前向网络。误差
反向传播的BP 算法简称BP 算法，其基本思想是梯度下降法。它采用梯度搜索技术，以期使网络的实际输出值与期望输出值的误差均方值为最小。
BP 网络特点:
（1）是一种多层网络，包括输入层、隐含层和输出层； （2）层与层之间采用全互连方式，同一层神经元之间不连接； （3）权值通过δ学习算法进行调节；
（4）神经元激发函数为S 函数；
（5）学习算法由正向传播和反向传播组成； （6）层与层的连接是单向的，信息的传播是双向的。 BP 网络结构:
含一个隐含层的BP 网络结构如图4-1所示，图中i 为输入层神经元，j 为隐层神经元，k 为输出层神经元。
u (k )
y n (k )
y (k )
图4-1 BP神经网络结构
BP 网络的逼近:
BP 网络逼近的结构如图4-2所示，图中变量k 为网络的迭代步骤。BP 为网络逼近器，y(k)为被控对象实际输出，y n (k)为BP 的输出。将系统输出y(k)及输入u(k)的值作为逼近器BP 的输入，将系统输出与网络输出的误差作为逼近器的调整信号。
图4-2 BP神经网络逼近
BP 算法的学习过程由正向传播和反向传播组成。在正向传播过程中，输入信息从输入层经隐层逐层处理，并传向输出层，每层神经元（节点）的状态只影响下一层神经元的状态。如果在输出层不能得到期望的输出，则转至反向传播，将误差信号（理想输出与实际输出之差）按联接通路反向计算，由梯度下降法调整各层神经元的权值，使误差信号减小。
（1）前向传播：计算网络的输出。
x =w ij x i
隐层神经元的输入为所有输入的加权之和： j
i
（都是向量形式 比如i=1,2,3 j=1,2,3 x 3=w13x 1+w23x 2+w33x 3） 1'
隐层神经元的输出采用S 函数激发： x j =f (x j ) =-x j
1+e
∂x ' j e -xj ' '
=x j (1-x j ) = 则 -xj 2∂x (1+e ) j
'
x =w x 输出层神经元的输出： k j 2j
j
网络输出与理想输出误差为： e (k ) =y (k ) -y (k ) n
1
误差性能指标函数为： E =e (k ) 2
2
∑
∑
（2）反向传播：采用δ学习算法，调整各层间的权值。
根据梯度下降法，权值的学习算法如下： 输出层及隐层的连接权值学习算法为：
∂x ∂E ∆w j 2=-η=-η⋅e (k ) ⋅=-η⋅e (k ) ⋅x ' j
∂w j 2∂w j 2
k+1时刻网络的权值为： w j 2(t +1) =w j 2(t ) +∆w j 2隐层及输入层连接权值学习算法为：
∂y ∂E
∆w ij =-η=-η⋅e (k ) ⋅
∂w ij ∂w ij
其中 ' '
∂x ∂x ∂x ∂y ∂y j j j =⋅ ⋅=w j 2⋅⋅x i =w j 2⋅x ' j (1-x ' j ) ⋅x i
∂w ij ∂x ' j ∂x j ∂w ij ∂x j
k+1时刻网络的权值为：
w ij (k +1) =
w ij (k ) +∆w ij
α，此时的权值如果考虑上次权值对本次权值变化的影响，需要加入动量因子
为：
w ij (k +1) =w ij (k ) +∆w ij +α(w ij (k ) -w ij (k
η为学习速率， 其中, α为动量因子。
∈0, 1∈0, 1
仿真图如示：
w j 2(k +1) =w j 2(k ) +∆w j 2+α(w j 2(k ) -w j 2(k -1))
-1))
η
[]
α
[]
程序如示：
%BP identification clear all; close all;
xite=0.50; alfa=0.05;
w2=rands(3,1); w2_1=w2;w2_2=w2_1;
w1=rands(2,3); w1_1=w1;w1_2=w1;
dw1=0*w1;
x=[0,0]';
u_1=0; y_1=0;
I=[0,0,0]'; Iout=[0,0,0]'; FI=[0,0,0]';
ts=0.001;
for k=1:1:1000
time(k)=k*ts;
u(k)=0.50*sin(3*2*pi*k*ts); y(k)=u_1^3+y_1/(1+y_1^2);
for j=1:1:3
I(j)=x'*w1(:,j);
Iout(j)=1/(1+exp(-I(j))); end
yn(k)=w2'*Iout; % Output of NNI networks
e(k)=y(k)-yn(k); % Error calculation
w2=w2_1+(xite*e(k))*Iout+alfa*(w2_1-w2_2);
for j=1:1:3
FI(j)=exp(-I(j))/(1+exp(-I(j)))^2; end
for i=1:1:2 for j=1:1:3
dw1(i,j)=e(k)*xite*FI(j)*w2(j)*x(i); end end
w1=w1_1+dw1+alfa*(w1_1-w1_2);
%%%%%%%%%%%%%%Jacobian%%%%%%%%%%%%%%%% yu=0;
for j=1:1:3
yu=yu+w2(j)*w1(1,j)*FI(j); end
dyu(k)=yu;
x(1)=u(k); x(2)=y(k);
w1_2=w1_1;w1_1=w1; w2_2=w2_1;w2_1=w2; u_1=u(k); y_1=y(k); end
figure(1);
plot(time,y,'r',time,yn,'b');
xlabel('times');ylabel('y and yn'); figure(2);
plot(time,y-yn,'r');
xlabel('times');ylabel('error');