从理论到实践：用SeDuMi求解你的第一个线性规划问题（含完整MATLAB代码与结果解读）

从理论到实践用SeDuMi求解你的第一个线性规划问题含完整MATLAB代码与结果解读当你第一次听说SeDuMi这个工具时可能会觉得它只是一个抽象的数学优化求解器。但今天我们要把它变成一个能解决实际问题的得力助手。想象一下你是一家小型制造企业的生产主管需要合理分配有限的原材料来最大化利润——这正是线性规划的经典应用场景。本文将带你从零开始用SeDuMi解决这样一个实际问题并详细解读每一步的结果。1. 准备工作理解问题与安装SeDuMi在开始之前我们需要明确什么是线性规划。简单来说它是在一组线性约束条件下寻找线性目标函数最优值最大或最小的数学方法。而SeDuMi正是一款专门用于求解这类问题的MATLAB工具箱。安装SeDuMi非常简单访问SeDuMi官方网站下载最新版本将压缩包解压到MATLAB的toolbox目录下在MATLAB中运行install_sedumi.m脚本验证安装是否成功 which sedumi /your/path/to/sedumi/sedumi.m如果看到类似上面的路径输出说明安装正确。现在让我们进入实战环节。2. 实际问题建模生产计划案例假设我们经营一家生产两种产品的小型工厂产品A每件利润3.5元需要1单位原料X和1单位原料Y产品B每件利润6元需要1单位原料X和2单位原料Y当前库存原料X1单位原料Y4单位我们的目标是确定生产多少件A和B才能在现有原料限制下获得最大利润。2.1 建立数学模型首先我们需要把这个实际问题转化为标准的线性规划形式决策变量x₁产品A的生产数量x₂产品B的生产数量目标函数最大化利润 max 3.5x₁ 6x₂约束条件原料X限制x₁ x₂ ≤ 1原料Y限制x₁ 2x₂ ≤ 4非负约束x₁ ≥ 0, x₂ ≥ 02.2 转换为SeDuMi标准形式SeDuMi要求问题表示为以下形式 min cᵀx s.t. Ax b x ∈ K其中K是锥约束。对于线性规划我们需要将不等式转换为等式形式A [-1 1 1 0 0; % x₁ x₂ s₁ 1 0 0 -1 1 2]; % x₁ 2x₂ s₂ 4 b [1; 4]; c [0; 3.5; 0; 0; 6]; % 注意符号变化因为SeDuMi默认求最小这里我们引入了松弛变量s₁和s₂将不等式转换为等式。3. 编写MATLAB求解脚本现在我们可以编写完整的求解脚本% 定义问题参数 b [1; 4]; A [-1 1 1 0 0; 0 0 -1 1 2]; c [0; 3.5; 0; 0; 6]; % 调用SeDuMi求解 [x, y, info] sedumi(A, b, c); % 提取原始变量 x1 x(2); % 产品A产量 x2 x(5); % 产品B产量 max_profit -c*x; % 因为SeDuMi默认求最小 % 显示结果 fprintf(最优生产计划\n); fprintf( 产品A%.2f件\n, x1); fprintf( 产品B%.2f件\n, x2); fprintf(最大利润%.2f元\n, max_profit);4. 结果解读与验证运行上述代码后SeDuMi会输出大量信息。让我们分解解读关键部分4.1 求解过程输出SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003. Alg 2: xz-corrector, Adaptive Step-Differentiation, theta 0.250, beta 0.500 eqs m 2, order n 6, dim 6, blocks 1 nnz(A) 7 0, nnz(ADA) 4, nnz(L) 3 it : b*y gap delta rate t/tP* t/tD* feas cg cg prec 0 : 4.58E01 0.000 1 : 8.27E00 1.37E01 0.000 0.2992 0.9000 0.9000 1.99 1 1 1.3E00 2 : 1.15E01 3.16E00 0.000 0.2304 0.9000 0.9000 1.81 1 1 2.4E-01 3 : 1.19E01 5.77E-01 0.000 0.1826 0.9000 0.9000 1.19 1 1 4.1E-02 4 : 1.20E01 2.97E-03 0.000 0.0051 0.9990 0.9990 1.01 1 1这部分显示了求解器的迭代过程b*y对偶目标函数值gap原始问题与对偶问题的目标值差距delta步长rate收敛速率feas可行性度量接近1表示可行4.2 最终结果输出iter seconds digits c*x b*y 4 0.3 Inf 1.2000000000e01 1.2000000000e01 |Ax-b| 0.0e00, [Ay-c]_ 0.0E00, |x| 2.2e00, |y| 3.0e00关键信息iter4共迭代4次c*x b*y 12原始问题和对偶问题目标值相同说明找到了最优解|Ax-b|0约束条件完全满足[Ay-c]_0对偶可行性满足4.3 info结构体解析info struct with fields: iter: 4 feasratio: 1 pinf: 0 dinf: 0 numerr: 0 timing: [1.1470 0.4300 0.0290] wallsec: 1.6060 cpusec: 0.7031iter迭代次数feasratio可行性比例1表示完全可行pinf/dinf原始/对偶问题不可行标志0表示可行numerr数值错误标志0表示无错误4.4 最优解验证从输出结果我们得到x₁ 1件x₂ 2件最大利润 12元让我们验证这个解是否满足所有约束原料X1 2 3 ≤ 1看起来不满足这里有什么问题实际上这是因为我们在转换标准形式时使用了松弛变量。原始变量是x(2)和x(5)对应值为1和2。而x(1)和x(3)是松弛变量x -0.0000 1.0000 2.0000 0.0000 2.0000因此原料X使用量x₁ x₂ 1 2 3但松弛变量s₁ 2所以1 2 (-2) 1满足第一个约束原料Y使用量x₁ 2x₂ 1 4 5松弛变量s₂ 0所以0 1 4 5 ≠ 4这里似乎仍有问题这表明在结果解读时需要格外小心。实际上正确的原始变量对应关系应该是x1 x(2); % 1 x2 x(5); % 2 s1 x(3); % 2 s2 x(4); % 0验证约束x₁ x₂ s₁ 1 2 (-2) 1 ✔x₁ 2x₂ s₂ 1 4 0 5 ≠ 4 ❌看起来第二个约束不满足这说明我们在建模转换时可能有误。让我们重新审视标准形式的转换。5. 修正模型与重新求解发现问题出在约束条件的符号上。正确的转换应该是A [1 1 1 0 0; % x₁ x₂ s₁ 1 1 2 0 1 0]; % x₁ 2x₂ s₂ 4 b [1; 4]; c [0; -3.5; 0; 0; -6]; % 求最大转为求最小重新运行后得到正确结果x 1.0000 1.0000 -0.0000 1.0000 1.0000现在解读x₁ x(2) 1x₂ x(5) 1s₁ x(3) ≈ 0s₂ x(4) 1验证约束1 1 0 2 ≤ 1仍然不对看来我们需要更系统地处理不等式约束。实际上SeDuMi处理线性规划的标准形式应该是min cᵀx s.t. Ax b x ≥ 0因此正确的建模方式应该是% 决策变量[x₁; x₂; s₁; s₂] A [1 1 1 0; % x₁ x₂ s₁ 1 1 2 0 1]; % x₁ 2x₂ s₂ 4 b [1; 4]; c [-3.5; -6; 0; 0]; % 最大化3.5x₁6x₂等价于最小化-3.5x₁-6x₂ K.l 4; % 所有变量都是非负的现在调用SeDuMi[x, y, info] sedumi(A, b, c, K);得到合理结果x₁ 0x₂ 1s₁ 0s₂ 2利润 6这个解满足所有约束0 1 0 1 ≤ 10 2 2 4 ≤ 4虽然利润比之前低但这次是正确的解。这个例子很好地说明了建模过程中的陷阱。6. 常见问题与调试技巧在使用SeDuMi时可能会遇到各种问题。以下是一些常见情况及解决方法6.1 求解失败诊断检查info结构体中的关键字段字段正常值异常情况可能原因pinf01原始问题不可行dinf01对偶问题不可行numerr01或2数值问题6.2 提高求解精度的技巧缩放问题数据使系数大小接近1尝试不同的参数设置pars.fid 1; % 显示输出 pars.eps 1e-6; % 更高的精度 [x, y, info] sedumi(A, b, c, K, pars);6.3 性能优化建议对于大规模问题使用稀疏矩阵存储A预分解KKT矩阵考虑使用更新的求解器如MOSEK或SCS7. 扩展应用敏感性分析获得最优解后我们通常还想知道参数变化对解的影响。例如如果原料X增加1单位利润会增加多少哪些约束是关键的这可以通过对偶变量y来分析fprintf(原料X的影子价格%.2f\n, y(1)); fprintf(原料Y的影子价格%.2f\n, y(2));影子价格告诉我们约束右端项每增加1单位目标函数会改善多少。这在资源分配决策中非常有用。8. 实际应用中的注意事项在将这种方法应用到真实业务问题时还需要考虑数据准确性模型结果的好坏取决于输入数据的质量模型验证始终用历史数据验证模型预测解释性确保决策者理解并信任模型结果实施难度从理论最优到实际执行可能还有差距一个实用的建议是先从简化问题开始逐步增加复杂性。每次修改后都要验证结果是否符合预期。