Self-Tuning PID Parameters with a Genetic Algorithm (Simulink Implementation)

In real-world engineering, PID parameters are mostly tuned by hand, which is time-consuming, labor-intensive, and heavily dependent on experience. Today I’ll go over the idea and implementation of automatically tuning PID parameters with a genetic algorithm (GA), carried out mainly on the MATLAB/Simulink platform.

A caveat up front: tuning PID with a genetic algorithm isn’t really practical in real engineering—the computational load is high, and an MCU often can’t handle it—but it still has some value as a research and learning method.

Building the Simulink Simulation Platform

First, build the Simulink model for the PID simulation.

Compared to an ordinary PID simulation model, the main difference is the introduction of a GA function to optimize the PID parameters. This requires an output port to pass the simulation results out, and an input port to bring in the population produced by the GA function. Both input and output are linked to the workspace, implemented with the simout and simin blocks respectively; the simin parameter also needs to be split into three channels by a bus block, connected to the P, I, and D ports. The constant block in the model (value 100) serves as the filter coefficient.

Designing the Optimization Objective Function

Drawing on the literature, the following optimization objective function is chosen:

Here e(t) is the system error, u(t) is the controller output (in a step simulation, this is the unit step 1(t)), tr is the rise time, and w1, w2, w3 are weights. To suppress overshoot, a penalty mechanism is introduced: as soon as overshoot occurs, the overshoot amount is counted separately as a term in the objective function.

The MATLAB function corresponding to the objective function is as follows; the simulation uses the parameters w1=1.2, w2=0.001, w3=4.0, w4=100:

function aim = fcn(u)

time=(0:0.001:3);  % Time series: sampling time 0.001 s, simulation time 3 s

output=u;

w1=1.2;
w2=0.001;
w3=4.0;
w4=100;

% Get the rise time
i1=1;
while(output(i1)<0.1)
    i1=i1+1;
end
t1=time(i1);

i2=1;
while(output(i2)<0.9)
    i2=i2+1;
end
t2=time(i2);

tr=t2-t1;

% Get the integral terms
int1=ones(1,3001);
int2=ones(1,3001);
int3=ones(1,3001);

for i=(1:3001)
    if output(i)<1
        int1(i)=0.001*abs(output(i)-1);
        int2(i)=0.001;
        int3(i)=0;
    else
        int1(i)=0;
        int2(i)=0.001;
        int3(i)=0.001*abs(output(i)-1);
    end
end

% Compute the objective function value
aim = w1*sum(int1)+w2*sum(int2)+w4*sum(int3)+w3*tr;

The rise time tr is computed using the standard 10%-to-90% definition. The integral part is handled in two segments depending on whether overshoot occurs: while the output has not exceeded the target value, only the integral of the absolute error is counted; once overshoot occurs, the overshoot amount is counted separately and multiplied by the larger weight w4=100, acting as a penalty.

The Interface Function Between the GA and Simulink

With the objective function in place, we also need an overall GA interface function whose input is the population of kp, ki, and kd, and whose output is the aim value corresponding to each individual. This function is essentially an implicit function: internally it calls the Simulink model to run the simulation, then passes the simulation result into fcn() to convert it into an aim value:

function z=ga(kp,ki,kd)

t=[0:0.001:3]';

% Define the simin struct and declare it as a global variable
field1='time';value1=ones(3001,1);
field3='values';value3=ones(3001,3);
field4='dimensions';value4=3;
field2='signals';value2=struct(field3,value3,field4,value4);

global simin;
simin=struct(field1,value1,field2,value2);

% Simulate each of the 40 individuals one by one and compute the corresponding aim value
for i=1:40
    gain=[kp(i),ki(i),kd(i)];
    gain_cal = [kp(i)*ones(3001,1),ki(i)*ones(3001,1),kd(i)*ones(3001,1)];

    simin.time = t;
    simin.signals.values = gain_cal;
    simin.signals.dimensions = 3;

    sim pid_test;

    z(i,1)=fcn(simout);
end

Each individual has to run a Simulink simulation on its own, so 40 individuals means 40 runs—this is the main reason it’s hard to use online in real engineering.

The Main Function: Genetic Algorithm Iteration

Finally, the main function, which handles population initialization, the evolutionary iterations (selection, crossover, mutation, reinsertion), and visualization of the results:

global simin;

global NIND;  % Define the population size as NIND; initial population of 40 individuals
NIND=40;

NVAR=3;      % Number of variable dimensions is 3
MAXGEN=100;  % Maximum number of generations is 100
PRECI=30;    % Number of binary bits in the gene encoding
GGAP=1;      % Generation gap of 1; total population size stays unchanged

trace=zeros(2,MAXGEN);  % Initialize the algorithm performance tracking matrix

Chrom=crtbp(NIND,PRECI*NVAR);  % Initialize the population

FieldD=[PRECI PRECI PRECI;10 0 0;100 10 10;1 1 1;0 0 0;1 1 1;1 1 1];
% Population evolution domain; here 10 0 0;100 10 10 means the lower and upper bounds of the three variables are [10,100], [0,10], [0,10] respectively
% The last four rows are: 1 - standard binary encoding, 0 - standard scaling, 1 - lower bound closed interval, 1 - upper bound closed interval

gen=0;  % Initial generation is 0

variable=bs2rv(Chrom,FieldD);  % Convert the binary population into decimal variables
ObjV=ga(variable(:,1),variable(:,2),variable(:,3));  % Compute the aim values of the initial population

while gen<MAXGEN

    FitnV=ranking(ObjV);                    % Fitness ranking; for a minimization problem use ObjV directly
    SelCh=select('sus',Chrom,FitnV,GGAP);  % Selection operation
    SelCh=recombin('xovsp',SelCh,0.7);     % Crossover operation
    SelCh=mut(SelCh);                       % Mutation operation

    variable=bs2rv(SelCh,FieldD);           % Convert binary to decimal variables
    ObjVSel=ga(variable(:,1),variable(:,2),variable(:,3));  % Compute the aim values after evolution

    [Chrom, ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);    % Reinsertion operation

    gen=gen+1;

    [Y,I]=min(ObjV);hold on;  % Find the minimum value in the evolved population

    trace(1,gen)=min(ObjV);               % Put the minimum value into the performance tracking matrix
    trace(2,gen)=sum(ObjV)/length(ObjV);  % Put the population's average aim value into the tracking matrix

end

variable=bs2rv(Chrom,FieldD);

figure(1);
plot(variable(:,1),'o');grid

figure(2);
plot(trace(1,:));
hold on;
plot(trace(2,:),'-.');grid
legend('Change in the solution','Change in the population mean')  % Visualize the change in the aim value and the change in the population mean

The search range set in FieldD is kp∈[10,100], ki∈[0,10], kd∈[0,10], which can be adjusted according to the specific plant. After the whole process finishes, the trace matrix records the optimum and the population mean for each generation, and the plot gives an intuitive view of how the algorithm converges.

Overall, the implementation approach is fairly clear: Simulink handles the simulation, the GA on the MATLAB side handles the optimization, and the two exchange data through the global variables simin/simout. The main limitation is that every fitness evaluation requires a full simulation run, which is slow; if needed later, one could consider replacing the simulation with an analytical model to speed up the computation.