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

雾雨霜星

已于 2024-01-29 09:38:59 修改

阅读量9.8k

点赞数 19

分类专栏：数学建模文章标签：算法贪心算法

于 2022-07-30 10:02:38 首次发布

本文链接：https://blog.csdn.net/weixin_44786238/article/details/126068290

版权

数学建模专栏收录该内容

1 篇文章 1 订阅

订阅专栏

ε-约束算法

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)\leq0$
使用ε约束算法转化问题为：
$\min f_{1}(x)\\ f_{2}(x)\leq \epsilon_{2} ,...,f_{n}(x)\leq \epsilon_{n}\\ h(x)=0\\ g(x)\leq0$
得到多个待优化的子问题，每个子问题得到一个解加入帕累托解集，最终从帕累托解集中选取一个最优解。

其中的参数 $\epsilon_{2}$ ，…， $\epsilon_{n}$ 通过计算payoff table后得到。

payoff table的计算：

求解出第 $i$ 个目标函数的最优值 $f_{i}(x_{i}^{*})$ ，得到其最优解 $x_{i}^{*}$ 。
将 $x_{i}^{*}$ 代入其他目标函数得到 ${f_{1}(x_{i}^{*}),f_{2}(x_{i}^{*}),...,f_{n}(x_{i}^{*})\}$ 。
对全部目标函数按照上述流程求解，得到payoff table矩阵如下：
$\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}$

优化过程：

此时可从payoff table矩阵中对每个目标函数得到最优解(U)和最劣解(SN)： $f_{i}^{U}=f_{i}(x_{i}^{*}),f_{i}^{SN}=f_{i}(x_{j}^{*})$ 。
根据决策者的偏好，选取一个主目标函数 $f_{k}(x)$ 。
对于主目标函数外的目标函数 $f_{i}(x)$ ，设置一个网格化分数 $q_{ij}\in\{1,2,...,q_{i,max}\}$ 。与网格搜索法本质相同。
由此计算除了主目标函数外的其余目标函数的ε约束如下：
$\epsilon_{ij}=f_{i}^{SN}-(\frac{f_{i}^{SN}-f_{i}^{U}}{q_{ij}})\cdot j\quad j=1,2,...,q_{i,max}$
因此每个目标函数 $f_{i}(x)$ 会存在 $q_{i,max}$ 个约束，由此得到优化子问题共计 $\prod_{i=0,i\neq k}^{n}q_{i,max}$ 个。
得到每个优化子问题如下：
$\min f_{k}(x)\\ \text{subject to}\quad f_{1}(x)\leq\epsilon_{1j},f_{2}(x)\leq\epsilon_{1l},...,f_{n}(x)\leq\epsilon_{1m},h(x)=0,g(x)\leq0$
其中， $j=1,2,..,q_{1,max};l=1,2,..,q_{2,max};...;m=1,2,..,q_{n,max};$
每次求出一个最优解，若在可行域内则加入帕累托解集，若不在可行域内则丢弃。

特点：

约束方法的一个理想特征是，我们可以通过适当地将值分配给 $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