一、问题描述

求解有容量限制的的设施位置问题，使用Benders分解。模型如下：

\[min\quad\sum^{locations}_{j=1}fixedCost_j//open_j+\sum^{locations}_{j=1}\sum^{clients}_{i=1}cost_{ij}×supply_{ij} \]

\(s.t.\)

\[\sum^{locations}_{j=1}supply_{ij}=1\quad\quad\forall{i\in{clients}} \]

\[\sum^{clients}_{i=1}supply_{ij}\leq{capacity}_{ij}×open_{j}\quad\quad\forall{j\in{locations}} \]

\[0\leq{supply}_{ij}\leq{1}\quad\quad\forall{i\in{clients}};\forall{j\in{locations}} \]

\[open_i=0或1\quad\quad\forall{i\in{clients}} \]

二、程序

1. sys.exit函数：优雅地退出程序

摘自博文： Python中的sys.exit函数：优雅地退出程序_Python 笔记_设计学院

在Python编程中，程序在运行过程中可能会遇到需要停止程序的情况，如果不加处理，程序运行到中途就被强制停止的话，可能会导致数据丢失，甚至可能会让程序异常崩溃。因此，对于Python程序退出的处理，我们可以使用Python的内置函数sys.exit()，进行优雅地退出程序。

（1）sys.exit函数的基本用法

Python内置函数 \(sys.exit()\) ，用于 退出程序 。如果该函数被调用时 不带任何参数 ，那么Python解释器将会以 状态码0 来退出程序， 表示程序运行成功，并返回控制台 。下面是 \(sys.exit()\) 基本用法的示例代码：

import sys
sys.exit()
如果该函数被调用时带有整数参数n，那么Python解释器将会以状态码n来退出程序，并返回控制台。下面是以状态码1退出程序的示例代码：
import sys
sys.exit(1)
注意，当状态码不等于0时，表示程序运行发生某种异常或错误，需要进一步处理。
（2）sys.exit程序终止时的清理工作
在程序结束之前，我们可能需要完成一些清理工作，比如关闭一些文件、释放一些内存等。如果程序直接使用强制停止的方式结束，就有可能使这些工作被忽略或未能完全完成。此时，我们可以利用try/finally语句来完善这部分清理工作：
import sys
import time
    # code block here
    time.sleep(5)
finally:
    # closing file or releasing resource
    print('clean up resources')
    sys.exit()
2. 程序实现
（1）函数简介
\(cpx = cplex.Cplex()\)
\(cpx.variables.add(obj,lb,ub,type)\)
简介：cpx.cariables.add()是用于向模型中添加变量的方法，在新版的Cplex中，为了节约内存空间，通常以range的格式存储数据。所以为了方便查看，一般还要将其转换为list。
参数详解：
obj(list): 变量的目标函数系数列表；
lb(list): 变量的下界列表。默认为0；
ub(list): 变量的上界列表。默认为正无穷。
types(list): 变量的类型。可以是以下值：'C': 连续变量;'I': 整数变量;'B': 二进制变量。默认值为 'C'。
简介：cplex.SparsePair()用于表示稀疏的线性表达式。它主要用于定义线性约束和目标函数，特别是在涉及大量变量且多数系数为0的情况下。使用稀疏表示可以显著提高效率和节省内存。
参数详解：
ind(list): 一个包含变量索引的列表。这些索引指向模型中的特定变量（也就是用于输入决策变量）。这个列表为int型数据列表时，代表的是决策变量对应索引列表；为字符串型数据时，代表的是决策变量。如\(x_1+2x_3=0\)，ind=["x1","x3"]或ind=[0,2]
val(list): 与 ind 中的变量索引相对应的系数列表（与决策变量一一对应的系数）
参数详解：
lin_expr: 这是线性表达式的列表，表示线性约束的左侧。每个线性表达式由变量的索引和相应的系数组成，通常使用cplex.SparsePair来表示。
senses: 这是一个字符列表，表示每个约束的符号（'L'表示“<=”，'E'表示“=”，'G'表示“>=”）
rhs: 这是一个数字列表，表示每个约束的右侧值。
注意：上面所有参数只支持一维列表的输入
\[x+y=5
\]
\[2x-y\leq{10}
  cpx.linear_constraints.add(
      lin_expr=[[["x", "y"], [1, 1]], [["x", "y"], [2, -1]]],
      senses=["E", "L"],
      rhs=[5, 10]
\(cpx.long\_annotations.add(name, defval)\)
简介：cpx.long_annotations.add()一般用于添加长注解(long annotations)。长注解是Cplex用来为决策变量、约束等模型元素添加元数据或“注解”的机制。这些注解可能会影响求解器如何解决模型，尤其是在高级策略和方法中。
参数详解：
name (string): 注解的名称。Cplex预定义了一些注解名称，如cpx.long_annotations.benders_annotation，用于Benders分解策略。
defval (int): 注解的默认值。这是当注解没有明确为某个模型元素设置值时使用的值。
简介：使用cpx.solution.get_status_string()可获得关于这个状态的描述性字符串。这个函数特别有用，因为它可以让你快速了解模型解的状态，并据此采取相应的决策。
"optimal": 表示找到了最优解。
"infeasible": 表示问题是不可行的。
"unbounded": 表示问题是无界的。
"feasible": 表示找到了一个可行解，但不一定是最优的。
"integer optimal solution": 表示在整数线性规划或混合整数线性规划问题中，已经找到了最优的整数解。
... （还有其他可能的状态）
整数容差：整数容差定义了一个变量距离其最近的整数值可以有多远，而仍然被认为是整数。例如，如果整数容差设置为0.1，那么一个值为0.9或1.1的变量仍然会被认为满足整数约束。但如果整数容差设置得更小，例如0.001，那么这两个值就不会被认为满足整数约束。这个参数的目的是为了处理数值误差和确保求解器在数值上是稳健的。
自定义整数容差：如设置为0.001：cpx.parameters.mip.tolerances.integrality.set(0.001)
输出所有决策变量的值：values = cpx.solution.get_values()，如\(values=[1,2,3]\)，那么代表\(x_1=1;x_2=2;x_3=3\)
指定决策变量输出：如果只想输出\(x_1\)和\(x_3\)的值，那么根据其索引，可以写成：values = cpx.solution.get_values([0, 2])
# 用于求解模型的Benders分解类型
# 定义Benders分解的三种类型：无Benders分解、自动Benders分解和带注释的Benders分解
NO_BENDERS = 1
AUTO_BENDERS = 2
ANNO_BENDERS = 3
# 输出如何使用该脚本的说明，并退出程序。
def usage():
    print("""\
    Usage: facility.py [options] [inputfile]
     where
       inputfile describes a capacitated facility location instance as in
       ../../../../examples/data/facility.dat. If no input file
       is specified read the file in example/data directory.
       Options are:
       -a solve problem with Benders letting CPLEX do the decomposition
       -b solve problem with Benders specifying a decomposition
       -d solve problem without using decomposition (default)
     Exiting...
    sys.exit(2)
# 解决有容量限制的设施位置问题(模型输入部分)
def facility(bendersopt):
    """输入参数（已知量）"""
    fixedcost=[ 480, 200, 320, 340, 300]
    cost=[[24, 74, 31, 51, 84],
          [57, 54, 86, 61, 68],
          [57, 67, 29, 91, 71],
          [54, 54, 65, 82, 94],
          [98, 81, 16, 61, 27],
          [13, 92, 34, 94, 87],
          [54, 72, 41, 12, 78],
          [54, 64, 65, 89, 89]]
    capacity=[3, 1, 2, 4, 1]
    num_locations = len(fixedcost)  #计算区域数量locations
    num_clients = len(cost)  #计算顾客数量clients
    # 创建Cplex模型实例
    cpx = cplex.Cplex()
    """输入目标函数"""
    #下面为输入决策变量为Open部分的目标函数
    #obj：变量的目标函数系数列表；lb:变量的下界列表。默认为0；ub:变量的上界列表。默认为正无穷。
    #types:变量的类型。可以是以下值：'C': 连续变量;'I': 整数变量;'B': 二进制变量。默认值为 'C'。
    #lb是元素为0，1×num_location的列表；ub是元素为1，1×num_location的列表；这两个对应的是决策变量open的上限和下限约束
    open_ = list(cpx.variables.add(obj=fixedcost,
                                   lb=[0] * num_locations,
                                   ub=[1] * num_locations,
                                   types=["B"] * num_locations))
    # 在Cplex中，当你向模型中添加变量或约束时，它会返回一个表示新添加对象的索引范围的序列，这个序列用list来存储
    # 比如上面的open_,代表的是目标函数系数由 fixedcost 定义，且其值范围在 0 和 1 之间，最后会生成open_=[0,1,2,3,4](就是5个open变量的索引)
    #输入决策变量为supply部分的目标函数
    supply = [None] * num_clients  #初始化supply[i]的部分，用于代表顾客i的索引
    for i in range(num_clients):
        # 目标：最小化使用某一位置的固定成本，以及从特定位置为客户提供服务的成本。
        # 因为cost是一个二维list，所以需要通过for去一维维输入，这里的type默认为‘C’(连续变量)
        supply[i] = list(cpx.variables.add(obj=cost[i],
                                           lb=[0.0] * num_locations,
                                           ub=[1.0] * num_locations))
        # 经过循环之后，cpx.variable.add()就会在open_的基础上，继续往后面加变量，所以supply索引从5开始，supply=[[5, 6, 7, 8, 9], [10, 11, 12, 13, 14],...]
    """输入约束条件"""
    # 定义每个客户必须被分配的约束
    #   sum(j in nbLocations) supply[i][j] == 1  for each i in nbClients
    for i in range(num_clients):
        cpx.linear_constraints.add(
            lin_expr=[cplex.SparsePair(
                ind=supply[i], val=[1.0] * num_locations)],
            senses=["E"],
            rhs=[1.0])
    # 定义每个位置的容量必须被遵守的约束
    #   sum(i in nbClients) supply[i][j] <= capacity[j] * open_[j]
    for j in range(num_locations):
        #ind: 先把supply中的每一列取出来组成一个列表，再在列表后面append(open_[j])
        ind = [supply[i][j] for i in range(num_clients)] + [open_[j]]
        val = [1.0] * num_clients + [-capacity[j]]
        cpx.linear_constraints.add(
            lin_expr=[cplex.SparsePair(ind=ind, val=val)],
            senses=["L"],
            rhs=[0.0])
    """设置Benders分解（如果需要的话）"""
    # 这个代码没有使用Benders分解
    # 带注释的Benders分解
    if bendersopt == ANNO_BENDERS:
        # 我们通过指定结构来执行Benders分解，
        # 通过使用注解告诉CPLEX哪些变量在主问题中。
        # 默认情况下，变量被分配值CPX_BENDERS_MASTERVALUE+1，因此进入工作区。
        # 变量 open_[j] 应该进入主问题，所以
        # 我们为它们分配值 CPX_BENDERS_MASTER_VALUE。
        mastervalue = cpx.long_annotations.benders_mastervalue
        idx = cpx.long_annotations.add(
            name=cpx.long_annotations.benders_annotation,
            defval=mastervalue + 1)
        objtype = cpx.long_annotations.object_type.variable
        cpx.long_annotations.set_values(idx, objtype,
                                        [(open_[x], mastervalue)
                                         for x in range(num_locations)])
        print("Solving with explicit Benders decomposition.")
    # 自动Benders分解
    elif bendersopt == AUTO_BENDERS:
        # 让CPLEX自动分解问题。在有容量的设施位置问题中，
        # 主问题的变量应该是整数变量。通过将Benders策略参数设置为Full，
        # CPLEX会将所有整数变量放入主问题，将所有连续变量放入一个子问题，
        # 并进一步分解那个子问题（如果可能的话）。
        cpx.parameters.benders.strategy.set(
            cpx.parameters.benders.strategy.values.full)
        print("Solving with automatic Benders decomposition.")
    # 无Benders分解
    elif bendersopt == NO_BENDERS:
        print("Solving without Benders decomposition.")
    # 否则直接报错
    else:
        raise ValueError("invalid bendersopt argument")
    """解决模型并显示解决方案"""
    cpx.solve()  #求解优化问题
    print("Solution status =", cpx.solution.get_status_string())  #返回解的状态
    print("Optimal value:", cpx.solution.get_objective_value())  #返回求得的目标函数值
    tol = cpx.parameters.mip.tolerances.integrality.get()  #返回整数容差参数的当前值
    values = cpx.solution.get_values()  #输出决策变量的值
    # 遍历所有可能的设施位置
    for j in [x for x in range(num_locations) if values[open_[x]] >= 1.0 - tol]:
        # 首先，open[j]为0-1变量，要认定open[j]=1，需满足open[j]∈[1-tol,1+tol]，所以这里用了values[open_[x]] >= 1.0 - tol]
        # 这里先找出open[j]=1的量，然后进行循环遍历
        # 同样的，这里先找出supply[x][j]为1的数值，并生成针对该区域j的服务顾客编号列表
        # 占位符{0}表示的是设施地点编号j,占位符{1}表示的是对应j的服务顾客列表
        print("Facility {0} is open, it serves clients {1}".format(
            j, " ".join([str(x) for x in range(num_clients)
                         if values[supply[x][j]] >= 1.0 - tol])))
def main():
    """处理命令行参数。"""
    # filename = "../../../examples/data/facility.dat"   # 默认的数据文件路径
    # 解析命令行参数，设置Benders分解选项或数据文件路径
    # 初始化Benders分解选项为NO_BENDERS
    benders = NO_BENDERS
    # 遍历命令行参数
    for arg in sys.argv[1:]:
        # 判断参数是否为选项（开始于"-"）
        if arg.startswith("-"):
            # 根据选项设定Benders分解的方式
            if arg == "-a":
                benders = AUTO_BENDERS
            elif arg == "-b":
                benders = ANNO_BENDERS
            elif arg == "-d":
                benders = NO_BENDERS
            else:  # 如果是未知选项则调用usage函数显示帮助信息
                usage()
        else:  # 如果参数不是选项，则认为它是文件路径
            filename = arg
    # 调用facility函数来解决问题
    facility(bendersopt=benders)
if __name__ == "__main__":
    main()