自适应PID算法学习（01）——单神经元PID控制

神经元PID

单神经元结构

  单神经元含有n个输入，仅1个输出，每个输入端可记作$x_i(i=1,2,...n)$，若该神经元为多元组成网络中某一层（输入层/输出层/隐含层）其中的一个单元，记该神经元输出为$o_j(j=1,2,...)$；不同输入端进入该单元的连接权值不同，记作${\omega }_{ij}$，表示连接第$j$个神经元的第$i$个输入端的连接权值。
  另外，神经元在接收多个输入时，应有累加-整合的过程，即先有：$ne{t}_j={\sum }_{i=1}^n{\omega }_{ij}\cdot x_{ij}$，再经历：$o_j=f(ne{t}_j)$才得到输出。累加运算也可表述为矩阵运算形式如：n维权向量${\omega }_j$的转置（列向量）与n维输入向量$x_j$（行向量）相乘，得到标量$ne{t}_j$，即有：
$$ne{t}_j={\omega }_j^T\cdot x_j\text{\hspace{1em}(1)}$$
  整合$ne{t}_j$的$f$即称为激活函数，可选用各种有阈值限制的非线性函数，如常用的有：切换函数$f(net)=sgn(net)$、反正切函数$f(net)=atan(net)$，等等。

学习规则

  神经元/神经网络的学习规则，即修正（更新）连接权值所选用的算法，可分类为有监督学习或无监督学习。

无监督Hebb学习

  连接权值的修正与输入$x_i$和输出$o_j$的乘积成正比，即有：
$$\Delta {\omega }_{ij}^{(k)}=\eta \cdot (x_i^{(k)}\cdot o_j^{(k)})\text{\hspace{1em}(2)}$$

有监督Delta学习

  在无监督Hebb的基础上，引入教师信号，将输出$o_j$替换为实际输出$o_j$相对期望输出$d_j$的误差，即有：
$$\Delta {\omega }_{ij}^{(k)}=\eta \cdot [x_i^{(k)}\cdot (d_j^{(k)}-o_j^{(k)})]\text{\hspace{1em}(3)}$$

有监督Hebb学习

  结合无监督Hebb和有监督Delta，即连接权更新正比于输入、输出、输出误差的乘积，即有：
$$\Delta {\omega }_{ij}^{(k)}=\eta \cdot [(d_j^{(k)}-o_j^{(k)})\cdot x_i^{(k)}\cdot o_j^{(k)}]\text{\hspace{1em}(4)}$$
上述各式中：
  上角标$(k)$表示当前的迭代轮次；$\Delta {\omega }_{ij}$为计算得的连接权值的修正量，即有：
$${\omega }_{ij}^{(k+1)}={\omega }_{ij}^{(k)}+\Delta {\omega }_{ij}^{(k+1)}\text{\hspace{1em}(5)}$$

神经元PID控制过程

输入更新

  对于单神经元PID控制而言，我们的目的是要自适应地更新三个控制器参数Kp，Ki，Kd，输入维数n=3；而PID控制器的输入量与传统PID一样，为当前被控对象的实际输出$y^{(k)}$相对目标输出$y_d^{(k)}$的误差，即有：
$$e^{(k)}=y_d^{(k)}-y^{(k)}\text{\hspace{1em}(6)}$$
  若PID控制器模式选用增量式PID，即$\Delta u^{(k)}=K_p \cdot (e^{(k)}-e{(k-1)})+K_i \cdot $对于神经元输入$x_{ij}$，应有：
$$\begin{array}{rl}& x_{1j}^{(k)}=e^{(k)}-e^{(k-1)}\\ & x_{2j}^{(k)}=e^{(k)}\\ & x_{3j}^{(k)}=e^{(k)}-2e^{(k-1)}+e^{(k-2)}\end{array}\text{\hspace{1em}(7)}$$

连接权值更新、归一化

  在已有更新的输入值$x_{ij}^{(k)}$和已有的输出值$o_j^{(k-1)}$后，即可更新本轮次神经元运算需用到的连接权值${\omega }_{ij}^{(k}$；设选用了有监督的Hebb学习规则，即同时调用了本轮更新的输出误差$e^{(k)}$，则有：
$$\begin{array}{rl}& {\omega }_{1j}^{(k)}={\omega }_{1j}^{(k-1)}+{\eta }_p\cdot [e^{(k)}\cdot x_1^{(k)}\cdot o_j^{(k-1)}]\\ & {\omega }_{2j}^{(k)}={\omega }_{2j}^{(k-1)}+{\eta }_i\cdot [e^{(k)}\cdot x_2^{(k)}\cdot o_j^{(k-1)}]\\ & {\omega }_{3j}^{(k)}={\omega }_{3j}^{(k-1)}+{\eta }_d\cdot [e^{(k)}\cdot x_3^{(k)}\cdot o_j^{(k-1)}]\end{array}\text{\hspace{1em}(11)}$$
式中，${\eta }_p$、${\eta }_i$、${\eta }_d$分别为神经元PID的比例项、积分项、微分项的学习速率。
  在得到上述更新的连接权值${\omega }_{ij}^{(k)}$后，还需要做类似于单位化的归一化操作，即有：
$${\omega }_{ij}^{(k)}={\omega }_{ij}^{(k)}/\sum _{i=1}^n|{\omega }_{ij}^{(k)}|,i=1,2,3\text{\hspace{1em}(10)}$$
  不难理解可知，一直在更新并用于与式(7)所示的$x_{1j}$、$x_{2j}$、$x_{3j}$相乘再累加的连接权值${\omega }_{1j}$、${\omega }_{2j}$、${\omega }_{3j}$，乘上神经元比例系数$ϵ$，即是我们神经元PID的自适应参数$K_p$、$K_i$、$K_d$。

输出更新

  单神经元的输出如同普通PID一样，为输入到被控对象的供给量$u$（控制器输出，即控制律），即经过神经元，输出得到控制律$u^{(k)}$变化量。另外，当前轮次k用到的是当前轮次k的输入、连接权，即有：
$$ne{t}_j^{(k)}=\sum _{i=1}^n{\omega }_{ij}^{(k)}\cdot x_{ij}^{(k)}\text{\hspace{1em}(8)}$$
  整合函数无非线性要求，直接成比例输出，则控制律更新有：
$$u^{(k)}=u^{(k-1)}+ϵ\cdot ne{t}_j^{(k)}\text{\hspace{1em}(9)}$$
式中，$ϵ$为神经元的比例系数。当然，神经元输出符号代换有：
$$o_j^{(k)}=u^{(k)}\text{\hspace{1em}(10)}$$
  上述各part表达为控制流图，即如下所示：
![]

代码实现和展示

主程序

  编写成Matlab代码，程序框架（主程序）如下所示：

close all;clear
% 迭代次数：
Num=2000;
% 迭代周期：
T=0.001;     % T=0.01s

% 连接权值初始化：
omega=[0.1,0.1,0.1];
% 输入变量初始化：
x=[0,0,0]';
% 输出变量初始化：
u(1)=0;     % u(k)
u_1=0;      % u(k-1)
u_2=0;      % u(k-2)
u_3=0;      % u(k-3)
net=0;      % net=omega'*x; u(k)=u(k-1)+eps*net
y(1)=0;     % 被控对象实际输出
y_1=0;      % y(k-1)
y_2=0;      % y(k-2)
% 误差变量初始化：
err_0(1)=0; % e(k)
err_1=0;    % e(k-1)
err_2=0;    % e(k-2)
% 控制律输出限幅:
u_max=1.2;
u_min=-1.2;
% 神经元参数设置：
eps=0.12;   % 神经元比例系数
eta_p=0.4; % 比例项学习速率
eta_i=0.35;  % 积分项学习速率
eta_d=0.4;  % 微分项学习速率
w1(1)=0;
w2(1)=0;
w3(1)=0;

%%% 进入控制-训练轮次：
for k=1:Num
   % 目标输出：
   yd(k)=Target_Oput(k,T);
   % 输入到被控对象，实际输出：
   y(k)=Controlled_Obj(y_2,y_1,u_2,u_1);
   % 计算输出误差：
   err_0(k)=yd(k)-y(k);

   % （神经元）输入更新：
   x(1)=err_0(k)-err_1;
   x(2)=err_0(k);
   x(3)=err_0(k)-2*err_1+err_2;
   
   % （神经元）连接权值更新
   omega=Hebb_Learning(omega,eta_p,eta_i,eta_d,...
       err_0(k),x,u_1);
   omega=Omega_Norm(omega);
   
   % （神经元）输出更新：
   u(k)=u_1+eps*omega*x;       % 输出控制律
   if(u(k)>u_max) u(k)=u_max;
   elseif(u(k)<u_min) u(k)=u_min;end
   
   % 误差序列更新：
   err_2=err_1;
   err_1=err_0(k);
   % 控制序列更新：
   u_3=u_2;
   u_2=u_1;
   u_1=u(k);
   % 输出序列更新：
   y_2=y_1;
   y_1=y(k);
   
   % 用于显示Kp,Ki,Kd：
   w1(k)=eps*omega(1);
   w2(k)=eps*omega(2);
   w3(k)=eps*omega(3);
end  

%%% 画图展示:
time=T:T:Num*T;
figure  % 输出响应
plot(time,yd,'r');hold on;
plot(time,y,'b');title("输出响应");
ylabel('yt');xlabel('time');legend('目标输出-yd','实际输出-y');
figure  % 误差响应+控制律变化
subplot(2,1,1);plot(time,err_0,'-r');title("误差响应");
ylabel('error');xlabel('time');
subplot(2,1,2);plot(time,u,'-g');title("控制律变化");
ylabel('ut');xlabel('time');
figure  % PID参数变化
subplot(3,1,1);plot(time,w1,'r');ylabel('Kp');title("PID参数");
subplot(3,1,2);plot(time,w2,'b');ylabel('Ki');
subplot(3,1,3);plot(time,w3,'k');ylabel('Kd');xlabel('time');

子函数——连接权值更新（使用有监督Hebb学习规则）

  在每轮神经元运算之前，更新连接权值${\omega }_{ij}$，这里示例用的是有监督Hebb学习规则，对应公式(11)，如下所示：

function [omega_new]=Hebb_Learning(omega,eta_p,eta_i,eta_d,...
    err,x,u)
    w1=omega(1);
    w2=omega(2);
    w3=omega(3);
    x1=x(1);
    x2=x(2);
    x3=x(3);
    if err>1;err=-err;end    % 随便加的，合法有效但有病
    w1=w1+eta_p*(err*x1*u);
    w2=w2+eta_i*(err*x2*u);
    w3=w3+eta_d*(err*x3*u);
    omega_new=[w1,w2,w3];
end

子函数——连接权值归一化

  对应公式(10)，完成连接权值${\omega }_{ij}$更新后的归一化（单位化）步骤，如下所示：

function [omega_new]=Omega_Norm(omega)
    w1=omega(1);
    w2=omega(2);
    w3=omega(3);
    wsum=abs(w1)+abs(w2)+abs(w3);
    w1=w1/wsum;
    w2=w2/wsum;
    w3=w3/wsum;
    omega_new=[w1,w2,w3];
end

子函数——目标输出设置

  设计跟踪目标输出为上下幅值为±1的正负阶跃信号（方波信号），如下示例有：

function [yd_k]=Target_Oput(k,T)
    yd_k=sign(sin(2*pi*k*T));
    % yd_k=sin(2*pi*k*T);
end

  注释中的是将目标输出切换为正弦函数$y_d=sin(2\pi t)$，你也可以试着复现下看看效果~

子函数——控制对象设置

  相当于我们在控制系统中所见到的传递函数G(s)的部分，只不过这里的对象是离散系统，并且将该数学描述换用了状态空间表达式的形式，具体怎么推导的可以看看我这段子程序开头的注释，如下有：

% 被控对象的状态空间表达式：
% dx=A*x+B*u；
% 令x1=y(k-2),x2=y(k-1)，则:dx1=y(k-1)=x2，dx2=y(k);
% y(k-1)=[0 1]*[y(k-2) y(k-1)]'+0*u，不用写;
% y(k)=[a21 a22]*[y(k-2) y(k-1)]'+[b21 b22]*[u2 u1]';
% u1=u(k-1),u2=u(k-2);
function [y_k]=Controlled_Obj(y_2,y_1,u_2,u_1)
    y_k=0.26*y_2+0.368*y_1+0.632*u_2+0.1*u_1;
end

运行结果展示

运行结果展示

  运行主程序代码（程序要“添加到路径”），输出该单神经元PID控制仿真结果（使用的是有监督Hebb学习规则），跟踪三个上下摆动周期，如下所示：
  目标输出–实际输出跟踪曲线：

  输出误差曲线+控制律输入曲线

  P-I-D系数迭代变化曲线

  

未完待续 ~ ~