Python：带约束的双目标优化问题的实现方法

一、优化问题的数学表述方法

一个单目标优化问题可以表示为如下的数学公式：

$$\begin{array}{c}\begin{array}{rlr}\mathrm{minimize}& f_m(x)& m=1,...,M\\ \mathrm{subject}\mathrm{to}& g_j(x)\le 0& j=1,...,J\\ & h_k(x)=0& k=1,...,K\\ & x_i^L\le x_i\le x_i^U& i=1,...,N\\ & x\in \Omega & \end{array}\end{array}$$

基于公式(1)，本文实现了一个具有两个约束的双目标优化问题，其数学定义如下所示：

$$\begin{array}{c}\begin{array}{rl}\mathrm{minimize}& f_1(x)=100(x_1^2+x_2^2)\\ \mathrm{maximize}& f_2(x)=-(x_1-1{)}^2-x_2^2\\ \mathrm{subject}\mathrm{to}& g_1(x)=2(x_1-0.1)(x_1-0.9)\le 0\\ & g_2(x)=20(x_1-0.4)(x_1-0.6)\ge 0\\ & -2\le x_1\le 2\\ & -2\le x_2\le 2\\ & x\in \mathbb{R}\end{array}\end{array}$$

问题(2)包含两个目标($M=2$)，其中$f_1(x)$为最小化问题，$f_2(x)$为最大化问题；并包括两个不等式约束($J=2$)，其中$g_1(x)$为小于等于不等式，$g_2(x)$为大于等于不等式。该问题包含两个变量($N=2$)：$x_1$与$x_2$，它们的区间均为$[-2,2]$。

二、求解最优化问题并可视化结果

• Pareto最优解的位置确定
  对于公式(2)表述的带约束双目标优化问题，这里首先分析Pareto最优解的位置。$f_1$的最小值位于$(0,0)$，而$f_2$的最大值位于$(1,0)$。因为$f_1$与$f_2$两个函数都是二次函数，所以通过两个最优值之间直线给出双目标的最优值。在不考虑约束的情况下，这意味着所有Pareto最优解都有一个共同点，即$x_2=0$及$x_1\in (0,1)$。

• 约束分析
  第一个约束仅仅决定于$x_1$并且满足问题的区间为$x_1\in (0.1,0.9)$。而满足第二个约束的区间为$x_1\in [-2,0.4)\cup (0.6,2]$。因此，最优的Pareto集合为$PS=\{(x_1,x_2)|(0.1\le x_1\le 0.4)\cup (0.6\le x_1\le 0.9)\cap x_2=0\}$。

• Python代码实现

import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', family='serif')

# 定义两个自变量 x1 与 x2，分别在[-2, 2]区间内等间距抽取500个点
X1, X2 = np.meshgrid(np.linspace(-2, 2, 500), np.linspace(-2, 2, 500))

# 目标函数 f1 与 f2 的定义
F1 = 100 * (X1**2 + X2**2)
F2 = (X1 - 1)**2 + X2**2

# 约束定义
G1 = 2 * (X1[0] - 0.1) * (X1[0] - 0.9)
G2 = 20 * (X1[0] - 0.4) * (X1[0] - 0.6)

levels = np.array([0.02, 0.1, 0.25, 0.5, 0.8])
plt.figure(figsize=(7, 5))

# 绘制 f1 的等高图
CS = plt.contour(X1, X2, F1, 10*levels, color='black', alpha=0.5)
CS.collections[0].set_label("$f_1(x)$")

# 绘制 f2 的等高图
CS = plt.contour(X1, X2, F2, levels, linestyles="dashed", colors='black', alpha=0.5)
CS.collections[0].set_label("$f_2(x)$")

# 绘制约束 g1，以绿色线条显示，并将满足Pareto最优集合以实线显示，其余用虚线显示
plt.plot(X1[0], G1, linewidth=2.0, color="green", linestyle='dotted')
plt.plot(X1[0][G1<0], G1[G1<0], label="$g_1(x)$", linewidth=2.0, color="green")

# 绘制约束 g2，以蓝色线条显示，并将满足Pareto最优集合以实线显示，其余用虚线显示
plt.plot(X1[0], G2, linewidth=2.0, color="blue", linestyle='dotted')
plt.plot(X1[0][X1[0]>0.6], G2[X1[0]>0.6], label="$g_2(x)$", linewidth=2.0, color="blue")
plt.plot(X1[0][X1[0]<0.4], G2[X1[0]<0.4], linewidth=2.0, color="blue")

# 绘制求解的最优化问题的Pareto前沿，以橘黄色显示
plt.plot(np.linspace(0.1, 0.4, 100), np.zeros(100), linewidth=3.0, color="orange")
plt.plot(np.linspace(0.6, 0.9, 100), np.zeros(100), linewidth=3.0, color="orange")

# 设置x、y轴的显示范围及其坐标名称
plt.xlim(-0.5, 1.5)
plt.ylim(-0.5, 1)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")

# 显示图例
plt.legend(loc="upper center", bbox_to_anchor=(0.5, 1.12), ncol=4, fancybox=True, shadow=False)

plt.tight_layout()

代码执行结果如下图所示：

接下来，我们将Pareto前沿集合映射到目标空间，得到问题(2)的Pareto前沿。对于问题(2)，Pareto前沿等式是基于目标$f_2$，而$f_2$是以$f_1$为自变量的。及经过上面的分析，我们可以得到Pareto的最优解集$x_2=0$，这样目标函数可以简化为$f_1(x)=100x_1^2$与$f_2(x)=-(x_1-1{)}^2$。进而，将$x_1$用$f_1$表示，并将其带入$f_2$得到如下所示的等式：

$$\begin{array}{r}{f}_2=-(\sqrt{\frac{f_1}{100}}-1{)}^2\end{array}$$

经过前面的分析Pareto集合为$PS=\{(x_1,x_2)|(0.1\le x_1\le 0.4)\cup (0.6\le x_1\le 0.9)\cap x_2=0\}$，然后结合公式(3)，我们可以得到Pareto集合在目标空间对应的区间为$[1,16]$与$[36,81]$，下面给出计算Pareto前沿的代码如下所示：

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(7, 5))

# 目标空间的函数关系，以第一个目标 f1 为自变量，第二个目标 f2 为因变量
f2 = lambda f1: -((f1/100) ** 0.5 - 1)**2
# 在目标空间上绘制Pareto前沿
F1_a, F1_b = np.linspace(1, 16, 300), np.linspace(36, 81, 3000)
F2_a, F2_b = f2(F1_a), f2(F1_b)

# 绘制结果
plt.rc('font', family='serif')
plt.plot(F1_a, F2_a, linewidth=2.0, color="green", label="Pareto-front")
plt.plot(F1_b, F2_b, linewidth=2.0, color="green")

plt.xlabel("$f_1$")
plt.ylabel("$f_2$")

plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.10), ncol=4, fancybox=True, shadow=False)

plt.tight_layout()
plt.show()

代码执行结果如下图所示：

总结

上面例子相对来说比较简单，我们可以得到Pareto集(Pareto Set)与Pareto前沿(Pareto Front)的数学表达式。然而对于复杂的问题，从数学上是很难推导出Pareto集与Pareto前沿，甚至是不可能的。此时，我们就需要优化算法，通过良好的基准算法找到(接近最优的)解。