new post: multi-objectives linear programming

docs/2022-12-17-多目标线性规划的模糊折衷解法.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]
+```
+

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])