pid算法的MATLAB仿真 – 用MATLAB进行pid算法仿真实验

算法理论知识

PID算法是工业应用中最广泛算法之一，在闭环系统的控制中，可自动对控制系统进行准确且迅速的校正。PID算法已经有100多年历史，在四轴飞行器，平衡小车、汽车定速巡航、温度控制器等场景均有应用。

PID算法：就是“比例（proportional）、积分（integral）、微分（derivative）”，是一种常见的“保持稳定”控制算法。

1、基础理论知识

常规的模拟PID控制系统原理框图如下所示：

因此可以得出e(t)和u(t)的关系：

其中：

Kp：比例增益，是调适参数；

Ki：积分增益，也是调适参数；

Kd：微分增益，也是调适参数；

e：误差=设定值（SP）- 回授值（PV）；

t：目前时间。

2、案例讲解

数学公式可能比较枯燥，通过以下例子，了解PID算法的应用。

例如，使用控制器使一锅水的温度保持在50℃，小于50℃就让它加热，大于50度就断电不就行了？

没错，在要求不高的情况下，确实可以这么干，如果换一种说法，你就知道问题出在哪里了。

如果控制对象是一辆汽车呢？要是希望汽车的车速保持在50km/h不动，这种方法就存在问题了。

设想一下，假如汽车的定速巡航电脑在某一时间测到车速是45km/h，它立刻命令发动机：加速！

结果，发动机那边突然来了个100%全油门，嗡的一下汽车急加速到了60km/h，这时电脑又发出命令：刹车！结果车子一会儿加速，一会儿又立马刹车减速，乘客吐了，肯定要骂人……

所以，在大多数场合中，用“开关量”来控制一个物理量就显得比较简单粗暴了，有时候是无法保持稳定的，因为单片机、传感器不是无限快的，采集、控制需要时间。

而且，控制对象具有惯性，比如将热水控制器拔掉，它的“余热”即热惯性可能还会使水温继续升高一小会。此时就需要使用PID控制算法了。

3、算法参数讲解

接着咱再来详细了解PID控制算法的三个最基本的参数：Kp比例增益、Ki积分增益、Kd微分增益。

（1）、Kp比例增益

Kp比例控制考虑当前误差，误差值和一个正值的常数Kp（表示比例）相乘。需要控制的量，比如水温，有它现在的当前值，也有我们期望的目标值。

当两者差距不大时，就让加热器“轻轻地”加热一下。

要是因为某些原因，温度降低了很多，就让加热器“稍稍用力”加热一下。

要是当前温度比目标温度低得多，就让加热器“开足马力”加热，尽快让水温到达目标附近。

这就是P的作用，跟开关控制方法相比，是不是“温文尔雅”了很多。

实际写程序时，就让偏差（目标减去当前）与调节装置的“调节力度”，建立一个一次函数的关系，就可以实现最基本的“比例”控制了~

Kp越大，调节作用越激进，Kp调小会让调节作用更保守。

若你正在制作一个平衡车，有了P的作用，你会发现，平衡车在平衡角度附近来回“狂抖”，比较难稳住。

（2）、Kd微分增益

Kd微分控制考虑将来误差，计算误差的一阶导，并和一个正值的常数Kd相乘。

有了P的作用，不难发现，只有P好像不能让平衡车站起来，水温也控制得晃晃悠悠，好像整个系统不是特别稳定，总是在“抖动”。

设想有一个弹簧：现在在平衡位置上，拉它一下，然后松手，这时它会震荡起来，因为阻力很小，它可能会震荡很长时间，才会重新停在平衡位置。

请想象一下：要是把上图所示的系统浸没在水里，同样拉它一下：这种情况下，重新停在平衡位置的时间就短得多。

此时需要一个控制作用，让被控制的物理量的“变化速度”趋于0，即类似于“阻尼”的作用。

因为，当比较接近目标时，P的控制作用就比较小了，越接近目标，P的作用越温柔，有很多内在的或者外部的因素，使控制量发生小范围的摆动。

D的作用就是让物理量的速度趋于0，只要什么时候，这个量具有了速度，D就向相反的方向用力，尽力刹住这个变化。

Kd参数越大，向速度相反方向刹车的力道就越强，如果是平衡小车，加上P和D两种控制作用，如果参数调节合适，它应该可以站起来了。

（3）、Ki积分增益

Ki积分控制考虑过去误差，将误差值过去一段时间和（误差和）乘以一个正值的常数Ki。

还是以热水为例，假如有个人把加热装置带到了非常冷的地方，开始烧水了，需要烧到50℃。

在P的作用下，水温慢慢升高，直到升高到45℃时，他发现了一个不好的事情：天气太冷，水散热的速度，和P控制的加热的速度相等了。

这可怎么办？

P兄这样想：我和目标已经很近了，只需要轻轻加热就可以了。

D兄这样想：加热和散热相等，温度没有波动，我好像不用调整什么。

于是，水温永远地停留在45℃，永远到不了50℃。

根据常识，我们知道，应该进一步增加加热的功率，可是增加多少该如何计算呢？

前辈科学家们想到的方法是真的巧妙，设置一个积分量，只要偏差存在，就不断地对偏差进行积分（累加），并反应在调节力度上。

这样一来，即使45℃和50℃相差不是太大，但是随着时间的推移，只要没达到目标温度，这个积分量就不断增加，系统就会慢慢意识到：还没有到达目标温度，该增加功率啦！

到了目标温度后，假设温度没有波动，积分值就不会再变动，这时，加热功率仍然等于散热功率，但是，温度是稳稳的50℃。

Ki的值越大，积分时乘的系数就越大，积分效果越明显，所以，I的作用就是，减小静态情况下的误差，让受控物理量尽可能接近目标值。

I在使用时还有个问题：需要设定积分限制，防止在刚开始加热时，就把积分量积得太大，难以控制。

4、PID算法参数调试

PID算法的参数调试是指通过调整控制参数（比例增益、积分增益/时间、微分增益/时间）让系统达到最佳的控制效果。

调试中稳定性（不会有发散性的震荡）是首要条件，此外，不同系统有不同的行为，不同的应用其需求也不同，而且这些需求还可能会互相冲突。

PID算法只有三个参数，在原理上容易说明，但PID算法参数调试是一个困难的工作，因为要符合一些特别的判据，而且PID控制有其限制存在。

（1）、稳定性

若PID算法控制器的参数未挑选妥当，其控制器输出可能是不稳定的，也就是其输出发散，过程中可能有震荡，也可能没有震荡，且其输出只受饱和或是机械损坏等原因所限制。不稳定一般是因为过大增益造成，特别是针对延迟时间很长的系统。

（2）、最佳性能

PID控制器的最佳性能可能和针对过程变化或是设定值变化有关，也会随应用而不同。

两个基本的需求是调整能力（regulation，干扰拒绝，使系统维持在设定值）及命令追随（设定值变化下，控制器输出追随设定值的反应速度）。有关命令追随的一些判据包括有上升时间及整定时间。有些应用可能因为安全考量，不允许输出超过设定值，也有些应用要求在到达设定值过程中的能量可以最小化。

（3）、各调试方法对比

方法	优点	缺点
人工调试	不需要数学，可以在线调试	需要有经验的工程师
齐格勒-尼克尔斯方法	被证实有效的方法，可以在线调试	会影响制程，需要试误，得到的参数可能使响应太快
软件工具	调试的一致性，可以在线调试或离线调试，可以配合计算机自动设计，包括阀及感测器的分析，可以下载前进行模拟，可以致远非稳态（NSS）的调试	需要成本或是训练
Cohen-Coon	好的程序模型	需要一些数学，需要离线调试，只对一阶系统有良好效果

（4）、调整PID参数对系统的影响

一、实验要求

学习PID控制器的基本原理，掌握PID参数的物理调节规律，通过仿真验证PID调节过程。选择例题中的任意两个实验进行仿真研究，改变参数后观察系统的动态性能和稳态性能的变化，确认是否和理论分析一致。

二、实验内容

1、实验1：具有变化积分率的PID控制器

A、代码

% PID Controller with changing integration rate
% 具有变化积分率的PID控制器

clear all;
close all;
clc

% Big time delay Plant 大延时装置

ts=20;
sys=tf([1],[60,1],'inputdelay',80);
dsys=c2d(sys,ts,'zoh');
[num,den]=tfdata(dsys,'v');
u_1=0;
u_2=0;
u_3=0;
u_4=0;
u_5=0;

y_1=0;y_2=0;y_3=0;
error_1=0;error_2=0;
ei=0;
for k = 1:1:200
    time(k)=k*ts;
    rin(k)=1.0;%Step Signal
   
    %Linear model
    yout(k)=-den(2)*y_1+num(2)*u_5;
    error(k)=rin(k)-yout(k);
    kp=0.45;kd=12;ki=0.0048;
    A=0.4;B=0.6;
    
    %T type integration
    ei=ei+(error(k)+error_1)/2*ts;
    M=1;
    if M==1%Changing integration rate
        if abs(error(k))<=B
            f(k)=1;
        elseif abs(error(k))>B&abs(error(k))<=A+B
            f(k)=(A-abs(error(k))+B)/A;
        else
            f(k)=0; 
        end
        
    elseif M==2 %Not changing integration rate
        f(k)=1; 
    end
    
    u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+ki*f(k)*ei;
    
    if u(k)>=10
        u(k)=10;
    end    
    if u(k)<=-10
        u(k)=-10;
    end
    
    %Return of PID parameters
    u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);
    y_3=y_2;y_2=y_1;y_1=yout(k);
    error_2=error_1;
    error_1=error(k);
    
end

figure(1);
plot(time,rin,'b',time,yout,'r');
xlabel('time(s)');ylabel('rin,yout');

figure(2);
plot(time,f,'r');
xlabel('time(s)');ylabel('Integration rate f');

B、实验结果

（1）默认效果

在没有更改代码的情况下，实验的效果如图一所示。

图一

（2）改变kp

①增大kp,理论上可以看到pid调节的更快更陡。实际效果如图2。Kp=3.

    kp=3;kd=12;ki=0.0048;

图2

实际上产生了自激震荡，说明kp不能调的过大，否则系统会不稳定。

②减小kp，理论上系统调节缓慢，会产生超调量，调节曲线不平滑，有毛刺。Kp=0.1

    kp=0.1;kd=12;ki=0.01;

图3

实际上也证明比例系数太小，系统变化时调节相对缓慢。

（3）改变kd

①增大kd，理论上系统稳定时波动很大，精度不高，且有很多毛刺与抖动。Kd=30.

kp=0.45;kd=30;ki=0.0048;

图4

实际上系统开始调节时会有些尖峰，调节曲线不平滑。

②减小kd，理论上系统调节稳定时变化缓慢，会有些超调量。Kd=3.

kp=0.45;kd=3;ki=0.0048;

图5

实际上也可以看出减小kd对系统影响不大，只是在系统刚进入稳定状态时有一个小小的超调量。

（4）改变ki

①增大ki,积分变大，系统系统变得不稳定，抖动波动很厉害。

kp=0.45;kd=3;ki=0.01;

图6

实际上可以其变化幅度大，超调量很高，系统很不稳定。

②减小ki,系统调节变慢，系统可能无法调节。

kp=0.45;kd=12;ki=0.001;

图7

系统无法达到设定的稳态值，无法调节。

2、实验2：串联系统控制

A、代码

%Series System Control 串联系统控制
clear all;
close all;
clc

ts=2;
sys1=tf(1,[10,1]);
dsys1=c2d(sys1,ts,'z');
[num1,den1]=tfdata(dsys1,'v');
sys2=tf(1,[10,1]);
dsys2=c2d(sys2,ts,'z');
[num2,den2]=tfdata(dsys2,'v');
dph=1/zpk('z',ts);
Gc2=dph/(dsys2*(1-dph));
[nump,denp]=tfdata(Gc2,'v');
u1_1=0.0;u2_1=0.0;
y1_1=0;y2_1=0;
e2_1=0;ei=0;
for k=1:1:2000
    time(k)=k*ts;
    r1(k)=1;
    %Linear model
    y1(k)=-den1(2)*y1_1+num1(2)*y2_1;%Main plant
    y2(k)=-den2(2)*y2_1+num2(2)*u2_1;%Assistant plant
    error(k)=r1(k)-y1(k);
    ei=ei+error(k);
    
    u1(k)=1.2*error(k)+0.02*ei;%Main Controller
    e2(k)=u1(k)-y2(k);%Assistant Controller
    u2(k)=-denp(2)*u2_1+nump(1)*e2(k)+nump(2)*e2_1;
    d2(k)=0.01*rands(1);
    u2(k)=u2(k)+d2(k);
    %----------Return of PID parameters------------
    u1_1=u1(k);
    u2_1=u2(k);
    e2_1=e2(k);
    y1_1=y1(k);
    y2_1=y2(k);
end

figure(1);%Assistant Control
plot(time,u1,'b',time,y2,'r');
xlabel('time(s)');ylabel('u1,y2');

figure(2);%Main Control
plot(time,r1,'b',time,y1,'r');
xlabel('time(s)');ylabel('r1,y1');

figure(3);
plot(time,d2,'r');
xlabel('time(s)');ylabel('disturbance');