diff --git "a/docs/2022-12-17-\345\244\232\347\233\256\346\240\207\347\272\277\346\200\247\350\247\204\345\210\222\347\232\204\346\250\241\347\263\212\346\212\230\350\241\267\350\247\243\346\263\225.md" "b/docs/2022-12-17-\345\244\232\347\233\256\346\240\207\347\272\277\346\200\247\350\247\204\345\210\222\347\232\204\346\250\241\347\263\212\346\212\230\350\241\267\350\247\243\346\263\225.md" new file mode 100644 index 0000000..3b00952 --- /dev/null +++ "b/docs/2022-12-17-\345\244\232\347\233\256\346\240\207\347\272\277\346\200\247\350\247\204\345\210\222\347\232\204\346\250\241\347\263\212\346\212\230\350\241\267\350\247\243\346\263\225.md" @@ -0,0 +1,192 @@ +--- +categories: [mathematics, optimization] +tags: [math, python] +mathjax: true +--- + +# 多目标线性规划的模糊折衷解法 + + +--- + +本文基于 Python `pulp` 建模框架实现求解、测试多目标线性规划的模糊折衷算法。 + + +## 多目标优化概述 + +多目标优化是多准则决策的一个领域,涉及同时优化多个目标函数的问题。各个目标之间通常相互制约,使得一个目标性能的改善往往是以损失其它目标性能为代价,即不可能存在一个使所有目标性能都达到最优的解。所以,对于多目标优化问题,其解通常是一个非劣解(非支配解、Pareto最优解、Pareto有效解)的集合——Pareto解集。研究人员从不同的角度研究多目标优化问题,从而在设置和解决多目标优化问题时存在不同的求解哲学和目标:可以是寻找Pareto解集,或者量化满足不同目标的折衷,或者找到满足人类决策者主观偏好的单一解决方案。 + +多目标优化算法可以归结为两类: + +- 传统优化算法,包括加权法、约束法和线性规划法等,本质是将多目标函数转化为单目标函数,通过采用单目标优化的方法达到对多目标函数的求解。 + +- 智能优化算法,包括进化算法(Evolutionary Algorithm)、粒子群算法(Particle Swarm Optimization)等。 + +> 本节参考:[多目标优化简述 ](https://imonce.github.io/2019/11/28/%E5%A4%9A%E7%9B%AE%E6%A0%87%E4%BC%98%E5%8C%96%E7%AE%80%E8%BF%B0/) + + +## 多目标线性规划的模糊折衷算法 + +多目标线性规划是多目标规划的子集,其中目标函数和约束都是线性函数形式。下面形式以最小化为例: + +$$ +\min \, \boldsymbol{Z} = \left[\boldsymbol{c}_1^T \boldsymbol{x}, ..., \boldsymbol{c}_n^T \boldsymbol{x} \right] \\ + +s.t. \quad \boldsymbol{A}\boldsymbol{x} \leq 0 \\ + +\quad\,\quad \boldsymbol{H}\boldsymbol{x} = 0 +$$ + +本文参考以下文章实现多目标线性规划的模糊折衷算法: + +> 李学全,李辉;多目标线性规划的模糊折衷算法[J];中南大学学报(自然科学版);2004年03期 + + +### 1. 算法 + +本质是转为一系列单目标线性规划问题。 +基本思想是构建目标函数的隶属度函数,然后 **弥补短板**(优化最差的分量),在此基础上 **提升整体**(优化所有分量之和)。 + +- 针对多目标规划中的各个目标,分别求取单目标情景下的最大值和最小值,依据最值构建各目标的隶属度函数; + + 例如对于目标函数分量 $z_i(\boldsymbol{x})$,分别求解两次单目标线性规划问题得到 $z_i^-=\min(z_i(\boldsymbol{x}))$ 和 $z_i^+=\max(z_i(\boldsymbol{x}))$,则相应隶属度函数: + + $$ + u_i(\boldsymbol{x}) = \frac{z_i(\boldsymbol{x})-z_i^-(\boldsymbol{x})}{z_i^+(\boldsymbol{x})-z_i^-(\boldsymbol{x})} + $$ + + +- 第一阶段,以最小化所有隶属度函数的最大值分量为规划目标,开展单目标线性规划; + + $$ + \min \, \lambda \\ + + s.t. \quad \lambda \geq u_i(\boldsymbol{x}) \\ + + \quad\,\quad \boldsymbol{A}\boldsymbol{x} \leq 0, \boldsymbol{H}\boldsymbol{x} = 0 + $$ + + 解得 $\boldsymbol{x}^*$,对应隶属度函数 $u^*(\boldsymbol{x})$。 + +- 第二阶段,以第一阶段求得的隶属度 $u^*(\boldsymbol{x})$ 为最大值约束,最小化所有隶属度函数之和,获得多目标折衷规划结果。 + + $$ + \min \, \lambda = \frac{\Sigma_i {\lambda_i}}{n} \\ + + s.t. \quad u_i(\boldsymbol{x}) \leq \lambda_i \leq u^*(\boldsymbol{x}) \\ + + \quad\,\quad \boldsymbol{A}\boldsymbol{x} \leq 0, \boldsymbol{H}\boldsymbol{x} = 0 + $$ + +综上,对于 $n$ 个目标的线性规划问题,需要进行 $2n+2$ 次单目标线性规划问题的求解。 + + +### 2. 实现 + +`pulp` 是一个基于 Python 的线性规划问题建模框架,自带 `cbc` 求解器,同时支持调用其他商用如 Gurobi、COPT 或开源求解器如 SCIP。接下来基于 `pulp` 实现以上算法。 + +首先设计一个多目标线性规划的求解类 `MOModel`,接收原问题的目标函数列表和约束列表。 + +```python +class MOModel: + def __init__(self, objectives:List[pulp.LpAffineExpression], + constraints:List[pulp.LpConstraint]=None) -> None: + '''Multiple objectives model based on `pulp`, including variables, objectives and constraints. + + Args: + objectives (list): objective functions list. NOTE: minimize each component. + constraints (list): constraints list. + kwargs (dict): solving parameters, e.g., solver name, rel_gap. + ''' + self.objectives = objectives or [] + self.constrains = constraints or [] +``` + +普通线性规划是求解原问题的基石,`pulp` 正好胜任: + +```python +def _solve_single_objective_problem(self, name:str, obj:pulp.LpAffineExpression, + cons:List[pulp.LpConstraint]) -> pulp.LpProblem: + # problem + prob = pulp.LpProblem(name=name, sense=pulp.LpMinimize) + + # objective and constraints + prob += obj + for c in cons: prob += c + + # solve + status = prob.solve() + assert status==1, f'Problem {name}: {pulp.LpStatus[status]} solution found.' + return prob +``` + +接下来只需按照前一小节的步骤构造目标函数和约束,依次求解即可。 + +```python +def solve(self): + '''Solve multiple objectives linear programming problem with two stages fuzzy compromise approach.''' + # stage 0: membership functions + logging.info('(FC-1) Solving the lower and upper bounds of each objective...') + min_range, max_range = self._solve_lower_and_upper_ranges() + memberships = self.membership_function(self.objectives, min_range, max_range) + + # stage 1: minimize the maximum component + logging.info('(FC-2) Minimizing the maximum component of objectives...') + self.solve_stage_1(memberships=memberships) + + # stage 2: minimize sum of all components by taking stage 1 results as a baseline + logging.info('(FC-3) Minimizing the sum of all objectives...') + m0 = self.membership_function(self.objective_values, min_range, max_range) + self.solve_stage_2(m0=m0, memberships=memberships) +``` + +完整代码参考: + +> https://github.com/dothinking/blog/blob/master/samples/fc.py + + + +## 算例 + +直接使用原论文中的例子: + +$$ +\max Z_1 = 2x_1 + 5x_2 + 7x_3 + x_4 \\ +\max Z_2 = 4x_1 + x_2 + 3x_3 + 11x_4 \\ +\max Z_3 = 9x_1 + 3x_2 + x_3 + 2x_4 \\ +\min W_1 = 1.5x_1 + 2x_2 + 0.3x_3 + 3x_4 \\ +\min W_2 = 0.5x_1 + x_2 + 0.7x_3 + 2x_4 \\ + +s.t. \quad 3x_1 + 4.5x_2 + 1.5x_3 + 7.5x_4 = 150 \\ +x_1 \geq 0, x_2 \geq 0, x_3 \geq 0, x_4 \geq 0 +$$ + + +注意将最大化目标函数取反统一为最小化,然后按照 `pulp` 语法建模即可。最终结果与论文一致。 + + +```python +# variables +x = [pulp.LpVariable(name=f'x_{i+1}', lowBound=0) for i in range(4)] + +# objectives +A = [ + [-2, -5, -7, -1], + [-4, -1, -3, -11], + [-9, -3, -1, -2], + [1.5, 2, 0.3, 3], + [0.5, 1, 0.7, 2] +] +objs = [pulp.lpDot(x, a) for a in A] + +# constraints +cons = [pulp.lpDot(x, [3, 4.5, 1.5, 7.5])==150] + +# solve +m = MOModel(objectives=objs, constraints=cons) +m.solve() + +# results +print([t.value() for t in x]) # [25.0, 0.0, 50.0, 0.0] +``` + diff --git a/samples/fc.py b/samples/fc.py new file mode 100644 index 0000000..d6ce012 --- /dev/null +++ b/samples/fc.py @@ -0,0 +1,136 @@ +''' +Two-stages fuzzy compromise approach for multi-objectives linear programming problems. + +李学全,李辉;多目标线性规划的模糊折衷算法[J];中南大学学报(自然科学版);2004年03期 +''' + +import logging +from typing import (List, Union) +import pulp + + +class MOModel: + def __init__(self, objectives:List[pulp.LpAffineExpression], + constraints:List[pulp.LpConstraint]=None) -> None: + '''Multiple objectives model based on `pulp`, including variables, objectives and constraints. + + Args: + objectives (list): objective functions list. NOTE: minimize each component. + constraints (list): constraints list. + kwargs (dict): solving parameters, e.g., solver name, rel_gap. + ''' + self.objectives = objectives or [] + self.constrains = constraints or [] + + @property + def objective_values(self): + return [obj.value() for obj in self.objectives] + + def membership_function(self, objs:List[Union[float, pulp.LpAffineExpression]], + min_range:list, max_range:list): + '''The membership function is interpreted as the fuzzy value of the decision maker, + which describes the behavior of indifference, preference or aversion toward uncertainty. + In solving fuzzy mathematical programming problems, a linear membership function defined + by two points, the upper and lower levels of acceptability is used because of its simplicity. + ''' + return [1/(u-l) * (x-l) for l,u,x in zip(min_range, max_range, objs)] + + def solve(self): + '''Solve multiple objectives linear programming problem with two stages fuzzy compromise approach.''' + # stage 0: membership functions + logging.info('(FC-1) Solving the lower and upper bounds of each objective...') + min_range, max_range = self._solve_lower_and_upper_ranges() + memberships = self.membership_function(self.objectives, min_range, max_range) + + # stage 1: minimize the maximum component + logging.info('(FC-2) Minimizing the maximum component of objectives...') + self.solve_stage_1(memberships=memberships) + + # stage 2: minimize sum of all components by taking stage 1 results as a baseline + logging.info('(FC-3) Minimizing the sum of all objectives...') + m0 = self.membership_function(self.objective_values, min_range, max_range) + self.solve_stage_2(m0=m0, memberships=memberships) + + def solve_stage_1(self, memberships:List[pulp.LpAffineExpression]): + '''Stage 1: minimize the maximum component of membership functions.''' + # new temp variable + x = pulp.LpVariable(name='lambda', lowBound=0, upBound=1) + + # the maximum component + lower_cons = [x>=m for m in memberships] + cons = self.constrains + lower_cons + + # solve problem + return self._solve_single_objective_problem(name='stage_1', obj=x, cons=cons) + + def solve_stage_2(self, m0:List[float], memberships:List[pulp.LpAffineExpression]): + '''Stage 2: minimize the sum of all membership function components.''' + # new temp variables + X = [pulp.LpVariable(name=f'lambda_{i}', lowBound=0) for i in range(len(self.objectives))] + + # optimize the bounds further + upper_cons = [x<=m for x,m in zip(X,m0)] # upper bound of each component + lower_cons = [x>=m for x,m in zip(X,memberships)] # lower bound of each component + cons = self.constrains + lower_cons + upper_cons + + # solve problem + return self._solve_single_objective_problem(name='stage_2', obj=pulp.lpSum(X), cons=cons) + + def _solve_lower_and_upper_ranges(self): + '''Calculate lower and upper ranges by solving single objective optimization for each objective component.''' + min_range, max_range = [], [] + for i,obj in enumerate(self.objectives): + # check or calculate lower range + res = obj.lowBound if hasattr(obj, 'lowBound') else None + if res is None: + p = self._solve_single_objective_problem(name=f'solve_lower_range_{i}', obj=obj, cons=self.constrains) + res = p.objective.value() + min_range.append(res) + + # check or calculate upper range + res = obj.upBound if hasattr(obj, 'upBound') else None + if res is None: + p = self._solve_single_objective_problem(name=f'solve_upper_range_{i}', obj=-obj, cons=self.constrains) + res = -p.objective.value() + max_range.append(res) + return min_range, max_range + + def _solve_single_objective_problem(self, name:str, obj:pulp.LpAffineExpression, + cons:List[pulp.LpConstraint]) -> pulp.LpProblem: + # problem + prob = pulp.LpProblem(name=name, sense=pulp.LpMinimize) + + # objective and constraints + prob += obj + for c in cons: prob += c + + # solve + status = prob.solve() + assert status==1, f'Problem {name}: {pulp.LpStatus[status]} solution found.' + return prob + + +if __name__=='__main__': + + # variables + x = [pulp.LpVariable(name=f'x_{i+1}', lowBound=0) for i in range(4)] + + # objectives + A = [ + [-2, -5, -7, -1], + [-4, -1, -3, -11], + [-9, -3, -1, -2], + [1.5, 2, 0.3, 3], + [0.5, 1, 0.7, 2] + ] + objs = [pulp.lpDot(x, a) for a in A] + + # constraints + cons = [pulp.lpDot(x, [3, 4.5, 1.5, 7.5])==150] + + # solve + m = MOModel(objectives=objs, constraints=cons) + m.solve() + + # results + print([t.value() for t in x])