基于某BP神经网络地自整定PID控制仿真

1、基于BP神经网络的自整定PID控制仿真实验六 基于BP神经网络的自整定PID控制仿真一、目的 1.熟悉神经网络的特征、结构以及学习算法 2.掌握神经网络自整定PID的工作原理3.了解神经网络的结构对控制效果的影响 4.掌握用MATLAB实现神经网络控制系统仿真的方法。 二、设备及条件 计算机系统 Matlab仿真软件 三、问题背景 在工业控制中，PID控制是工业控制中最常用的方法。这是因为PID控制器结构简单、实现简单，控制效果良好，已得到广泛应用。但是，PID具有一定的局限性：被控制对象参数随时间变化时，控制器的参数难以自动调整以适应外界环境的变化。为了使控制器具有较好的自适应性，实现控制器

2、参数的自动调整，可以采用神经网络控制的方法。利用人工神经网络的自学习这一特性，并结合传统的PID控制理论，构造神经网络PID控制器，实现控制器参数的自动调整。 基于BP神经网络的PID控制器结构如图1所示。控制器由两部分组成：一是常规PID控制器，用以直接对对象进行闭环控制，且三个参数在线整定；二是神经网络NN，根据系统的运行状态，学习调整权系数，从而调整PID参数，达到某种性能指标的最优化。 图1 基于神经网络的PID控制器结构四、基于BP神经网络的PID设计1设计原理 神经网络采用结构为4-5-3型的BP网络，如图2。图2 BP网络结构其中，输出层激励函数取非负的Sigmoid函数，隐层取

3、正负对称的Sigmoid函数。被控对象为一时变非线性对象，数学模型可表示为： 式中，系数a(k)是慢时变的，。为保证控制器有一定的动态跟踪能力，选定神经网络的输入层输入为网络的学习过程由正向和反向传播两部分组成。如果输出层不能得到期望输出，那么转入反向传播过程，通过修改各层神经元的权值，使得输出误差信号最小。 输出层节点分别对应三个可调参数 取性能指标函数为： 设其中 若PID控制器采用采用增量式数字PID控制算法，则有 2.网络权系数调整网络权系数的修正采用梯度下降法。根据相关数学知识，针对指定因变量的梯度代表了使因变量增速最大的自变量变化方向，故而其反方向代表了因变量下降最快的自变量变化方

4、向，如果我们选取性能指标为因变量，网络各层权系数为自变量，则对应梯度的负方向就是权系数调整的最佳方向，因此，实际上构成了一个有目标的搜索算法，对最终结果的收敛性提供了有力的保证。对应于本题采用的4-5-3型的BP网络，梯度负方向的计算：（1） 隐含层-输出层：其中：为指标函数 为隐含层-输出层权系数矩阵元素 为被控对象输出 为PID控制器输出 为输出层输出为输出层输入根据所选用神经网络的数学模型，易知：其中: 为隐含层输出，为输出层激励函数，为其偏导数。另外，直接的数学表达不容易获得，但我们可以使用它的符号函数来近似，仍可以保证参数修正方向的正确性，而由此造成其模的误差只影响参数调整的速度，它

5、可以通过调整学习速率来得以补偿。故而最终有：（2） 输入层-隐含层：其中： ，为输出层激励函数，为其偏导数。故而最终有：（3） 网络权值修正量梯度的负方向给出了网络权值修正的方向，乘以系数，即为权值修正量。其中代表了算法每次在梯度负方向搜索的步长，称为网络的学习速率，另外，如果考虑上次权值对本次权值变化的影响，需要加入动量（平滑）因子 ，此时新的权值为:3.程序流程 步骤1：设定初始状态和参数初始值，包括随机产生初始BP神经网络权值系数，设定初始输入输出值为零，设定学习速率和惯性系数，计数器设为k=1，并设定计数上限等。步骤2：计算产生BP神经网络隐含层输入。本程序为采样获得e（k），并结合储

6、存的e（k-1），e（k-2），及常数1作为隐含层输入。前两次的e（k-1），e（k-2）并未真实产生直接取0.步骤3：前向传播计算。包括:（1）BP神经网络前向传播计算，得到输出层输出，, .；（2）增量式PID控制器计算控制器输出；（3）被控对象模型计算输出值步骤4：反向传播计算。包括：（1） 修正输出层的权系数；（2） 修正隐含层的权系数；步骤5：参数更新步骤6：如果k达到设定的次数上限，则结束；否则，k=k+1，并返回步骤2.总的程序流程图如下：图3 程序流程图五、仿真程序及分析学习速率取0.25，平滑因子取0.05。经过多次运行，最终得到一个较好的结果，并将此结果的稳态权重作为初始权

7、重。在调试中发现，即使稍微改变学习速率和平滑因子，还是需要经过多次运行，才能得到一个较好的结果。如果用得到的较好结果的稳态权重作为初始权值，学习速率和平滑因子的取值可以更加随意。（1）输入为阶跃信号(2)输入为正弦信号程序代码：%基于BP神经网络的自整定PID控制仿真clear all;close all;xite = 0.25; %学习速率alfa = 0.05; %平滑因子S=2; %选择输入信号的类型（1：阶跃信号；2：正弦信号）IN = 4; %输入层、隐含层、输出层节点数H = 5;OUT= 3;if S=1 wi= -0.4129 -0.2553 -0.7973 -0.1004 -

8、0.2771 0.2676 0.4234 0.3484 -0.6914 0.2740 -0.1590 -0.8642 -0.3915 0.1627 -0.6956 -0.7668 -0.4133 0.2296 -0.5542 -0.2671;% wi = 0.5*rands(H,IN); %权值系数初始值 wi_1 = wi; wi_2 = wi; wo= 0.5661 0.2004 0.9433 0.1832 0.5971 0.4185 0.2750 0.6734 0.9408 0.4597 0.4348 0.0402 0.9523 0.8143 0.2773;% wo = 0.5*rand

9、s(OUT,H); wo_1 = wo; wo_2 = wo;else wi= -0.4257 0.3975 -0.2219 0.0629 -0.3548 -0.4002 -0.3585 -0.2581 -0.4162 -0.4559 -0.1764 -0.5627 -0.0390 0.0586 -0.3188 -0.3945 -0.3042 0.2744 0.0130 -0.3788;% wi = 0.5*rands(H,IN); %权值系数初始值 wi_1 = wi; wi_2 = wi; wo= 0.2279 0.5382 -0.1358 0.6441 0.0699 0.1584 0.2

10、123 0.1166 0.4402 0.6627 0.7242 0.6211 0.6540 0.3450 -0.1486;% wo = 0.5*rands(OUT,H); wo_1 = wo; wo_2 = wo;endu_1 = 0; y_1 = 0;Oh = zeros(H,1);error_2 = 0;error_1 = 0;ts = 0.001;for k = 1:6000 time(k) = k*ts; if S = 1 rin(k) = 1.0; else rin(k) = sin(2*pi*k*ts); end %被控对象为非线性模型 a(k) = 1.2*(1-0.8*exp(

11、-0.1*k); yout(k) = a(k)*y_1/(1+y_12)+u_1; error(k) = rin(k) - yout(k); xi = error_2 error_1 error(k) 1; %输入层输入 epid = error(k)-error_1;error(k);error(k)-2*error_1+error_2; %计算隐含层输出 I = xi*wi; for j = 1:H Oh(j) = (exp(I(j)-exp(-I(j)/(exp(I(j)+exp(-I(j); end %计算输出层输出 K = wo*Oh; for l = 1:OUT O(l) = ex

12、p(K(l)/(exp(K(l)+exp(-K(l); end %得到kp,ki,kd的值 kp(k) = O(1); ki(k) = O(2); kd(k) = O(3); Kpid = kp(k) ki(k) kd(k); %计算控制器的输出 du(k)= Kpid * epid; u(k) = u_1 + du(k); if u(k) 10 u(k) = 10; end dyu(k) = sign(yout(k)-y_1)/(u(k)-u_1+0.); %输出层权系数修正 for l = 1:OUT dK(l) = 2/(exp(K(l)+exp(-K(l)2; end for l =

13、1:OUT delta3(l) = error(k)*dyu(k)*epid(l)*dK(l); end for l = 1:OUT for j = 1:H d_wo = xite*delta3(l)*Oh(j)+alfa*(wo_1-wo_2); end end wo = wo_1+d_wo+alfa*(wo_1-wo_2); %隐含层权系数修正 for j = 1:H dO(j) = 4/(exp(I(j)+exp(-I(j)2; end sigma = delta3*wo; for j = 1:H delta2(j) = dO(j)*sigma(j); end d_wi = xite *

14、 delta2 * xi; wi = wi_1 + d_wi + alfa * (wi_1-wi_2); %参数更新 u_1 = u(k); y_1 = yout(k); wo_2 = wo_1; wo_1 = wo; wi_2 = wi_1; wi_1 = wi; error_2 = error_1; error_1 = error(k);endfigure(1);plot(time,rin,r,time,yout,b);xlabel(time(s);ylabel(rin,yout);figure(2);plot(time,error,r);xlabel(time(s);ylabel(error);figu