线性规划

线性规划（Linear Programming, LP）是指目标函数与约束条件均是关于变量的线性函数的优化问题。线性规划是最简单的优化问题类型，本节只讨论连续变量的线性规划。

关于线性规划，我们有很多工具可以使用。本节将会使用：

[2]:

import scipy, cvxpy, pyomo

for lib in (scipy, cvxpy, pyomo):
    print(lib.__name__, lib.__version__, sep='\t')

scipy   1.10.1
cvxpy   1.4.2
pyomo   6.7.1

例题

考虑这样一个问题：

工厂现有的某机器能生产甲、乙、丙三种产品，它们的具体信息如下所示：

每种产品的单件利润
每种产品在生产单件时所需的人工工时
每种产品在生产单件时所需的机器工时
每种产品的数量

产品	利润	人工工时	机器工时	产品数量
A	10	2	4
B	7	0.5	3
C	3	1	1
约束		= 45	<= 100	<= 60

假设该机器平均每周固定分配人工 45 小时、机器工作不超过 100 小时，且不应生产超过 60 件产品；生产的产品全部可以卖出。求如何分配三种产品的生产数量（可以是小数），使得平均每周的利润最大。

我们列出优化模型如下，它是一个线性规划问题：

\begin{align*} \min\ & f(x_1, x_2, x_3) = -10x_1 - 7x_2 - 3x_3 \\ \text{s.t.}\quad & 4x_1 + 3x_2 + x_3 \leq 100 \\ & x_1 + x_2 + x_3 \leq 60 \\ & 2x_1 + 0.5x_2 + x_3 = 45 \\ & x_1, x_2, x_3 \geq 0 \end{align*}

scipy.optimize.linprog

Scipy 可能是许多人第一次在 Python 上解线性规划问题时使用的工具。scipy.optimize.linprog 需要我们使用矩阵形式书写代码：

\begin{align*} \min\ & f(x) = \boldsymbol{C}x \\ \text{s.t.}\quad & \boldsymbol{A}_{leq} x \leq b_{leq} \\ & \boldsymbol{A}_{eq} x = b_{eq} \\ & \ell \leq x \leq u \end{align*}

[3]:

import numpy as np
from scipy.optimize import linprog

c = np.array([-10, -7, -3])
A_leq = np.array([[4, 3, 1],
                  [1, 1, 1]])
b_leq = np.array([100, 60])
A_eq = np.array([[2, 0.5, 1]])
b_eq = np.array([45])
x_bounds = [(0, None) for _ in range(3)]
result = linprog(c, A_leq, b_leq, A_eq, b_eq, bounds=x_bounds)
print(result.message)

Optimization terminated successfully. (HiGHS Status 7: Optimal)

以上 Optimal 状态信息表示 LP 问题成功解得最优值。我们利用 result.x 查看取得最优值时，变量的取值：

[4]:

x = np.array(result.x)
x

[4]:

array([ 0., 22., 34.])

查看目标函数最优值 result.fun ，或者使用 c @ x 手动计算：

[5]:

print("Optimal:\t", result.fun)
print("Manual compute:\t", c @ x)

Optimal:         -256.0
Manual compute:  -256.0

还可以用 result.slack 查看不等式约束的满足情况，或用 result.con 查看等式约束的满足情况：

[6]:

print(f"Leq constraints: {b_leq - (A_leq @ x).flatten()}",
      f"\t{result.slack = }",
      f"Eq constraints: {b_eq - (A_eq @ x).flatten()}",
      f"\t{result.con = }",
      sep='\n')

Leq constraints: [0. 4.]
        result.slack = array([0., 4.])
Eq constraints: [0.]
        result.con = array([0.])

cvxpy

cvxpy 的用法也比较类似，不过需要用 cvxpy.Variable() 方法声明变量。

cvxpy 可以使用 <= 或者 >= 来表示带等的不等于号；但是，不能像 1<= x <=2 这样连写。

[7]:

import numpy as np
import cvxpy as cp

c = np.array([-10, -7, -3])
A_leq = np.array([[4, 3, 1],
                  [1, 1, 1]])
b_leq = np.array([100, 60])
A_eq = np.array([[2, 0.5, 1]])
b_eq = np.array([45])

x = cp.Variable(3)
x_lower = np.zeros(3)
problem = cp.Problem(cp.Minimize(c @ x),
                     [A_leq @ x <= b_leq, A_eq @ x <= b_eq, x >= x_lower])
problem.solve()

[7]:

-255.99999919804102

LP 问题的解结果（例如，求得最优 cvxpy.OPTIMAL 、不可行 cvxpy.INFEASIBLE 或者未约束 cvxpy.UNBOUNDED 等）：

[8]:

print(f"Status: {problem.status}")
print(problem.status == cp.OPTIMAL)

Status: optimal
True

LP 问题的解存放在 x.value 与 problem.value 中：

[9]:

print("A feasible solution:\t", x.value)
print("Optimal value:\t\t", problem.value)

A feasible solution:     [9.66046589e-07 2.19999993e+01 3.39999981e+01]
Optimal value:           -255.99999919804102

对偶解（Dual Solution）中符号为正的第 i 个值表示在最优解时约束 \(a_i x \leq b_i\) 恰好取得等号。这也意味着改变 \(b_i\) 会引起最优解的改变。

[10]:

constraint_labels = ["A_leq", "A_eq", "x_bounds"]
print("Dual solution:")
for cgroup, label in zip(problem.constraints, constraint_labels):
     print(f"    {label:<10}", cgroup.dual_value)

Dual solution:
    A_leq      [2.19999999e+00 1.01474852e-07]
    A_eq       [0.79999995]
    x_bounds   [3.99999883e-01 9.66102358e-09 6.00201485e-09]

我们传入了三组约束。在上述对偶解中，小于号约束组 A_leq 的第二值为 0，表示原问题的 \(x_1 + x_2 + x_3 \leq 60\) 是松弛的（在最优解时未取得等号）。同理，变量约束组 x_bounds 的第二、三值也是松弛的，表示 \(x_2, x_3\) 的最优解值不在边界上（而显然 \(x_1^* = 0\) 在边界上）。

pyomo (AbstractModel)

通过 pyomo 提供的接口，我们调用需要的求解器来解优化问题。我们将前述的线性规划问题

\begin{align*} \min\ & f(x_1, x_2, x_3) = -10x_1 - 7x_2 - 3x_3 \\ \text{s.t.}\quad & 4x_1 + 3x_2 + x_3 \leq 100 \\ & x_1 + x_2 + x_3 \leq 60 \\ & 2x_1 + 0.5x_2 + x_3 = 45 \\ & x_1, x_2, x_3 \geq 0 \end{align*}

抽象成规范的优化问题表述的形式，以便使用 AbstractModel 来建立问题（如果愿意，也可以使用 ConcreteModel ）。下面是原问题的表述（目标函数写为最大化形式）：

集合
- \(J\)：产品下标集
参数
- \(c_j\)：产品 \(j\) 的利润，\(j\in J\)
- \(a^{(machine)}_j\)：产品 \(j\) 的机器工时（正数），\(j\in J\)
- \(a^{(product)}_j\)：产品 \(j\) 的生产数量（非负数），\(j\in J\)
- \(a^{(human)}_j\)：产品 \(j\) 的人工工时（正数），\(j\in J\)
- \(ub^{(machine)}\): 总机器工时上限（正数），标量
- \(ub^{(product)}\)：总产品数量上限（正数），标量
- \(eb^{(human)}\)：固定的总人工工时（正数），标量
变量
- \(x_j\)：产品 \(j\) 的生产数量，\(j\in J\)
目标函数
- 最大化产品的总利润： \(\max \sum_{j=1}^n c_j x_j,\quad j\in J\)
约束条件：以下均有 \(\forall j\in J\)
- 机器工时上限：\(\sum_{j=1}^n a^{(machine)}_{j} x_j \leq ub^{(machine)}\)
- 产品生产数量上限：\(\sum_{j=1}^n a^{(product)}_{j} x_j \leq ub^{(product)}\)
- 固定人工工时：\(\sum_{j=1}^n a^{(human)}_{j} x_j = eb^{(human)}\)
- 产品数量非负：\(x_j \geq 0\)

本例用 pyomo 调用 CPLEX（需要预先安装）来解上述线性规划问题。

约束条件中的变量非负由 within=pyo.NonNegativeReals 给出，不用再显示地写出。

[11]:

import pyomo.environ as pyo
infinity = float('inf')

model = pyo.AbstractModel()

# --- 集合 ---
model.J = pyo.Set()

# --- 参数 ---
model.c = pyo.Param(model.J, within=pyo.NonNegativeReals)

model.a_machine = pyo.Param(model.J, within=pyo.PositiveReals)
model.a_product = pyo.Param(model.J, within=pyo.NonNegativeReals)
model.a_human = pyo.Param(model.J, within=pyo.PositiveReals)

model.ub_machine = pyo.Param(within=pyo.PositiveReals)
model.ub_product = pyo.Param(within=pyo.PositiveReals)
model.eb_human = pyo.Param(within=pyo.PositiveReals)

# --- 变量 ---
model.x = pyo.Var(model.J, within=pyo.NonNegativeReals)

# --- 目标函数 ---
def sale_rule(model):
    return sum(model.c[j] * model.x[j] for j in model.J)

model.objectives = pyo.Objective(rule=sale_rule, sense=pyo.maximize)  # 最大化

# --- 约束条件
def machine_constraint_rule(model):
    machine_hours = sum(model.a_machine[j] * model.x[j] for j in model.J)
    return machine_hours <= model.ub_machine

model.constraint_machine = pyo.Constraint(rule=machine_constraint_rule)

def product_constraint_rule(model):
    total_product = sum(model.a_product[j] * model.x[j] for j in model.J)
      # 另一种写法: (None, total_product, model.ub_product)
    return total_product <= model.ub_product

model.constraint_product = pyo.Constraint(rule=product_constraint_rule)

def human_constraint_rule(model):
    human_hours = sum(model.a_human[j] * model.x[j] for j in model.J)
    return human_hours == model.eb_human

model.constraint_human = pyo.Constraint(rule=human_constraint_rule)

导入模型参数值

我们可以直接使用 Python 字典变量来传入数据，也可以使用外部文件。先看使用字典，如果参数不含键名，它使用 None 字段给出参数的值：

[12]:

datadict = {
    None: {
        "J": {None: ['A', 'B', 'C']},
        "c": {'A': 10, 'B': 7, 'C': 3},
        "a_machine": {'A': 4, 'B': 3, 'C': 1},
        "a_product": {'A': 1, 'B': 1, 'C': 1},
        "a_human": {'A': 2, 'B': 0.5, 'C': 1},
        "ub_machine": {None: 100},
        "ub_product": {None: 60},
        "eb_human": {None: 45}
    }
}

如果使用外部文件，我们可以准备一个 YAML 文件（需要 pyyaml 库）或者 JSON 文件来指定参数的值。使用 YAML 文件也是我喜欢的做法（因为相比 JSON 不需要输入很多引号），例如：

[13]:

data_file = "dataset/LP_machine.yaml"
with open(data_file, "r") as f:
    data_rawstr = f.read()
print(data_rawstr)

J: [A, B, C]
c:
  A: 10
  B: 7
  C: 3
a_machine:
  A: 4
  B: 3
  C: 1
a_product:
  A: 1
  B: 1
  C: 1
a_human:
  A: 2
  B: 0.5
  C: 1
ub_machine: 100.0
ub_product: 60.0
eb_human: 45

又或者等价地，我们也可以准备一个 pyomo Data Command Files 格式（类似 AMPL 格式，但不完全相同）的 .dat 文件（具体参考 pyomo 文档 - Working with Abstract Models: Data Command Files）。这样的好处是直接写起来简单，缺点是需要练习这种（对许多人来说都是）新的语法、且难以通过脚本来自动化地更改参数的值。例如：

[14]:

data_file = "dataset/LP_machine.dat"
with open(data_file, "r") as f:
    data_rawstr = f.read()
print(data_rawstr)

param:  J:      c   a_machine   a_product   a_human :=
        "A"     10  4           1           2
        "B"     7   3           1           0.5
        "C"     3   1           1           1       ;
param ub_machine    :=  100;
param ub_product    :=  60;
param eb_human      :=  45;

接下来，我们用 create_instance() 实例化这个模型，其中选项 filename= 指定外部数据（或者用 data= 指定 Python 字典数据）。这里使用 YAML 文件：

[15]:

data_file = "dataset/LP_machine.yaml"
# 或者使用：model_instance = model.create_instance(data=datadict)
model_instance = model.create_instance(filename=data_file)
model_instance.is_constructed()

[15]:

True

[16]:

model_instance.pprint()

1 Set Declarations
    J : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {'A', 'B', 'C'}

7 Param Declarations
    a_human : Size=3, Index=J, Domain=PositiveReals, Default=None, Mutable=False
        Key : Value
          A :     2
          B :   0.5
          C :     1
    a_machine : Size=3, Index=J, Domain=PositiveReals, Default=None, Mutable=False
        Key : Value
          A :     4
          B :     3
          C :     1
    a_product : Size=3, Index=J, Domain=NonNegativeReals, Default=None, Mutable=False
        Key : Value
          A :     1
          B :     1
          C :     1
    c : Size=3, Index=J, Domain=NonNegativeReals, Default=None, Mutable=False
        Key : Value
          A :    10
          B :     7
          C :     3
    eb_human : Size=1, Index=None, Domain=PositiveReals, Default=None, Mutable=False
        Key  : Value
        None :    45
    ub_machine : Size=1, Index=None, Domain=PositiveReals, Default=None, Mutable=False
        Key  : Value
        None : 100.0
    ub_product : Size=1, Index=None, Domain=PositiveReals, Default=None, Mutable=False
        Key  : Value
        None :  60.0

1 Var Declarations
    x : Size=3, Index=J
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          A :     0 :  None :  None : False :  True : NonNegativeReals
          B :     0 :  None :  None : False :  True : NonNegativeReals
          C :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    objectives : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 10*x[A] + 7*x[B] + 3*x[C]

3 Constraint Declarations
    constraint_human : Size=1, Index=None, Active=True
        Key  : Lower : Body                     : Upper : Active
        None :  45.0 : 2*x[A] + 0.5*x[B] + x[C] :  45.0 :   True
    constraint_machine : Size=1, Index=None, Active=True
        Key  : Lower : Body                   : Upper : Active
        None :  -Inf : 4*x[A] + 3*x[B] + x[C] : 100.0 :   True
    constraint_product : Size=1, Index=None, Active=True
        Key  : Lower : Body               : Upper : Active
        None :  -Inf : x[A] + x[B] + x[C] :  60.0 :   True

13 Declarations: J c a_machine a_product a_human ub_machine ub_product eb_human x objectives constraint_machine constraint_product constraint_human

求解

在求解之前，先确保当前计算机上安装了求解器。可以使用 pyomo help --solvers 来查看可用的求解器。前带 + 号的是当前可用的求解器。

本例通过 pyomo.SolverFactory() 函数指定了 Cplex_direct 求解器（请确认安装），然后进行求解（并使用 tee=True 选项打印信息）。求解器的安装路径应当在 PATH 环境变量中，或者可以用 SolverFactory() 的选项 executable='path/to/dir' 指定位置。

[17]:

solver = pyo.SolverFactory('cplex_direct')
result = solver.solve(model_instance, tee=True)

Version identifier: 22.1.1.0 | 2023-06-15 | d64d5bd77
CPXPARAM_Read_DataCheck                          1
Tried aggregator 1 time.
No LP presolve or aggregator reductions.
Presolve time = 0.00 sec. (0.00 ticks)

Iteration log . . .
Iteration:     1   Dual infeasibility =             0.000000
Iteration:     2   Dual objective     =           256.000000

检查求解状态 .solver.status：

[18]:

print(f"Solver status: {result.solver.status}")

Solver status: ok

打印最优解结果：

[19]:

print(result)


Problem:
- Name:
  Lower bound: 256.0
  Upper bound: 256.0
  Number of objectives: 1
  Number of constraints: 3
  Number of variables: 3
  Number of binary variables: 0
  Number of integer variables: 0
  Number of continuous variables: 3
  Number of nonzeros: None
  Sense: -1
Solver:
- Name: CPLEX 22.1.1.0
  Status: ok
  Wallclock time: 0.0009968280792236328
  Termination condition: optimal
Solution:
- number of solutions: 0
  number of solutions displayed: 0

查看最优解的对应的变量值。你也可以把 pyo.Var 换成 pyo.Constraint 或 pyo.Param 来查看约束或者参数的情况。

[20]:

def get_opt_values(inst, component=pyo.Var):
    data = dict()
    for group in inst.component_objects(component, active=True):
        data[str(group)] = {k: pyo.value(v) for k, v in group.items()}
    return data

print(f"Variable:\n\t{get_opt_values(model_instance)}")

Variable:
        {'x': {'A': 0.0, 'B': 21.999999999999996, 'C': 34.0}}

查看约束的松弛情况。分别使用 lslack() 与 uslack() 来获取约束的下界与上界处的松弛信息：

[21]:

def get_slacks(inst):
    data = dict()
    for group in inst.component_objects(pyo.Constraint, active=True):
        data[str(group)] = (group.lslack(), group.uslack())
    return data

get_slacks(model_instance)

[21]:

{'constraint_machine': (inf, 1.4210854715202004e-14),
 'constraint_product': (inf, 4.0),
 'constraint_human': (0.0, 0.0)}

由于 0 值表示未松弛，我们可以看出：

machine 与 product 没有下界约束（因为它们是小于等于约束）
machine 的上界约束取到等号、human 的上下界约束均取到等号（因为是等式约束）
product 的上界约束尚有 4.0 的松弛量，即最优解中生产的产品数量比约束要少 4 个