在我的博文 PLC的PID控制器代码 中介绍了用于PLC的PID控制算法，既然代码出来了，那就要验证其正确性了，从本文开始，使用两篇文章介绍通过不同的方式对PID进行仿真，本来从《先进PID控制MATLAB仿真》这篇文章中，我们知道有三种仿真形式，但纯M代码的仿真会把PLC编写的PID代码部分打乱，为了尽量减少代码的改动，这里采用Simulink+M函数以及Simulink+S函数两种方式进行。

     本文讲述使用M函数的形式。

     1>输入参数的规划

    t:考虑到M函数需要初始化，这里通过传递仿真时间给M函数，通过判断其t==0来初始化变量

    SP:来自仿真系统给定的设定值

    PV:来自仿真系统反馈的值

    Flag:用于切换不同的输入函数波形

   2>输出参数的规划

   MV：这里只输出一个操作值

好了，M函数的编写很简单，直接上代码：

function [ MV ] = PID_PosExpMFun1( t,SP,PV,Flag)
   persistent ek1;      %上一次偏差
   persistent integral;   %积分量
   persistent lastMV; %上一次的输出
   persistent lastUd; %上一次的微分项值
   persistent lastPV; %上一次的PV值
   IsHSPID=0; %是否使用不完全微分
   IsDSpeed=0; %是否使用变速积分(=1,使用)
   IsDGama=0;  %是否使用微分先行算法(=1,采用)
   MVH=100;    %输出上限
   MVL=-100;      %输出下限
   Deadband=0.001;     %死区
   Ts=1;         %采样周期
   Tf=0.01;         %不完全微分滤波系数(a=Tf/(Ts+Tf))

   IUse=60;     %积分分离法中确定偏差值在所属范围内启用积分项
   IVar=60;     %变积分控制参数
   DGama=10000;    %微分先行系数
  
   switch Flag %选择调节参数,针对不同波形
       case 1,
            Kp=1.8;       %比例系数
            Ti=220;     %积分实现
            Td=50;         %微分时间
       case 2,
            Kp=1.5;Ti=280;Td=50;
            
        case 3,
            Kp=1.5;       %比例系数
            Ti=220;     %积分实现
            Td=50;         %微分时间
        case 4,
            Kp=1.5;       %比例系数
            Ti=220;     %积分实现
            Td=60;         %微分时间
            IsHSPID=1;
            Tf=6;
   end
   if t==0 %初始化函数
        ek1=0.0;
        integral=0.0;
        lastMV=0.0;
        lastUd=0.0;
        lastPV=0.0;
        MV=0.0;
        return;
   end
    ek=SP-PV; %当前偏差
    if ek < -Deadband   %死区
          ek = ek + Deadband ;
    elseif ek > Deadband
          ek = ek - Deadband ;
    else
          ek = 0.0 ;  
    end
    %计算比例项
     pterm=Kp*ek;
     %计算积分项               
      if abs(ek)
           %抗积分饱和:当上一次的输出达到饱和后,且偏差朝反方向变化,则可累加积分
            IsISum=(lastMVMVL+0.01)...
                 ||(lastMV>=MVH-0.01 && ek<=0.0)...
                 ||(lastMV<=MVL+0.01 && ek>=0.0);
            if IsISum %抗积分饱和
                 if IsDSpeed  %使用变速积分
                     integral=integral+Kp*Ts*ek*(1-abs(ek)/IVar)/Ti;
                 else
                     integral=integral+Kp*Ts*ek/Ti;
                 end
            end
       else
             integral=0.0;
      end   
       %微分项
     if IsHSPID  %使用不完全微分
          dterm =(ek-ek1)*Kp*Td/(Ts+Tf)+Tf*lastUd/(Ts+Tf);
     elseif IsDGama %采用微分先行算法
         dterm=lastUd*DGama*Td/(DGama*Td+Ts)+PV*(Td+Ts)/(DGama*Td+Ts)...
                        -lastPV*Td/(DGama*Td+Ts);
     else
          dterm=(ek-ek1)*Kp*Td/Ts;
     end
      %计算输出值
     MV=pterm+integral+dterm;
     if MV>MVH  %限幅
           MV=MVH;
     elseif  MV           MV=MVL;
     end
      ek1=ek;  %记录上一次的偏差
      lastMV=MV;%记录上一次的输出值
      lastUd=dterm;%储存上一次的微分项值
      lastPV=PV;%记录上一次的测量值
end

检查上述代码，大家应该看到，后半部分代码几乎与PLC中的PID代码一样(只是格式稍微有点区别而已),这样我们可以以最接近原PLC的方式进行模拟。

下面来搭建Simulink仿真平台，先上图：

http://s16/mw690/002wlscIzy6XsrrWBOfbf&690
通过上图中的Constant模块可以切换不同的输入波形，Transport Delay设置为延时40秒，仿真参数设定如下：

在示波器中设定Y轴上下限-200-200：

http://s14/mw690/002wlscIzy6XsrMfrqZ9d&690

最后设定一个正玄波的参数设置(其他波形的设置类似)

http://s1/mw690/002wlscIzy6Xss5Ye3K90&690
通过上述的配置，可以看到仿真效果了：

       问题说明：由上述代码我们可以发现，在M函数中，若要保持前一状态的结果，使用的是persistent关键字创建变量，它类似其他高级语言中的静态变量，可以保存变量结果，但它有一个缺点，那就是这种变量具有共享特性，当在同一个Simulink模型空间中多个位置对这个M函数进行调用时，由于静态变量的共享特性，不同位置的调用都会修改这些静态变量，从而导致不可预料的结果，虽然在一个空间中M函数被多次调用的场合不多，但也可能存在，此时要么再创建一个M函数的拷贝，换一个名称调用，或者使用S函数来实现，下一次博客会讲到S函数的仿真。