基于模糊PID的精准温控算法

最近根据实际的项目进行了模拟仿真验证，以下为仿真程序

clear all;
a=readfis('sanjiao');%读取模糊推理的参数，采三角形隶属度，重心法解模糊

%% 传递函数定义部分
ts=1;  %采样周期
s=tf('s');
sys=tf(0.2048,[200.9,1],'inputdelay',12.97);  %传递函数表达式
dsys=c2d(sys,ts,'z');  %将连续的时间模型转换成离散的时间模型
[num,den]=tfdata(dsys,'v');  %获得离散还建模型的分子分母矩阵

%% 模糊PID参数设定
u_1=0.0;u_2=0.0;u_3=0.0;u_5=0;u_4=0;u_6=0;u_7=0;u_8=0;u_9=0;u_10=0;u_11=0;u_12=0;u_13=0;u_14=0;  
y_1=15;  %初始环境温度设定
y_2=0;y_3=0;
x=[0,0,0]';  %误差/累加误差/误差积分初始化
error_1=0;  %上一时刻温差初始化
e_1=0.0;  %误差初始化
ec_1=0.0;  %误差变化率初始化

%% PID参数初值初始值
kp0=33;
ki0=0.0000001;
kd0=150;

%% 升温过程
for k=1:1:4000
time(k)=k*ts; 
rin(k)=20;  %设定温度值
%解模糊部分
k_pid=evalfis([e_1,ec_1],a);  %计算模糊推理输出KP.KI.KD参数增量
    kp(k)=kp0+k_pid(1); %Kp解模糊结果
    ki(k)=ki0+k_pid(2); %Ki解模糊结果
    kd(k)=kd0+k_pid(3); %Kd解模糊结果

yout(k)=den(1)*y_1+num(2)*u_13+num(1)*u_14;  %z变换后的离散对象

%输入参数及使用参数计算
error(k)=rin(k)-yout(k);  %输入温度差
x(1)=error(k);   %x为输入温度差
x(2)=error(k)-error_1;  %误差变化率计算
x(3)=x(3)+error(k);  %累计误差计算

%% 积分分离PID开关选择，有利于减小超调量
%M=1不使用积分分离，M=2使用分段式积分分离PID，M=3使用开关式积分分离
M=1;
if M==3  %开关式积分分离
    if abs(error(k)>=1)
        beta=0;
    else beta=0.9;
    end
end

if M==2;  %分段式积分分离PID
    if abs(error(k)>=30&abs(error(k)<=40))beta=0.000000001;
    elseif abs(error(k)>=20&abs(error(k)<=30))beta=0.000001;
        elseif abs(error(k)>=5&abs(error(k)<=20))beta=0.00002;
            elseif abs(error(k)>=0&abs(error(k)<=5))beta=0.00005;
    else beta=0;
    end
    
else if M==1;  %无积分分离PID
         beta=1;
    end
end

%% 实时输出功率计算
u(k)=kp(k)*x(1)+kd(k)*x(2)+ki(k)*x(3)*beta;%实时输出功率功率

%% 迭代过程        
u_14=u_13;
u_13=u_12;
u_12=u_11;
u_11=u_10;
u_10=u_9;
u_9=u_8;
u_8=u_7;
u_7=u_6;
u_6=u_5;
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);
e_1=x(1);
ec_1=x(2);
error_2=error_1;
error_1=error(k); 
end

%% 画图部分
showrule(a);
%time,rin,'b',
figure(1);plot(time,rin,'b',time,yout,'r');
xlabel('时间/s');ylabel( '实时温度曲线/℃');
legend('设定温度','模糊PID');
title('温度变化曲线图');
grid on;
hold on;

figure(2);plot(time,kp,'r');
xlabel('时间/s');ylabel( 'Kp值 ');
title('Kp变化曲线图');
grid on;

figure(3);plot(time,ki,'g');
xlabel('时间/s');ylabel( 'Ki值 ');
title('Ki变化曲线图')
grid on;

figure(4);plot(time,kd,'b');
xlabel('时间/s');ylabel( 'Kd值 ');
title('Kd变化曲线图')
grid on;

figure(5);plot(time,u,'m');
xlabel('时间/s');ylabel( '输出PWM ');
title('实时输出功率')
grid on;

figure(6);
plotfis(a);%绘制模糊推理系统的推理过程结构框图
fuzzy sanjiao.fis

%% 以下为响应性能指标评价
%误差泛函数指标，4项指标越小越好
k=1:1:1000;
IAE=trapz(k,abs(error(k)))         % IAE 
ISE=trapz(k,error(k).^2)           % ISE 
ITAE=trapz(k,k.*abs(error(k)))     % ITAE
ITSE=trapz(k,k.*(error(k).^2))     % ITSE

%% 传统指标计算
n=length(k);
[ymax,ind]=max(yout);  %y是原系统的阶跃响应，是一个二维矩阵，返回本身曲线的最大值和对应的时间矩阵变量的第几位
yss=yout(n);  %我们这里就设我们时间的最终值对应的就是稳态值
mp=(ymax-yss)/yss;  %由超调量定义可以计算超调量
max_temp=ymax-yss;

for j=1:n  %该层循环是计算上升时间即第一次达到稳态值的时间
if yout(j)<=yss&yout(j+1)>=yss  %第一次达到稳态值
tr=k(j+1);  %用tr存储稳态值
break
end
end

for i=n:-1:1  %该层循环用于计算调整时间，我们是从最终时间（无穷远的时间）往前开始遍历查找，否则从前开始查找的话会出错
    if yout(i)>=1.02*yss|yout(i)<=0.98*yss  %当得到的值与稳态值差的绝对值与稳态值之比为2%时说明到达了调整时间
        ts=k(i);  %ts存储调整时间
        break
    end
end

%% 误差性能指标输出
disp('峰值时间为:')  
disp(k(ind))
disp('上升时间为:')
disp(tr)
disp('调整时间为:')
disp(ts)
disp('超调量为:')
disp(mp)
disp('超调温度为:')
disp(max_temp)

%%以上为本人根据实际项目进行的仿真验证实验部分，如有错误欢迎指正！
%%如果想要整个项目文件的话可以留邮箱，有空看到会发的！欢迎交流
%%PS额最近也比较忙，我把全部资料放在下载了，包括论文和源码，额默认
好像系统自己设置的积分。链接如下所示：https://download.csdn.net/download/qq_39727487/12550181