PID controller optimization in simulink using parallel simulations in objective function

I am tunning a PID controller in simulink using optimization tool (Genetic Algorithim). It is working fine but iterative simulations are in serial and taking too much time. I saw an answer to run simulations in parallel by using Simulink.SimulationInput.

I tried this solution and my simulation is running fast but the optimization values are not correct ,also the Fitness value graph is abnormal. I guess the problem is in the new objetive function for parallel simulations " cost(i2,1)=Simout.ITAE(end)". I would appriciate if an expert like you could revise the objective fucntion as i am stucked on it for many days. I am attaching my objective function.

Please note that simulation output block name ITAE (to workspace) , higlighted in green, is the cost function

Matlab version: R2021b

Simulink Model and PID controller varriables

Live Editor to optimize the PID gains

% Set nondefault solver options

options2 = optimoptions("ga","UseVectorized",true,"CreationFcn",...

"gacreationuniform","SelectionFcn","selectiontournament","MutationFcn",...

"mutationadaptfeasible","Display","iter","PlotFcn","gaplotbestf");

% Solve

[solution2,objectiveValue] = ga(@tunning,nVariables,[],[],[],[],...

zeros(nVariables,1),repmat(5,nVariables,1),[],[],options2);

The serial objective function which is Working fine but slow is as follows:

function cost = tunning(k)

assignin('base','k',k);

sim('test.slx');

cost= ITAE(length(ITAE));

end

Parallel simulation Objective function providing Wrong output is as follows:

Answer 1

[Sam Chak]

Hi @Irfan Khan

I'm unfamiliar with the Simulink approach, but in MATLAB, the script to find the PID gains subject to ITAE criterion looks like this.

Due to time limit to run the simulation, smaller population size, bounds and shorter simulation duration are used.

fitfun  = @costfun;         % call cost function
PopSize = 3;                % Population size
MaxGen  = 6;                % Maximum generation
nvars   = 3;                % number of design variables: P,I,D
A       = -eye(nvars);      % linear constraints subject to A*k <= b 
b       = zeros(nvars, 1);  % (to search in the positive direction)
Aeq     = [];               
beq     = [];               % linear constraints subject to Aeq*k = beq
lb      = [30.1 51.3 13.7]; % lower bounds for lb < k < ub
ub      = [30.3 51.5 13.9]; % upper bounds for lb < k < ub
nonlcon = [];               % nonlinear constraints subject to C(k) <= 0 and Ceq(k) = 0

[k, fval, exitflag, output] = ga(fitfun, nvars, A, b, Aeq, beq, lb, ub, nonlcon)

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

k = 1×3

30.1343 51.3026 13.7872

fval = 0.0353

exitflag = 1

output = struct with fields:
      problemtype: 'linearconstraints'
         rngstate: [1×1 struct]
      generations: 51
        funccount: 2450
          message: 'Optimization terminated: average change in the fitness value less than options.FunctionTolerance.'
    maxconstraint: 0
       hybridflag: []

% Objective function to be optimized
function J = costfun(k)

    % Transfer functions
    G1 = tf(1, [0.2 1]);
    G2 = tf(1, [0.5 1]);
    G3 = tf(1, [10 0.8]);
    Gp = G1*G2*G3;
    H  = 20;
    
    % Compensator
    kp = k(1);                 % Proportional gain
    ki = k(2);                 % Integral gain
    kd = k(3);                 % Derivative gain
    Tf = 1/165.277936355792;   % 1st-order Derivative filter coefficient (1/N)
    Gc = pid(kp, ki, kd, Tf);
    
    % Closed-loop system (disturbance path)
    Gd = minreal(feedback(G3, H*Gc*G1*G2));
    
    % Cost function
    dt     = 0.01;
    ty     = 0:dt:8;
    [y, t] = step(Gd, ty);
    itae   = t.*abs(0 - H*y(:));
    J      = trapz(t, itae);
    
end

[Sam Chak]

Using the PID gains from the GA search, the zero steady-state output due to the step disturbance input is obtained.

G1 = tf(1, [0.2 1]);
G2 = tf(1, [0.5 1]);
G3 = tf(1, [10 0.8]);
Gp = G1*G2*G3

Gp =
 
                1
  ------------------------------
  s^3 + 7.08 s^2 + 10.56 s + 0.8
 
Continuous-time transfer function.

H  = 20;

% Compensator
kp = 30.1343;               % P gain from GA
ki = 51.3026;               % I gain from GA
kd = 13.7872;               % D gain from GA
Tf = 1/165.277936355792;    % 1/N (fixed)
Gc = pid(kp, ki, kd, Tf)

Gc =
 
             1            s    
  Kp + Ki * --- + Kd * --------
             s          Tf*s+1 

  with Kp = 30.1, Ki = 51.3, Kd = 13.8, Tf = 0.00605
 
Continuous-time PIDF controller in parallel form.

% Step Disturbance Response
Gd = minreal(feedback(G3, H*Gc*G1*G2))

Gd =
 
              0.1 s^4 + 17.23 s^3 + 116.7 s^2 + 165.3 s
  -----------------------------------------------------------------
  s^5 + 172.4 s^4 + 1181 s^3 + 4.792e04 s^2 + 1.008e05 s + 1.696e05
 
Continuous-time transfer function.

step(Gd, 10);

% Cost function
dt     = 0.001;
ty     = 0:dt:10;
[y, t] = step(Gd, ty);
itae   = t.*abs(0 - H*y(:));
J      = trapz(t, itae)

J = 0.0354

[Farshid R]

Hi @Sam Chak

Could I require a question about optimization?

Many thanks

[Sam Chak]

Hi @Farshid R

Yes, you can. And you are encouraged to post a new question, if it is not closely or directly related the OP's question or problem.

[Farshid R]

Thank you, I created a new question, I would be happy if you could help me. Here is the link:

https://www.mathworks.com/matlabcentral/answers/1812615-optimization-with-fmincon-command-in-simulink