关于BP神经网络PID控制的一些问题

小弟最近在做一个毕业设计，是基于BP神经网络PID的温度控制，最近要中期检查了，要把仿真程序拿出来，使用matlab仿真程序仿真，我在网上找了很多，都是基本那个程序，但是程序中我有点不懂，还有仿真图像也不是我想要的那种，比如设置温度在20度，那么阶跃响应后的曲线最后徘徊在20左右，但是那个仿真程序我又改不来，不是很会，所以求助哈各位网友，希望网友能够提供一些方法和基于温度控制那种改好的程序，小弟不胜感谢，或者教我这个程序怎么改，温度传递函数应该是个惯性加延迟环节，比如G(s)=1/60S+1*exp(-200s),
%基于BP神经网络的PID控制器设计
clc;
clear;
NI=3;NH=10;NO=3;                          %选取网络结构为：3-10-3            
%输入层到隐层权值初始化
Wih=0.01*rands(NI,NH);Wih_1=Wih;Wih_2=Wih;
%隐层到输出层权值初始化
Who=0.01*rands(NH,NO);Who_1=Who;Who_2=Who;
%变量初始化
u=0;u_1=0;u_2=0;ly=0;err_1=0;err_2=0;
eta=0.5;                                   %给定学习率
a=0.1;                                     %给定动量项因子
ts=0.002;                                  %给定采样时间
%信号源选择：1—阶跃；2—正弦；3—方波；
souce=1;                                   %选择阶跃输入
%开始参数自学习实时控制
for k=1:1000
    t(k)=k*ts;
    %判断信号源
    switch souce
        case 1
            rin(k)=1;
        case 2
            rin(k)=sin(2*pi*t(k));
        case 3
            rin(k)=sign(sin(2*pi*t(k)));
        otherwise
            error('信号源错误!');         
    end
    %计算被控对象输出
    y(k)=1.2*(1-0.8*exp(-0.1*k))*ly/(1+ly^2)+u_1;   
    err(k)=rin(k)-y(k);                      %计算误差
    %-------------计算BP网络各层输入及输出-------------
    x1=err(k)-err_1; x2=err(k); x3=err(k)-2*err_1+err_2;
    %输入层输入输出，注意：必要时输入要做归一化
    N1i=[x1,x2,x3];
    N1o=N1i;
    %隐层输入输出
    N2i=N1o*Wih;
    N2o=tanh(N2i);                          
    %输出层输入输出
    N3i=Who'*N2o';                           
    N3o=1./(1+exp(-N3i));                    %注意：一定是点除！
    Kp(k)=N3o(1); Ki(k)=N3o(2); Kd(k)=N3o(3);%记录每次参数
    Kpid=[Kp(k),Ki(k),Kd(k)];                %PID参数
    %-----------计算控制量（使用增量式PID控制规律）---------
    epid=[x1;x2;x3];   
    du=Kpid*epid;                             
    u(k)=u_1+0.4*du;
    %------------------------修正权值------------------------
    %修正输出层权值
    dg=exp(-N3o)./((1+exp(-N3o)).^2);         %激活函数g(x)一阶导
    %加0.0001为防止开始时除数为零
    dlta3=err(k)*sign((y(k)-ly)/((u_1-u_2)+0.0001))*epid.*dg;
    for i=1:NO
        DWho(:,i)=eta*dlta3(i)*N2o';
    end
    Who=Who+DWho+a*(Who_1-Who_2);
    %修正隐层权值
    df=2*exp(-N2o)./((1+exp(-N2o)).^2);       %激活函数f(x)一阶导f'(x)   
    for i=1:NH
        dlta2(i)=Who(i,:)*dlta3*df(i);
    end
    for i=1:NH
        DWih(:,i)=eta*dlta2(i)*N1o';
    end
    Wih=Wih+DWih+a*(Wih_1-Wih_2);
    %----------更新数据-----------
    ly=y(k);
    u_2=u_1;        u_1=u(k);
    err_2=err_1;    err_1=err(k);
    Wih_2=Wih_1;    Wih_1=Wih;
    Who_2=Who_1;    Who_1=Who;
end
%绘图
figure(1)
plot(t,rin,t,y,t,u,t,err,'linewidth',1.5);
title('BP网络PID控制效果');
legend('期望曲线','跟踪曲线','控制量','跟踪误差');
xlabel('时间');ylabel('幅值');
figure(2)
subplot(3,1,1);plot(t,Kp,'linewidth',1.5);
title('PID参数变化曲线');ylabel('Kp');
subplot(3,1,2);plot(t,Ki,'linewidth',1.5);
ylabel('Ki');
subplot(3,1,3);plot(t,Kd,'linewidth',1.5);
ylabel('Kd');xlabel('时间');