Example of Designing a PID Controller Using Particle Swarm Optimization in MATLAB

First, fill in the parameters for the particle swarm method, such as the number of particles, inertia factor, learning factors, etc. Then, input the coefficients of the numerator and denominator of the system’s forward transfer function, determine the cost function, and set the number of iterations. Next, loop through the calculations to obtain the optimal PID parameters. These parameters are then used to construct the closed-loop transfer function, which is tested with a step signal to evaluate the system’s output results.

%% Particle Swarm Optimization Method Design PID Controller %%

G(s)=(s+3)/(s^3+7*s^2+8*s+1.2);

clear all;close all;clc;

% Parameter settings
param1=0.5;  % Inertia weight
param2=1.8;  % Cognitive learning factor
param3=1.6;  % Social learning factor
Dim1=3;  % Number of parameters
swarmNum=100;  % Number of particles
minFit1=0.1;  % Minimum fitness value
LoopMax=999;  % Maximum iterations
maxV1=1.2;  % Maximum velocity
minV1=-1.2;  % Minimum velocity
uper1=[600, 20, 20];  % Upper bounds
Low1=[0.01, 0.01, 0.01];  % Lower bounds
range1=ones(swarmNum, 1)*(uper1-Low1);  % Range of parameters

% Initialize particle swarm
swarmGrp=rand(swarmNum,Dim1).*range1+ones(swarmNum,1)*Low1;
swarmGrp(1,:)=[3.0, 0.2, 0.1];  % Initialize velocity
stepV1=rand(swarmNum,Dim1)*(maxV1-minV1)+minV1;
fitSwarm1=zeros(swarmNum,1);  % Fitness values

% Evaluate fitness
for k1=1:1:swarmNum
    fitSwarm1(k1)=feval(@myFun1,swarmGrp(k1,:));
end
[bestFit,bestID]=min(fitSwarm1);
totBestLizi=swarmGrp(bestID,:);  % Best particle

% Store PID parameters
fitFind=zeros(1,LoopMax);
Kp1=zeros(1,LoopMax);
Ki1=zeros(1,LoopMax);
Kd1=zeros(1,LoopMax);

for Loop=1:1:LoopMax
    for k1=1:1:swarmNum
        stepV1(k1,:)=param1*stepV1(k1,:)+param2*rand*(gloBestGrp(k1,:)-swarmGrp(k1,:))+...
            param3*rand*(totBestLizi-swarmGrp(k1,:));
        if stepV1(k1,:) > maxV1
            stepV1(k1,:) = maxV1;
        end
        if stepV1(k1,:) < minV1
            stepV1(k1,:) = minV1;
        end
        % Update position
        swarmGrp(k1,:)=swarmGrp(k1,:)+stepV1(k1,:);
        for k2=1:1:Dim1
            if swarmGrp(k1,k2) > uper1(k2)
                swarmGrp(k1,k2) = uper1(k2);
            end
            if swarmGrp(k1,k2) < Low1(k2)
                swarmGrp(k1,k2) = Low1(k2);
            end
        end  % for k2=1:1:Dim1
        fitSwarm1(k1)=feval(@myFun1, swarmGrp(k1,:));
        % Update best particle
        if fitSwarm1(k1) < fitBest(k1)
            fitBest(k1)=fitSwarm1(k1);
            gloBestGrp(k1,:)=swarmGrp(k1,:);
        end
        % Update global best
        if fitSwarm1(k1) < totFitBest
            totFitBest=fitSwarm1(k1);
            totBestLizi=swarmGrp(k1,:);
        end
    end  % for k1=1:1:swarmNum
    fitFind(Loop)=totFitBest;
    Kp1(Loop)=totBestLizi(1);
    Ki1(Loop)=totBestLizi(2);
    Kd1(Loop)=totBestLizi(3);
end  % for Loop=1:1:LoopMax

figure(1);plot(fitFind);title('Fitness');
figure(2);plot(Kp1);title('Kp');
figure(3);plot(Ki1);title('Ki');
figure(4);plot(Kd1);title('Kd');

Kp=Kp1(end);
Ki=Ki1(end);
Kd=Kd1(end);
num1=[Kp, Ki, Kd];
den1=[1, 0];
sys1=tf(num1,den1);
num2=[1,3];
den2=[1, 7, 8, 1.2];
sys2=tf(num2,den2);
sys3=series(sys1,sys2);
sys4=feedback(sys3,1,-1);
den4=sys4.Denominator;
den5=den4{1};
disp('Characteristic equation roots:');
roots(den5);
figure(5);step(sys4,30);title('Step Response');
disp('System');

function vaL1=myFun1(PID1)
    num1=[PID1(3), PID1(1), PID1(2)];
    den1=[1, 0];
    sys1=tf(num1,den1);  % G(s)=(s+3)/(s^3+7*s^2+8*s+1.2);
    num2=[1,3];
    den2=[1, 7, 8, 1.2];
    sys2=tf(num2,den2);
    sys3=series(sys1,sys2);
    sys4=feedback(sys3,1,-1);
    t1=0:0.01:20;
    Len=length(t1);
    u1=ones(1,Len);
    [outY,t2]=lsim(sys4,u1,t1);
    vaL1=0;
    for k1=1:1:Len
        vaL1=vaL1+t2(k1)*abs(outY(k1)-1);
    end
end