课程笔记神经网络原理及应用
,

单神经元PID控制器的设计与仿真

文章目录
  1. 1. 绪论
    1. 1.1. 神经网络
    2. 1.2. 单神经元原理
    3. 1.3. 单神经元PID控制器原理
    4. 1.4. 单神经元PID控制器改进
  2. 2. 单神经元PID控制器仿真分析
    1. 2.1. 仿真模型搭建
    2. 2.2. S-function模块
    3. 2.3. S-function编写
  3. 3. 参考资料

预习神经网络原理及应用课程大作业所整理的笔记。
本文将介绍单神经元PID控制器原理、有监督Hebb学习规则、单神经元PID控制器的改进、Simlink仿真具体步骤、S-function的使用方法。依照上述原理与方法完成单神经元PID控制器的设计与仿真。
下一篇文章将介绍基于BP神经网络的PID控制器的设计与仿真。

绪论

神经网络

人工神经网络(Artificial Neural Networks，ANN)，简称神经网络，是指一类受生物学和神经科学启发得来的的数学模型。主要通过对人脑的神经元网络进行抽象，构建人工神经元，并按照一定拓扑结构来建立人工神经元之间的连接，来模拟生物的神经网络。
人工神经网络由许许多多信息处理单位（对应于人脑中的神经元）根据不同的相互作用强度（即连接权值）互连而成，其实质是对生物学上的人脑神经体系的一种近似与模拟，人工神经网络能够像人体神经系统一样对外界知识进行检索、处理与存储。

单神经元原理

单神经元，又被称作MP神经元，由心理学家McCulloch和数学家Pitts于1943年根据生物神经元的结构提出。
图1 单神经元模型

单神经元模型如上图所示，其中\(x_1,x_2,\cdots,x_n\)为接收到的信息，\(\omega_{i1},\omega_{i2},\cdots,\omega_{in}\)为连接的权值。将各个信息值按照权重相加并做一定处理，所得到的结果被称为“净输入”，其表达式为：
\(net_i=\sum\limits_{j = 1}^n \omega_{ij}x_j-\theta_i\)
式中\(\theta_i\)为神经元\(i\)的阈值。
最终的输出值\(y_i\)与净输入\(net_i\)的关系式为：\(y_i=g(net_i)\)，这里的\(g\)被称为激活函数。

单神经元PID控制器原理

图2 单神经元PID控制框图

单神经元PID控制框图如上图所示，图中\(r_{in}\)为系统输入值（设定值），\(y_{out}\)为系统实际输出值，输入转换模块将\(r_{in}\)与\(y_{out}\)变换为神经元的输入状态变量\(x_i(k)\)，这三个状态变量与PID控制器的三个参数（\(K_p,K_i,K_d\)）相关，PID控制器采用增量式算法，其中：
\(x_1(k)=e(k)\)，对应PID控制的积分项，与误差量相等；
\(x_2(k)=\Delta e(k)=e(k)-e(k-1)\)，对应PID控制的比例项，反映误差变化的方向；
\(x_3(k)=\Delta^2 e(k)=e(k)-2e(k-1)+e(k-2)\)，对应PID控制的微分项。

取\(\omega_i(k)\)为\(x_i(k)\)的神经元连接权值，取\(K\)为神经元的比例系数（要求\(K>0\)），则有：
\(u(k)=u(k-1)+K\sum\limits_{i = 1}^3 \omega_i(k)x_i(k)\)
\(\Delta u(k)=K[\omega_1(k)e(k)+\omega_2(k)\Delta e(k)+\omega_3(k)\Delta^2 e(k)]\)
由此可知，\(K\)值越大，系统整体速度越快，但同时超调量也会增多，而且影响系统整体稳定性。

神经元通过有监督型Hebb学习，对神经元连接权值进行实时修正，以实现对PID控制器参数的自适应调整。
这里取教师信号（系统性能指标）\(z(k)=e(k)\)，得到\(k\)时刻的连接权值为：
\(\omega_i(k)=\omega_i(k-1)+\eta z(k)u(k)x_i(k)\)
式中\(\eta\)为学习速率，取大于零的较小值。

单神经元PID控制器改进

为解决当\(\eta\)的值较大时，单神经元PID的学习速度增大的同时可能造成的收敛速度降低问题，可以让输入状态变量的连接权值以不同学习速率分开调整。
这里对不同状态变量采用不同的连接权值学习速率\(\eta_P,\eta_I,\eta_D\)，可以得到：
\(\omega_1(k)=\omega_1(k-1)+\eta_I z(k)u(k-1)x_1(k)\)
\(\omega_2(k)=\omega_2(k-1)+\eta_P z(k)u(k-1)x_2(k)\)
\(\omega_3(k)=\omega_3(k-1)+\eta_D z(k)u(k-1)x_3(k)\)
\(u(k)=u(k-1)+K\sum\limits_{i=1}^3 \omega_i^{'}(k)x_i(k)\)
\(\omega_i^{'}(k)=\frac{\omega_i(k)}{\sum\limits_{j=1}^3 |\omega_j(k)|}\)（归一化处理）

而实际经验表明，PID参数不仅与\(e(k)\)有关，还与其变化值\(\Delta e(k)\)有关，将后者也加入到连接权值调整公式中，可以得到：
\(\omega_1(k)=\omega_1(k-1)+\eta_I z(k)u(k-1)[e(k)+\Delta e(k)]\)
\(\omega_2(k)=\omega_2(k-1)+\eta_P z(k)u(k-1)[e(k)+\Delta e(k)]\)
\(\omega_3(k)=\omega_3(k-1)+\eta_D z(k)u(k-1)[e(k)+\Delta e(k)]\)

单神经元PID控制器仿真分析

仿真模型搭建

使用Simulink构建一个简单的连续系统仿真模型，如下图所示：
图3 单神经元PID控制系统仿真模型
其中输入为阶跃信号，传递函数为\(G(s)=\frac{1}{20s+1}\)，单神经元PID控制器如下图所示：
图4 单神经元PID控制器系统仿真模型

仿真模型最中央的“exp_pidf”是一个S函数(S-function)模块，它是一个Matlab脚本文件，在仿真模型中实现单神经元PID的相关运算。
S函数模块左侧的Mux总线与该函数的五个输入值相连，它们分别是\(e(k),e(k-1),e(k-2),u(k-1),u(k)\)。
这里\(e(k)\)是误差变量，也是整个PID控制器的输入变量。在\(e(k)\)的基础上串联一个延时环节(Unit Delay)就会得到\(e(k-1)\)，他表示的是上一次的误差。同理串联两个延时环节得到\(e(k-2)\)。
\(u(k)\)是控制器输出量，同理可知串联一个延时环节得到的\(u(k-1)\)是上一次的输出量。
S函数模块右侧的Demux总线与该函数的四个输出值相连，它们是控制器输出值\(u(k)\)与PID的三个参数\(K_i,K_p,K_d\)。
这里\(u(k)\)是必须输出的值，而输出PID的三个参数值，是为了观察仿真运行过程中三个参数的变化情况。

S-function模块

S函数(S-function)是System-function系统函数的缩写，是指采用非图形化的方式（计算机语言）描述的一个功能块。
首先要在Simulink中，找到S函数模块，打开库浏览器，找到Simulink->User-Defined Functions->S Function。
图5 S函数模块参数

打开S函数模块，可以看到“模块参数”界面，有三个空。
第一个是S-Function名称，这个名称就是自己创建的matlab脚本的名称，脚本与仿真文件需要在同一文件夹内，点击“编辑”可以打开该脚本文件。
第二个是S-Function参数，它既不是输入参数也不是输出参数，而是外部参数，这些参数将被导入到脚本文件中参与计算。本次仿真外部参数为不同的连接权值学习速率\(\eta_P,\eta_I,\eta_D\)和神经元比例系数\(K\)，填写为[ni;np;nd],K。
最后一个是S-Function模块，在S-function是用C语言编写并用MEX工具编译的C-MEX文件时，需要填写该参数，这里不做修改。

而如果导入了外部参数，还需要进行另外的一步，那就是封装。
右键S函数模块->封装->编辑封装，就能打开封装编辑器(Mask Editor)。
图6 封装编辑器
将之前的四个外部参数(np,ni,nd,K)按如图所示添加进去，这里的参数名称和前面“模块参数”中的名称保持一致。
完成这一步后，再打开S函数模块，会直接弹出设置外部参数的对话框，更方便调整参数。
图7 参数设置

S-function编写

下面介绍S函数的编写原则，并基于此写出单神经元PID控制器的S函数。

1
function [sys,x0,str,ts,simStateCompliance] = exp_pidf(t,x,u,flag,deltak,K)
这是主函数的第一行，函数的定义，exp_pidf是函数名，其右侧为输入参数（包括外部参数），等号左侧是输出参数。

t是时间变量；
x是状态变量，对应前面的\(x_i(k)\)或\(\omega_i(k)\)；
u是输入向量，对应仿真文件中S函数模块与左侧Mux总线相连的五个输入值\(e(k),e(k-1),e(k-2),u(k-1),u(k)\)，它们在S函数程序中被表示为\(u(1),u(2),u(3),u(4),u(5)\)；
flag是标志位，在仿真程序运行的不同阶段，会根据标志位运行不同的子函数。
deltak是外部参数，这里deltak=[ni;np;nd]，表示控制器积分、比例、微分的学习速率，它是一个列向量。
K是外部参数，神经元比例系数。这里的外部参数与前面“模块参数”界面中设置的参数要一一对应。

sys是程序返回值，它不一定是控制器输出量，返回的值取决于flag值，也就是不同的子函数；
x0是状态变量初始值；
str默认为空，无需设置；
ts是采样时间，包含采样时间和偏移量；
simStateComplicance是附加变量；

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
switch flag,
  case 0,
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
%调用初始化函数
  case 2,
    sys=mdlUpdate(t,x,u,deltak);
%更新状态变量，此算法有三个状态变量，三个权值系数w1(k),w2(k),w3(k)
  case 3,
    sys=mdlOutputs(t,x,u,K);
%调用输出函数，输出控制值信号
  case {1,4,9},
    sys=[];
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end

接下来完成这个主函数，它是一个switch-case结构，选择不同的标志位以进入不同的子函数。
其中flag=0时进入初始化函数，flag=2时更新状态变量，flag=3时计算控制器输出值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
%初始化函数
sizes = simsizes;
sizes.NumContStates  = 0;
sizes.NumDiscStates  = 3;
%定义三个离散的状态变量，w1(k),w2(k),w3(k)
sizes.NumOutputs     = 4;
%定义四个输出，算法中的控制变量u（k)和三个PID参数Kp,Ki,Kd
sizes.NumInputs      = 5;
%定义该系统有五个输入
sizes.DirFeedthrough = 1;
%直接馈通，输出和输入有关存在直接馈通参数为1
sizes.NumSampleTimes = 1;  
%采样时间个数，默认值
sys = simsizes(sizes);
x0  = [0.1,0.1,0.1];
%定义三个状态变量的初始值，即w1(0)=0.1,w2(0)=0.1,w3(0)=0.1
str = [];
ts  = [-1 0];
%S-function设置为[-1 0]为按照其所连接块的速率来运行
simStateCompliance = 'UnknownSimState';

初始化函数如上所示，主要关注离散的状态变量数量(NumDiscStates)、输出量数量(NumOutputs)、输入量数量(NumInputs)。
状态变量初始值x0设成不为零的较小数，可以避免可能出现的程序报错。

1
2
function sys=mdlUpdate(t,x,u,deltak)
sys = x+deltak*u(1)*u(5)*(2*u(1)-u(2));

状态变量更新函数如上所示。
其中x为上一次状态变量，\(x=[\omega_1(k-1);\omega_2(k-1);\omega_3(k-1)]\)
sys为更新后状态变量，\(sys=[\omega_1(k);\omega_2(k);\omega_3(k)]\)
其对应的公式为：
\(\omega(k)=\omega(k-1)+\eta z(k)u(k-1)[e(k)+\Delta e(k)]\)

1
2
3
function sys=mdlOutputs(t,x,u,K)
xx = [u(1) u(1)-u(2) u(1)+u(3)-2*u(2)];
sys = [u(4)+K*xx*x/sum(abs(x));K*x/sum(abs(x))];

控制器输出值计算函数如上所示。 其中xx对应\(x_i(k)\)
sys的第一个参数是控制器输出值，对应公式：
\(u(k)=u(k-1)+K\sum\limits_{i=1}^3 \omega_i^{'}(k)x_i(k)\)
sys的后面三个输出值为PID控制器参数\(K_i,K_p,K_d\)，对应公式为：
\(K_i=K \omega_1^{'}(k)\)
\(K_p=K \omega_2^{'}(k)\)
\(K_d=K \omega_3^{'}(k)\)
注：这里的xx是行向量，x是列向量，只有这样向量相乘才会得到一个数。

参考资料