多目标优化（三）简单的 MOEA/D

多目标优化（三）简单的 MOEA/D

写在前面:

① MOEA/D提供了一个简单但是有效的方法，将分解的方法引入到多目标进化计算中。

② 分解策略主要使用三种聚合函数：1）权重聚合方法；2）切比雪夫方法；3）基于惩罚的边界交叉方法

算法概要：

step1：获取算法参数（包括种群大小、迭代终止条件、函数问题参数、生成权重向量、生成邻居矩阵）

step2：初始化种群（生成NP个初始解，即每个子问题的初始解，每个子问题对应一个权重向量）

step3：生成候选解，更新各目标的最小值（for 最小值问题）

step4：按照MOEA/D的分解策略更新各个子问题的最优解

step5：重复上述过程，直至满足迭代终止条件

clc; clear; close all;
format long
global NP x_num f_num nei_num iter_max fun x_min x_max
global weight_vec f_min  nei pop

%% 获取函数参数
fun = 'ZDT1'; test              % 选择测试函数 DTLZ1
weight_vec = gen_weight_vec();  % 产生权重向量
nei_num = NP / 5;               % 邻居规模大小
nei = gen_nei();                % 生成邻居矩阵
iter_max =  250;                % 最大迭代次数
NP = 100;                      % 种群大小，也是子问题 subproblem 的个数

%% 初始种群
pop = x_min + (x_max - x_min) .* rand(NP, x_num);
for i = 1:NP
    pop(i, x_num+1:x_num+f_num) = object_fun(pop);
end
f_min = min(pop(:,x_num+1:x_num+f_num));

%% 迭代循环
for iter = 1:iter_max
    
    % 生成候选种群
    cpop = gen_candidate();
    f_min_y = min(cpop(:, x_num+1:x_num+f_num));
    f_min = min(f_min, f_min_y);
    
    % MOEA/D 分解策略
    % pop = weighted(cpop);
  	% pop = tchebycheff(cpop);
    pop = penalty_boundary_intersection(cpop);
    
    % 可视化
    h = scatter(pop(:,x_num+1),pop(:,x_num+2),'*');
    title(num2str(iter))
    drawnow
    if iter < iter_max
        delete(h);
    end
    
end

子程序：

1、gen_weight_vec：与NSGA-III中产生参考点的方法一样，这里产生的权重向量个数与种群大小相同，意思是每个个体看成是一个子问题，每个子问题对应一个权重向量。

https://blog.csdn.net/weixin_41942547/article/details/114082875?spm=1001.2014.3001.5502

2、gen_nei：对于某个子问题的权重向量，找到离它最近的nei_num个权重向量（欧氏距离）

function Nei_pick = gen_nei()

    global weight_vec nei_num NP
    dis = pdist2(weight_vec, weight_vec);
    
    for i = 1:NP
        
        [~, idx] = sort(dis(i, :));
        Nei_pick(i, :) = idx(1: nei_num); 
        
    end
    
end

3、gen_candidate：个人认为能生成候选解就行了。这里用的方法是从该子问题的邻居中选择两个进行交叉变异（这样生成的候选解被“认为”是在子问题附近的，子问题更有机会得到更新）

% 个体更新
function cpop = gen_candidate()

    global x_num pop NP nei nei_num
    pc = 0.7; pm = 1 / x_num; yita1 = 2; yita2 = 5; cpop = pop*0;
    
    for i = 1:NP
        
        pick = randperm(nei_num, 2); % 从邻居中随机挑选两个个体
        x1 = pop(nei(i, pick(1)), :); x2 = pop(nei(i, pick(2)), :); cpop(i, :) = x1;
        % crossover
        for j = 1:x_num
            if rand < pc 
                u = rand;
                if u < 0.5
                    gama = (2 * u)^(1 / (yita1 + 1));
                else
                    gama = (1 / (2 * (1 - u)))^(1 / (yita1+1));
                end
                cpop(i, j) = 0.5 * ((1 + gama) * x1(j) + (1 - gama) * x2(j));
            end
        end
        
        % mutation
        for j = 1:x_num
            if rand < pm
                u = rand;
                if u < 0.5
                    delta = (2 * u)^(1 / (yita2 + 1)) - 1;
                else
                    delta = 1 - (2 * (1 - u))^(1 / (yita2 + 1));
                end
                cpop(i, j) = cpop(i, j) + delta;
            end
        end
        
    end
    
    % 判界
    cpop = judge_bound(cpop);
    
end

4、MOEA/D 分解策略：公式参考https://blog.csdn.net/sinat_33231573/article/details/80271801

1）、Weighted Sum Approach

加权求和法，是那个人为设定了各个目标所占的比重，然后将多目标转化为单目标的方法吧！

这里应该更高级一点，设置的比重就是权重向量，对于每个子问题，候选解加权求和小于原子问题解的加权求和的话，就能更新这个子问题的解了。

function rpop = weighted(cpop)
    
    global weight_vec x_num f_num f_min nei pop NP
    rpop = pop;
    
    % 计算加权求和后的值 
    for i = 1:NP
        
        y = cpop(i, :); 
        for m = 1:size(nei, 2)

            weighted_x = sum((weight_vec(nei(i, m), :)) .* (pop(nei(i, m), x_num+1:x_num+f_num) - f_min));
            weighted_y = sum((weight_vec(nei(i, m), :)) .* (y(x_num+1:x_num+f_num) - f_min));
     	
            if weighted_y <= weighted_x
                rpop(nei(i, m), :) = y;
            end
            
        end
        
    end
    
end
        

2）、Tchebycheff Approach

https://blog.csdn.net/jinjiahao5299/article/details/76045936

该方法大致思想就是减少最大差距从而将个体逼近pareto前言面。

function rpop = tchebycheff(cpop)
    
    global weight_vec x_num f_num f_min nei pop NP
    rpop = pop;
    
    % 计算切比雪夫值 
    for i = 1:NP
        
        y = cpop(i, :); 
        for m = 1:size(nei, 2)

            tche_cheff_x = max(weight_vec(nei(i, m), :) .* abs(pop(nei(i, m), x_num+1:x_num+f_num) - f_min));
            tche_cheff_y = max(weight_vec(nei(i, m), :) .* abs(y(x_num+1:x_num+f_num) - f_min));
     	
            if tche_cheff_y <= tche_cheff_x
                rpop(nei(i, m), :) = y;
            end
            
        end
        
    end
    
end

3）、Boundary Intersection Approach

子问题的解可以在空间上的任意位置，那么它就有两个因素来评价解的好坏：1）解到权重向量的距离越小越好；2）解在权重向量上的映射点离最优点（所有目标都最小）越小越好。

这里解到权重向量的惩罚因子设置为 2）的值，效果更好点儿！！

function rpop = penalty_boundary_intersection(cpop)
    
    global weight_vec x_num f_num f_min nei pop NP
    rpop = pop;
    
    % 计算两个距离值 
    for i = 1:NP
        
        y = cpop(i, :); 
        % penalty = mean(pop(i, x_num+1:x_num+f_num));
        for m = 1:size(nei, 2)
            
            fx = pop(nei(i, m), x_num+1:x_num+f_num); w = weight_vec(nei(i, m), :);
            boundary_x = norm((f_min - fx) * w') / norm(w);
            boundary_x = boundary_x + boundary_x * norm(fx - (f_min - boundary_x * w));
            boundary_y = norm((f_min - y(x_num+1:x_num+f_num)) * w') / norm(w);
            boundary_y = boundary_y + boundary_y * norm(y(x_num+1:x_num+f_num) - (f_min - boundary_y * w));
            
            if boundary_y <= boundary_x
                rpop(nei(i, m), :) = y;
            end
            
        end
        
    end
    
end
        

完整可执行代码：https://download.csdn.net/download/weixin_41942547/15479604