算法-多目标优化-遗传算法

多目标优化

这段时间接触了不少多目标优化方面的知识，在此做一些粗浅的总结。

多目标优化（MOO，即Multi-Objective Optimization）是指在一个优化问题中，需要同时考虑多个目标函数，即多个目标变量。

在众多目标函数中，某些目标函数之间会存在“互斥关系”。假设所有的目标函数都是最小化问题，且目前存在两种解$x_1，x_2$，$f_1(x_1)$的值低于$f_1(x_2)$，然而$f_2(x_1)$的值高于$f_2(x_2)$，我们期望存在一种解$x$，能够同时做到$f_1(x)$最小化，且$f_2(x)$最小化，然后在现实世界中呈现的例子里大多数都是$x_1$和$x_2$这样的互斥解。

例如，制作某一商品，成本较低的原料制作的成品价格较低，而成本较高的原料制作的成品价格则较高，成本与成品价格之间就存在这样的互斥关系，想要获得一个成本最低且成品价格颇高的商品几乎是很难的。那么在这样一个互斥关系中，如何寻求一个相对合适的解呢？这就是多目标优化问题想要解决的核心问题。

NSGA

NSGA（Non-dominated Sorting Genetic Algorithm）即非支配排序遗传算法。简单来说，遗传算法就是遵循了生物种群的自然选择，也就是优胜劣汰。我们可以把想要得到的解形象化为一个生物种群，对于适应生存环境的种群个体而言，他们显然具有生存繁衍的能力，反之，不适应生存环境的个体就会被淘汰。对于适者，他们的后代通过两两交配遗传了父母双方的优良基因（基因交叉），同时在基因交叉后，后代的基因也会产生变异，从而成为新的个体，经过一轮又一轮的生存繁衍，基因通过选择、淘汰、交叉、变异等操作，最终形成了能够稳定适应生存环境的强大种群，种群中每一个个体都具备了生存繁衍的能力，我们称此时的种群中的每一个个体为有效解。

NSGA-II

NSGA-II是NSGA的改进版本。

为什么要改进NSGA呢？

• NSGA的复杂度较高：$O(MN^3)$。其中，M为目标函数数量，N为种群内的个体数量。
• NSGA没有精英策略
• 需要指定共享半径

随便找了一篇利用多目标优化的论文来进行了复现，A novel multi-objective optimization method for the pressurized reservoir in hydraulic robotics，虽然没有源码，但是应该复现起来也不算很难。这里使用的是pymoo框架。

从论文里可以发现，共有六个变量、六个目标函数、七个约束。

六个变量分别为：$d_r(60~90),l_r(140~200),l_s(180~400),d_c(30~84),d_w(2.5~8.5),n_a(8~25，论文后续又改为了6~12)$

六个目标函数为：

• $k_s * (l_s - l_r - l_p) / a_r$（原论文中公式错误）
• $p_v = k_s * v_n / a_r^2$（原论文中公式错误）
• $v_m = a_r * (l_r - l_p - d_w*(2 + n_a)) / 1000$（原论文中公式错误）
• $m_r = (m_h + m_p + m_s + m_c) / 1000$（原论文中公式错误）
• $v_r = l_r * d_r^2 / 1000$（原论文中公式错误）
• f_a = p_i * a_r

其中:

• $k_s = (G_ * d_w^4) / (8 * d_c^3 * n_a)$
• $a_r = \pi * d_r^2 / 4$
• $m_h = \pi * d_r * l_r * l_w * rho_h * 1000$（原论文中公式错误）
• $m_p = \pi * d_r^2 * l_t * rho_p * 1000 / 4$（原论文中公式错误）
• $m_s = \pi * d_w^2 * rho_s * \sqrt{((\pi * d_c * n_a)^2 + l_s^2)} * 1000 / 4$（原论文中公式错误）
• $m_c = \pi * d_r^2 * l_c * rho_c * 1000 / 2$。（原论文中公式错误）

七个约束分别为：

• $p_i \geq 0.1$
• $p_v \leq 0.5$
• $v_r \leq 1200$（论文在约束提出后经过论证又改为了1000）
• $m_r \leq 1000$
• $v_m \geq 200$（论文在约束提出后经过论证又改为了450）
• $f_a \leq 1000$
• $r_h \leq 2.4$

其中，$r_h = l_r / d_r$

世界就是一个巨大的草台班子，这论文怎么发的出来的。

算了，把复现的源码也一并贴出来吧。

import numpy as np
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

class HydraulicTankProblem(Problem):
    def __init__(self):
        n_var = 6
        xl = np.array([60, 140, 180, 30, 2.5, 6])
        xu = np.array([90, 200, 400, 84, 8.5, 12])
        n_obj = 6
        n_constr = 7
        super().__init__(n_var=n_var, n_obj=n_obj, n_constr=n_constr, xl=xl, xu=xu)

    def _evaluate(self, X, out, *args, **kwargs):
        G_ = 79000
        rho_c = 7.85 * 0.001  # 将密度单位转换为 g/mm^3
        rho_s = 7.85 * 0.001
        rho_h = 2.78 * 0.001
        rho_p = 2.78 * 0.001
        l_w = 3.0
        l_c = 3.0
        l_p = 15.0
        l_t = 6.0
        v_n = 350 * 1000  # 将体积单位转换为 mm^3

        n_pop = X.shape[0]
        F = np.zeros((n_pop, self.n_obj))
        G = np.zeros((n_pop, self.n_constr))

        for i in range(n_pop):
            d_r = X[i, 0]
            l_r = X[i, 1]
            l_s = X[i, 2]
            d_c = X[i, 3]
            d_w = X[i, 4]
            n_a = X[i, 5]

            k_s = (G_ * d_w**4) / (8 * d_c**3 * n_a)
            a_r = np.pi * d_r**2 / 4
            p_i = k_s * (l_s - l_r - l_p) / a_r
            p_v = k_s * v_n / a_r**2
            v_r = l_r * d_r ** 2 / 1000  # 将空间占用单位转换为 ml
            m_h = np.pi * d_r * l_r * l_w * rho_h * 1000  # 将质量单位转换为 g
            m_p = np.pi * d_r**2 * l_t * rho_p * 1000 / 4
            m_s = np.pi * d_w**2 * rho_s * np.sqrt((np.pi * d_c * n_a)**2 + l_s**2) * 1000 / 4
            m_c = np.pi * d_r**2 * l_c * rho_c * 1000 / 2
            m_r = (m_h + m_p + m_s + m_c) / 1000
            v_m = a_r * (l_r - l_p - d_w*(2 + n_a)) / 1000  # 将最大体积单位转换为 ml
            f_a = p_i * a_r

            r_h = l_r / d_r

            F[i, 0] = -p_i
            F[i, 1] = p_v
            F[i, 2] = v_r
            F[i, 3] = m_r
            F[i, 4] = -v_m
            F[i, 5] = f_a

            G[i, 0] = 0.1 - p_i if p_i < 0.1 else 0
            G[i, 1] = p_v - 0.5 if p_v > 0.5 else 0
            G[i, 2] = v_r - 1000 if v_r > 1000 else 0
            G[i, 3] = m_r - 1000 if m_r > 1000 else 0
            G[i, 4] = 450 - v_m if v_m < 450 else 0
            G[i, 5] = f_a - 1000 if f_a > 1000 else 0
            G[i, 6] = r_h - 2.4 if r_h > 2.4 else 0

        out["F"] = F
        out["G"] = G

problem = HydraulicTankProblem()
algorithm = NSGA2(
    pop_size=120,
    eliminate_duplicates=True
)
res = minimize(
    problem,
    algorithm,
    ('n_gen', 60),
    seed=1,
    verbose=False
)

print(res.F)
plot = Scatter()
plot.add(res.F, marker='.', facecolor='none', edgecolor='red')
plot.show()