多目标遗传算法代码

% function nsga_2(pro)% Main Function% Main program to run the NSGA-II MOEA.% Read the corresponding documentation to learn more about multiobjective% optimization using evolutionary algorithms.% initialize_variables has two arguments; First being the population size% and the second the problem number. 1 corresponds to MOP1 and 2% corresponds to MOP2.%inp_para_definition=input_parameters_definition;% Initialize the variables% Declare the variables and initialize their values% pop - population% gen - generations% pro - problem number%clear;clc;tic;pop = 100; %每一代的种群数gen = 100; %总共的代数pro = 2;%问题选择1或者2,见switchswitch procase 1% M is the number of objectives.M = 2;% V is the number of decision variables. In this case it is% difficult to visualize the decision variables space while the% objective space is just two dimensional.V = 6;case 2M = 3;V = 12;case 3% case 1和case 2用来对整个算法进行常规验证，作为调试之用；case 3为本工程所需；M = 2; %(output parameters 个数)V = 8; %(input parameters 个数)K = 10;end% Initialize the populationchromosome = initialize_variables(pop,pro);% Sort the initialized population% Sort the population using non-domination-sort. This returns two columns% for each individual which are the rank and the crowding distance% corresponding to their position in the front they belong. 真是牛 X 了。chromosome = non_domination_sort_mod(chromosome,pro);% Start the evolution process% The following are performed in each generation% Select the parents% Perfrom crossover and Mutation operator% Perform Selectionfor i = 1 : gen% Select the parents% Parents are selected for reproduction to generate offspring. The% original NSGA-II uses a binary tournament selection based on the% crowded-comparision operator. The arguments are% pool - size of the mating pool. It is common to have this to be half the%population size.% tour - Tournament size. Original NSGA-II uses a binary tournament%selection, but to see the effect of tournament size this is kept%arbitary, to be choosen by the user.pool = round(pop/2);tour = 2;%下面进行二人锦标赛配对，新的群体规模是原来群体的一半parent_chromosome = tournament_selection(chromosome,pool,tour);% Perfrom crossover and Mutation operator% The original NSGA-II algorithm uses Simulated Binary Crossover (SBX) and% Polynomial crossover. Crossover probability pc = 0.9 and mutation% probability is pm = 1/n, where n is the number of decision variables.% Both real-coded GA and binary-coded GA are implemented in the original% algorithm, while in this program only the real-coded GA is considered.% The distribution indeices for crossover and mutation operators as mu = 20% and mum = 20 respectively.mu = 20;mum = 20;%针对对象是上一步产生的新的个体parent_chromosome%对parent_chromosome每次操作以较大的概率进行交叉(产生两个新的候选人)，或 者较小的概率变异(一个新的候选人)操作，这样%就会产生较多的新个体offspring_chromosome = genetic_operator(parent_chromosome,pro,mu,mum);% Intermediate population% Intermediate population is the combined population of parents and% offsprings of the current generation. The population size is almost 1 and % half times the initial population.main_pop,temp=size(chromosome); offspring_pop,temp=size(offspring_chromosome);intermediate_chromosome(1:main_pop,:)=chromosome;intermediate_chromosome(main_pop+l:main_pop+offspring_pop,l:M+V)=offspring_chromosome；%intermediate_chromosome=inter_chromo(chromosome,offspring_chromosome,pro);% Non-domination-sort of intermediate population% The intermediate population is sorted again based on non-domination sort% before the replacement operator is performed on the intermediate% population.intermediate_chromosome =.non_domination_sort_mod(intermediate_chromosome,pro);% Perform Selection% Once the intermediate population is sorted only the best solution is% selected based on it rank and crowding distance. Each front is filled in% ascending order until the addition of population size is reached. The% last front is included in the population based on the individuals with% least crowding distancechromosome = replace_chromosome(intermediate_chromosome,pro,pop);if mod(i,10)fprintf(%dn,i);endend% Result% Save the result in ASCII text format.save solution.txt chromosome -ASCII% Visualize% The following is used to visualize the result for the given problem.switch procase 1plot(chromosome(:,V + 1),chromosome(:,V + 2),y+);title(MOP1 using NSGA-II);xlabel(f(x_1);ylabel(f(x_2);case 2plot3(chromosome(:,V + 1),chromosome(:,V + 2),chromosome(:,V + 3),*);title(MOP2 using NSGA-II);xlabel(f(x_1);ylabel(f(x_2);zlabel(f(x_3);end%disp(run time is:)%toc;%function f = initialize_variables(N,problem) % function f = initialize_variables(N,problem)% N - Population size% problem - takes integer values 1 and 2 where,%T for M0P1%2 for MOP2% This function initializes the population with N individuals and each% individual having M decision variables based on the selected problem.% M = 6 for problem MOP1 and M = 12 for problem MOP2. The objective space% for MOP1 is 2 dimensional while for MOP2 is 3 dimensional.% Both the MOPs has 0 to 1 as its range for all the decision variables. min = 0;max = 1;switch problemcase 1M = 6;K = 8;% k=决策变量(M=6)+目标变量(K-M=2)=8case 2M = 12; K = 15;case 3% case 1和case 2用来对整个算法进行常规验证，作为调试之用；case 3为本工程所需；M = 8; %(input parameters 个数)K = 10;endfor i = 1 : N% Initialize the decision variablesforj = 1 : Mf(i,j) = rand(1); % i.e f(i,j) = min + (max - min)*rand(1);end% Evaluate the objective functionf(i,M + 1: K) = evaluate_objective(f(i,:),problem);end% %function f = evaluate_objective(x,problem)% Function to evaluate the objective functions for the given input vector% x. x has the decision variablesswitch problemcase 1f =;% Objective function onef(l) = 1 - exp(-4*x(l)*(sin(6*pi*x(l)人6;sum = 0;for i = 2 : 6sum = sum + x(i)/4;end% Intermediate functiong_x = 1 + 9*(sum)人(0.25);% Objective function onef(2) = g_x*(1 - (f(1)/(g_x)人2);case 2f =;% Intermediate functiong_x = 0;for i = 3 : 12g_x = g_x + (x(i) - 0.5)A2;end% Objective function onef(1) = (1 + g_x)*cos(0.5*pi*x(1)*cos(0.5*pi*x(2);% Objective function twof(2) = (1 + g_x)*cos(0.5*pi*x(1)*sin(0.5*pi*x(2);% Objective function threef(3) = (1 + g_x)*sin(0.5*pi*x(1);case 3f =;% Objective function onef(1) = 0;% Objective function onef(2) = 0;end% Non-Donimation Sort%按照目标函数最小了好% This function sort the current popultion based on non-domination. All the% individuals in the first front are given a rank of 1, the second front% individuals are assigned rank 2 and so on. After assigning the rank the% crowding in each front is calculated.function f = non_domination_sort_mod(x,problem)N,M = size(x);switch problemcase 1M = 2;V = 6;case 2M = 3;V = 12;case 3% case 1和case 2用来对整个算法进行常规验证，作为调试之用；case 3为本工程所需；M = 2; %(output parameters 个数)V = 8; %(input parameters 个数)K = 10;endfront = 1;% There is nothing to this assignment, used only to manipulate easily in% MATLAB.F(front).f =; individual =;for i = 1 : N% Number of individuals that dominate this individual 支配 i 的解的个数individual(i).n = 0;% Individuals which this individual dominate被 i 支配的解individual(i).p =;forj = 1 : Ndom_less = 0;dom_equal = 0; dom_more = 0;for k = 1 : Mif (x(i,V + k) return a increasing orer data sequencetemp,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),:); % 对解(染色体)进行排序，按照 front分层end%到这里分层结束，下面就是计算距离了current_index = 0;% Find the crowding distance for each individual in each front%，计算不同层的拥挤距离是没有意义的for front = 1 : (length(F) - 1)objective =;distance = 0;y =;previous_index = current_index + 1;for i = 1 : length(F(front).f)y(i,:) = sorted_based_on_front(current_index + i,:);endcurrent_index = current_index + i;% Sort each individual based on the objective sorted_based_on_objective =;for i = 1 : Msorted_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),:);endf_max =.sorted_based_on_objective(length(index_of_objectives), V + i);% 最大值f_min = sorted_based_on_objective(l, V + i);%最小值y(index_of_objectives(length(index_of_objectives),M + V + 1 + i).=Inf; % 最大值放lnfy(index_of_objectives(1),M + V + 1 + i) = Inf; % 最小值放 lnffor j = 2 : length(index_of_objectives) - 1next_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;elsey(index_of_objectives(j),M + V + 1 + i)=.(next_obj - previous_obj)/(f_max - f_min);endendenddistance =;distance(:,1) = zeros(length(F(front).f),1);for i = 1 : Mdistance(:,1) = distance(:,1) + y(:,M + V + 1 + i);endy(:,M + V + 2) = distance;y = y(:,1 : M + V + 2); z(previous_index:current_index,:) = y;endf = z;% function f = replace_chromosome(intermediate_chromosome,pro,pop)% replace_chromosome(intermediate_chromosome,pro,pop)% This function replaces the chromosomes based on rank and crowding% distance. Initially until the population size is reached each front is% added one by one until addition of a complete front which results in% exceeding the population size. At this point the chromosomes in that% front is added subsequently to the population based on crowding distance.N,V = size(intermediate_chromosome);switch procase 1M = 2;V = 6;case 2M = 3;V = 12;case 3% case 1和case 2用来对整个算法进行常规验证，作为调试之用；case 3为本工程所需；M = 2; %(output parameters 个数)V = 8; %(input parameters 个数)K = 10;end% Get the index for the population sort based on the rank temp,index = sort(intermediate_chromosome(:,M + V + 1);% Now sort the individuals based on the indexfor i = 1 : Nsorted_chromosome(i,:) = intermediate_chromosome(index(i),:);end% Find the maximum rank in the current populationmax_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_rankcurrent_index = max(find(sorted_chromosome(:,M + V + 1) = i);if current_index popremaining = pop - previous_index;temp_pop =.sorted_chromosome(previous_index + 1 : current_index,:); temp_pop_1=temp_pop*(-1);temp_sort_1,temp_sort_index=.sort(temp_pop_1(:, M + V + 2); %在编译的时候出现的问题找到了，主要是对 sort(a,descend)并不支持。temp_sort=temp_sort_1*(-1);for j = 1 : remainingf(previous_index + j,:) = temp_pop(temp_sort_index(j),:);endreturn;elseif current_index 1while isempty(find(candidate(1 : j - 1) = candidate(j)candidate(j) = round(pop*rand(1);if candidate(j) = 0candidate(j) = 1;endendendendfor j = 1 : tour_sizec_obj_rank(j) = chromosome(candidate(j),rank); %提取出阶数 c_obj_distance(j) = chromosome(candidate(j),distance); %提取出距离 endmin_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(l); %层数，拥挤距离都相同，随便选个就可以了endf(i,:) = chromosome(candidate(min_candidate(max_candidate),:);elsef(i,:) = chromosome(candidate(min_candidate(1),:);endend%function f = genetic_operator(parent_chromosome,pro,mu,mum)% This function is utilized to produce offsprings from parent chromosomes.% The genetic operators corssover and mutation which are carried out with% slight modifications from the original design. For more information read% the document enclosed.%input_parameters=inp_para_definition; %定义了一个输入变量的结构，N,M = size(parent_chromosome); % 注意这里的M其实是个垃圾，后面语句会对其重新 赋值利用 switch procase 1M = 2;V = 6;case 2M = 3;V = 12;case 3% case 1和case 2用来对整个算法进行常规验证，作为调试之用；case 3为本工程所需；M = 2; %(output parameters 个数)V = 8; %(input parameters 个数)K = 10;endp = 1;was_crossover = 0;was_mutation = 0;l_limit=0;u_limit=1;for i = 1 : Nif rand(1) 0.9child_1 =;child_2 =;parent_1 = round(N*rand(1);if parent_1 1parent_1 = 1;endparent_2 = round(N*rand(1);if parent_2 1 parent_2 = 1;endwhile isequal(parent_chromosome(parent_l,:),parent_chromosome(parent_2,:) parent_2 = round(N*rand(1);if parent_2 1parent_2 = 1;endend%目的是随机找出两个不同的候选人parent=parent_chromosome(parent_1,:);parent_2 = parent_chromosome(parent_2,:);forj = 1 : V% SBX (Simulated Binary Crossover)% Generate a random numberu(j) = rand(1);if u(j) uimitchild_1(j) = uimit;elseif child_1(j) uimitchild_2(j) = uimit;elseif child_2(j) l_limitchild_2(j) = l_limit;endendchild_1(:,V + 1: M + V) = evaluate_objective(child_1,pro);child_2(:,V + 1: M + V) = evaluate_objective(child_2,pro); was_crossover = 1;was_mutation = 0;elseparent_3 = round(N*rand(l);if parent_3 1 parent_3 = 1;end% Make sure that the mutation does not result in variables out of% the search space. For both the MOPs the range for decision space% is 0,1. In case different variables have different decision% space each variable can be assigned a range.child_3 = parent_chromosome(parent_3,:);forj = 1 : Vr(j) = rand(1);if r(j) u_limitchild_3(j) = u_limit;elseif child_3(j) l_limitchild_3(j) =l_limit;endendchild_3(:,V + 1: M + V) = evaluate_objective(child_3,pro); was_mutation = 1;was_crossover = 0;endif was_crossoverchild(p,:) = child_1;child(p+1,:) = child_2; was_cossover = 0;p = p + 2;elseif was_mutationchild(p,:) = child_3(1,1 : M + V); was_mutation = 0;p = p + 1;endend%以较大概率进行交叉操作，较小的概率进行变异操作f = child;