1. 线性规划与单纯形算法基础我第一次接触线性规划问题时是在大学运筹学课上。教授在黑板上画出一个二维可行域用等高线演示最优解的寻找过程。这种用数学方法解决资源分配问题的思路让我着迷而单纯形算法就是解决这类问题的经典方法。线性规划问题的标准形式包含三个关键要素决策变量、约束条件和目标函数。举个例子假设你经营一家小工厂生产两种产品A和B。产品A每件利润100元B每件利润150元。但生产过程中受到原料、工时等限制。如何安排生产计划使利润最大化这就是典型的线性规划问题。单纯形算法的核心思想是在可行域的顶点之间移动逐步逼近最优解。想象你在一片多面体山脉中寻找最高点每次沿着最陡峭的边向上爬直到找不到更高的位置为止。这个几何直观对应到代数操作就是通过基变换来迭代改进解的质量。两阶段单纯形算法的出现是为了处理更复杂的情况。当约束条件包含等式或大于等于形式时初始基本可行解可能不明显。第一阶段通过引入人工变量构造辅助问题第二阶段再用找到的可行解优化原问题。这就好比先修路再开车——第一阶段修建通往可行域的道路第二阶段沿着道路寻找最优解。2. 算法框架设计与数据结构当我第一次尝试实现单纯形算法时最头疼的就是数据结构的设计。经过多次重构最终采用了面向对象的方式封装核心组件。下面是我们需要准备的主要数据结构class LinearProgram { private: int m, n; // 约束总数和变量数 double **a; // 单纯形表二维数组 int *basic; // 基本变量索引 int *nonbasic; // 非基本变量索引 double minmax; // 1为max-1为min };单纯形表是这个实现的核心我们用二维数组a来存储。这里有个编程细节通常数学上的单纯形表从1开始索引而C数组从0开始。我的处理方式是第0行存储目标函数第0列存储右侧常数项这样既符合数学习惯又适应编程需求。初始化阶段需要特别注意约束条件的处理。对于不等式约束我们需要添加松弛变量对于等式约束可能需要人工变量。例如x1 2x3 ≤ 18 → 添加松弛变量s1x1 2x3 s1 18x1 x2 x3 x4 9 → 需要人工变量a1x1 x2 x3 x4 a1 9第一阶段的目标是消除这些人工变量为此我们构造辅助目标函数最小化所有人工变量之和。只有当这个最小值能为0时原问题才有可行解。3. 核心算法步骤实现3.1 入基变量选择选择入基变量就像在股市中选择最有潜力的股票——我们要找能让目标函数增长最快的方向。在代码中这体现为在z行目标函数行寻找最大正系数int LinearProgram::enter(int objrow) { double temp DBL_EPSILON; int col 0; for(int j1; jn1; j) { if(nonbasic[j]n2 a[objrow][j]temp) { col j; temp a[objrow][j]; break; // Bland法则避免循环 } } return col; }这里有几个实现细节值得注意使用DBL_EPSILON作为阈值避免浮点误差通过nonbasic[j]n2确保只考虑真正的决策变量和松弛变量采用Bland法则选择下标最小的候选变量防止算法陷入循环3.2 离基变量确定选定入基变量后我们需要确定哪个基本变量应该离开基。这就像确定登山时哪个装备需要留下——受限制最多的变量必须离基int LinearProgram::leave(int col) { double temp DBL_MAX; int row 0; for(int i1; im; i) { double val a[i][col]; if(val DBL_EPSILON) { val a[i][0]/val; // 计算theta比值 if(val temp) { row i; temp val; } } } return row; }比值测试(theta规则)确保转换后所有变量仍保持非负。如果发现入基变量所在列所有元素都非正说明问题无界——就像登山时发现可以无限向上没有最高点。3.3 转轴变换实现转轴变换是单纯形法的核心操作相当于在代数层面从一个顶点移动到相邻顶点。这个过程需要精确的矩阵行变换void LinearProgram::pivot(int row, int col) { // 标准化主元行 for(int j0; jn1; j) { if(j ! col) a[row][j] / a[row][col]; } a[row][col] 1.0/a[row][col]; // 更新其他行 for(int i0; im1; i) { if(i ! row) { for(int j0; jn1; j) { if(j ! col) { a[i][j] - a[i][col]*a[row][j]; if(fabs(a[i][j]) DBL_EPSILON) a[i][j] 0.0; } } a[i][col] -a[i][col]*a[row][col]; } } swapbasic(row, col); }这个实现中我特别注意了数值稳定性问题。当元素绝对值小于DBL_EPSILON时直接设为0避免浮点误差累积。转轴后记得更新基本变量和非基本变量的索引关系。4. 边界情况处理与优化4.1 两阶段法实现实际项目中遇到的线性规划问题往往没有明显的初始可行解。两阶段法就像分两步走的战略int LinearProgram::compute() { if(error 0) return error; // 第一阶段寻找初始可行解 if(m ! m1) { error phase1(); if(error 0) return error; } // 第二阶段优化原问题 return phase2(); }第一阶段构造辅助问题时需要将所有人工变量的和作为新目标函数// 在构造函数中构造辅助目标函数 for(int j1; jn; j) { double value 0.0; for(int im11; im; i) { value a[i][j]; } a[m1][j] value; // 辅助目标函数行 }4.2 退化与循环处理单纯形法可能遇到退化情况——某个theta值为0导致目标函数没有改进。我采用Bland法则来避免循环在入基变量选择时遇到第一个符合条件的变量就选择它在离基变量选择时遇到相同theta值选择下标更小的变量这种方法虽然可能稍微降低效率但能保证算法必定终止。在实际测试中对于规模不大的问题性能影响可以忽略。4.3 无解与无界判断完善的算法实现必须能够识别特殊状况void LinearProgram::solve() { error compute(); switch(error) { case 0: output(); break; case 1: cout 输入数据错误 endl; break; case 2: cout 无界解 endl; break; case 3: cout 无可行解 endl; } }无界解的判断发生在离基变量选择时如果入基变量列没有正元素无可行解的判断在第一阶段结束时如果人工变量之和不能降到0。5. 实际应用与性能考量将理论算法转化为实际可用的代码还需要考虑许多工程细节。我在实现过程中总结了几个关键点输入输出处理设计友好的数据输入格式很重要。我的实现支持以下结构1 // 1max, -1min 4 4 // 约束数m4, 变量数n4 2 1 1 // ≤约束2个约束1个≥约束1个 ... // 约束系数 ... // 目标函数系数数值稳定性单纯形法涉及大量浮点运算必须处理舍入误差。我采取的措施包括使用相对误差比较fabs(x) DBL_EPSILON而非x 0定期对接近0的元素进行归零处理在转轴运算中使用高精度计算顺序内存管理动态二维数组的创建和释放需要特别注意void LinearProgram::Make2DArray(double ** a, int m, int n) { a new double*[m]; for(int i0; im; i) a[i] new double[n]; } void LinearProgram::Delet2DArray(double ** a, int m, int n) { for(int i0; im; i) delete[] a[i]; delete[] a; }扩展性考虑虽然这个实现使用标准两阶段法但架构设计允许轻松扩展可以添加对偶单纯形法的实现支持灵敏度分析功能集成稀疏矩阵存储优化在实际应用中这个算法可以解决生产计划、资源分配、投资组合等多种优化问题。我曾用它优化过一个车间的排产问题将生产效率提升了约15%。当然对于超大规模问题可能需要考虑更高级的内点法或商业求解器但对于中等规模问题这个实现已经足够强大。