惩罚函数方法介绍和简单代码实现

avatar 2024年06月17日20:34:57 0 1294 views
博主分享免费Java教学视频，B站账号：Java刘哥 ，长期提供技术问题解决、项目定制：本站商品点此

惩罚函数方法通过对不可行解进行惩罚，将约束优化问题转化为无约束优化问题，从而简化约束处理过程。

惩罚函数方法可以通过静态、动态和自适应策略来平衡目标函数和约束之间的关系。

一、简单原理

简单来说，将原始目标函数和惩罚项组合成适应度函数 F(x)，形式如下： 

1）就是在求目标值(适应度)的时候，同时求一下这个个体是否违反约束，如果违反约束，求一下违反多少

  1. function penalty = evaluatePenalty(x, penalty_factor)
  2.     % 计算罚函数值
  3.     % 初始化罚函数值
  4.     penalty = 0;
  5.     % 约束1：g1(x) = x1^2 + x2^2 - 1 <= 0（圆内）
  6.     g1 = x(1)^2 + x(2)^2 - 1;
  7.     % 约束2：g2(x) = x1 + x2 - 1.5 <= 0（直线下方）
  8.     g2 = x(1) + x(2) - 1.5;
  9.     % 如果违反约束，则累加罚函数值。penalty_factor这里设置的是1000
  10.     if g1 > 0
  11.         penalty = penalty + penalty_factor * g1;
  12.     end
  13.     if g2 > 0
  14.         penalty = penalty + penalty_factor * g2;
  15.     end
  16. end

这里我们的罚函数比较简单，就是 罚函数值 =  约束违反程度 * 罚系数

 

2）求下一代种群时，通常是取目标值(适应度)最小的前pop_size个

这里我们需要对【目标值+罚函数值】一起排序

  1. function [next_pop, next_obj, next_penalty] = selectNextGeneration(
  2.     combined_pop, combined_obj, combined_penalty, pop_size)
  3.     % 选择下一代种群
  4.     % 初始化下一代种群和目标函数值、罚函数值矩阵
  5.     next_pop = [];
  6.     next_obj = [];
  7.     next_penalty = [];
  8.     
  9.     % 计算每个个体的总目标值（目标函数值加上罚函数值）
  10.     total_obj = sum(combined_obj, 2) + combined_penalty;
  11.     
  12.     % 对总目标值进行排序，选择前pop_size个个体
  13.     [~, sorted_indices] = sort(total_obj); % 升序排序
  14.     
  15.     % 选择前pop_size个个体作为下一代
  16.     selected_indices = sorted_indices(1:pop_size);      % 选择排序后，目标值值最小的前pop_size个
  17.     next_pop = combined_pop(selected_indices, :);       % 获得新种群
  18.     next_obj = combined_obj(selected_indices, :);       % 求目标值(直接从合并里取，下面)
  19.     next_penalty = combined_penalty(selected_indices);  % 求罚函数值
  20. end

 

二、罚函数类型

惩罚函数方法主要有以下几种类型：

  1. 固定惩罚函数（Static Penalty Function）：

    • 惩罚系数在优化过程中保持不变。
    • 优点：实现简单，计算量小。
    • 缺点：惩罚系数难以选择，容易导致搜索过程陷入局部最优解或无法收敛。

  2. 动态惩罚函数（Dynamic Penalty Function）：

    • 惩罚系数在优化过程中逐渐增大或减小。
    • 优点：通过动态调整惩罚系数，能够更好地平衡探索与开发，提高搜索效果。
    • 缺点：需要设定合适的动态调整策略，增加了参数选择的复杂性。

  3. 自适应惩罚函数（Adaptive Penalty Function）：

    • 惩罚系数根据个体的适应度或种群的进化状态自适应调整。
    • 优点：能够根据优化过程中的情况自适应调整，效果较好。
    • 缺点：实现复杂，需要根据具体问题设定合适的自适应策略。

  4. 非线性惩罚函数（Nonlinear Penalty Function）：

    • 惩罚项为非线性函数，如平方惩罚或指数惩罚。
    • 优点：能够更灵活地处理违反约束的解，效果较好。
    • 缺点：需要选择合适的非线性函数形式，增加了参数选择的复杂性。

 

三、完整代码

  1. function PenaltyFunction_MOO()
  2.     % 参数设置
  3.     pop_size = 100; % 种群规模
  4.     num_generations = 250; % 进化代数
  5.     num_vars = 2; % 决策变量数量
  6.     num_obj = 2; % 目标函数数量
  7.     penalty_factor = 1000; % 罚系数
  8.     lb = [-1, -1]; % 决策变量下界
  9.     ub = [2, 2]; % 决策变量上界
  11.     % 初始化种群
  12.     Population = initPopulation(pop_size, num_vars, lb, ub);
  14.     % 评估种群目标函数值和罚函数值
  15.     [Obj, Penalty] = evaluatePopulationWithPenalty(Population, num_obj, penalty_factor);
  17.     % 进化过程
  18.     for gen = 1:num_generations
  19.         % 选择父代
  20.         Parents = tournamentSelection(Population, Obj, Penalty, pop_size);
  22.         % 交叉和变异生成子代
  23.         Offspring = crossoverAndMutation(Parents, lb, ub);
  25.         % 评估子代目标函数值和罚函数值
  26.         [OffspringObj, OffspringPenalty] = evaluatePopulationWithPenalty(Offspring, num_obj, penalty_factor);
  28.         % 合并种群和子代
  29.         CombinedPop = [Population; Offspring];
  30.         CombinedObj = [Obj; OffspringObj];
  31.         CombinedPenalty = [Penalty; OffspringPenalty];
  33.         % 选择前 N 个个体作为下一代种群
  34.         [Population, Obj, Penalty] = selectNextGeneration(CombinedPop, CombinedObj, CombinedPenalty, pop_size);
  35.     end
  37.     % 绘制Pareto前沿
  38.     plotParetoFront(Obj);
  39. end
  41. function pop = initPopulation(pop_size, num_vars, lb, ub)
  42.     % 初始化种群
  43.     % 使用随机数生成初始种群，每个个体的变量值在上下界之间
  44.     pop = rand(pop_size, num_vars) .* (ub - lb) + lb;
  45. end
  47. function [obj, penalty] = evaluatePopulationWithPenalty(pop, num_obj, penalty_factor)
  48.     % 评估种群目标函数值和罚函数值
  49.     % 初始化目标函数值和罚函数值矩阵
  50.     obj = zeros(size(pop, 1), num_obj);
  51.     penalty = zeros(size(pop, 1), 1);
  52.     for i = 1:size(pop, 1)
  53.         % 计算每个个体的目标函数值
  54.         obj(i, :) = evaluateIndividual(pop(i, :), num_obj);
  55.         % 计算每个个体的罚函数值
  56.         penalty(i) = evaluatePenalty(pop(i, :), penalty_factor);
  57.     end
  58. end
  60. function obj = evaluateIndividual(x, num_obj)
  61.     % 目标函数
  62.     % 初始化目标函数值
  63.     obj = zeros(1, num_obj);
  64.     % 目标函数1：f1(x) = x1^2 + x2^2
  65.     obj(1) = x(1)^2 + x(2)^2;
  66.     % 目标函数2：f2(x) = (x1 - 1)^2 + x2^2
  67.     obj(2) = (x(1) - 1)^2 + x(2)^2;
  68. end
  70. function penalty = evaluatePenalty(x, penalty_factor)
  71.     % 计算罚函数值
  72.     % 初始化罚函数值
  73.     penalty = 0;
  74.     % 约束1：g1(x) = x1^2 + x2^2 - 1 <= 0（圆内）
  75.     g1 = x(1)^2 + x(2)^2 - 1;
  76.     % 约束2：g2(x) = x1 + x2 - 1.5 <= 0（直线下方）
  77.     g2 = x(1) + x(2) - 1.5;
  78.     % 如果违反约束，则累加罚函数值
  79.     if g1 > 0
  80.         penalty = penalty + penalty_factor * g1;
  81.     end
  82.     if g2 > 0
  83.         penalty = penalty + penalty_factor * g2;
  84.     end
  85. end
  87. function parents = tournamentSelection(pop, obj, penalty, pool_size)
  88.     % 锦标赛选择
  89.     % 初始化父代矩阵
  90.     parents = zeros(pool_size, size(pop, 2));
  91.     for i = 1:pool_size
  92.         % 随机选择两个个体
  93.         a = randi(size(pop, 1));
  94.         b = randi(size(pop, 1));
  95.         % 比较两个个体，选择较优者
  96.         if better(a, b, obj, penalty)
  97.             winner = a;
  98.         else
  99.             winner = b;
  100.         end
  101.         % 将较优个体加入父代
  102.         parents(i, :) = pop(winner, :);
  103.     end
  104. end
  106. function better_flag = better(a, b, obj, penalty)
  107.     % 比较两个个体，判断哪个更好
  108.     % 如果两个个体的罚函数值相同，比较目标函数值
  109.     if penalty(a) == penalty(b)
  110.         better_flag = dominates(obj(a, :), obj(b, :));
  111.     else
  112.         % 否则，选择罚函数值较小的个体
  113.         better_flag = penalty(a) < penalty(b);
  114.     end
  115. end
  117. function d = dominates(obj1, obj2)
  118.     % 支配关系判断
  119.     % 如果obj1在所有目标上不劣于obj2，并且在至少一个目标上优于obj2
  120.     d = all(obj1 <= obj2) && any(obj1 < obj2);
  121. end
  123. function offspring = crossoverAndMutation(parents, lb, ub)
  124.     % 交叉和变异操作
  125.     % 初始化子代矩阵
  126.     offspring = zeros(size(parents));
  127.     for i = 1:size(parents, 1)
  128.         % 选择两个父代个体
  129.         parent1 = parents(i, :);
  130.         parent2 = parents(mod(i, size(parents, 1)) + 1, :);
  131.         % 均匀交叉，生成子代
  132.         child = 0.5 * (parent1 + parent2);
  133.         % 变异操作
  134.         mutation_rate = 1 / length(lb);
  135.         for j = 1:length(lb)
  136.             if rand() < mutation_rate
  137.                 % 对子代的变量进行变异
  138.                 child(j) = lb(j) + (ub(j) - lb(j)) * rand();
  139.             end
  140.         end
  141.         % 将生成的子代加入子代矩阵
  142.         offspring(i, :) = child;
  143.     end
  144. end
  146. function [next_pop, next_obj, next_penalty] = selectNextGeneration(combined_pop, combined_obj, combined_penalty, pop_size)
  147.     % 选择下一代种群
  148.     % 初始化下一代种群和目标函数值、罚函数值矩阵
  149.     next_pop = [];
  150.     next_obj = [];
  151.     next_penalty = [];
  152.     
  153.     % 计算每个个体的总目标值（目标函数值加上罚函数值）
  154.     total_obj = sum(combined_obj, 2) + combined_penalty;
  155.     
  156.     % 对总目标值进行排序，选择前pop_size个个体
  157.     [~, sorted_indices] = sort(total_obj); % 升序排序
  158.     
  159.     % 选择前pop_size个个体作为下一代
  160.     selected_indices = sorted_indices(1:pop_size);
  161.     next_pop = combined_pop(selected_indices, :);
  162.     next_obj = combined_obj(selected_indices, :);
  163.     next_penalty = combined_penalty(selected_indices);
  165.     % % 选择前pop_size个个体作为下一代
  166.     % for i = 1:pop_size
  167.     %     next_pop = [next_pop; combined_pop(sorted_indices(i), :)];
  168.     %     next_obj = [next_obj; combined_obj(sorted_indices(i), :)];
  169.     %     next_penalty = [next_penalty; combined_penalty(sorted_indices(i))];
  170.     % end
  171. end
  173. function plotParetoFront(obj)
  174.     % 绘制Pareto前沿
  175.     figure;
  176.     scatter(obj(:, 1), obj(:, 2), 'filled');
  177.     title('Pareto Front');
  178.     xlabel('Objective 1');
  179.     ylabel('Objective 2');
  180.     grid on;
  181. end

 

历史上的今天
六月
17日
  • 微信
  • 交流学习，资料分享
  • weinxin
  • 个人淘宝
  • 店铺名：言曌博客咨询部

  • (部分商品未及时上架淘宝)
avatar

发表评论

avatar 登录者：匿名
匿名评论，评论回复后会有邮件通知

  

已通过评论：0   待审核评论数：0