资讯动态

Java线性规划求解器实现:单纯形法与工程化数值稳定

发布时间:2026/11/6 6:46:40 来源:尧图企业网站定制
简介面向运筹学课程学习者与 Java 开发者这份线性规划算法实现以两阶段法求解线性规划问题重点处理模型标准化、人工变量引入以及最优解计算适合课程实验、算法复现与小型优化问题验证。压缩包共 11 个文件约 13KB主要由 Java 源码与编译后的 class 文件构成辅以 XML 配置、项目描述文件和忽略规则目录精简但保留了完整工程组织方式便于导入开发环境阅读和调试。已有 3481 人学习说明其在教学与自学场景中具有一定参考价值。读者可从中梳理目标函数、约束条件与决策变量的建模思路理解标准化、松弛变量、人工变量及两阶段求解流程并参考控制台参数录入、矩阵与列表存储、结果输出等实现细节为后续补充测试用例、优化数据结构或扩展求解能力提供可操作的代码基础。1. 排产配比算不明白的时候我选择自己写一个求解器先说一个具体场景饲料配方排产原料七八种约束包括蛋白上下限、库存、产能、成本四十几条不等式决策变量三十多个。业务侧一直用表格软件手动凑方案改一次参数就得重试半天还经常凑出个不满足约束的“看起来还行”的方案。想让业务系统在报价时实时算一遍把计算塞进 Java 服务里随手一搜才发现线性规划在运筹学教材里讲得清清楚楚能直接嵌进工程、能 new 出来调用的 Java 实现却得自己拼。这份 Java 版线性规划算法实现干的就是这件事把单纯形法、两阶段法、大 M 法从课本公式变成一组类输入一组约束和目标函数输出最优解、状态和迭代次数。它适合两类人一类是要把求解放进后端做实时测算的开发者一类是想照着能跑的代码把单纯形表真正吃透的人。2. 工程骨架一次 solve 调用背后的四个包算法代码最容易写成一坨几百行的solve()里又建表又算比值又判状态跑得通但改不动。这个项目一开始就把边界划清楚了模型只描述问题标准型只负责变换单纯形表只负责一次消元动作求解器负责串流程。分层不是为了好看是为了后面换算法的时候不用动模型加约束类型的时候不用动表。判断一个求解器写得行不行看它的调用链能不能一眼读到底LinearProgram进来建成标准型灌进SimplexTableau循环 pivot 到最优最后把解回填到原始变量上。中间任何一步出问题状态能定位到具体阶段而不是一句“求解失败”。2.1 目录结构与模块职责包结构基本对应上面那条调用链linprog/ ├── model/ # 问题描述纯数据 │ ├── LinearProgram.java # 目标函数 约束集合 max/min │ ├── Constraint.java # 一行约束系数、关系符、右端项 │ ├── Relation.java # LE / GE / EQ │ └── ObjectiveSense.java # MAX / MIN ├── core/ # 算法内核不依赖业务 │ ├── StandardForm.java # 标准型A、b、c、变量映射 │ ├── StandardFormBuilder.java # 模型 - 标准型的变换 │ ├── SimplexTableau.java # 单纯形表 pivot │ └── Numeric.java # EPS、比较、相对容差工具 ├── solver/ │ ├── LpSolver.java # 求解器接口 │ ├── SimplexSolver.java # 有初始可行基时的主流程 │ ├── TwoPhaseSolver.java # 两阶段法 │ ├── BigMSolver.java # 大 M 法 │ └── SolverOptions.java # 迭代上限、EPS、入基规则 ├── result/ │ ├── SolveStatus.java # OPTIMAL / INFEASIBLE / UNBOUNDED / ... │ └── SolveResult.java # 目标值、变量值、对偶值、迭代数 └── demo/ └── ProductionPlanDemo.java # 一个可运行的配比例子包职责改动频率model只放问题描述不含算法逻辑低加变量类型时才动core标准型变换、单纯形表、数值工具中改精度策略时动solver流程编排、状态判定中换算法时新增实现result输出结构低这个划分有个好处core完全不 import 任何业务类可以单独拿去跑单元测试solver依赖core但不依赖具体模型字段换一套约束类型时影响面可控。血泪经验是很多人图省事把标准型变换写在求解器里后面想加一个“变量有上下界”的支持得把整段 pivot 逻辑翻出来改非常容易改出玄学 bug。2.2 数据模型业务问题怎么翻译成标准型模型层要足够“笨”笨到只描述数学对象不掺任何求解假设。目标函数一个系数数组约束一条一条存关系符用枚举而不是字符串避免和≤混用导致的坑。public final class Constraint { private final double[] coeffs; // 长度必须等于变量个数 private final Relation relation; // LE / GE / EQ private final double rhs; // 右端项允许为负由标准型阶段归一化 public Constraint(double[] coeffs, Relation relation, double rhs) { this.coeffs coeffs.clone(); // 防御性拷贝防止外部数组被改 this.relation relation; this.rhs rhs; } public double[] coeffs() { return coeffs.clone(); } public Relation relation() { return relation; } public double rhs() { return rhs; } }LinearProgram持有ObjectiveSense、目标系数c和约束列表构造时校验所有约束的系数长度一致public final class LinearProgram { private final ObjectiveSense sense; private final double[] objective; private final ListConstraint constraints; public LinearProgram(ObjectiveSense sense, double[] objective, ListConstraint constraints) { for (Constraint c : constraints) { if (c.coeffs().length ! objective.length) { throw new IllegalArgumentException(约束系数维度与目标函数不一致); } } this.sense sense; this.objective objective.clone(); this.constraints List.copyOf(constraints); } // getter 省略 }参数说明objective的长度就是决策变量个数constraints里每条的长度必须跟它对齐这是整个项目唯一的硬性维度约束。校验放在构造器里而不是求解时是为了把“模型本身写错”和“算法算不出来”这两类问题分开——前者是调用方的 bug后者才是求解器的 bug。这个区分在排查阶段极其省事。2.3 求解器接口与调用链接口只暴露一个方法把所有可调项收进SolverOptions避免方法签名越来越长public interface LpSolver { SolveResult solve(LinearProgram lp, SolverOptions options); } public record SolverOptions( double eps, // 数值容差默认 1e-9 int maxIterations, // 迭代上限默认 10000 PivotRule pivotRule, // DANTZIG / BLAND / HARRIS boolean collectDual // 是否回填对偶值 ) { public static SolverOptions defaults() { return new SolverOptions(1e-9, 10_000, PivotRule.DANTZIG, true); } }一次完整调用是这样的LinearProgram lp new LinearProgram( ObjectiveSense.MAX, new double[]{12.0, 9.5, 7.0}, // 三种产品的单位利润 List.of( new Constraint(new double[]{2, 1, 1}, Relation.LE, 100), // 原料 A new Constraint(new double[]{1, 2, 3}, Relation.LE, 120), // 原料 B new Constraint(new double[]{1, 1, 1}, Relation.GE, 10) // 最低产量 )); SolveResult r new TwoPhaseSolver().solve(lp, SolverOptions.defaults()); System.out.println(r.status()); // OPTIMAL System.out.println(r.objective()); // 目标值 System.out.println(Arrays.toString(r.variableValues()));SolveResult用 record 承载字段包括status、objective、variableValues长度对齐原始变量个数不是标准型变量数、dualValues、iterations、phase1Objective两阶段法第一阶段的最优值排查不可行时非常有用。这里有个设计取舍值得说一句variableValues只回填原始变量松弛变量和人工变量的值不往外交。原因是松弛变量本身有业务含义剩余产能、富余库存但不同建模方式下位置不一样交给调用方自己从StandardForm里映射更不容易错。常见做法是在SolveResult里额外挂一个MapString, Double的业务字段映射这个可以按项目需要加。3. 单纯形法核心一张表怎么迭代出最优解单纯形法的直觉不复杂从一个可行基出发每次找一个能让目标变好的方向走一步走到不能再好为止。落到代码上所有的“走一步”都归一成一次主元消去所有的“还能不能更好”都归一成看检验数的符号。真正难的不是理解这个思路是把从业务模型到标准型、从标准型到初始基、从初始基到最优解的每一步变换都写对错一步结果就是错的而且错得悄无声息。这一章按求解顺序拆先看模型怎么变标准型再看入基出基怎么选然后是 pivot 的实现最后是无初始可行基时的两阶段法。3.1 标准型构造松弛、剩余、人工三种附加变量求解器内部只认一种形式min c^T y约束全为等式A y b所有变量y 0右端项b 0。业务模型到这一步要过四道变换任何一道漏了都可能导致结果错误原形式变换动作新增变量a·x b加松弛变量s 01 个a·x b减剩余变量s 0再加人工变量2 个a·x b直接加人工变量1 个b 0整行乘 -1同时翻转关系符0 个max c^T x目标系数取负结果再取负0 个最后那条b 0的归一化最容易被漏掉它必须发生在加松弛/人工变量之前因为关系符翻转会改变该加松弛还是剩余。public static StandardForm build(LinearProgram lp) { int nOrig lp.variableCount(); int m lp.constraints().size(); // 先做一遍归一化b 0 的行整行取反并翻转关系符 Listdouble[] rows new ArrayList(); ListRelation rels new ArrayList(); ListDouble rhss new ArrayList(); for (Constraint c : lp.constraints()) { double[] coeffs c.coeffs(); double rhs c.rhs(); Relation rel c.relation(); if (rhs 0) { for (int j 0; j coeffs.length; j) coeffs[j] -coeffs[j]; rhs -rhs; rel rel.flip(); // LE - GEEQ 不变 } rows.add(coeffs); rels.add(rel); rhss.add(rhs); } // 再统计附加列松弛/剩余占一列人工变量各占一列 int nTotal nOrig; int[] slackCol new int[m]; int[] artCol new int[m]; Arrays.fill(slackCol, -1); Arrays.fill(artCol, -1); for (int i 0; i m; i) { switch (rels.get(i)) { case LE - slackCol[i] nTotal; // s case GE - { slackCol[i] nTotal; // -s artCol[i] nTotal; // a } case EQ - artCol[i] nTotal; // a } } // 之后按列号填 A 矩阵、目标向量并记录每个原始变量的列下标 return assemble(rows, rels, rhss, lp, nOrig, nTotal, slackCol, artCol); }逻辑说明slackCol[i]记的是第 i 条约束的松弛/剩余变量列号artCol[i]是人工变量列号-1 表示没有。这两个数组是后面构造初始基、以及第二阶段删人工变量的关键索引。参数上nOrig是原始变量数nTotal是标准型总列数assemble里会按nOrig之后依次追加的原则填系数松弛变量在对应行填 1LE或 -1GE人工变量在对应行填 1。目标向量也要跟着扩原始变量位置填sense MAX ? -c[j] : c[j]松弛/剩余变量位置填 0人工变量位置在第一阶段填 1最小化人工变量之和第二阶段填 0或在塑成大 M 法时填一个大数 M。提示sense MAX时目标取负这一步务必在标准型构造里做一次、在结果回填时再做一次两个地方成对出现。只在构造时取负忘了回填结果符号就是反的而且最优解变量值是对的非常迷惑人。3.2 检验数与入基怎么判断还能不能更优单纯形表用一个二维数组存最后一行是目标行最后一列是右端项final class SimplexTableau { final int m; // 约束行数 final int n; // 标准型变量列数 final double[][] a; // (m 1) x (n 1)a[m] 是目标行a[*][n] 是 rhs final int[] basis; // 长度 mbasis[i] 是第 i 行的基变量列号 }最小化问题里目标行非基变量位置上的值就是检验数reduced cost。检验数为负说明把这个变量从 0 往上抬能让目标继续下降全非负就到最优了。入基列的选择规则直接决定收敛速度也决定会不会陷入循环。/** Dantzig 规则挑检验数最负的一列入基收敛快但可能循环 */ int enteringByDantzig() { int best -1; double bestVal -Numeric.EPS; for (int j 0; j n; j) { if (isBasic(j)) continue; // 基变量的检验数恒为 0跳过 double rc a[m][j]; if (rc bestVal) { bestVal rc; best j; } } return best; // 返回 -1 表示已达最优 }参数说明a[m][j]是第 j 列的检验数阈值用-Numeric.EPS而不是 0是为了避免极小负数被当成“还能优化”而导致几十次空转。isBasic(j)遍历basis数组判断O(m) 开销规模小的时候无所谓约束上千条时建议额外维护一个boolean[] inBasis做 O(1) 判断pivot 时同步更新。Bland 规则是另一个极端挑列号最小的负检验数列理论上保证不循环代价是迭代次数可能翻几倍。工程上常见的折中是先用 Dantzig 跑检测到目标值连续若干次下降幅度小于阈值时切到 Bland兼顾速度和稳定。3.3 最小比值与主元消去pivot 那一圈循环入基列定了接下来要决定这一列往上抬多高。限制来自那些系数为正的行抬得太高会让某个基变量变成负数破坏可行性。最小比值规则就是算出每行的“天花板”取最小的那个。/** 返回出基行号-1 表示该列无正系数问题无界 */ int leaving(int col, boolean useBlandTieBreak) { int best -1; double bestRatio Double.POSITIVE_INFINITY; for (int i 0; i m; i) { double aic a[i][col]; if (aic Numeric.EPS) continue; // 非正系数不构成上界 double ratio a[i][n] / aic; if (ratio bestRatio - Numeric.EPS) { bestRatio ratio; best i; } else if (useBlandTieBreak Math.abs(ratio - bestRatio) Numeric.EPS best 0 basis[i] basis[best]) { best i; // 比值相同取基变量列号更小的 } } return best; }逻辑说明aic EPS的行直接跳过因为这一列在这一行上的系数非正抬高变量不会让该行的基变量变负。如果所有行都被跳过best保持 -1说明这个方向可以无限走问题无界。useBlandTieBreak打开时比值相等的情况下选基变量列号更小的行这是 Bland 规则的另一半配合 Bland 入基能严格防循环。主元消去本身是整段代码里最“机械”的部分但也是最容易埋浮点残渣的地方void pivot(int r, int col) { double p a[r][col]; for (int j 0; j n; j) a[r][j] / p; // 主元行归一化 for (int i 0; i m; i) { // 注意 m目标行也要消 if (i r) continue; double f a[i][col]; if (Math.abs(f) Numeric.EPS) continue; for (int j 0; j n; j) a[i][j] - f * a[r][j]; a[i][col] 0.0; // 显式清零压住 1e-17 这类残渣 } basis[r] col; }参数说明r是出基行col是入基列p是主元。外层循环上界是m含目标行漏掉目标行是最常见的新手错误——表看着在动检验数却永远不更新跑满迭代上限返回“未收敛”。内层循环结束后把a[i][col]强行置 0是因为理论上它应该精确等于 0但浮点算出来可能是1e-17累积几百次迭代后这些残渣会让检验数出现假阴性误判为已最优。主循环public SolveResult solve() { for (int iter 0; iter options.maxIterations(); iter) { int col tableau.entering(); if (col 0) return buildOptimal(); // 无负检验数最优 int row tableau.leaving(col, options.pivotRule() PivotRule.BLAND); if (row 0) return SolveResult.unbounded(); tableau.pivot(row, col); } return SolveResult.iterationLimit(); }3.4 两阶段法第一阶段把人工变量赶出去上一节的流程有个前提初始基是可行的。但原始约束里只要出现或松弛变量凑不出单位阵就得靠人工变量搭一个“人造”的单位基。问题是这个基对应的解不一定满足原约束所以第一阶段先最小化人工变量之和如果最优值大于 0说明原问题无可行解等于 0说明人工变量已经全部退基可以进第二阶段。这里有个关键动作人工变量进基之后目标行上它们对应的位置不是 0必须先把目标行消成规范形否则检验数是错的。/** 把目标行里基变量对应的系数消成 0构造规范形 */ private void canonicalize(double[] phaseCost) { System.arraycopy(phaseCost, 0, a[m], 0, n); // 先铺目标行 a[m][n] 0.0; // 目标值初值 for (int i 0; i m; i) { int b basis[i]; if (Math.abs(phaseCost[b]) Numeric.EPS) continue; double f phaseCost[b]; for (int j 0; j n; j) a[m][j] - f * a[i][j]; } }逻辑说明phaseCost是当前阶段的目标系数向量。第一阶段里人工变量位置填 1、其余填 0第二阶段里原始变量位置填变换后的目标系数、人工变量位置填 0。遍历每一行如果该行基变量在目标向量里有非零系数就用整行去减把这个系数消掉。不做这一步人工变量明明在基里、检验数却显示成 1算法会误以为还能优化反复在人工变量上打转。第一阶段的收尾有两件事判断可行性以及把人工变量列从表里彻底清出去。if (phase1Optimum Numeric.EPS) { return SolveResult.infeasible(phase1Optimum); // 人工变量消不掉原问题不可行 } // 人工变量可能还留在基里取值为 0 的退化情况把这些行替换掉 for (int i 0; i m; i) { if (basis[i] firstArtificial) { int replace findNonArtificialColumn(i); // 找一个非人工列系数非零的 if (replace 0) continue; // 整行都是人工变量冗余约束可删 pivot(i, replace); } } canonicalize(phase2Cost);参数说明firstArtificial是人工变量的起始列号标准型构造时按nOrig 松弛/剩余 人工的顺序排所以人工变量总是连续的一段。findNonArtificialColumn在某行所有非人工列系数都为 0 时返回 -1意味着这行是冗余约束多个约束线性相关导致的直接忽略不影响结果但建议记进日志因为它通常说明建模时有重复约束。4. 数值稳定与退化double 撑不住的那几个地方纯理论推导里检验数是精确的分数比较大小没有歧义。换成 double 就全变了两个本该相等的数差 1e-16一个本该是 0 的检验数算出来 -3e-17。这些误差不会立刻爆出来而是让算法在几个几乎一样的解之间来回倒腾或者提前宣布最优。数值稳定性听起来像玄学实际上是这类算法最需要拿捏的地方。4.1 EPS 阈值与相对容差绝对阈值是最省事的做法一个大类常量管全部比较public final class Numeric { public static final double EPS 1e-9; /** 相对比较a 是否“大于” b */ public static boolean gt(double a, double b) { return a - b EPS * Math.max(1.0, Math.max(Math.abs(a), Math.abs(b))); } /** 判断是否为零兼顾量级 */ public static boolean isZero(double a) { return Math.abs(a) EPS; } }为什么不用纯绝对阈值目标值可能只有 0.5也可能有 1e7同一个1e-9对前者太松、对后者太严。乘上max(|a|, |b|, 1)做相对缩放能把不同量级的问题拉到同一标准上。实际项目里我一般会准备两套EPS_FEAS用于可行性判定约束是否满足EPS_OPT用于最优性判定检验数是否够负前者通常可以放宽到 1e-7因为输入数据本身往往只有几位有效数字。注意EPS不要设得比输入数据的精度还小。输入成本是 12.34两位小数你把EPS设成 1e-12算法会认真地区分 12.3400000001 和 12.34纯属给自己找麻烦还会让迭代次数暴涨。4.2 Harris 比值检验与 Bland 规则防循环退化是另一个经典麻烦多个约束在同一点交汇最小比值出现并列出基行选哪个都能走但选错就可能回到之前的状态形成循环。理论上 Bland 规则能保证不循环实践上还有 Harris 两阶段比值检验能同时改善数值稳定性和减少退化。/** Harris 两阶段比值检验先按放宽的右端项定阈值再在其中挑主元最大的 */ int leavingHarris(int col) { double maxRatio Double.NEGATIVE_INFINITY; for (int i 0; i m; i) { if (a[i][col] Numeric.EPS) continue; double ratio (a[i][n] Numeric.EPS_FEAS) / a[i][col]; if (ratio maxRatio) maxRatio ratio; } int best -1; double bestPivot -1; for (int i 0; i m; i) { if (a[i][col] Numeric.EPS) continue; double ratio a[i][n] / a[i][col]; if (ratio maxRatio a[i][col] bestPivot) { bestPivot a[i][col]; best i; } } return best; }逻辑说明第一轮用一个略微放大的右端项算出一个“安全上界”maxRatio第二轮只在比值不超过这个上界的行里挑主元绝对值最大的那个。主元越大消元时的误差放大越小得到的表数值质量越好。代价是要多扫一遍扫描成本只有 O(m)相对于 pivot 的 O(mn) 可以忽略。什么时候用哪套规则我的习惯是小规模约束几十条默认 Dantzig快一旦检测到目标值连续 20 次迭代改善小于 1e-9切 Bland 或 Harris 兜底对精度要求极高的场景直接上 Harris。4.3 有理数复算给结果加一道后悔药double 版的求解器再稳也有一种情况说不清结果到底对不对。尤其是当输入全是整数或简单分数时理论上最优解应该是个漂亮的分数跑出来却是 3.9999999997这时候你没法判断是浮点误差还是算法真的错了。我的做法是准备一个精简版的有理数实现只在调试和测试时用public final class Rational implements ComparableRational { private final long num; // 分子带符号 private final long den; // 分母恒为正 public static Rational of(long n, long d) { if (d 0) throw new ArithmeticException(分母为零); if (d 0) { n -n; d -d; } long g gcd(Math.abs(n), d); return new Rational(n / g, d / g); // 始终保持最简 } public Rational minus(Rational o) { return of(num * o.den - o.num * den, den * o.den); } public int signum() { return Long.signum(num); } Override public int compareTo(Rational o) { return Long.compare(num * o.den, o.num * den); } // plus / times / divide / gcd 省略 }参数说明num、den用long存储构造时约分并保证分母为正这样比较大小只需要交叉相乘。这套实现不适合跑大问题——分子分母会迅速膨胀到溢出——但用来复算一个 5 变量 8 约束的教学例子绰绰有余能精确判断最优解应该是13/3还是4。实际用法是写一个ExactSimplexSolver只实现 pivot 和检验数比较跑同一份模型把 double 版的结果按 1e-6 精度取整后跟有理数版对比。对不上就说明 double 版有 bug对得上就放心。这个对拍过程我一般会固化进单元测试成本不高收益是后面改任何数值策略都有一道兜底。5. 避坑排查五类跑不通的现场与处理求解器出问题时最难的是它不报错只给你一个看起来合理的错答案。下面这五条是实际调试中反复遇到的按“现象 → 原因 → 解决”整理遇到新问题可以先往这五类上套。第一类结果全 0或者目标值明显偏离预期。现象变量值全是 0status却返回OPTIMAL或者目标值是业务预期的相反数。原因max转min的取负做了两次或一次都没做。标准型构造时目标系数取负、结果回填时再取负这是成对操作如果构造时取了负、回填时又直接拿a[m][n]当目标值就只对了一半。另外a[m][n]存的始终是“最小化问题”的目标值对原始max问题必须再取一次负。解决把目标值回填单独抽一个方法mapObjective(double internal)里面只写一次符号逻辑所有返回路径都走它。加一条断言内部目标值取负后如果跟手工算的预期差太远直接抛异常而不是返回。第二类报无界但业务上明显有界。现象status UNBOUNDEDiterations很小就返回了。原因通常是某条约束的松弛变量符号搞反了。约束减剩余变量如果误写成加那一行的系数全为正等于给了变量一个可以无限增长的通道还有一种情况是关系符归一化时只翻了系数没翻关系符或者反过来。解决在标准型构造完成后打印一次 A 矩阵和 b 向量手工核对每条约束的附加变量符号。更稳的做法是加一个validateStandardForm()检查b是否全非负、每个松弛/剩余变量是否在对应行有且仅有一个 ±1 系数、每条GE约束是否都配了人工变量。这个校验在调试期常开上线后可以关掉。第三类迭代几百次不收敛或者日志里目标值反复横跳。现象跑满maxIterations返回ITERATION_LIMIT或者目标值在两个数之间来回跳。原因退化加 Dantzig 规则导致的循环。经典例子是 Beale 构造的那个循环模型约束数量不多但最小比值每次都并列Dantzig 规则会稳稳地转圈。另一种可能是EPS设得太小导致本该相等的比值被判断成不等选了一个数值上更差的行。解决切 Bland 规则验证一下——如果切过去能收敛基本可以确认是循环。长期方案是默认开启 Harris 比值检验或者在检测到目标改善停滞时自动切规则。顺带把maxIterations从 10000 降到 2000循环时会更快暴露而不是白白跑几十秒。第四类约束明明满足却判不可行。现象手工代入一组明显可行的解status却返回INFEASIBLE。原因第一阶段判定阈值问题。理论上第一阶段最优值应该精确为 0实际算出来是1e-8这个量级如果直接用 0判断就误判了。还有一种可能是人工变量残留在基里值为 0 的退化情况第二阶段没清理干净导致行数错位。解决第一阶段判定改用phase1Optimum EPS_FEAS并把这个值记进SolveResult不可行时打日志看它到底多大——如果只有 1e-8那就是阈值问题如果到了 0.5 这个量级那就是真的不可行回去查约束。人工变量残留那部分canonicalize之前一定要把基里的人工变量 pivot 出去。第五类小模型很快规模一上来就卡住。现象20 个变量 30 条约束毫秒级200 个变量 300 条约束跑了几十秒还没完。原因单纯形表是稠密的每次 pivot 的复杂度是 O(mn)迭代次数又大致随 m 增长。300 行 400 列的表一次 pivot 就是 12 万次乘加几百次迭代下来量级就上去了。这是算法本身的特性不是代码写错了。解决先降低迭代次数用对偶单纯形或热启动拿上一次的解当初始基应对“参数微调后重算”的场景往往几次迭代就收敛。真要上大规模稠密表这条路走到头了换稀疏存储或者直接上内点法——内点法的迭代次数对规模不敏感一般在几十次量级只是实现复杂度高不少。6. 二次开发从能跑通到能接业务拿到能跑通的结果只是起点。真要用起来很快会冒出“变量必须是整数”“要输出影子价格”“求解太慢要加缓存”这类需求。整个项目的分层决定了这些改造的难度不一样下面按改动成本从低到高排。6.1 三个稳定的扩展点扩展点一对偶值与灵敏度分析。最优单纯形表里目标行上松弛变量那一列的值取负就是对应约束的影子价格对偶变量。这个不需要改任何算法只要在buildOptimal()里按slackCol[i]索引取值回填即可。业务含义是“这条约束放松一个单位目标能改善多少”做资源分配时非常有用。扩展点二新增约束类型。比如变量有上下界。最小改动是在StandardFormBuilder里为每个有上界的变量预置一行x_j ub走现成的LE分支不动core。稍微优雅一点的做法是在Constraint上挂一个isBound()标记比价时特殊处理。前一种十几行代码就能跑后一种要动SimplexTableau的界处理逻辑。扩展点三换算法。LpSolver只有一个方法新增一个InteriorPointSolver或者DualSimplexSolver实现它调用方一行不改。这里的前提是SolveResult的字段足够通用——如果某个新算法想返回额外信息加字段而不是改语义否则老调用方会受影响。改造方向主要改动点大致工作量风险点对偶值回填buildOptimal()半天松弛列索引映射变量上下界StandardFormBuilder12 天列数膨胀后性能热启动新增solve(lp, warmStart)23 天初始基可行性校验整数变量新增分支定界外壳12 周节点管理与剪枝策略内点法新增求解器实现24 周数值稳定性比单纯形更敏感整数规划那个方向要特别注意分支定界不是改单纯形是在外层反复构造子问题再调它。所以架构上应该是BranchAndBoundSolver依赖LpSolver而不是继承它别把树搜索塞进 pivot 循环里那样后面想换底层求解器会非常难受。6.2 对拍与回归怎么证明改动没把算法改坏数值算法的回归测试不能只断言“结果等于 X”因为浮点误差会让相等判断时对时错。我一般用三层验证Test void productionPlanShouldMatchReference() { LinearProgram lp DemoModels.productionPlan(); SolveResult r new TwoPhaseSolver().solve(lp, SolverOptions.defaults()); assertEquals(SolveStatus.OPTIMAL, r.status()); assertEquals(437.5, r.objective(), 1e-6); // 目标值按容差比 // 关键一步把解代回原始约束逐条验证可行性 for (Constraint c : lp.constraints()) { double lhs dot(c.coeffs(), r.variableValues()); assertTrue(c.relation().holds(lhs, c.rhs(), 1e-6), () - 约束被违反: lhs lhs , rhs c.rhs()); } }第一层是状态和目标值按容差比较第二层是可行性校验把解代回每一条原始约束这一步比目标值更可靠——目标值算错往往只是精度问题约束违反了那就是实打实的 bug第三层是跟有理数版或暴力枚举顶点版对拍只在小规模模型上跑作为开发期的额外保险。改数值参数EPS、pivot 规则的时候我会把这三层全跑一遍。经验是只改EPS通常前两层就够动了 pivot 规则或者标准型构造第三层必须过因为这类改动影响的是路径而不只是终点小模型上最容易暴露逻辑错误。从那以后我每次改求解器核心哪怕只是把Math.abs(f) EPS换成Numeric.isZero(f)都会强制走一遍可行性校验和对拍——数值算法的 bug 不会当场炸它会安安静静地给你一个看起来合理的错答案等到业务侧拿着错误的排产方案去下单那就不是改几行代码能解决的事了。希望这些拆解能帮到你少走几年弯路。本文还有配套的精品资源点击获取

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑