浅谈BP神经网络PID控制算法及matlab仿真

本文是对BP神经网络PID控制算法的数学描述及仿真实验，若有错误之处，欢迎指正！

BP神经网络

简述

老规矩不废话，直接上链接 BP神经网络维基百科
BP神经网络是人工神经网络中的一种常用结构，其由输入层(input)-隐含层(hidding)-输出层三层构成(output)。

上图中，$B1$是输入层，$B2$是隐含层，$B3$是输出层，$W12$是输入-隐含层之间的权重系数矩阵，$W23$是隐含-输出层的权重系数矩阵，$H$是隐含层神经元函数的激励函数，$b2$是隐含层神经元的偏置阈值，$O$是输出层神经元的激励函数，$b3$是输出层神经元的偏置阈值。
假设输入向量为$X=[1\times 3]$，输入-隐含层系数矩阵为$W_{12}=[3\times N]$，则隐含层的中每个神经元的输出为：
$$Hiddin{g}_{O}utput=H((X\cdot W_{12}{)}_i-b_2),i=1,2,3,...N$$
对于偏置阈值来说，同样也会影响到损失函数$J$的大小，因此也要将其和权重系数一样作为可变参数，应用梯度下降求解最优值。
利用矩阵运算的技巧，假设输入层的输入为1x3的向量$X=[x_1,x_2,x_3]$，对应的权重系数为$W=[3\times N]$，则为了方便计算，可以将输入向量变为$X^{\prime }=[x_1,x_2,x_3,1]$，相应的权重系数会多出一行变为$W^{\prime }=[4\times N]$，最后这一行就对应每一个隐含层神经元的偏置阈值。

流程

BP神经网络的运作机理主要包含两个流程

1. 前向激励传播：将输入送入网络进行运算获得结果。
2. 反向传播更新权重：将运算结果得到的反馈信息计算损失函数，以最优化目标更新各神经元连接权重值。

BP神经网络PID控制算法(BPNN-PID)

与单神经元PID控制算法的对比

从输入输出关系上来看，单神经元PID控制算法输入三项误差项，输出控制律。而BP神经网络输入的是系统输入$r(k)$，系统输出$y(k)$，系统误差$e(k)$，输出的是PID参数。
可以理解成BP神经网络根据系统输入输出和误差项的信息和损失函数的反馈来不断迭代优化各位置的权重系数，以实现损失函数达到最小的目的。

为什么输入的三项分别是$r(k),y(k),e(k)$而不是其他的变量。

• 首先排除妈妈生的（狗头）
• 个人认为，在进行权值更新的损失函数就是以系统误差作为基础进行计算的——$J=\frac{1}{2}(r(k)-y(k){)}^2$，因此对于整个BP神经网络部分来说，输出的PID系数值的优劣(以损失函数大小评判)是直接和输入$r(k)y(k)$以及二者的差$e(k)$决定的。因此才将这三个作为BP神经网络的输入变量。
• 从神经网络的角度来说，输入的数据是服务于损失函数的训练样本，样本中就需要包含训练结果($y$)及参考答案($r$)以及他们的差距。当然输入的信息不是一成不变的，根据不同的代价函数和系统结构，可以提供更多数据类型来达到最终的收敛目的。

前向激励

在BPNN-PID控制算法中，前向激励是指从输入层信号输入通过隐含层最终从输出层输出的通道，也是算法中的主要流程之一。

• 假设隐含层中的神经元个数为$N$，隐含层和输出层的激励函数分别为$f(x)，g(x)$：
  $$\begin{gathered} f(x)=\frac{e^x-e^{-x}}{e^x+e^{-x}} \\ g(x)=\frac{e^x}{e^x+e^{-x}} \end{gathered}$$
• 假设从输入层-隐含层的权重系数为$w_{ij}(i=1,2,3;j=1,2,3,...,N)$，隐含层-输出层的权重系数为$v_{ij}(i=1,2,3,...,N;j=1,2,3)$。对应的系数矩阵分别为
  W = [ w 11 w 12 . . . w 1 N w 21 w 22 . . . w 2 N w 31 w 32 . . . w 3 N ] V = [ v 11 v 12 v 13 v 21 v 22 v 23 . . . . . . . . . v N 1 v N 2 v N 3 ] W = \begin{bmatrix} w_{11}& w_{12}& ...&w_{1N} \\ w_{21}& w_{22}& ...&w_{2N}\\ w_{31}& w_{32}& ...&w_{3N} \end{bmatrix}\\ V = \begin{bmatrix} v_{11}& v_{12}& v_{13} \\ v_{21}& v_{22}& v_{23}\\ ...& ...& ...&\\ v_{N1}& v_{N2}& v_{N3} \end{bmatrix}\\ W= ​w11​w21​w31​​w12​w22​w32​​.........​w1N​w2N​w3N​​ ​V= ​v11​v21​...vN1​​v12​v22​...vN2​​v13​v23​...vN3​​​ ​

假设三个输入分别为$x_1，x_2，x_3$，则隐含层第一个神经元的输入为
$$h{i}_1=x_1\ast w_{11}+x_2\ast w_{21}+x_3\ast w_{31}$$

假设隐含层N个输出为$h{o}_1,h{o}_2,h{o}_3,...,h{o}_N$，则输出层第一个神经元的输入为

$$o{i}_1=h{o}_1\ast v_{11}+h{o}_2\ast v_{21}+...+h{o}_N\ast v_{N1}$$
此处内容用于辅助理解权重系数下标数字的含义，为后续公式推导提供帮助

• 为了简化公式，采用矩阵形式来描述前向激励过程。输入向量为$X=[x_1,x_2,x_3]$1x3的行向量，则隐含层的输入向量为
  $$Hi=X\cdot W$$
  显然$Hi$为1xN的行向量。
  隐含层的激励函数为$f(x)=\frac{e^x-e^{-x}}{e^x+e^{-x}}$，则隐含层的激励输出为
  $$Ho=f(Hi)=\frac{e^{Hi}-e^{-Hi}}{e^{Hi}+e^{-Hi}}$$
  显然，对向量$Hi$的每一个元素进行函数映射转换后，不改变其向量的维度，$Ho$仍然为1xN的行向量。

  输出层的神经元输入向量为$Oi$1x3的行向量为
  $$Oi=Ho\cdot V$$
  输出层的神经元激励输出向量为$Oo$1x3的行向量为
  $$Oo=g(Oi)=g(x)=\frac{e^{Oi}}{e^{Oi}+e^{-Oi}}$$

  在matlab中，对一个向量或者矩阵进行函数运算会得到维度相同的新向量或新矩阵，其元素是原先向量或矩阵中每个元素经过函数之后的结果。

• $BPNN-PID$控制算法输出的是$Kp,Ki,Kd$三个参数，即$Oo=[Kp,Ki,Kd]$。
  设误差向量$E$3x1列向量为
  E ⃗ = [ e r r o r k − e r r o r k − 1 e r r o r k e r r o r k − 2 e r r o r k − 1 + e r r o r k − 2 ] \vec{E}=\begin{bmatrix} error_{k}-error_{k-1} \\ error_{k}\\ error_{k}-2error_{k-1}+error_{k-2} \end{bmatrix} E = ​errork​−errork−1​errork​errork​−2errork−1​+errork−2​​ ​
  根据增量PID公式计算控制器增量$\Delta U$
  $$\Delta U=Oo\cdot E$$
  则控制律公式为
  $$U_k=U_{k-1}+\Delta U$$
  假设系统函数满足以下关系
  $$\begin{gathered} {\mathbf{x}}^{{}^{\prime }}=\mathbf{h}(\mathbf{x},\mathbf{U}) \\ \mathbf{y}=\mathbf{x} \end{gathered}$$
  其中$\mathbf{x}$为状态量，$\mathbf{y}$为系统输出，$\mathbf{h}$可以为任意形式的函数形式，包括各类非线性函数。

反向传播

反向传播的目的是根据系统响应结果，来判断当前系统参数的优化方向，使得系统响应按照期望的方向演变。

• 如何评价响应结果
  对于机器来说，它判断当前系统响应状态好坏的方式只能是通过数值，因此需要定义能够描述系统响应标准的数学表达式，以表达式的计算结果衡量响应结果。这种数学表达式就是损失函数cost function，也可以称作代价函数。常见的表达式为：
  $$\mathbf{J}=\frac{1}{2}({\mathbf{y}}_{\mathbf{d}}-\mathbf{y}{)}^2$$
  其中，${\mathbf{y}}_{\mathbf{d}},\mathbf{y}$分别系统期望和系统状态。显然，$\mathbf{J}$越小，系统的响应性能越好。

• 如何判断系统参数的优化方向
  从数值上很容易看出，我们优化的目的就是让损失函数尽可能的变小甚至达到并保持在0。也就是说，损失函数$\mathbf{J}$下降的方向就是优化的方向，在工程上一般用梯度下降法来实现。而能够影响到损失函数$\mathbf{J}$数值大小的系统参数只有各层之间的权重系数（若神经元中的偏置阈值不恒为0的话同样会作为系统参数影响损失函数数值大小，此文中不讨论）。

  以下为梯度下降法的基本公式基本公式（式1）和动量梯度下降法（式2）
  $$\begin{gathered} \Delta X_k=-\eta \frac{\partial \mathbf{J}}{\partial X} \\ \Delta X_k=-\eta \frac{\partial \mathbf{J}}{\partial X}+\alpha \Delta X_{k-1} \end{gathered}$$
  其中，$\eta$为学习速率或者叫下降步长，$\alpha$为惯性因子。
  在整个神经网络中，存在着总共$6N$个权重系数，为了找到最优的方向，则令$\mathbf{J}$对每个权重系数求偏导，根据梯度下降法得到变化值，最终完成整个网络参数的优化迭代。以输入层-隐含层权重系数 $W$为例，它是第一个输入连接隐含层第一个神经元的权重系数，则它的梯度下降过程为：
  $$\begin{gathered} \Delta W=-\eta \frac{\partial \mathbf{J}}{\partial W} \\ \frac{\partial \mathbf{J}}{\partial W}=\frac{\partial \mathbf{J}}{\partial \mathbf{y}}\frac{\partial \mathbf{y}}{\partial \mathbf{u}}\frac{\partial \mathbf{u}}{\partial \mathbf{Oo}}\frac{\partial \mathbf{Oo}}{\partial \mathbf{Oi}}\frac{\partial \mathbf{Oi}}{\partial \mathbf{Ho}}\frac{\partial \mathbf{Ho}}{\partial \mathbf{Hi}}\frac{\partial \mathbf{Hi}}{\partial \mathbf{W}} \end{gathered}$$

这里面设计到对向量求偏导，对矩阵求偏导，其本质是复合函数的求偏导，因为将矩阵运算拆开来就是线性方程组，其本质依然是链式求偏导法则。
例如，$u$ = [$u_1$,$u_2$,$u_3$]，$v$=[$v_1$,$v_2$,$v_3$]${}^T$。并且满足$K=u\cdot v$，则可以将$K$等效成：
$$K=u_1{v}_1+u_2{v}_2+u_3{v}_3=f(u_1,u_2,u_3,v_1,v_2,v_3)$$
则当我们令$K$对$u$求导时，得到：
$$\frac{\partial K}{\partial u}=\left[\begin{array}{c}\frac{\partial K}{\partial u_1}\frac{\partial K}{\partial u_2}\frac{\partial K}{\partial u_3}\end{array}\right]$$
其中，$\frac{\partial K}{\partial u_1},\frac{\partial K}{\partial u_2},\frac{\partial K}{\partial u_3}$分别为：
$$\begin{gathered} \frac{\partial K}{\partial u_1}=\frac{\partial f}{\partial u_1}=v_1 \\ \frac{\partial K}{\partial u_2}=\frac{\partial f}{\partial u_2}=v_2 \\ \frac{\partial K}{\partial u_3}=\frac{\partial f}{\partial u_3}=v_3 \end{gathered}$$
即$\frac{\partial K}{\partial u}=v^T$。同理可得$\frac{\partial K}{\partial v}=u^T$
当$K$变成一个向量或者是矩阵时，情况稍微复杂一些，但本质是一样的。
假设 l ⃗ = [ l 1 , l 2 , l 3 ] , M = [ m 11 m 12 m 21 m 22 m 31 m 32 ] \vec{l}=[l_1,l_2,l_3],\mathbf{M}=\begin{bmatrix} m_{11}& m_{12}\\ m_{21}& m_{22}\\ m_{31}&m_{32}\end{bmatrix} l =[l1​,l2​,l3​],M= ​m11​m21​m31​​m12​m22​m32​​ ​，令$n=l\cdot \mathbf{M}=[l_1{m}_{11}+l_2{m}_{21}+l_3{m}_{31},l_1{m}_{12}+l_1{m}_{22}+l_1{m}_{32}]$,可以等效成一个方程组：
$$\{\begin{array}{c}{n}_1=l_1{m}_{11}+l_2{m}_{21}+l_3{m}_{31}=f_1(l_1,l_2,l_3,m_{11},m_{21},m_{31})\\ n_2=l_1{m}_{12}+l_1{m}_{22}+l_1{m}_{32}=f_2(l_1,l_2,l_3,m_{12},m_{22},m_{32})\end{array}$$
则当令$n$对向量$l$求导时，得到：
$$\frac{\partial n}{\partial l}=[\frac{\partial n_1}{\partial l};\frac{\partial n_2}{\partial l}]=\left[\begin{array}{ccc}\frac{\partial n_1}{\partial l_1}& \frac{\partial n_1}{\partial l_2}& \frac{\partial n_1}{\partial l_3}\\ \frac{\partial n_2}{\partial l_1}& \frac{\partial n_2}{\partial l_2}& \frac{\partial n_2}{\partial l_3}\end{array}\right]=M^T$$
同理可得
$$\frac{\partial n}{\partial \mathbf{M}}=l^{{}^{\prime }}$$
但事实上这个梯度下降过程的本质是损失函数对每个权值参数的直接求偏导，再将每个结果放回权值矩阵中。

$$\begin{gathered} \frac{\partial \mathbf{J}}{\partial \mathbf{y}}=\frac{\partial \frac{1}{2}(r-y{)}^2}{\partial \mathbf{y}}=-\mathbf{error} \\ \frac{\partial \mathbf{y}}{\partial \mathbf{u}}=\mathbf{sign}(\frac{\mathbf{\Delta y}}{\mathbf{\Delta u}}) \\ \frac{\partial \mathbf{u}}{\partial \mathbf{Oo}}=\frac{\partial \mathbf{\Delta u}}{\partial \mathbf{Oo}}=\frac{\partial \mathbf{Oo}\cdot \mathbf{E}}{\partial \mathbf{Oo}}=\mathbf{E} \\ \frac{\partial \mathbf{Oo}}{\partial \mathbf{Oi}}=\frac{\partial \mathbf{g}(\mathbf{Oi})}{\partial \mathbf{Oi}}={\mathbf{g}}^{\prime }(\mathbf{Oi}) \\ \frac{\partial \mathbf{Oi}}{\partial \mathbf{Ho}}=\frac{\partial \mathbf{Ho}\mathbf{V}}{\partial \mathbf{Ho}}=\mathbf{V} \\ \frac{\partial \mathbf{Ho}}{\partial \mathbf{Hi}}=\frac{\partial \mathbf{f}(\mathbf{Hi})}{\partial \mathbf{Hi}}={\mathbf{f}}^{\prime }(\mathbf{Hi}) \\ \frac{\partial \mathbf{Hi}}{\partial \mathbf{W}}=\frac{\partial \mathbf{X}\cdot \mathbf{W}}{\partial \mathbf{W}}=\mathbf{X} \end{gathered}$$
整理可以得到对于输入-隐含层的权重系数的梯度下降公式为：
$$\Delta W=-\eta \frac{\partial \mathbf{J}}{\partial W}=\eta \cdot \mathbf{error}\cdot \mathbf{sign}(\frac{\mathbf{\Delta y}}{\mathbf{\Delta u}})\cdot \mathbf{E}\cdot {\mathbf{g}}^{\prime }(\mathbf{Oi})\cdot \mathbf{V}\cdot {\mathbf{f}}^{\prime }(\mathbf{Hi})\cdot \mathbf{X}$$

同理可得，损失函数对隐含-输出层权重系数的梯度下降过程为：
$$\begin{gathered} \Delta V=-\eta \frac{\partial \mathbf{J}}{\partial V} \\ \frac{\partial \mathbf{J}}{\partial V}=\frac{\partial \mathbf{J}}{\partial \mathbf{y}}\frac{\partial \mathbf{y}}{\partial \mathbf{u}}\frac{\partial \mathbf{u}}{\partial \mathbf{Oo}}\frac{\partial \mathbf{Oo}}{\partial \mathbf{Oi}}\frac{\partial \mathbf{Oi}}{\partial \mathbf{V}} \end{gathered}$$
每一项的公式为：
$$\begin{gathered} \frac{\partial \mathbf{J}}{\partial \mathbf{y}}=\frac{\partial \frac{1}{2}(r-y{)}^2}{\partial \mathbf{y}}=-\mathbf{error} \\ \frac{\partial \mathbf{y}}{\partial \mathbf{u}}=sign(\frac{\mathbf{\Delta y}}{\mathbf{\Delta u}}) \\ \frac{\partial \mathbf{u}}{\partial \mathbf{Oo}}=\frac{\partial \mathbf{\Delta u}}{\partial \mathbf{Oo}}=\frac{\partial \mathbf{Oo}\cdot \mathbf{E}}{\partial \mathbf{Oo}}=\mathbf{E} \\ \frac{\partial \mathbf{Oo}}{\partial \mathbf{Oi}}=\frac{\partial \mathbf{g}(\mathbf{Oi})}{\partial \mathbf{Oi}}={\mathbf{g}}^{\prime }(\mathbf{Oi}) \\ \frac{\partial \mathbf{Oi}}{\partial \mathbf{V}}=\frac{\partial \mathbf{Ho}\mathbf{V}}{\partial \mathbf{V}}=\mathbf{Ho} \end{gathered}$$
整理可以得到对于隐含-输出层的权重系数的梯度下降公式为：
$$\Delta V=-\eta \frac{\partial \mathbf{J}}{\partial V}=\eta \cdot \mathbf{error}\cdot \mathbf{sign}(\frac{\mathbf{\Delta y}}{\mathbf{\Delta u}})\cdot \mathbf{E}\cdot {\mathbf{g}}^{\prime }(\mathbf{Oi})\cdot \mathbf{Ho}$$

这里面的矩阵向量直接相乘是不对的，会出现维度不满足矩阵相乘条件的情况。但其实也只需要做一些转置或者元素直接相乘就可以完成。具体的原理用下面的两个图来描述。

这里就是求偏导链式法则的体现，只不过以矩阵形式展示可以简化表达式，这里就不再赘述了，自行推导！

matlab仿真

废话不多说，直接上代码

close,clear all;
warning off
xita = 1; %学习因子
alfa = 0.05; %惯量因子
hl_number = 5; %隐藏层神经元个数

hidding_wi = 0.5*rands(4,hl_number); %4xN
output_wi = 0.5*rands(hl_number+1,3);  %Nx3
%抽奖

%初始化隐含层和输出层的两组权重系数
hidding_wi_1 = hidding_wi;
hidding_wi_2 = hidding_wi;
output_wi_1 = output_wi;
output_wi_2 = output_wi;
%权重系数记录

error_1 = 0;error_2 = 0;
u_1 = 0;u_2 = 0;u_3 = 0;
y_1 = 0;y_2 = 0;y_3 = 0;
%控制律及系统状态序列
last_error_ventor = [0;0;0];
ts = 0.001;
for t = 1:1:5000
    time(t) = t*ts; 
    rin(t) = 1*sin(5*2*pi*t*ts);
    %输入信号
    a(t) = 1.2*(1-0.8*exp(-0.1*t));
    y(t) = a(t)*y_1/(1+y_1^2)+u_1; 
    %系统模型输出
    J(t) = 0.5*(rin(t)-y(t))^2;
    
    error(t) = rin(t) - y(t);
    %计算误差
    x_input = [rin(t),y(t),error(t),1];
    %神经网络输入向量  1x4
    error_ventor = [error(t)-error_1;error(t);error(t)-2*error_1+error_2];
    %计算误差向量  3x1

    hiddingLayer_input = x_input * hidding_wi;
    %计算隐含层输入  1xN
    
    hiddingLayer_output = (exp(hiddingLayer_input) - exp(-hiddingLayer_input))./(exp(hiddingLayer_input) + exp(-hiddingLayer_input));
    %计算隐含层输出 激励函数为 (exp(x) - exp(-x))/(exp(x) + exp(-x) ) 
    %1xN
   
    outputLayer_input_offset = [hiddingLayer_output(:);1]';
    %加入输出层神经元的偏置  1x(N+1)
    outputLayer_input = outputLayer_input_offset * output_wi;
    %计算输出层输入 1x3

    outputLayer_output = exp(outputLayer_input) ./ (exp(outputLayer_input) + exp(-outputLayer_input));
    %计算输出层输出 1x3 
    %激励函数为 exp(x) / (exp(x) + exp(-x))
    
    delta_u = outputLayer_output * error_ventor;
    %计算控制律增量
  
    u(t) = u_1 + delta_u;
    %计算控制律

    dJdy = error_1;
    dydu = sign((y_1-y_2)/(u_1-u_2+0.0000001));
    dudoutput3 = last_error_ventor;
    %3x1
    doutput3dnet3 = 2 ./((exp(outputLayer_input)+exp(-outputLayer_input)).^2);
    %1x3  g'(x)
    dnet3dv = outputLayer_input_offset;
    %1x(N+1)
    delta_output_wi = (xita * dJdy * dydu) .* dudoutput3 .* doutput3dnet3' * dnet3dv + alfa*(output_wi_1-output_wi_2)';
    %计算输出层权重增量
    output_wi_2 = output_wi_1;
    output_wi_1 = output_wi;
    output_wi = output_wi_1 + delta_output_wi';
    %更新输出层权重系数
    
    dJdy = error_1;
    dydu = sign((y_1-y_2)/(u_1-u_2+0.0000001));
    dudoutput3 = last_error_ventor;
    %3x1
    doutput3dnet3 = 2 ./((exp(outputLayer_input)+exp(-outputLayer_input)).^2);
    %1x3
    dnet3dv = hiddingLayer_output;
    %1xN
    dnet3doutput2 = output_wi_1(1:hl_number,:);
    % Nx3
    doutput2dnet2 =  4 ./(exp(hiddingLayer_input)+exp(-hiddingLayer_input)).^2;
    % 1x(N)  针对的是N个隐含层神经元
    dnet2dw = x_input;
    % 1x4  针对三个输入和一个偏置   
    delta_hidding_wi = ( ( dnet3doutput2 * ( (xita * dJdy *dydu ) .* dudoutput3 .* doutput3dnet3' ) ) .* doutput2dnet2' ) * dnet2dw + alfa * (hidding_wi_1 - hidding_wi_2)';
    %计算隐含层权重增量
    hidding_wi_2 = hidding_wi_1;
    hidding_wi_1 = hidding_wi;
    hidding_wi = hidding_wi_1 + delta_hidding_wi';
    %更新隐含层权重系数

    error_2 = error_1;
    error_1 = error(t);
    y_3 = y_2;y_2 = y_1;y_1 = y(t);    
    u_3 = u_2;u_2 = u_1;u_1 = u(t);
    last_error_ventor = error_ventor;
    %更新序列
end

%绘图
figure(1);
plot(time,rin,'r',time,y,'b','linewidth',2);
xlabel('time(s)');ylabel('rin,yout');   % 输入和输出

运行结果如下图

总结

从测试的规律来看，学习因子对于整个神经网络的影响较大。过大的学习因子会让系统出现振荡，过小的学习因子则会使系统在初始阶段发散，无法收敛（系统参数的初始值是随机产生的）。当然也可以考虑训练出一组较为理想的参数作为系统初始值，再以较小的学习因子来提高控制器的精度，使系统的动态性能提高。
另外，神经网络的迭代周期对于控制器性能的影响也较大，理想上周期越小，系统状态拟合的更好，收敛的越快。文章来源地址https://www.toymoban.com/news/detail-454488.html

到了这里，关于浅谈BP神经网络PID控制算法及matlab仿真的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！