灰狼算法整定PID参数，传递函数可以任意修改，附MATLAB代码

今天的主题是：采用灰狼算法整定各种常见传递函数的PID参数。

PID参数整定的方法有很多，其中包括正交试验法、极点配置法，另一类为智能优化算法，包括遗传算法、粒子群算法、神经网络算法。今天给大家带来一期，采用灰狼算法整定PID参数，并比较整定前后的时域指标。

传递函数包含以下五种，基本上包含了所有常见的传递函数。

测试在阶跃响应下，采用灰狼算法整定PID的结果图如下：

01

取T = 50; K = 5; n = 3，（这几个参数在代码中可自行修改）

GWO优化前后结果如下：

控制台输入结果如下：

02

取Tz = [40,20,60,100,20];  K = 6;  n = 5，（这几个参数在代码中可自行修改）

GWO优化前后结果如下：

控制台输入结果如下：

03

取T = 15; K = 10; n= 5;  tt = 5; （这几个参数在代码中可自行修改）

GWO优化前后结果如下：

控制台输入结果如下：

04

取T1 = 10; K1 = 2; n1=3;  tt1 = 0.02;  T2 = 2; K2 = 0.1; n2=4;  tt2 = 0.001; （这几个参数在代码中可自行修改）

GWO优化前后结果如下：

控制台输入结果如下：

05

取K = 10;T = 50; n = 4;  tt = 5; （这几个参数在代码中可自行修改）

GWO优化前后结果如下：

控制台输入结果如下：

代码：

clear
clc
%% 五种适应度函数分别 对应图片1-5
flag = 1;   %1-5分别对应图片1-5，图片中对应的参数请自行修改即可
KP0 = 0.1;  %优化前先手动给一个P , I , D指标，方便与优化后做对比，这三个参数也可以自行更改
Ti0 = 100;
Kd0 = 0.5;
if flag == 1
    T = 50; K = 5; n = 3;                           %这几个参数可以自行修改
    PI_1;                                           %PI_1-5分别对应5个传递函数。采用手动给的PID指标算响应曲线
elseif flag==2
    Tz = [40,20,60,100,20];  K = 6;  n = 5;         %这里注意，n和Tz变量的长度一致，这几个参数可以自行修改
    PI_2;                                           %采用手动给的PID指标算响应曲线
elseif flag==3
    T = 15; K = 10; n= 5;  tt = 5;                  %这几个参数可以自行修改
    PI_3;                                           %采用手动给的PID指标算响应曲线
elseif flag==4
    T1 = 10; K1 = 2; n1=3;  tt1 = 0.02;                %这几个参数可以自行修改
    T2 = 2; K2 = 0.1; n2=4;  tt2 = 0.001;  
    PI_4;                                           %采用手动给的PID指标算响应曲线
elseif flag==5
    K = 10;T = 50; n = 4;  tt = 5;                  %这几个参数可以自行修改
    PI_5;                                           %采用手动给的PID指标算响应曲线
else
    T = 50; K = 5; n = 3;                           %这几个参数可以自行修改 
    PI_1;                                           %采用手动给的PID指标算响应曲线
end


% Grey Wolf Optimizer
SearchAgents_no = 50;
dim = 3;  %优化参数个数
ub = [100,200,10];%优化变量的上限
lb = [0.01,20,0];%优化变量的下限
Max_iter = 100;
% initialize alpha, beta, and delta_pos
Alpha_pos=zeros(1,dim);
Alpha_score=inf; %change this to -inf for maximization problems
Beta_pos=zeros(1,dim);
Beta_score=inf; %change this to -inf for maximization problems
Delta_pos=zeros(1,dim);
Delta_score=inf; %change this to -inf for maximization problems


%Initialize the positions of search agents
Positions=initialization(SearchAgents_no,dim,ub,lb);


Convergence_curve=zeros(1,Max_iter);


l=0;% Loop counter


% Main loop
while l<Max_iter
    for i=1:size(Positions,1)  
        
       % Return back the search agents that go beyond the boundaries of the search space
        Flag4ub=Positions(i,:)>ub;
        Flag4lb=Positions(i,:)<lb;
        Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
            
        % Update Alpha, Beta, and Delta
        if fitness<Alpha_score 
            Alpha_score=fitness; % Update alpha
            Alpha_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness<Beta_score 
            Beta_score=fitness; % Update beta
            Beta_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score 
            Delta_score=fitness; % Update delta
            Delta_pos=Positions(i,:);
        end
    end
    
    
    a=2-l*((2)/Max_iter); % a decreases linearly fron 2 to 0
    
    % Update the Position of search agents including omegas
    for i=1:size(Positions,1)
        for j=1:size(Positions,2)     
                       
            r1=rand(); % r1 is a random number in [0,1]
            r2=rand(); % r2 is a random number in [0,1]
            
            A1=2*a*r1-a; % Equation (3.3)
            C1=2*r2; % Equation (3.4)
            
            D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1
            X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
                       
            r1=rand();
            r2=rand();
            
            A2=2*a*r1-a; % Equation (3.3)
            C2=2*r2; % Equation (3.4)
            
            D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); % Equation (3.5)-part 2
            X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2       
            
            r1=rand();
            r2=rand(); 
            
            A3=2*a*r1-a; % Equation (3.3)
            C3=2*r2; % Equation (3.4)
            
            D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); % Equation (3.5)-part 3
            X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3             
            
            Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)
            
        end
    end
    l=l+1;    
    Convergence_curve(l)=Alpha_score;
end


Best_pos = Alpha_pos;


if flag==1
    [Qbg,y_gwo,u,t]=Object1(Best_pos,T,K,n);
elseif flag==2
    [Qbg,y_gwo,u,t]=Object2(Best_pos,Tz,K,n);
elseif flag==3
    [Qbg,y_gwo,u,t]=Object3(Best_pos,T,K,n,tt);
elseif flag==4
    [Qbg,y_gwo,u,t]=Object4(Best_pos,T1,K1,tt1,n1,T2,K2,tt2,n2);
elseif flag==5
    [Qbg,y_gwo,u,t]=Object5(Best_pos,T,K,tt,n);
else
    [Qbg,y_gwo,u,t]=Object1(Best_pos,T,K,n);
end


%计算优化后的时域性能指标
[Overshoot,tp,tr,ts]=performance_index(y_gwo,t);
disp(['优化后的阶跃响应曲线参数：   超调量：',num2str(Overshoot),'，峰值时间：',num2str(tp),'，上升时间：',num2str(tr),'，调节时间：',num2str(ts)])
disp(['优化后的PID参数为：',num2str(Best_pos)])
[Overshoot0,tp0,tr0,ts0]=performance_index(y0,t0);
disp(['优化后的阶跃响应曲线参数：   超调量：',num2str(Overshoot0),'，峰值时间：',num2str(tp0),'，上升时间：',num2str(tr0),'，调节时间：',num2str(ts0)])
disp(['优化前的PID参数为：',num2str([KP0,Ti0,Kd0])])


%% 作图
figure
plot(t0,y0,"Color",'b','LineWidth',1);
hold on
plot(t,y_gwo,"Color",'r','LineWidth',1);
hold on
plot(t,ones(length(t),1),'LineStyle','--',"Color",'k','LineWidth',1);
legend('优化前','GWO优化后');
xlabel('时间/s');
ylabel('输出');

观察代码第4行，可以通过修改参数flag，任意切换五种传递函数，除此之外，大家可以自行将传递函数改成自己想要的任意参数。

代码注释很详细。获取完整代码方式：点击下方卡片回复关键词获取。

关键词：GWOPID