2019.7.17

很意外本人這篇文章受到很多人的關(guān)注，在此把源碼貼出來(lái)供大家更好的理解學(xué)習(xí)。

https://download.csdn.net/download/joekepler/10590751

========================分割=====================================

本人最近研究NSGA2算法，網(wǎng)上有很多示例代碼，但是基本沒(méi)有注釋，代碼看起來(lái)很頭疼，因此我最近把整個(gè)代碼研讀了一遍，并做上中文注釋，希望可以幫助到一些和我一樣的初學(xué)者們。貼出代碼之前，首先介紹一下NSGA2遺傳算法的流程圖：流程圖中我把每個(gè)詳細(xì)的步驟用號(hào)碼標(biāo)出來(lái)，對(duì)應(yīng)下文的代碼部分。

?

?

?

首先貼出主函數(shù)代碼，對(duì)應(yīng)整個(gè)流程圖：

function nsga_2_optimization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %此處可以更改 %更多機(jī)器學(xué)習(xí)內(nèi)容請(qǐng)?jiān)L問(wèn)omegaxyz.com pop = 200; %種群數(shù)量 gen = 500; %迭代次數(shù) M = 2; %目標(biāo)函數(shù)數(shù)量 V = 30; %維度（決策變量的個(gè)數(shù)） min_range = zeros(1, V); %下界 生成1*30的個(gè)體向量 全為0 max_range = ones(1,V); %上界 生成1*30的個(gè)體向量 全為1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% chromosome = initialize_variables(pop, M, V, min_range, max_range);%初始化種群 chromosome = non_domination_sort_mod(chromosome, M, V);%對(duì)初始化種群進(jìn)行非支配快速排序和擁擠度計(jì)算for i = 1 : genpool = round(pop/2);%round() 四舍五入取整 交配池大小tour = 2;%競(jìng)標(biāo)賽? 參賽選手個(gè)數(shù)parent_chromosome = tournament_selection(chromosome, pool, tour);%競(jìng)標(biāo)賽選擇適合繁殖的父代mu = 20;%交叉和變異算法的分布指數(shù)mum = 20;offspring_chromosome = genetic_operator(parent_chromosome,M, V, mu, mum, min_range, max_range);%進(jìn)行交叉變異產(chǎn)生子代 該代碼中使用模擬二進(jìn)制交叉和多項(xiàng)式變異 采用實(shí)數(shù)編碼[main_pop,~] = size(chromosome);%父代種群的大小[offspring_pop,~] = size(offspring_chromosome);%子代種群的大小clear tempintermediate_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);%對(duì)新的種群進(jìn)行快速非支配排序chromosome = replace_chromosome(intermediate_chromosome, M, V, pop);%選擇合并種群中前N個(gè)優(yōu)先的個(gè)體組成新種群if ~mod(i,100)clc;fprintf('%d generations completed\n',i);end endif M == 2plot(chromosome(:,V + 1),chromosome(:,V + 2),'*');xlabel('f_1'); ylabel('f_2');title('Pareto Optimal Front'); elseif M == 3plot3(chromosome(:,V + 1),chromosome(:,V + 2),chromosome(:,V + 3),'*');xlabel('f_1'); ylabel('f_2'); zlabel('f_3');title('Pareto Optimal Surface'); end ?

1 初始化代碼

function f = initialize_variables(N, M, V, min_range, max_range)%f是一個(gè)由種群個(gè)體組成的矩陣 min = min_range; max = max_range; K = M + V;%%K是數(shù)組的總元素個(gè)數(shù)。為了便于計(jì)算，決策變量和目標(biāo)函數(shù)串在一起形成一個(gè)數(shù)組。 %對(duì)于交叉和變異，利用目標(biāo)變量對(duì)決策變量進(jìn)行選擇 for i = 1 : Nfor j = 1 : Vf(i,j) = min(j) + (max(j) - min(j))*rand(1);%f(i j)表示的是種群中第i個(gè)個(gè)體中的第j個(gè)決策變量，%這行代碼為每個(gè)個(gè)體的所有決策變量在約束條件內(nèi)隨機(jī)取值endf(i,V + 1: K) = evaluate_objective(f(i,:), M, V); % M是目標(biāo)函數(shù)數(shù)量 V是決策變量個(gè)數(shù)%為了簡(jiǎn)化計(jì)算將對(duì)應(yīng)的目標(biāo)函數(shù)值儲(chǔ)存在染色體的V + 1 到 K的位置。 end

2 快速非支配排序和擁擠度計(jì)算代碼

%% 對(duì)初始種群開(kāi)始排序 快速非支配排序 % 使用非支配排序?qū)ΨN群進(jìn)行排序。該函數(shù)返回每個(gè)個(gè)體對(duì)應(yīng)的排序值和擁擠距離，是一個(gè)兩列的矩陣。 % 并將排序值和擁擠距離添加到染色體矩陣中 function f = non_domination_sort_mod(x, M, V) [N, ~] = size(x);% N為矩陣x的行數(shù)，也是種群的數(shù)量 clear m front = 1; F(front).f = []; individual = [];for i = 1 : Nindividual(i).n = 0;%n是個(gè)體i被支配的個(gè)體數(shù)量individual(i).p = [];%p是被個(gè)體i支配的個(gè)體集合for j = 1 : Ndom_less = 0;dom_equal = 0;dom_more = 0;for k = 1 : M %判斷個(gè)體i和個(gè)體j的支配關(guān)系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;elsedom_more = dom_more + 1;endendif dom_less == 0 && dom_equal ~= M % 說(shuō)明i受j支配，相應(yīng)的n加1individual(i).n = individual(i).n + 1;elseif dom_more == 0 && dom_equal ~= M % 說(shuō)明i支配j,把j加入i的支配合集中individual(i).p = [individual(i).p j];endend if individual(i).n == 0 %個(gè)體i非支配等級(jí)排序最高，屬于當(dāng)前最優(yōu)解集，相應(yīng)的染色體中攜帶代表排序數(shù)的信息x(i,M + V + 1) = 1;F(front).f = [F(front).f i];%等級(jí)為1的非支配解集end end %上面的代碼是為了找出等級(jí)最高的非支配解集 %下面的代碼是為了給其他個(gè)體進(jìn)行分級(jí) while ~isempty(F(front).f)Q = []; %存放下一個(gè)front集合for i = 1 : length(F(front).f)%循環(huán)當(dāng)前支配解集中的個(gè)體if ~isempty(individual(F(front).f(i)).p)%個(gè)體i有自己所支配的解集for j = 1 : length(individual(F(front).f(i)).p)%循環(huán)個(gè)體i所支配解集中的個(gè)體individual(individual(F(front).f(i)).p(j)).n = ...%...表示的是與下一行代碼是相連的， 這里表示個(gè)體j的被支配個(gè)數(shù)減1individual(individual(F(front).f(i)).p(j)).n - 1;if individual(individual(F(front).f(i)).p(j)).n == 0% 如果q是非支配解集，則放入集合Q中x(individual(F(front).f(i)).p(j),M + V + 1) = ...%個(gè)體染色體中加入分級(jí)信息front + 1;Q = [Q individual(F(front).f(i)).p(j)];endendendendfront = front + 1;F(front).f = Q; end[temp,index_of_fronts] = sort(x(:,M + V + 1));%對(duì)個(gè)體的代表排序等級(jí)的列向量進(jìn)行升序排序 index_of_fronts表示排序后的值對(duì)應(yīng)原來(lái)的索引 for i = 1 : length(index_of_fronts)sorted_based_on_front(i,:) = x(index_of_fronts(i),:);%sorted_based_on_front中存放的是x矩陣按照排序等級(jí)升序排序后的矩陣 end current_index = 0;%% Crowding distance 計(jì)算每個(gè)個(gè)體的擁擠度f(wàn)or front = 1 : (length(F) - 1)%這里減1是因?yàn)榇a55行這里，F的最后一個(gè)元素為空，這樣才能跳出循環(huán)。所以一共有l(wèi)ength-1個(gè)排序等級(jí)distance = 0;y = [];previous_index = current_index + 1;for i = 1 : length(F(front).f)y(i,:) = sorted_based_on_front(current_index + i,:);%y中存放的是排序等級(jí)為front的集合矩陣endcurrent_index = current_index + i;%current_index =isorted_based_on_objective = [];%存放基于擁擠距離排序的矩陣for i = 1 : M[sorted_based_on_objective, index_of_objectives] = ...sort(y(:,V + i));%按照目標(biāo)函數(shù)值排序sorted_based_on_objective = [];for j = 1 : length(index_of_objectives)sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);% sorted_based_on_objective存放按照目標(biāo)函數(shù)值排序后的x矩陣endf_max = ...sorted_based_on_objective(length(index_of_objectives), V + i);%fmax為目標(biāo)函數(shù)最大值 fmin為目標(biāo)函數(shù)最小值f_min = sorted_based_on_objective(1, V + i);y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i)...%對(duì)排序后的第一個(gè)個(gè)體和最后一個(gè)個(gè)體的距離設(shè)為無(wú)窮大= Inf;y(index_of_objectives(1),M + V + 1 + i) = Inf;for j = 2 : length(index_of_objectives) - 1%循環(huán)集合中除了第一個(gè)和最后一個(gè)的個(gè)體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;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; end f = z();%得到的是已經(jīng)包含等級(jí)和擁擠度的種群矩陣 并且已經(jīng)按等級(jí)排序排序

3 競(jìng)標(biāo)賽選擇代碼

function f = tournament_selection(chromosome, pool_size, tour_size) [pop, variables] = size(chromosome);%獲得種群的個(gè)體數(shù)量和決策變量數(shù)量 rank = variables - 1;%個(gè)體向量中排序值所在位置 distance = variables;%個(gè)體向量中擁擠度所在位置 %競(jìng)標(biāo)賽選擇法，每次隨機(jī)選擇兩個(gè)個(gè)體，優(yōu)先選擇排序等級(jí)高的個(gè)體，如果排序等級(jí)一樣，優(yōu)選選擇擁擠度大的個(gè)體 for i = 1 : pool_sizefor j = 1 : tour_sizecandidate(j) = round(pop*rand(1));%隨機(jī)選擇參賽個(gè)體if candidate(j) == 0candidate(j) = 1;endif j > 1while ~isempty(find(candidate(1 : j - 1) == candidate(j)))%防止兩個(gè)參賽個(gè)體是同一個(gè)candidate(j) = round(pop*rand(1));if candidate(j) == 0candidate(j) = 1;endendendendfor j = 1 : tour_size% 記錄每個(gè)參賽者的排序等級(jí) 擁擠度c_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));%選擇排序等級(jí)較小的參賽者，find返回該參賽者的索引if length(min_candidate) ~= 1%如果兩個(gè)參賽者的排序等級(jí)相等 則繼續(xù)比較擁擠度 優(yōu)先選擇擁擠度大的個(gè)體max_candidate = ...find(c_obj_distance(min_candidate) == max(c_obj_distance(min_candidate)));if length(max_candidate) ~= 1max_candidate = max_candidate(1);endf(i,:) = chromosome(candidate(min_candidate(max_candidate)),:);elsef(i,:) = chromosome(candidate(min_candidate(1)),:);end end

4、5 交叉 變異代碼

function f = genetic_operator(parent_chromosome, M, V, mu, mum, l_limit, u_limit) [N,m] = size(parent_chromosome);%N是交配池中的個(gè)體數(shù)量clear m p = 1; was_crossover = 0;%是否交叉標(biāo)志位 was_mutation = 0;%是否變異標(biāo)志位for i = 1 : N%這里雖然循環(huán)N次，但是每次循環(huán)都會(huì)有概率產(chǎn)生2個(gè)或者1個(gè)子代，所以最終產(chǎn)生的子代個(gè)體數(shù)量大約是2N個(gè)if rand(1) < 0.9%交叉概率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 < 1parent_2 = 1;endwhile isequal(parent_chromosome(parent_1,:),parent_chromosome(parent_2,:))parent_2 = round(N*rand(1));if parent_2 < 1parent_2 = 1;endendparent_1 = parent_chromosome(parent_1,:);parent_2 = parent_chromosome(parent_2,:);for j = 1 : Vu(j) = rand(1);if u(j) <= 0.5bq(j) = (2*u(j))^(1/(mu+1));elsebq(j) = (1/(2*(1 - u(j))))^(1/(mu+1));endchild_1(j) = ...0.5*(((1 + bq(j))*parent_1(j)) + (1 - bq(j))*parent_2(j));child_2(j) = ...0.5*(((1 - bq(j))*parent_1(j)) + (1 + bq(j))*parent_2(j));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);endif 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);endendchild_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;else%if >0.9parent_3 = round(N*rand(1));if parent_3 < 1parent_3 = 1;endchild_3 = parent_chromosome(parent_3,:);for j = 1 : Vr(j) = rand(1);if r(j) < 0.5delta(j) = (2*r(j))^(1/(mum+1)) - 1;elsedelta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1));endchild_3(j) = child_3(j) + delta(j);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);endend child_3(:,V + 1: M + V) = evaluate_objective(child_3, M, V);was_mutation = 1;was_crossover = 0;end% if <0.9if 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;end endf = child;

交叉算法選擇的是模擬二進(jìn)制交叉，變異算法選擇的是多項(xiàng)式變異， 算法的具體過(guò)程大家可以在網(wǎng)上查閱一下

8 生成新的種群（精英策略）

function f = replace_chromosome(intermediate_chromosome, M, V,pop)%精英選擇策略[N, m] = size(intermediate_chromosome); [temp,index] = sort(intermediate_chromosome(:,M + V + 1));clear temp m for i = 1 : Nsorted_chromosome(i,:) = intermediate_chromosome(index(i),:); endmax_rank = max(intermediate_chromosome(:,M + V + 1));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_sort,temp_sort_index] = ...sort(temp_pop(:, M + V + 2),'descend');for j = 1 : remainingf(previous_index + j,:) = temp_pop(temp_sort_index(j),:);endreturn;elseif current_index < popf(previous_index + 1 : current_index, :) = ...sorted_chromosome(previous_index + 1 : current_index, :);elsef(previous_index + 1 : current_index, :) = ...sorted_chromosome(previous_index + 1 : current_index, :);return;endprevious_index = current_index; end ?

本例子選擇的測(cè)試函數(shù)是ZDT1， 目標(biāo)評(píng)價(jià)函數(shù)如下：

function f = evaluate_objective(x, M, V)%%計(jì)算每個(gè)個(gè)體的M個(gè)目標(biāo)函數(shù)值 f = []; f(1) = x(1); g = 1; sum = 0; for i = 2:Vsum = sum + x(i); end sum = 9*(sum / (V-1)); g = g + sum; f(2) = g * (1 - sqrt(x(1) / g)); end ?

經(jīng)歷500次迭代后的pareto最優(yōu)解集：

總結(jié)

以上是生活随笔為你收集整理的NSGA2 算法MATLAB完整代码 中文注释详解的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。