【Python】【原创】线性规划与单纯形法求解
本帖最后由 白冥 于 2025-3-6 17:49 编辑概述
线性规划（Linear Programming, LP）是数学优化领域的核心方法，用于在线性约束下寻找目标函数的极值（最大化或最小化）。
本贴旨在实现线性规划问题的求解，采用的是两阶段单纯形法（Two-Phase Simplex Method）。单纯形法是求解线性规划问题最为经典和高效的算法之一，适用于各类具有线性约束与线性目标函数的优化问题。本模块不仅支持传统的非负变量（positive variables），还支持变量取负或取任意实数值的情况，并对不等式约束（≤、≥）与等式约束（=）进行了相应的预处理和转化，以便将所有问题转换为标准的单纯形表形式进行求解。
通过它，你可以方便地将目标函数系数、约束矩阵和右侧常数项输入，同时指定各变量的取值范围以及各约束的类型，再指定求解问题的类型（最大化或最小化），代码将返回一个最优解、最优目标函数值以及求解状态（例如求解成功、无可行解或无界解）。
系统依赖与环境要求
本代码基于 Python 编写，建议使用 Python 3.8+（需支持 enum 和 match-case 语法）及以上版本。模块中使用了以下 Python 内置模块和标准库：
● enum：用于定义枚举类型，帮助用户在编程时对约束类型、变量取值范围、求解问题种类和求解状态进行标识。
● typing：主要用于定义类型提示（Type Hints），提高代码可读性和可维护性。
● math：主要用于使用无限大（inf）这一常量，用于判断无界解的情况。
目前代码没有使用任何第三方依赖，所有功能均基于 Python 标准库实现。但在后续扩展中，若需要实现更多高级功能，可考虑引入诸如 numpy、pandas、matplotlib 等常用库。
代码可在各种操作系统（如 Windows、Linux、macOS）上运行，只要安装了合适版本的 Python 环境即可。同时，由于算法涉及矩阵操作，若求解大规模线性规划问题时，建议在硬件性能较好的设备上运行。
线性规划与单纯形法基础
线性规划的标准形式：
目标函数：最大化或最小化 cᵀ·x
约束条件：Ax≤b （或≥, =）
变量约束：x≥0
标准形式要求所有约束为等式（Ax=b）且变量非负。代码的核心任务是将用户输入的任意约束（包括不等式、自由变量）转换为标准形式，从而应用单纯形法。
单纯形法概述：
单纯形法是一种迭代算法，通过在多维凸多面体的顶点间移动来寻找最优解。其核心步骤包括：
初始化：构造初始可行基解。
选择进入变量（Pivot Column）：根据目标函数系数选择能改进目标的方向。
选择离开变量（Pivot Row）：通过最小比值测试保持解的可行性。
高斯消元：更新单纯形表。
终止条件判断：所有目标系数非负（最小化问题）时达到最优。
两阶段法：
当原问题缺乏明显可行解时，需引入两阶段法：
阶段一：构造辅助问题，最小化人工变量之和，寻找初始可行基。
阶段二：若阶段一成功，用原目标函数继续求解。
主要类和枚举说明
在代码中，共定义了四个枚举类，分别用于描述约束条件、变量取值范围、求解问题类型以及求解状态。使用这些枚举类型的目的是使代码逻辑更加清晰和具备更高的可读性。
EqConstraint：
class EqConstraint(Enum):
LE = "le"
GE = "ge"
EQ = "eq" LE (≤)：表示约束为“小于等于”型约束
GE (≥)：表示约束为“大于等于”型约束
EQ (=)：表示约束为“等于”型约束
使用该枚举，可以在定义约束条件时明确标识各个约束的类型，在后续处理中能够根据不同约束类型分别添加松弛变量或人工变量。
VarConstraint：
class VarConstraint(Enum):
P = "positive"
N = "negative"
A = "all" P (positive)：表示变量为非负变量
N (negative)：表示变量为非正变量
A (all)：表示变量可取任意实数
由于标准单纯形法要求变量非负，因此对于取负或任意实数值的变量，代码通过变量转换的方式，将其转化为多个非负变量的组合来处理。
Problem：
class Problem(Enum):
MAX = "max"
MIN = "min" MAX：目标函数为最大化问题。
MIN：目标函数为最小化问题。
为了统一求解过程，代码将所有问题均转换为最大化问题。如果原问题为最小化，则通过对目标函数系数取负的方式转化为等价的最大化问题，求解结束后再调整结果符号。
State：
class State(Enum):
OK = "ok"
NFS = "no feasible solution"
UR = "unbounded range"
OK：表示求解成功，找到了最优解。
NFS (No Feasible Solution)：表示约束条件不满足，问题无可行解。
UR (Unbounded Range)：表示目标函数无限制增大（或减小），问题无界。
求解状态在函数返回结果中起到了指示求解是否成功的作用，你可以据此判断返回的最优解和目标值是否有效。
主函数
模块的核心函数为 lp，其定义如下：
def lp(c, A, b, x_card: List, A_card: List, problem: Problem): 参数
c：目标函数的系数列表，表示线性规划中目标函数的权重。例如，对于最大化 z=3x₁+5x₂ ，则c=
A：约束系数矩阵，每一行代表一个约束条件中各个变量的系数。其结构为二维列表，例如：
A = [
,
] 表示第一个约束为 1x₁+2x₂（与对应 b 值比较），第二个约束为 3x₁+1x₂。
b：约束右侧常数项的列表，表示各约束条件的常数项。与 A 中的每一行对应，例如 b = 表示第一个约束 1x₁+2x₂ ≤ 10（或者其他关系），以及第二个约束 3x₁+1x₂ ≤ 15。
x_card：一个列表，其中的每个元素均为 VarConstraint 枚举类型，用来指定每个变量的取值范围。其长度应与目标函数中变量个数一致。
若某个变量为非负，则应为 VarConstraint.P；
若为非正，则为 VarConstraint.N；
若变量允许取任意值，则为 VarConstraint.A。
A_card：一个列表，其中的每个元素均为 EqConstraint 枚举类型，用来指定每一条约束的类型。其长度应与约束数量（即 A 的行数）一致。
VarConstraint.LE 表示该约束为 “小于等于” 型；
VarConstraint.GE 表示 “大于等于” 型；
VarConstraint.EQ 表示 “等于” 型。
problem：Problem 枚举类型，表示当前问题为最大化问题（Problem.MAX）还是最小化问题（Problem.MIN）。
返回值说明：
函数返回一个元组，包含以下三部分内容：
solve：求得的最优解，类型为列表，包含原始变量的解。对于经过变量转换后的变量，会在最后转换回原始变量形式。
target：目标函数的最优值。若问题为最小化，函数会在计算结束后取负还原目标函数值。
State：求解状态，标记为 State.OK（求解成功）、State.NFS（无可行解）或 State.UR（无界解）。
初始化逻辑
主要分为以下几个步骤，对输入数据进行复制、目标函数调整、变量约束处理、引入松弛和剩余变量与右侧常数项调整，具体流程如下：
数据复制：
在函数开始时，为避免对原始数据造成修改，代码会将输入的目标函数系数 c、约束矩阵 A 以及常数项 b 进行复制：
A_eq = deepcopy(A)
b_eq = b[:]
c_original = c[:]
这样可以确保在后续处理中对 A、b、c 的修改不会影响调用者提供的原始数据。
目标函数的调整：
为了适应单纯形法的求解，代码在目标函数处理上做了以下工作：
若问题类型为 最大化（MAX），则将目标函数系数取负，转化为最大化问题求解。这是因为传统单纯形法主要针对最大化问题设计。
在最终求解结束后，再对目标函数值取负还原为最小化问题的真实目标值。
if problem == Problem.MAX:
c_original = [-c_i for c_i in c_original] 在最后返回结果前，对 target 值也进行相应处理。
变量约束处理（x_card）
由于单纯形法要求所有变量均为非负，因此对于输入中可能出现的非负以外的变量（包括负变量和可取任意值的变量），代码采取了如下处理策略：
对于非正变量（VarConstraint.N）：
将变量取负后，再作为非负变量处理。具体来说，对每个涉及该变量的约束系数和目标函数系数取负，等价于令y =-x，从而保证y≥0。
对于任意变量（VarConstraint.A）：
任意实数变量通常通过引入两个非负变量表示，即令x = x⁺ - x⁻，且x⁺,x⁻≥0。代码中通过在 A_copy 中插入新的变量列，并在 c_copy 中添加相应系数来实现这一转换。这样既保持了变量非负的要求，又能表示原始变量的全取值范围。
代码段如下：
p = 0
x_match = dict()
for col in range(cols):
x_match = p
if x_card == VarConstraint.N:
for a_k in A_copy:
a_k = -a_k
c_copy = -c_copy
elif x_card == VarConstraint.A:
for a_k in A_copy:
a_k.insert(p+1, -a_k)
c_copy.insert(p+1, -c_copy)
p += 1
p += 1 其中，变量 p 用来跟踪经过转换后新变量所在的索引位置，x_match 则用于记录原始变量与转换后变量之间的映射关系，便于最后恢复原始解。
约束条件处理（A_card）
松弛变量的引入：
对于不等式约束（非等式约束），需要引入额外的变量使得所有约束均转化为等式。
对于 “小于等于” 型约束（LE），引入松弛变量 s≥0 ，使得 aₖx + s = bₖ。
对于 “大于等于” 型约束（GE），引入剩余变量s≥0，并调整约束为aₖx - s = bₖ。
le_match = {}
for row in range(rows):
row_card = A_card
match row_card:
case EqConstraint.LE:
A_eq = for a_k in A_eq]
A_eq[-1] = 1
c_original +=
le_match = len(A_eq)
case EqConstraint.GE:
A_eq = for a_k in A_eq]
A_eq[-1] = -1
c_original +=
le_match 记录哪些行添加了松弛变量，这些变量可能在阶段一作为初始基变量。
右侧常数项的调整：为保证所有常数项 bₖ 均为非负数（这是单纯形法要求的一个前提），代码会遍历每个约束，如果对应的 bₖ 值为负，则对该行所有系数取负，并将 bₖ 的值也取负，从而使得不等式方向保持一致。
两阶段法的实现
构造辅助问题：
若原问题无可行初始基，需引入人工变量并构造辅助问题：
目标函数：最小化人工变量之和。
约束：原约束 + 人工变量构成的单位矩阵。
代码实现：
arti_rows = list(set(range(rows)) - set(le_match.keys()))
arti_T = [*row + + *(rows - row -1) for row in arti_rows]
arti = [ for i in range(len(arti_T)) ]for j in range(len(rows))]
A_aug =
c_auxiliary = *len(A_eq) + *len(arti_rows)
base = list(le_match.values()) + list(range(len(A_eq), len(A_aug))) arti_T 生成单位矩阵的转置，确保每行对应一个人工变量。
base 初始基包含松弛变量和人工变量。
构造单纯形表：
单纯形表是一个增广矩阵，包含约束系数、右侧常数和目标函数行。
z = c_auxiliary +
tableau = + z 最后一行是目标函数（辅助问题中为人工变量之和）。
关键算法细节
进入变量的选择：
代码中选择目标行中最小负系数：
min_neg = min([:-1] if z_i < 0])
entering = tableau[-1][:-1].index(min_neg) 若所有系数非负，达到最优解。
离开变量的选择（最小比值测试）：
计算每一行的θ=bₖ/entringₖ（当且仅当entringₖ>0）
theta = [ for row in range(rows)]
for row in range(rows):
if pivot_col > 0:
val = tableau[-1]/pivot_col
theta=
min_the = min(theta, key = lambda x: x) 若所有θ为无穷大，问题无界
枢轴变换：
以主元（pivot）为中心，归一化该行并消去其他行的对应列：
pivot = tableau
tableau = ]
for row in range(rows + 1):
if row != leaving:
factor = tableau
tableau = , tableau)] 此操作确保进入变量成为基变量，同时保持解的可行性。
终止
第一阶段终止条件：
若辅助问题最优值 > 0，原问题无可行解（返回 State.NFS）。
若最优值 = 0，准备阶段二。
第二阶段准备工作：
移除人工变量，将目标函数替换为原问题的系数。
tableau = + t_k[-1] for t_k in tableau]
tableau[-1] = c_original + 若基变量中仍含人工变量，说明存在冗余约束，可能返回无解。
第二阶段结束工作：
根据扩展后的变量还原原始解：
for row in range(rows):
p = x_match
match x_card:
case VarConstraint.P:
solve = p_solve
case VarConstraint.N:
solve = -p_solve
case VarConstraint.A:
solve = p_solve - p_solve 自由变量需合并其拆分后的正负部分。
使用示例
以下通过一个具体示例，说明如何调用 lp 函数求解一个线性规划问题。
假设需要求解如下线性规划问题：
max z = 3x₁+5x₂
满足约束：
1x₁+2x₂≤10
3x₁+1x₂≤15
x₁,x₂≥0
对应参数设置如下：
c =
A = [, ]
b =
x_card = （两变量均为非负）
A_card = （两条约束均为“小于等于”）
problem = Problem.MAX
调用示例代码：
# 导入必要的枚举类型和 lp 函数
from your_module import lp, VarConstraint, EqConstraint, Problem, State
c =
A = [
,
]
b =
x_card =
A_card =
problem = Problem.MAX
solution, optimal_value, status = lp(c, A, b, x_card, A_card, problem)
if status == State.OK:
print("最优解：", solution)
print("最优目标值：", optimal_value)
elif status == State.UR:
print("问题无界。")
elif status == State.NFS:
print("无可行解。") 运行结果应输出最优解及对应目标函数值。
代码结构图
lp()
├── 初始化参数与数据副本
├── 处理变量约束（正数/负数/自由）
├── 处理约束类型（添加松弛变量）
├── 构造辅助问题（两阶段法）
├── 单纯形法迭代
│ ├── 选择进入变量
│ ├── 选择离开变量
│ └── 更新单纯形表
└── 返回结果与状态码
源代码
from enum import Enum
from typing import List
from math import inf, isinf
from copy import deepcopy
class EqConstraint(Enum):
LE = "le"
GE = "ge"
EQ = "eq"
class VarConstraint(Enum):
P = "positive"
N = "negative"
A = "all"
class Problem(Enum):
MAX = "max"
MIN = "min"
class State(Enum):
OK = "ok"
NFS = "no feasible solution"
UR = "unbounded range"
def lp(c, A, b, x_card:List, A_card:List, problem:Problem):
rows = len(A)
cols = len(A)
A_eq = deepcopy(A)
b_eq = b[:]
c_original = c[:]
if problem == Problem.MAX:
c_original = [-c_i for c_i in c_original]
p = 0
x_match = dict()
for col in range(cols):
x_match = p
if x_card == VarConstraint.N:
for a_k in A_eq:
a_k = -a_k
c_original = -c_original
elif x_card == VarConstraint.A:
for a_k in A_eq:
a_k.insert(p+1, -a_k)
c_original.insert(p+1, -c_original)
p += 1
p += 1
le_match = {}
for row in range(rows):
row_card = A_card
match row_card:
case EqConstraint.LE:
A_eq = for a_k in A_eq]
A_eq[-1] = 1
c_original +=
le_match = len(A_eq)
case EqConstraint.GE:
A_eq = for a_k in A_eq]
A_eq[-1] = -1
c_original +=
for row in range(rows):
a_k = A_eq
if b < 0:
a_k = [-a_i for a_i in a_k]
b = -b
arti_rows = list(set(range(rows)) - set(le_match.keys()))
arti_T = [*row + + *(rows - row -1) for row in arti_rows]
arti = [ for i in range(len(arti_T)) ]for j in range(len(rows))]
A_aug =
c_auxiliary = *len(A_eq) + *len(arti_rows)
base = list(le_match.values()) + list(range(len(A_eq), len(A_aug)))
z = c_auxiliary +
tableau = + z
solve = * cols
for time in :
while True:
if all(z_i > 0 for z_i in tableau[-1][:-1]):
break
min_neg = min([:-1] if z_i < 0])
entering = tableau[-1][:-1].index(min_neg)
pivot_col = for t_k in tableau[:-1]]
if all(p_k <= 0 for p_k in pivot_col):
match problem:
case Problem.MAX:
return (solve, inf, State.UR)
case Problem.MIN:
return (solve, -inf, State.UR)
theta = [ for row in range(rows)]
for row in range(rows):
if pivot_col > 0:
val = tableau[-1]/pivot_col
theta=
min_the = min(theta, key = lambda x: x)
if not isinf(min_the):
match (time, problem):
case (0, Problem.MAX) :
return (solve, None, State.NFS)
case (0, Problem.MIN) :
return (solve, None, State.NFS)
case (1, Problem.MAX):
return (solve, inf, State.UR)
case (1, Problem.MIN):
return (solve, -inf, State.UR)
leaving = min_the
pivot = tableau
tableau = ]
for row in range(rows + 1):
if row != leaving:
factor = tableau
tableau = , tableau)]
base = entering
match time:
case 0:
if tableau[-1][-1] != 0:
return (solve, None, State.NFS)
tableau = + t_k[-1] for t_k in tableau]
tableau[-1] = c_original +
for row in range(rows):
var = base
if var < len(A_eq):
coeff = tableau[-1]
tableau[-1] = , tableau)]
for row in range(rows):
if base >= len(A_eq):
return (solve, None, State.NFS)
case 1:
p_solve = *len(A_eq)
for row in range(rows):
var = base
p_solve = tableau[-1]
match problem:
case Problem.MAX:
target = -tableau[-1][-1]
case Problem.MIN:
target = tableau[-1][-1]
for row in range(rows):
p = x_match
match x_card:
case VarConstraint.P:
solve = p_solve
case VarConstraint.N:
solve = -p_solve
case VarConstraint.A:
solve = p_solve - p_solve
return (solve, target, State.OK)
调试与错误处理
常见输入错误：
维度不匹配：确保 len(c) == cols 且 len(A_card) == rows。
无效枚举值：检查 x_card 和 A_card 是否均为合法枚举。
局限性
计算复杂度：最坏情况下为指数时间（但对实际问题通常高效）。
数值稳定性：浮点运算误差可能导致错误解
啊，本可来这是为了色色，为什么会有py强健本可的脑子啊{:6_198:} 大佬简直恐怖如斯。。在一片色色中简直了:curse: 看着有种知识从大脑滑过的感觉 好像是组合数学学过这些内容，不过那个课程不讲python的实现 hhh连运筹学都来了，代码还挺简略的，把中间过程加个print理论上就可以抄作业了{:6_167:} 这大佬也太厉害了呀，在一堆带颜色的内容里头，显得格外突出呢。 {:3_59:} 好详尽的知识贴，Python简直是我的一生之敌，看着脑子都疼了:dizzy:这几天报名了NECCS，希望到考试的时候可以过吧（保佑） 这是什么，{:4_88:}我在哪里？ 这素在？这一下子直接给我干回大三了，真的会被吓晕 我在gamemale上学习高数。嗯，一下回到了大学课堂上了。 我就围观一下楼主写代码，感觉老年人跟不上时代了啊，当年学的时候完全没有代码的事 非常感谢楼主提供的数学攻击{:4_92:}有种可以立刻睡着的美感 线性规划算法在算法题中见到的比较少，个人见到的多的还是DP动态规划，DP算法简直就是一生之敌，小心翼翼地递推也不一定对:) 原来泥潭是学术论坛{:4_88:} 相当棒的代码呢，让人确实感觉到数学与计算机代码之间的紧密联系呢 这是给我干哪来了，这还是泥潭吗（） 单纯形法确实很经典。学优化理论的时候作业就是用python写单纯形法的求解器 试图汲取知识，，，，，汲取失败{:4_100:} 在色色中显得清新脱俗{:4_94:}
页:
[1]
2