## 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); end g = g + 9*sum / (V-1)); f(2) = g * (1 - sqrt(x(1) / g)); end

The NSGA2 code is implemented as follows:

function NSGAII() clc; format compact; tic; 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 %-------overlapping 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 end %-------variation 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 g=1+9*sum(newpopulation(:,2:poplength),2)./(poplength-1); 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 fnum=fnum+1; d=cz; for i=1:size(functionvalue,1) if ~d(i) for j=i+1:size(functionvalue,1) if ~d(j) k=1; for m=2:size(functionvalue,2) if functionvalue_sorted(i,m)>functionvalue_sorted(j,m) k=0; break end end if k d(j)=true; end end end frontvalue(newsite(i))=fnum; cz(i)=true; end end end %-------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 fnum=fnum+1; end 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 [~,newsite]=sortrows(functionvalue(popu,i)); distancevalue(newsite(1))=inf; distancevalue(newsite(end))=inf; for j=2:length(popu)-1 distancevalue(newsite(j))=distancevalue(newsite(j))+(functionvalue(popu(newsite(j+1)),i)-functionvalue(popu(newsite(j-1)),i))/(fmax(i)-fmin(i)); end end 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 end %---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 axis([0,1,0,1]);xlabel('F_1');ylabel('F_2');title('ZDT1') end

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) y(1)=x(1)^4-10*x(1)^2+x(1)*x(2)+x(2)^4-x(1)^2*x(2)^2 y(2)=x(2)^4-x(1)^2*x(2)^2+x(1)^4+x(1)*x(2)

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. clc; 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 A=[]; b=[]; %Linear inequality constraint Aeq=[]; beq=[]; %linear equality constraint options=gaoptimset('paretoFraction',0.3,'populationsize',500,'generations',200,... 'stallGenLimit',200,'TolFun',1e-100,'PlotFcns',@gaplotpareto); %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) clc; fprintf('%d generations completed\n',i); end end figure(1); 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'); end figure(2); plot(timef,rin,'r',timef,yout,'b'); xlabel('Time(s)');ylabel('rin,yout'); %% //*******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); end Kpidi=f(i,:); %Single feasible solution f(i,V + 1: K) = evaluate_objective(Kpidi, M, V); end end %% //*******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; else dom_more = dom_more + 1; end end 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]; end end if individual(i).n == 0 x(i,M + V + 1) = 1; F(front).f = [F(front).f i]; end end 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)]; end end end end front = front + 1; F(front).f = Q; end [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),:); end 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,:); end 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),:); end 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; else y(index_of_objectives(j),M + V + 1 + i) = ... (next_obj - previous_obj)/(f_max - f_min); end end end distance = []; distance(:,1) = zeros(length(F(front).f),1); for i = 1 : M distance(:,1) = distance(:,1) + y(:,M + V + 1 + i); end y(:,M + V + 2) = distance; y = y(:,1 : M + V + 2); z(previous_index:current_index,:) = y; end f = z(); end %% //*******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; end 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; end end end end for j = 1 : tour_size c_obj_rank(j) = chromosome(candidate(j),rank); c_obj_distance(j) = chromosome(candidate(j),distance); end 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); end f(i,:) = chromosome(candidate(min_candidate(max_candidate)),:); else f(i,:) = chromosome(candidate(min_candidate(1)),:); end end end %% //*******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; end % Select the second parent parent_2 = round(N*rand(1)); if parent_2 < 1 parent_2 = 1; end % 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; end end % 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)); else bq(j) = (1/(2*(1 - u(j))))^(1/(mu+1)); end % 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); end 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); end end 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. else % Select at random the parent. parent_3 = round(N*rand(1)); if parent_3 < 1 parent_3 = 1; end % 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; else delta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1)); end % 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); end end child_3(:,V + 1: M + V) = evaluate_objective(child_3, M, V); % Set the mutation flag was_mutation = 1; was_crossover = 0; end 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; end end f = child; end %% //*******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),:); end % 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),:); end return; 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, :); else % Add all the individuals with rank i into the population. f(previous_index + 1 : current_index, :) = ... sorted_chromosome(previous_index + 1 : current_index, :); return; end % Get the index for the last added individual. previous_index = current_index; end end

### 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.