[Multi-Objective Optimization]: Conventional e-Constraint Method / Augmented e-Constraint Method

Introduction

最近在Pyomo上继续学习怎么使用多目标优化，基于文章

Mavrotas G . Effective implementation of theε-constraint method in Multi-Objective Mathematical Programming problems[J]. Applied Mathematics & Computation, 2009, 213(2):455-465.

真的很佩服这种单枪匹马写出高水平论文的大佬。

Pareto optimal

The Pareto optimal (or efficient, non-dominated, non-inferior) solutions are the solutions that cannot be improved in one objective function without deteriorating their performance in at least one of the rest.

Conventional e-constraint method

example_1：

传统 e-constraint method 的思想：将f2(x)、f3(x)...fp(x)单独求最优后的值作为约束条件，这样求解的过程中，所有可行解都不会比该约束对应的解要低，满足了Pareto optimal 所要求的cannot be improved in one objective function 性质。

example_1变为：

传统 e-constraint method 的步骤（以example_1为例）：

1. 求解 max f1(x)为目标函数下，f2(x), f3(x)...fp(x)的值，这些值为对应function的可行区间端点值,记为f2_min, f3_min...fp_min。
2. 求解余下各个objective function的值，max f2(x), max f3(x)...max fp(x)，这些值为对应function的可行区间的端点值，记为f2_max, f3_max...fp_max。至此，得到了p-1个可行区间[(f2_min, f2_max),(f3_min, f3_max)...(fp_min, fp_max)]。
3. 在p-1个可行区间内设定步长，得到Pareto optimal解集的个数。若设置每个区间内步长统一为n，则一个区间内可以得到(f2_max-f2_min)/n个解，p-1个区间内就有(p-1)*(f2_max-f2_min)/n。
4. 进行p-1次循环，每次的步长值作为约束条件添加进优化模型中：

   for i in p-1 
       for j in range（fi_min, fi_max, yourstep): 
           model.e = i
           model.f2 <= model.e
           solve(model)
   

注意事项：

1. 在求解过程中，求出端点值fi_min和fi_max后需要有一个比较过程，因为根据不同的算例，二者求得的顺序有可能相反，在规定步长的代码中最大值和最小值传入顺序（先大后小）要正确否则会报错。
2. 在相应求解最小值的过程中，应改为

   model.f2 <= model.e

Augmented e-constraint method

增广e-constraint方法主要是为了避免传统e-constraint法可能会遇到的支配解也进入解集中的情况，并对寻优进行加速：

一个小李子：

# min [f1(x)=x2-x1, f2(x)=x1+x2]
# st  x1^2 - 2*x1 + 1 <= x2
#     0 <= x1 <= 1
#     0 <= x2 <= 1

一个垃圾实现：

POT = ConcreteModel(name='Pareto Optimization Test')
POT.x1 = Var(bounds=(0, 1), within=NonNegativeReals)
POT.x2 = Var(bounds=(0, 1), within=NonNegativeReals)

POT.c1 = Constraint(expr=POT.x1 * POT.x1 - POT.x1 * 2 + 1 <= POT.x2)

POT.f1 = Var()
POT.f2 = Var()
POT.e = Param(initialize=0, mutable=True)

POT.c2 = Constraint(expr=POT.f1 == POT.x2 - POT.x1)
POT.c3 = Constraint(expr=POT.f2 == POT.x2 + POT.x1)

POT.obj1 = Objective(expr=POT.f1, sense=minimize)
POT.obj2 = Objective(expr=POT.f2, sense=minimize)

POT.obj2.deactivate()
solver = SolverFactory('ipopt')
print(solver.solve(POT))
# POT.display()
f1_min = value(POT.f1)
f2_max = value(POT.f2)

POT.obj1.deactivate()
POT.obj2.activate()

solver = SolverFactory('ipopt')
print(solver.solve(POT))
# POT.display()
f1_max = value(POT.f1)
f2_min = value(POT.f2)

POT.eCons = Constraint(expr=POT.f2 <= POT.e)
POT.obj1.activate()
POT.obj2.deactivate()

n = 10
list1 = np.linspace(f2_min, f2_max, n)
# step = abs(f2_max-f2_min)/n
# steps = list(range(f2_max, f2_min, step)) + f2_max
# steps.sort()

a = []
b = []


POT.del_component(POT.obj1)
POT.del_component(POT.obj2)
POT.del_component(POT.eCons)

POT.eps = Param(initialize=0.001)
POT.s = Var(within=NonNegativeReals)

POT.obj1 = Objective(expr=POT.f1 - POT.s * POT.eps, sense=minimize)
POT.C_e = Constraint(expr=POT.f2 + POT.s == POT.e)
#
#
for i in list1:
    POT.e = i

    solver = SolverFactory('ipopt')
    solver.solve(POT)
    a.append(value(POT.x1))
    b.append(value(POT.x2))

plt.plot(a, b, 'o-.')
plt.title('efficient Pareto-front')
plt.grid(True)
plt.show()

注意事项：

需要注意求解最大值和最小值时obj fun对于slack variable的操作：  

初中数学都忘了，搞了我一天，两行泪。