【多目标规划问题求解】ε-约束算法_约束法多目标规划问题求解

ε-约束算法

TIME:2022/07/30
Author:雾雨霜星
Web:雾雨霜星的小站
小站原文链接：https://www.shuangxing.top/#/post?id=27

参考

[1]Nima Amjady and Jamshid Aghaei and Heidar Ali Shayanfar. Stochastic Multiobjective Market Clearing of Joint Energy and Reserves Auctions Ensuring Power System Security[J]. IEEE Transactions on Power Systems, 2009, 24(4) : 1841-1854.

[2]J. Aghaei and N. Amjady and H.A. Shayanfar. Multi-objective electricity market clearing considering dynamic security by lexicographic optimization and augmented epsilon constraint method[J]. Applied Soft Computing Journal, 2011, 11(4) : 3846-3858.

简介

用于计算多目标数学问题的一种计算方法。

选取一个主目标函数，将其余目标函数转化为约束，从而计算每个子优化目标，得到帕累托解集。

方法

原多目标优化问题如下：
min ⁡ { f 1 ( x ) , f 2 ( x ) , . . . , f n ( x ) } h ( x ) = 0 g ( x ) ≤ 0 \min \{f_{1}(x),f_{2}(x),...,f_{n}(x)\}\\ h(x)=0\\ g(x)\leq0 min{f1​(x),f2​(x),...,fn​(x)}h(x)=0g(x)≤0
使用ε约束算法转化问题为：
$$\begin{gathered} \min f_1(x) \\ {f}_2(x)\le {ϵ}_2,...,f_n(x)\le {ϵ}_n \\ h(x)=0 \\ g(x)\le 0 \end{gathered}$$
得到多个待优化的子问题，每个子问题得到一个解加入帕累托解集，最终从帕累托解集中选取一个最优解。

其中的参数${ϵ}_2$，…，${ϵ}_n$通过计算payoff table后得到。

payoff table的计算：

1. 求解出第$i$个目标函数的最优值$f_i(x_i^{\ast })$，得到其最优解$x_i^{\ast }$。

2. 将$x_i^{\ast }$代入其他目标函数得到 { f 1 ( x i ∗ ) , f 2 ( x i ∗ ) , . . . , f n ( x i ∗ ) } \{f_{1}(x_{i}^{*}),f_{2}(x_{i}^{*}),...,f_{n}(x_{i}^{*})\} {f1​(xi∗​),f2​(xi∗​),...,fn​(xi∗​)}。

3. 对全部目标函数按照上述流程求解，得到payoff table矩阵如下：
   ( f 1 ( x 1 ∗ ) . . . f i ( x 1 ∗ ) . . . f n ( x 1 ∗ ) ⋮ ⋱ ⋮ f 1 ( x i ∗ ) . . . f i ( x i ∗ ) . . . f n ( x i ∗ ) ⋮ ⋱ ⋮ f 1 ( x n ∗ ) . . . f i ( x n ∗ ) . . . f n ( x n ∗ ) )

   $$\begin{pmatrix} f_{1}(x_{1}^{*})&...&f_{i}(x_{1}^{*})...&f_{n}(x_{1}^{*})\\ \vdots&\ddots&&\vdots\\ f_{1}(x_{i}^{*})&...&f_{i}(x_{i}^{*})...&f_{n}(x_{i}^{*})\\ \vdots&\ddots&&\vdots\\ f_{1}(x_{n}^{*})&...&f_{i}(x_{n}^{*})...&f_{n}(x_{n}^{*}) \end{pmatrix}$$ ​f1​(x1∗​)⋮f1​(xi∗​)⋮f1​(xn∗​)​...⋱...⋱...​fi​(x1∗​)...fi​(xi∗​)...fi​(xn∗​)...​fn​(x1∗​)⋮fn​(xi∗​)⋮fn​(xn∗​)​ ​

优化过程：

1. 此时可从payoff table矩阵中对每个目标函数得到最优解(U)和最劣解(SN)：$f_i^U=f_i(x_i^{\ast }),f_i^{SN}=f_i(x_j^{\ast })$。

2. 根据决策者的偏好，选取一个主目标函数$f_k(x)$。

3. 对于主目标函数外的目标函数$f_i(x)$，设置一个网格化分数 q i j ∈ { 1 , 2 , . . . , q i , m a x } q_{ij}\in\{1,2,...,q_{i,max}\} qij​∈{1,2,...,qi,max​}。与网格搜索法本质相同。

4. 由此计算除了主目标函数外的其余目标函数的ε约束如下：
   $${ϵ}_{ij}=f_i^{SN}-(\frac{f_i^{SN}-f_i^U}{q_{ij}})\cdot jj=1,2,...,q_{i,max}$$
   因此每个目标函数$f_i(x)$会存在$q_{i,max}$个约束，由此得到优化子问题共计${\prod }_{i=0,i\ne k}^n{q}_{i,max}$个。

5. 得到每个优化子问题如下：
   $$\begin{gathered} \min f_k(x) \\ \text{subject to}{f}_1(x)\le {ϵ}_{1j},f_2(x)\le {ϵ}_{1l},...,f_n(x)\le {ϵ}_{1m},h(x)=0,g(x)\le 0 \end{gathered}$$
   其中，$j=1,2,..,q_{1,max};l=1,2,..,q_{2,max};...;m=1,2,..,q_{n,max};$

6. 每次求出一个最优解，若在可行域内则加入帕累托解集，若不在可行域内则丢弃。

特点：

约束方法的一个理想特征是，我们可以通过适当地将值分配给$q_{ij}$来控制有效集表示的密度。网格点的数量越高，有效集的表示越密集，但计算时间越长。在有效集的密度和计算时间之间进行权衡总是可取的。

============================================================================
2023-08-29更新：
由于该文较为热门，很多人私信询问我的代码，因此我将自己的代码加入了文章中，方便各位学习，如有引用和转载此文章，请在文中进行说明标注此文出处。

# 说明：
# 1. 进行多目标优化，目标函数为：综合投资成本(ainv)、投资回收期(pb)、系统不确定适应度(a)。
# 2. 此处ε-约束法所得的每个子优化问采用Gurobi进行求解。
# 3. 默认系统不确定适应度a目标函数最劣为0，最优为1。
# 4. 由于此处我的Gurobi进行计算需要提供系统不确定适应度，因此默认以a=0的最劣来求其余目标函数的最优情况。
# 5. 选取不确定适应度a作为ε-约束算法的主目标函数。
# ----------------- 主程序 ----------------
# 记录开始时间
time_start = time.time()

# 非主目标函数的目标函数网格化参数（即上文所说的网格化分数qij）
ainv_Q = 10
pb_Q = 10

# 所求解的集合
solution_pool = []

# 初始化 ξ约束法的payoff矩阵
ainv_targetValue = []
pb_targetValue = []
a_targetValue = []

# ε-约束算法首先需要得到各个目标函数的最优和最劣，最劣从其他目标函数最优时将变量代入计算该目标得到
# ainv最优
td = ainv_gbm(a=0)  # ainv_gbm是计算ainv单目标优化的Gurobi函数
ainv_targetValue.append(td['annual_inv'])  # td['annual_inv']是ainv_gbm函数Gurobi优化结果中ainv的值
pb_targetValue.append(td['pb'])  # td['pb']是ainv_gbm函数Gurobi优化结果中pb的值
a_targetValue.append(0)

# pb最优
td = pb_gbm(a=0) # pb_gbm是计算pb单目标优化的Gurobi函数
ainv_targetValue.append(td['annual_inv'])
pb_targetValue.append(td['pb'])
a_targetValue.append(0)

# a最优（此处是为了得到a最优下的ainv和pb的值用于构建payoff矩阵）
td, o_flag = igdt_gbm()   # igdt_gbm是计算a单目标优化的Gurobi函数
ainv_targetValue.append(td['annual_inv'])
pb_targetValue.append(td['pb'])
a_targetValue.append(td['a'])  # td['a']是igdt_gbm函数Gurobi优化结果中a的值


# 获取最优值与最劣值
annualInv_U = min(ainv_targetValue)
annualInv_SN = max(ainv_targetValue)
paybackPeriod_U = min(pb_targetValue)
paybackPeriod_SN = max(pb_targetValue)
print("=======================================")
print("annualInv_U:", annualInv_U)
print("annualInv_SN:", annualInv_SN)
print("paybackPeriod_U:", paybackPeriod_U)
print("paybackPeriod_SN:", paybackPeriod_SN)
print("=======================================")

# ξ约束模型求解开始
for q1 in range(ainv_Q):
    for q2 in range(pb_Q):
        # 指示程序进展
        print('the ξ gens:', q1, q2)
        # 获取本轮ξ约束
        annualInv_target_c = annualInv_SN - (annualInv_SN - annualInv_U) / ainv_Q * q1
        paybackPeriod_target_c = paybackPeriod_SN - (paybackPeriod_SN - paybackPeriod_U) / pb_Q * q2
        # 子优化问题的求解
        # 通过传入annualInv_target_c与paybackPeriod_target_c在Gurobi内构建约束
        td_e, flag_e = igdt_gbm(ainv_c=annualInv_target_c, pb_c=paybackPeriod_target_c)
        # flag_e 是标志Gurobi求解是否有解的变量
        if not flag_e: continue
        # 更新所求解集合
        else: solution_pool.append([td_e['annual_inv'], td_e['pb'], 1 - td_e['a']])

# 记录结束时间
time_end = time.time()
time_sum = time_end - time_start

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

后续对所得解集合进行非支配排序即可得到帕累托解集。
具体帕累托解集中选取最优解，可参考模糊决策方法：
https://www.shuangxing.top/#/post?id=29