[Simulink] NSGA-2 optimization algorithm tuning PID controller parameters -- first-order controlled object with time delay

Posted by Styles2304 on Fri, 04 Mar 2022 06:17:30 +0100

0 research background

Write before:
 1. This code is based on MATLAB2019a. Errors may be reported in the lower version or different versions, and the opening of mdl file or slx file may fail;
 2. If the running time is too long, observe whether the number of iterations changes.
 3. This blog is attached with the code and detailed introduction. If you reprint it, please indicate the source;
 4. If this blog happens to be related to your research, you are welcome to consult qq1366196286;

Reference blog
[Simulink] PSO optimization algorithm for tuning PID controller parameters (I)
[Simulink] GA optimization algorithm for tuning PID controller parameters (III) -- first-order controlled object with time delay

   continue to explain this blog: [Simulink] NSGA-2 optimization algorithm tuning PID controller parameters (IV) -- first-order controlled object with time delay

1. Brief introduction of multi-objective optimization algorithm

   for the basic concept, please refer to the literature:
[1] Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2):182-197.
[2] Zhang Q, Li H. MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[J]. IEEE Transactions on Evolutionary Computation, 2007, 11(6):712-731.

1.1 brief introduction of nsga-2 algorithm

   NSGA - Ⅱ is one of the most popular multi-objective genetic algorithms at present. It reduces the complexity of non inferior sorting genetic algorithm. It has the advantages of fast running speed and good convergence of solution set. It has become the benchmark of the performance of other multi-objective optimization algorithms.
   NSGA - Ⅱ is improved on the basis of the first generation non dominated sorting genetic algorithm. Its improvement is mainly aimed at the three aspects mentioned above:
① A fast non dominated sorting algorithm is proposed. On the one hand, it reduces the computational complexity. On the other hand, it combines the parent population with the child population, so that the next generation population can be selected from double space, so as to retain all the best individuals;
② The elite strategy is introduced to ensure that some excellent population individuals will not be discarded in the process of evolution, so as to improve the accuracy of optimization results;
③ The use of crowding degree and crowding degree comparison operator not only overcomes the defect of manually specifying shared parameters in NSGA, but also takes it as the comparison standard between individuals in the population, so that the individuals in the quasi Pareto domain can be evenly extended to the whole Pareto domain, and the diversity of the population is guaranteed.

Figure 4-1 calculation of congestion degree and congestion distance

Algorithm purpose:
For the current m individuals, select N individuals (M > N).

NSGA-II key algorithm (steps):
1. First find the pareto solution for M individuals. Then we get the set of these pareto, such as F1, F2.
2. Put all individuals of F1 into N. if N is not full, continue to put F2 until Fk cannot be put into N (space) that has been put into F1, F2,..., F(k-1). Fk is solved at this time.
3. For individuals in FK, find the crowding distance Lk[i](crowding distance) of each individual in FK, sort it in descending order according to Lk[i] in FK, and put it into N until N is full.

Figure 4-2 generation of parents and offspring of NSGA-II

Among them, NSGA-II key subroutine algorithm
1) Fast non dominated sorting algorithm
   the key of multi-objective optimization problem is to find the Pareto optimal solution set. NSGA-II fast non dominated sorting is to obtain Fi by layering the population m according to the individual non inferior solution level, which is used to make the solution close to the Pareto optimal solution. This is a cyclic fitness grading process. First, find the non dominated solution set in the group, record it as F1, assign all its individuals to the non dominated order irank=1 (where irank is the non dominated order value of individual i), remove it from the whole group M, and then continue to find the non dominated solution set in the remaining group, record it as F2, and the individuals in F2 are given irank=2, It is known that the whole population is stratified, and the non dominated order value in Fi layer is the same.

2) Individual crowding distance
   selective sorting is required in the same layer Fk, which is sorted according to the size of individual crowding distance. Individual crowding distance is the distance between individuals i+1 and i-1 adjacent to i on Fk. Its calculation steps are as follows:
① Initialize the individual distance of the same layer so that L[i]d=0 (represents the congestion distance of any individual I).
② The individuals in the same layer are arranged in ascending order according to the value of the m-th objective function.
③ Individuals on the edge of ranking should be given a choice advantage.
④ For the individuals in the middle of the ranking, calculate the crowding distance:
⑤ For different objective functions, repeat steps ② to ④ to obtain the crowding distance L[i]d of individual I. limited selection of individuals with large crowding distance can make the calculation results evenly distributed in the target space and maintain the diversity of the population.

3) Elite strategy selection algorithm
  keep the excellent individuals in the parent generation directly into the offspring to prevent the loss of Pareto optimal solution. Select the index to optimize the population Ri synthesized by parent Ci and offspring Di to form a new parent Ci+1 First eliminate the scheme whose scheme test flag is infeasible in the parent generation, and then put the whole layer of population into Ci+1 from low to high according to the non dominated order value irank until Fk in a certain layer exceeds the limit of N. finally, fill Ci+1 according to the crowding distance until the population number is n.

   among them, Figure 4-3 shows the Pareto front of the double objective function, in which the smaller the objective functions f1 and f2, the better. The solid line and the dotted line constitute the feasible region. The solid line represents the Pareto frontier, that is, the surface composed of the target vectors corresponding to all Pareto optimal solutions. A multi-objective optimization problem corresponds to a Pareto frontier. In addition, points a, B and C are located on the Pareto front. The solutions of these three points are Pareto optimal solutions, and there is no dominant or dominant relationship among them. D. The solutions of points E and F are feasible solutions and non Pareto optimal solutions. The solution of point a dominates the solution of point F, or the solution of point a is Pareto dominant compared with the solution of point F. Similarly, there is a dominant or dominant relationship between points B and C and points D, e and F.

Figure 4-3 Pareto front of non dominated solution under double objective function

Figure 4-4 general steps of NSGA-II algorithm

Pseudo code of NSGA-2 algorithm

Pseudo code of sorting algorithm

def  fast_nondominated_sort( P ):
    F  =  [ ]
     for  p  in  P:
        Sp  =  [ ]
         np  =  0
          for  q  in  P:
              if  p  >  q:                 # If p dominates q, add q to the Sp list
                 Sp.append( q )
              else   if  p  <  q:         # If p is dominated by q, add np to 1
                 np  +=   1

         if  np  ==  0:
            p_rank  =   1          # If the individual's np is 0, the individual is Pareto level 1
      F1.append( p )
    F.append( F1 )
    i  =  0
     while  F[i]:
        Q  =  [ ]
         for  p  in  F[i]:
             for  q  in  Sp:         # Sort all individuals in the Sp set
                nq  -=   1
                 if  nq  ==  0:      # If the dominant number of the individual is 0, the individual is a non dominant individual
                    q_rank  =  i + 2      # The individual Pareto level is the current highest level plus 1. At this time, the initial value of i is 0, so add 2
                    Q.append( q )
        F.append( Q )
        i  +=   1

Pseudocode of individual crowding distance operator

def  crowding_distance_assignment( I )
        nLen  =  len( I )         # Number of individuals in I
     for  i  in  I:
                i.distance  =  0     # Initialize the crowding distance of all individuals
     for  objFun  in  M:         # M is the list of all objective functions
                I  =  sort( I, objFun )     # Sort in ascending order according to the objective function objFun
                I[0]  =  I[ len[I] - 1  ]  =  ∞     # Set the distance between the first and last individuals to infinity
                 for  i  in  xrange(  1 , len(I)  -   2  ):
                        I[i].distance  =  I[i].distance  +  ( objFun( I[i + 1 ] )  -  objFun( I[i - 1 ] ) ) / (Max(objFun())  -  Min(objFun()) )

1.3 case implementation of nsga-2 algorithm

   this blog posts the Matlab implementation of NSGA - Ⅱ algorithm (the test function is ZDT1), the original link: Matlab implementation of NSGA - Ⅱ algorithm (the test function is ZDT1)

The first test function ZDT1

Figure 4-5 test function ZDT1

The code is as follows:

function f = evaluate_objective(x, M, V)
f = [];
f(1) = x(1);
g = 1;
sum = 0;
for i = 2:V
    sum = sum + x(i);
g = g + 9*sum / (V-1));
f(2) = g * (1 - sqrt(x(1) / g));

The NSGA2 code is implemented as follows:

function NSGAII()
format compact;
hold on
%---initialization/Parameter setting
    generations=100;                                %Number of iterations
    popnum=100;                                     %Population size(Must be even)
    poplength=30;                                   %Individual length
    minvalue=repmat(zeros(1,poplength),popnum,1);   %Individual minimum
    maxvalue=repmat(ones(1,poplength),popnum,1);    %Individual maximum    
    population=rand(popnum,poplength).*(maxvalue-minvalue)+minvalue;    %Generate a new initial population
%---Start iterative evolution
    for gene=1:generations                      %Start iteration
        newpopulation=zeros(popnum,poplength);  %Offspring population
        for i=1:popnum/2                        %Cross generation offspring
            k=randperm(popnum);                 %Two parents were randomly selected from the population,Binary League method is not used
            beta=(-1).^round(rand(1,poplength)).*abs(randn(1,poplength))*1.481; %The normal distribution is used to cross produce two offspring
            newpopulation(i*2-1,:)=(population(k(1),:)+population(k(2),:))/2+beta.*(population(k(1),:)-population(k(2),:))./2;    %Produce the first generation           
            newpopulation(i*2,:)=(population(k(1),:)+population(k(2),:))/2-beta.*(population(k(1),:)-population(k(2),:))./2;      %Produce the second generation
        k=rand(size(newpopulation));    %Randomly select the loci to be mutated
        miu=rand(size(newpopulation));  %Polynomial variation
        temp=k<1/poplength & miu<0.5;   %Loci to be mutated
        newpopulation(temp)=newpopulation(temp)+(maxvalue(temp)-minvalue(temp)).*((2.*miu(temp)+(1-2.*miu(temp)).*(1-(newpopulation(temp)-minvalue(temp))./(maxvalue(temp)-minvalue(temp))).^21).^(1/21)-1);        %Variation I
        newpopulation(temp)=newpopulation(temp)+(maxvalue(temp)-minvalue(temp)).*(1-(2.*(1-miu(temp))+2.*(miu(temp)-0.5).*(1-(maxvalue(temp)-newpopulation(temp))./(maxvalue(temp)-minvalue(temp))).^21).^(1/21));  %Variation II
%-------Cross border processing/Population merging        
        newpopulation(newpopulation>maxvalue)=maxvalue(newpopulation>maxvalue); %Offspring over upper bound processing
        newpopulation(newpopulation<minvalue)=minvalue(newpopulation<minvalue); %Offspring cross lower bound processing
        newpopulation=[population;newpopulation];                               %Combined parent-child population
%-------Calculate the objective function value        
        functionvalue=zeros(size(newpopulation,1),2);           %Objective function values of the combined population,The problem here is ZDT1
        functionvalue(:,1)=newpopulation(:,1);                  %Calculate the value of the first dimensional objective function
        functionvalue(:,2)=g.*(1-(newpopulation(:,1)./g).^0.5); %Calculate the value of the second objective function
%-------Non dominated sorting        
        fnum=0;                                             %Currently assigned leading edge number
        cz=false(1,size(functionvalue,1));                  %Record whether the individual has been assigned a number
        frontvalue=zeros(size(cz));                         %Leading edge number of each individual
        [functionvalue_sorted,newsite]=sortrows(functionvalue);    %Sort the population according to the target value of the first dimension
        while ~all(cz)                                      %Start iteratively judging the frontiers of each individual,Adopt improved deductive sort
            for i=1:size(functionvalue,1)
                if ~d(i)
                    for j=i+1:size(functionvalue,1)
                        if ~d(j)
                            for m=2:size(functionvalue,2)
                                if functionvalue_sorted(i,m)>functionvalue_sorted(j,m)
                            if k
%-------Calculate congestion distance/Select the next generation of individuals        
        fnum=0;                                                                 %Current frontier
        while numel(frontvalue,frontvalue<=fnum+1)<=popnum                      %Judge how many individuals of the first face can be completely put into the next generation population
        newnum=numel(frontvalue,frontvalue<=fnum);                              %front fnum Number of individuals per face
        population(1:newnum,:)=newpopulation(frontvalue<=fnum,:);               %Before fnum Individual faces are copied into the next generation                       
        popu=find(frontvalue==fnum+1);                                          %popu Record No fnum+1 Individual number on one face
        distancevalue=zeros(size(popu));                                        %popu Crowding distance of each individual
        fmax=max(functionvalue(popu,:),[],1);                                   %popu Maximum value per dimension
        fmin=min(functionvalue(popu,:),[],1);                                   %popu Minimum value per dimension
        for i=1:size(functionvalue,2)                                           %Sub target calculation on each target popu Crowding distance of each individual
            for j=2:length(popu)-1
        popu=-sortrows(-[distancevalue;popu]')';                                %Ranked in descending order of crowding distance fnum+1 Individual on one face
        population(newnum+1:popnum,:)=newpopulation(popu(2,1:popnum-newnum),:);	%Will be the first fnum+1 Front with large crowding distance popnum-newnum Individuals are copied into the next generation        
%---Program output    
    fprintf('Completed,time consuming%4s second\n',num2str(toc));          %The final time of the program
    output=sortrows(functionvalue(frontvalue==1,:));    %final result:Function value of non dominated solution in population
    plot(output(:,1),output(:,2),'*b');                 %Mapping

The running result is shown in the following figure:

Figure 4-6 example implementation results of nsga-2 algorithm

Second test function

The code implementation of the test function is as follows:

function y=f(x)

The code implementation of NSGA-2 MATLAB built-in function is as follows:

% At present, there are many multi-objective optimization algorithms,?Kalyanmoy?Deb Fast non dominated sorting genetic algorithm with elite strategy(NSGA-II)
%It is undoubtedly the most widely used and successful one. The algorithm used in this paper is MATLAB Self contained function gamultiobj,Based on this function NSGA-II An improved multi-objective optimization algorithm.

clear all;
fitnessfcn=@object_fun;  %Fitness function handle,objective function 
nvars=2;                 %Number of variables x1 x2
lb=[-5,-5];              %lower limit
ub=[5,5];                %upper limit
b=[];                    %Linear inequality constraint
beq=[];                  %linear equality constraint 
%Optimal individual coefficient paretoFaction Is 0.3;Population size populationsize Is 100,
%Maximum evolutionary algebra generations Stop algebra for 200 stallGenlimit 200,
%Fitness function deviation TolFun Set to 1 e-100 function gaplotpareto: draw Pareto front end
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)   %MATLAB Self contained function gamultiobj

Figure 4-7 implementation result of nsga-2 matlab function

2 design of controlled object and fitness function

Introduction of controlled object

   taking the transfer function of the second-order type I time-delay system as an example, the NSGA-2 algorithm is used to optimize the PID parameters, in which the system is set as the sampling time of 1 ms, the instruction is the unit step signal, and the simulation running time is 1.0 s. Among them, the performance optimization function is Best_J adopts the Integral of Time Multiplied by the Absolute Value of Error (ITAE). At the same time, in order to avoid overshoot caused by too large control quantity, in the performance optimization function best_ Add the square term of PID controller input to J. This is only the fitness function of one of the objectives. In addition, it is necessary to design the fitness function of other objectives, such as the design of fitness function in 2.2.
   the transfer function of second-order type I delay system is as follows. It can be designed according to its own actual system m document

2.2 design of fitness function

(1) Objective fitness function 1
   in order to obtain a satisfactory transition process, the absolute value of error time integration performance index is used as the fitness evaluation function Best_J. At the same time, in order to prevent the control input from being too large, in Best_J add the scoring item of control input, as shown in formula (1-2).

In equation (1-2), e(t) is the system output error, and u(t) is the input of PID controller, ρ 1, ρ 2 is the weight value.
   in order to avoid overshoot, the penalty function is used to give priority to the overshoot, as shown in equation (1-3).

In formula (1-3) ρ 3>>max( ρ 1, ρ 2 and ρ 4) , and y(t) is the output of the controlled object, ey(t)=y(t)-y(t-1).
This is the design of the first objective fitness function.
(2) Objective fitness function 2
The rise time, adjustment time and delay time of the system are superimposed according to the principle of time minimization to construct the design of the second objective fitness function.

3PID parameter setting

3.1 PID parameter tuning of nsga-2 algorithm

global rin yout timef
pop = 500; %Population number
gen = 100; %Number of iterations
M = 2;  %Target quantity
V = 3;   %dimension

MinX(1) = 0*ones(1); %Kp
MaxX(1) = 30*ones(1);

MinX(2) = 0*ones(1); %Kp
MaxX(2) = 1.0*ones(1);

MinX(3) = 0*ones(1); %Kp
MaxX(3) = 1.0*ones(1);

min_range = [MinX(1) MinX(2) MinX(3)]; %Lower bound
max_range = [MaxX(1) MaxX(2) MaxX(3)];   %upper bound

%% initialization
chromosome = initialize_variables(pop, M, V, min_range, max_range);
%Non dominated quick sequencing and congestion calculation are carried out for the initial population

chromosome = non_domination_sort_mod(chromosome, M, V);

%% Main program loop
for i = 1 : gen
pool = round(pop/2); %Mating pool size
tour = 2;            %Number of contestants in the bidding competition
%Selection of parents suitable for breeding in the competition
parent_chromosome = tournament_selection(chromosome, pool, tour);
mu = 20;
mum = 20;
offspring_chromosome = genetic_operator(parent_chromosome,M, V, mu, mum, min_range, max_range);
[main_pop,~] = size(chromosome);
[offspring_pop,~] = size(offspring_chromosome);
clear temp
intermediate_chromosome(1:main_pop,:) = chromosome;
intermediate_chromosome(main_pop + 1 : main_pop + offspring_pop,1 : M+V) = offspring_chromosome;
intermediate_chromosome = non_domination_sort_mod(intermediate_chromosome, M, V);
chromosome = replace_chromosome(intermediate_chromosome, M, V, pop);
    if ~mod(i,100)
    fprintf('%d generations completed\n',i);
if M == 2
plot(chromosome(:,V + 1),chromosome(:,V + 2),'*');
xlabel('f_1'); ylabel('f_2');
title('Pareto Optimal Front');
elseif M == 3
plot3(chromosome(:,V + 1),chromosome(:,V + 2),chromosome(:,V + 3),'*');
xlabel('f_1'); ylabel('f_2'); zlabel('f_3');
title('Pareto Optimal Surface');

%% //*******Subroutine initialize_variables**********//
function f = initialize_variables(N, M, V, min_range, max_range)
min = min_range;
max = max_range;

K = M + V;

    for i = 1 : N
        for j = 1 : V
        f(i,j) = min(j) + (max(j) - min(j))*rand(1);
        Kpidi=f(i,:);    %Single feasible solution
    f(i,V + 1: K) = evaluate_objective(Kpidi, M, V);


%% //*******Subroutine non_domination_sort_mod**********//
function f = non_domination_sort_mod(x, M, V)

[N, ~] = size(x);
clear m
front = 1;
F(front).f = [];
individual = [];
for i = 1 : N
individual(i).n = 0;
individual(i).p = [];
for j = 1 : N
dom_less = 0;
dom_equal = 0;
dom_more = 0;
for k = 1 : M
if (x(i,V + k) < x(j,V + k))
dom_less = dom_less + 1;
elseif (x(i,V + k) == x(j,V + k))
dom_equal = dom_equal + 1;
dom_more = dom_more + 1;
if dom_less == 0 && dom_equal ~= M
individual(i).n = individual(i).n + 1;
elseif dom_more == 0 && dom_equal ~= M
individual(i).p = [individual(i).p j];
if individual(i).n == 0
x(i,M + V + 1) = 1;
F(front).f = [F(front).f i];
while ~isempty(F(front).f)
Q = [];
for i = 1 : length(F(front).f)
if ~isempty(individual(F(front).f(i)).p)
for j = 1 : length(individual(F(front).f(i)).p)
individual(individual(F(front).f(i)).p(j)).n = ...
individual(individual(F(front).f(i)).p(j)).n - 1;
if individual(individual(F(front).f(i)).p(j)).n == 0
x(individual(F(front).f(i)).p(j),M + V + 1) = ...
front + 1;
Q = [Q individual(F(front).f(i)).p(j)];
front =  front + 1;
F(front).f = Q;
[temp,index_of_fronts] = sort(x(:,M + V + 1));
for i = 1 : length(index_of_fronts)
sorted_based_on_front(i,:) = x(index_of_fronts(i),:);
current_index = 0;
%% Crowding distance
for front = 1 : (length(F) - 1)
distance = 0;
y = [];
previous_index = current_index + 1;
for i = 1 : length(F(front).f)
y(i,:) = sorted_based_on_front(current_index + i,:);
current_index = current_index + i;
sorted_based_on_objective = [];
for i = 1 : M
[sorted_based_on_objective, index_of_objectives] = ...
sort(y(:,V + i));
sorted_based_on_objective = [];
for j = 1 : length(index_of_objectives)
sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);
f_max = ...
sorted_based_on_objective(length(index_of_objectives), V + i);
f_min = sorted_based_on_objective(1, V + i);
y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i)...
= Inf;
y(index_of_objectives(1),M + V + 1 + i) = Inf;
for j = 2 : length(index_of_objectives) - 1
next_obj  = sorted_based_on_objective(j + 1,V + i);
previous_obj  = sorted_based_on_objective(j - 1,V + i);
if (f_max - f_min == 0)
y(index_of_objectives(j),M + V + 1 + i) = Inf;
y(index_of_objectives(j),M + V + 1 + i) = ...
(next_obj - previous_obj)/(f_max - f_min);
distance = [];
distance(:,1) = zeros(length(F(front).f),1);
for i = 1 : M
distance(:,1) = distance(:,1) + y(:,M + V + 1 + i);
y(:,M + V + 2) = distance;
y = y(:,1 : M + V + 2);
z(previous_index:current_index,:) = y;
f = z();

%% //*******Subroutine tournament_selection**********//
function f = tournament_selection(chromosome, pool_size, tour_size)
[pop, variables] = size(chromosome);
rank = variables - 1;
distance = variables;
for i = 1 : pool_size
for j = 1 : tour_size
candidate(j) = round(pop*rand(1));
if candidate(j) == 0
candidate(j) = 1;
if j > 1
while ~isempty(find(candidate(1 : j - 1) == candidate(j)))
candidate(j) = round(pop*rand(1));
if candidate(j) == 0
candidate(j) = 1;
for j = 1 : tour_size
c_obj_rank(j) = chromosome(candidate(j),rank);
c_obj_distance(j) = chromosome(candidate(j),distance);
min_candidate = ...
find(c_obj_rank == min(c_obj_rank));
if length(min_candidate) ~= 1
max_candidate = ...
find(c_obj_distance(min_candidate) == max(c_obj_distance(min_candidate)));
if length(max_candidate) ~= 1
max_candidate = max_candidate(1);
f(i,:) = chromosome(candidate(min_candidate(max_candidate)),:);
f(i,:) = chromosome(candidate(min_candidate(1)),:);

%% //*******Subroutine genetic_operator**********//
function f  = genetic_operator(parent_chromosome, M, V, mu, mum, l_limit, u_limit)
[N,m] = size(parent_chromosome);
clear m
p = 1;
was_crossover = 0;
was_mutation = 0;
for i = 1 : N
% With 90 % probability perform crossover
if rand(1) < 0.9
% Initialize the children to be null vector.
child_1 = [];
child_2 = [];
% Select the first parent
parent_1 = round(N*rand(1));
if parent_1 < 1
parent_1 = 1;
% Select the second parent
parent_2 = round(N*rand(1));
if parent_2 < 1
parent_2 = 1;
% Make sure both the parents are not the same. 
while isequal(parent_chromosome(parent_1,:),parent_chromosome(parent_2,:))
parent_2 = round(N*rand(1));
if parent_2 < 1
parent_2 = 1;
% Get the chromosome information for each randomnly selected
% parents
parent_1 = parent_chromosome(parent_1,:);
parent_2 = parent_chromosome(parent_2,:);
% Perform corssover for each decision variable in the chromosome.
for j = 1 : V
% SBX (Simulated Binary Crossover).
% For more information about SBX refer the enclosed pdf file.
% Generate a random number
u(j) = rand(1);
if u(j) <= 0.5
bq(j) = (2*u(j))^(1/(mu+1));
bq(j) = (1/(2*(1 - u(j))))^(1/(mu+1));
% Generate the jth element of first child
child_1(j) = ...
0.5*(((1 + bq(j))*parent_1(j)) + (1 - bq(j))*parent_2(j));
% Generate the jth element of second child
child_2(j) = ...
0.5*(((1 - bq(j))*parent_1(j)) + (1 + bq(j))*parent_2(j));
% Make sure that the generated element is within the specified
% decision space else set it to the appropriate extrema.
if child_1(j) > u_limit(j)
child_1(j) = u_limit(j);
elseif child_1(j) < l_limit(j)
child_1(j) = l_limit(j);
if child_2(j) > u_limit(j)
child_2(j) = u_limit(j);
elseif child_2(j) < l_limit(j)
child_2(j) = l_limit(j);
child_1(:,V + 1: M + V) = evaluate_objective(child_1, M, V);
child_2(:,V + 1: M + V) = evaluate_objective(child_2, M, V);
was_crossover = 1;
was_mutation = 0;
% With 10 % probability perform mutation. Mutation is based on
% polynomial mutation. 
% Select at random the parent.
parent_3 = round(N*rand(1));
if parent_3 < 1
parent_3 = 1;
% Get the chromosome information for the randomnly selected parent.
child_3 = parent_chromosome(parent_3,:);
% Perform mutation on eact element of the selected parent.
for j = 1 : V
r(j) = rand(1);
if r(j) < 0.5
delta(j) = (2*r(j))^(1/(mum+1)) - 1;
delta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1));
% Generate the corresponding child element.
child_3(j) = child_3(j) + delta(j);
% Make sure that the generated element is within the decision
% space.
if child_3(j) > u_limit(j)
child_3(j) = u_limit(j);
elseif child_3(j) < l_limit(j)
child_3(j) = l_limit(j);
child_3(:,V + 1: M + V) = evaluate_objective(child_3, M, V);
% Set the mutation flag
was_mutation = 1;
was_crossover = 0;
if was_crossover
child(p,:) = child_1;
child(p+1,:) = child_2;
was_cossover = 0;
p = p + 2;
elseif was_mutation
child(p,:) = child_3(1,1 : M + V);
was_mutation = 0;
p = p + 1;
f = child;

%% //*******Subroutine replace_chromosome**********//
function f  = replace_chromosome(intermediate_chromosome, M, V,pop)
[N, m] = size(intermediate_chromosome);
% Get the index for the population sort based on the rank
[temp,index] = sort(intermediate_chromosome(:,M + V + 1));
clear temp m
% Now sort the individuals based on the index
for i = 1 : N
sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));
% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
% Get the index for current rank i.e the last the last element in the
% sorted_chromosome with rank i. 
current_index = max(find(sorted_chromosome(:,M + V + 1) == i));
% Check to see if the population is filled if all the individuals with
% rank i is added to the population. 
if current_index > pop
% If so then find the number of individuals with in with current
% rank i.
remaining = pop - previous_index;
% Get information about the individuals in the current rank i.
temp_pop = ...
sorted_chromosome(previous_index + 1 : current_index, :);
% Sort the individuals with rank i in the descending order based on
% the crowding distance.
[temp_sort,temp_sort_index] = ...
sort(temp_pop(:, M + V + 2),'descend');
% Start filling individuals into the population in descending order
% until the population is filled.
for j = 1 : remaining
f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
elseif current_index < pop
% Add all the individuals with rank i into the population.
f(previous_index + 1 : current_index, :) = ...
sorted_chromosome(previous_index + 1 : current_index, :);
% Add all the individuals with rank i into the population.
f(previous_index + 1 : current_index, :) = ...
sorted_chromosome(previous_index + 1 : current_index, :);
% Get the index for the last added individual.
previous_index = current_index;

3.3 simulation diagram of PID parameter setting of nsga-2 algorithm

Figure 4-8 fitness function of nsga-2 algorithm

Figure 4-9 simulation results of PID parameter setting of nsga-2 algorithm

4 references

[1] Yuan Lei, Hu Bingxin, Wei keyin, et al Modern permanent magnet synchronous motor control principle and MATLAB simulation [M] Beijing: Beijing University of Aeronautics and Astronautics Press, 2016 (03): 12-18
[2]Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2):182-197.
[3]Zhang Q, Li H. MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[J]. IEEE Transactions on Evolutionary Computation, 2007, 11(6):712-731.
[4]M. Andersson, S. Bandaru and A. H. C. Ng, "Towards optimal algorithmic parameters for simulation-based multi-objective optimization," 2016 IEEE Congress on Evolutionary Computation (CEC), Vancouver, BC, 2016, pp. 5162-5169.