1. 从一道简单的生产计划题说起如果你正在准备数学建模竞赛或者你的专业课程里涉及到运筹优化那么“线性规划”这个词你肯定不陌生。它听起来有点学术但说白了就是一种在给定限制条件下寻找最优方案的方法。比如一个工厂生产两种产品每种产品需要不同的原料和工时原料和工时是有限的每种产品利润不同。厂长想知道怎么安排生产计划才能让总利润最高这就是一个典型的线性规划问题。我第一次在数学建模国赛里用上线性规划是处理一个资源调度问题。当时手忙脚乱对着理论公式推导了半天代码写得磕磕绊绊。后来才发现有了Matlab这个工具很多复杂的计算和求解过程都能被极大地简化。Matlab里的linprog函数就像是一个封装好的“规划求解器”你只需要把问题“翻译”成它认识的语言——也就是目标函数和约束条件的系数矩阵它就能给你算出最优解。这不仅仅是省去了手算的麻烦更重要的是它让你能把精力从繁琐的计算中解放出来专注于问题本身的分析和建模。所以这篇内容我想从一个建模者的角度而不是纯数学的角度来聊聊怎么用Matlab搞定线性规划。我们会从最标准的形式入手理解每个参数的意义然后通过几个从易到难的例子看看它如何解决实际问题比如资源分配、投资组合甚至一些简单的机器学习问题。最后我还会分享一些我踩过的坑比如解不出来、解无界或者结果不符合常识时该怎么去排查和调整你的模型。无论你是数学建模的新手还是想更高效地使用Matlab进行优化计算希望这些经验能帮到你。2. 线性规划的标准形式与Matlab的“输入语言”在让Matlab帮你求解之前你们俩得先对一下“暗号”。Matlab的linprog函数只认一种标准形式的线性规划问题。如果你手里的问题长得不一样那就必须先给它“化化妆”变成标准样子。2.1 线性规划的标准形式Matlab的linprog函数要求问题必须是如下形式最小化一个线性目标函数服从于线性等式和不等式约束。 具体写出来是这样的求x最小化f^T * x满足A * x bAeq * x beqlb x ub这里每一个符号都代表一个矩阵或向量x这是我们要求解的决策变量向量。比如生产计划里x1和x2就代表两种产品的产量。f目标函数的系数向量。f^T * x展开就是f1*x1 f2*x2 ...。注意linprog默认是求最小值。如果你的原始问题是求最大利润那很简单把利润系数都乘以-1求这个新目标的最小值等价于求原目标的最大值。A和b线性不等式约束的系数矩阵和右端向量。A * x b可以表示很多限制比如“原料消耗总量不能超过库存”、“工作时间不能超过8小时”。Aeq和beq线性等式约束的系数矩阵和右端向量。比如“两种产品的产出比例必须为1:1”或者“所有投资资金必须全部用完”。lb和ub决策变量的下界和上界向量。这通常代表了物理意义比如产量不能为负lb [0, 0]或者某个设备的产能有上限。注意这是最基本、最完整的形式。在实际使用时A,Aeq,lb,ub这些参数如果不存可以用空矩阵[]来代替。但f是必须的。2.2 一个简单的翻译例子假设我们有一个经典的生产计划问题 一家工厂生产桌子和椅子。生产一张桌子需要4单位木材和2单位工时利润为50元。生产一把椅子需要2单位木材和1单位工时利润为30元。现有木材100单位工时60单位。 问如何安排生产生产多少桌子和椅子使得总利润最大第一步定义决策变量设x1为桌子产量x2为椅子产量。第二步建立数学模型目标最大化总利润Z 50*x1 30*x2约束木材约束4*x1 2*x2 100工时约束2*x1 1*x2 60非负约束x1 0,x2 0第三步翻译成Matlab标准形式因为linprog求最小我们将最大化问题转化为最小化minimize -Z -50*x1 -30*x2。所以f [-50; -30]。不等式约束A [4, 2; 2, 1]每一行是一个约束的系数b [100; 60]没有等式约束所以Aeq [], beq []。变量非负所以下界lb [0; 0]上界ub不设限用[]表示。现在这个生产计划问题就变成了Matlab能听懂的“语言”linprog(f, A, b, Aeq, beq, lb, ub)。接下来我们就可以在Matlab里求解了。3. 核心函数linprog详解与基础求解掌握了“翻译”规则我们就可以召唤Matlab的求解器了。linprog函数是优化工具箱里的核心它的基本调用语法就是我们上面提到的。但要想用好它还得了解它返回什么以及有哪些选项可以调整。3.1 linprog的基本调用与输出我们接着用上面的生产计划例子在Matlab中实现求解。% 定义参数 f [-50; -30]; % 目标函数系数求最小所以利润取负 A [4, 2; 2, 1]; % 不等式约束系数矩阵 b [100; 60]; % 不等式约束右端项 Aeq []; % 无等式约束 beq []; lb [0; 0]; % 变量下界 ub []; % 变量上界无限制 % 调用linprog求解 [x, fval, exitflag, output] linprog(f, A, b, Aeq, beq, lb, ub); % 输出结果 disp(最优生产计划); disp([桌子数量 x1 , num2str(x(1))]); disp([椅子数量 x2 , num2str(x(2))]); disp([最大利润为, num2str(-fval)]); % 注意fval是最小化目标函数的值取负得到最大利润 disp([求解器退出状态 exitflag , num2str(exitflag)]); disp(output.message);运行这段代码你会得到类似下面的结果最优生产计划 桌子数量 x1 20 椅子数量 x2 10 最大利润为1300 求解器退出状态 exitflag 1 Optimization terminated.exitflag 1表示求解器成功找到了最优解。输出参数的含义如下x最优解向量即决策变量的值。这里x [20; 10]表示生产20张桌子和10把椅子利润最大。fval在最优解x处目标函数的值。因为我们求的是min -50x1-30x2所以fval -1300。原问题的最大利润就是-fval 1300。exitflag求解器的退出状态。这是非常重要的诊断信息。1函数收敛到最优解x。0迭代次数超过options.MaxIterations或函数计算次数超过options.MaxFunctionEvaluations。-2无可行解。这意味着你给出的约束条件互相矛盾找不到任何一个点能同时满足所有约束。比如你要求x1 x2 5同时又要求x1 10这显然不可能。-3问题无界。这意味着在满足约束的条件下目标函数值可以无限地减小对于最小化问题。比如你的约束只有x1 0目标函数是min -x1那么x1越大目标函数值就越小没有下限。其他负值代表求解过程中遇到了其他错误。output一个结构体包含关于优化过程的详细信息比如迭代次数、算法、收敛信息等。output.message通常包含了可读的退出原因。3.2 求解选项options设置默认情况下linprog会使用内点法求解。但对于不同规模、不同特性的问题你可能需要调整求解选项。这通过optimoptions函数来设置。% 创建优化选项 options optimoptions(linprog, ... Display, iter, ... % 显示每次迭代信息 Algorithm, dual-simplex, ... % 使用对偶单纯形法 MaxIterations, 1000, ... % 最大迭代次数 OptimalityTolerance, 1e-8); % 最优性容差 % 使用选项重新求解 [x, fval, exitflag, output] linprog(f, A, b, Aeq, beq, lb, ub, [], options);常用的选项包括Display: 控制输出信息的详细程度。off不显示final默认只显示最终结果iter显示每次迭代信息用于调试。Algorithm: 选择求解算法。对于线性规划主要有dual-simplex默认对偶单纯形法。通常对中等规模问题非常高效尤其在重新求解一个类似问题时。interior-point-legacy内点法旧版。对于大规模稀疏问题可能更快。interior-point内点法新版R2015b以后。MaxIterations/MaxFunctionEvaluations: 设置迭代或函数计算的上限防止程序陷入无限循环。OptimalityTolerance/ConstraintTolerance: 优化终止的容差。如果结果对精度要求极高可以适当调小这些值如1e-10但可能会增加计算时间。对于初学者大部分问题使用默认选项即可。当你遇到求解失败、速度慢或者对结果精度有疑问时才需要考虑调整这些选项。4. 从理论到实战三类典型建模案例拆解懂了基本操作我们来看几个更贴近实际数学建模赛题的例子。这些例子会涉及到如何将文字描述转化为数学模型以及处理一些常见的变体。4.1 案例一资源分配与成本最小化不等式约束为主问题描述一个饲料公司需要生产一种混合饲料要求每公斤饲料中营养成分A至少含18单位B至少含12单位C至少含16单位。现有三种原料可供选择其每公斤成本及营养成分含量如下表原料成本元/公斤营养A单位营养B单位营养C单位甲4321乙2213丙3132问如何配比这三种原料才能在满足营养要求的前提下使每公斤饲料的成本最低建模与求解决策变量设每公斤饲料中使用甲、乙、丙三种原料的量分别为x1, x2, x3公斤。目标函数最小化总成本min Z 4*x1 2*x2 3*x3。约束条件营养A要求3*x1 2*x2 1*x3 18营养B要求2*x1 1*x2 3*x3 12营养C要求1*x1 3*x2 2*x3 16非负约束x1, x2, x3 0隐含约束这里没有要求总重量必须等于1公斤因为x1, x2, x3本身就是重量目标是最小化单位重量的成本。如果要求“配制恰好1公斤饲料”则需要增加等式约束x1 x2 x3 1此时x1, x2, x3代表的是比例。本题没有这个要求所以不添加。Matlab标准化注意约束是“大于等于”而linprog标准形式是“小于等于”。我们需要两边乘以-1来转换。例如3x12x2x3 18等价于-3x1 -2x2 -x3 -18。% 饲料配比问题 - 成本最小化 f [4; 2; 3]; % 成本系数直接求最小 % 不等式约束 ( 转化为 ) A [-3, -2, -1; % 对应营养A: -3x1-2x2-x3 -18 -2, -1, -3; % 对应营养B: -2x1-x2-3x3 -12 -1, -3, -2]; % 对应营养C: -x1-3x2-2x3 -16 b [-18; -12; -16]; % 无等式约束 Aeq []; beq []; % 变量非负 lb [0; 0; 0]; ub []; [x_opt, fval_opt, exitflag] linprog(f, A, b, Aeq, beq, lb, ub); if exitflag 0 disp(最优原料配比公斤/每公斤饲料); disp([甲: , num2str(x_opt(1))]); disp([乙: , num2str(x_opt(2))]); disp([丙: , num2str(x_opt(3))]); disp([最低成本元/每公斤饲料: , num2str(fval_opt)]); else disp(求解失败请检查模型和约束。); end运行后你可能会得到类似x_opt [2; 0; 6]的结果表示每公斤饲料使用2公斤甲和6公斤丙总成本为4*23*626元。注意这里x1x2x38意味着为了满足高营养要求需要8公斤原料混合出“1公斤”符合要求的饲料这在逻辑上是“浓缩”饲料模型本身是合理的。如果要求总重为1则需要添加等式约束结果和意义会完全不同。4.2 案例二投资组合优化含等式约束与边界问题描述你有100万资金打算投资到三种资产股票A、债券B、基金C。已知期望年收益率股票A为10%债券B为5%基金C为8%。风险系数假设股票A为8债券B为2基金C为5。你希望总投资风险不超过500风险系数加权和并且为了分散风险规定对基金C的投资比例不低于20%。同时你与银行有一个协议必须将总资金的10%用于购买该银行发行的一款特殊债券D收益率3%这款债券必须单独考虑。问如何分配资金包括债券D在满足风险、比例和全额投资的要求下使期望年收益最大建模与求解决策变量设投资于股票A、债券B、基金C、特殊债券D的资金分别为x1, x2, x3, x4万元。目标函数最大化总期望收益max Z 0.10*x1 0.05*x2 0.08*x3 0.03*x4。约束条件总投资额x1 x2 x3 x4 100全部资金用完等式约束风险约束8*x1 2*x2 5*x3 500特殊债券D风险假设为0基金C比例约束x3 0.2 * 100即x3 20特殊债券D约束x4 0.1 * 100即x4 10这也是一个等式约束非负约束x1, x2, x3, x4 0Matlab标准化目标函数取负求最小。注意现在我们有两个等式约束总资金约束和债券D约束。% 投资组合问题 f -[0.10; 0.05; 0.08; 0.03]; % 收益系数求最大转为求最小 % 不等式约束风险约束和基金C下限 A [8, 2, 5, 0; % 风险约束系数 0, 0, -1, 0]; % x3 20 等价于 -x3 -20 b [500; -20]; % 等式约束总资金和债券D Aeq [1, 1, 1, 1; % x1x2x3x4 100 0, 0, 0, 1]; % x4 10 beq [100; 10]; % 变量非负 lb [0; 0; 0; 0]; ub []; % 无上限 [x_opt, fval_opt] linprog(f, A, b, Aeq, beq, lb, ub); disp(最优投资方案万元); disp([股票A: , num2str(x_opt(1))]); disp([债券B: , num2str(x_opt(2))]); disp([基金C: , num2str(x_opt(3))]); disp([特殊债券D: , num2str(x_opt(4))]); disp([最大年期望收益万元: , num2str(-fval_opt)]);这个例子展示了如何同时处理等式和不等式约束以及如何用矩阵Aeq和beq来表达多个等式关系。4.3 案例三数据拟合与简单机器学习线性回归作为线性规划线性规划不仅能解决资源分配问题还可以用于某些类型的数据拟合。考虑一个简单的中位数回归L1回归问题。与普通最小二乘法L2最小化误差平方和不同L1回归最小化误差的绝对值之和对异常值更不敏感。问题描述给定一组数据点(x_i, y_i)我们想用一条直线y a*x b来拟合。最小二乘法的目标是min Σ|y_i - (a*x_i b)|^2而L1回归的目标是min Σ|y_i - (a*x_i b)|。这个绝对值目标函数不是线性的但可以通过一个经典技巧转化为线性规划。转化技巧对于每个数据点i我们引入两个非负辅助变量u_i和v_i令y_i - (a*x_i b) u_i - v_i 且u_i 0, v_i 0。 那么绝对值|y_i - (a*x_i b)|就可以用u_i v_i来等价表示。因为对于任意实数它都可以表示成两个非负数的差而其绝对值就是这两个非负数的和。这样原问题就转化为决策变量a, b, u_1, v_1, u_2, v_2, ..., u_n, v_n目标函数min (u_1v_1 u_2v_2 ... u_nv_n)约束条件对于每个i有a*x_i b u_i - v_i y_i且u_i, v_i 0。a, b无约束可为负。% 生成带异常值的示例数据 rng(0); % 固定随机种子使结果可重复 x (1:20); y_true 2*x 5; y y_true randn(20,1)*3; % 加入正态噪声 y(10) y(10) 30; % 在第10个点加入一个异常值 % 使用线性规划进行L1回归 n length(x); % 决策变量顺序: [a; b; u1; v1; u2; v2; ...; un; vn] % 目标函数系数: a和b的系数为0u_i和v_i的系数为1 f [0; 0; ones(2*n, 1)]; % 等式约束: a*x_i b u_i - v_i y_i Aeq zeros(n, 2 2*n); for i 1:n Aeq(i, 1) x(i); % a 的系数 Aeq(i, 2) 1; % b 的系数 Aeq(i, 22*(i-1)1) 1; % u_i 的系数 Aeq(i, 22*(i-1)2) -1; % v_i 的系数 end beq y; % 不等式约束: 无 A []; b []; % 边界约束: a, b 无界u_i, v_i 0 lb [-inf; -inf; zeros(2*n, 1)]; ub []; % 求解线性规划 options optimoptions(linprog, Display, off); [z, fval] linprog(f, A, b, Aeq, beq, lb, ub, [], options); a_l1 z(1); b_l1 z(2); % 对比普通最小二乘法L2回归 X [x, ones(size(x))]; coeff_l2 X \ y; a_l2 coeff_l2(1); b_l2 coeff_l2(2); % 绘图对比 figure; scatter(x, y, bo, DisplayName, 数据点含异常值); hold on; plot(x, a_l1*x b_l1, r-, LineWidth, 2, DisplayName, L1回归线性规划); plot(x, a_l2*x b_l2, g--, LineWidth, 2, DisplayName, L2回归最小二乘); plot(x, y_true, k:, LineWidth, 1.5, DisplayName, 真实关系); legend(Location, best); xlabel(x); ylabel(y); title(L1回归 vs L2回归对异常值的鲁棒性); grid on; hold off; disp([L1回归结果: y , num2str(a_l1), * x , num2str(b_l1)]); disp([L2回归结果: y , num2str(a_l2), * x , num2str(b_l2)]);运行这段代码你会直观地看到由于异常值的存在绿色的L2回归线被明显拉偏而红色的L1回归线则更接近真实的黑色虚线。这个例子展示了线性规划在稳健统计中的应用也体现了将非线性的绝对值问题转化为线性规划的巧妙思路。在数学建模中这种转化能力非常重要。5. 求解失败问题无界或无解常见错误排查指南在实际使用中你可能会遇到linprog返回负的exitflag比如-2或-3这意味着求解失败。别慌这通常是你的模型出了问题而不是Matlab的错。下面是一些排查思路。5.1 问题无界exitflag -3表现求解器告诉你目标函数值可以无限小对于最小化问题。可能原因与排查缺少关键约束这是最常见的原因。检查你是否漏掉了对决策变量的必要限制。例如在最大化利润的生产问题中如果你忘记了原料或工时的约束产量就可以无限大利润也就无限大对应最小化问题的负利润无限小。约束方向错误检查你的不等式约束符号。如果你本意是消耗 资源但不小心写成了消耗 资源并且资源量很小这可能导致变量为了满足这个“宽松”的约束而变得非常大。变量无下界对于求最小化的问题如果目标函数系数f中存在负数且对应的变量没有下界lb为-inf则该变量趋向负无穷时目标函数值也会趋向负无穷。调试方法首先简化问题。尝试去掉一些约束或者固定一些变量的值看问题是否变得有解。这能帮你定位是哪个约束缺失或错误。其次检查模型假设。回到问题描述确保每一个现实中的限制条件都在模型中有对应的数学表达式。最后可视化对于二维或三维问题。你可以手动绘制约束条件围成的区域可行域看看它是否在一个方向上是不封闭的。例如只有x0和y0两个约束可行域就是第一象限对于目标函数min -x-y解就是无界的。5.2 无可行解exitflag -2表现求解器找不到任何一个点能同时满足所有约束条件。可能原因与排查约束条件互相矛盾这是根本原因。比如你要求x1 x2 5同时又要求x1 10且x2 0。这两个条件不可能同时成立。等式约束过于严格多个等式约束可能共同定义了一个空集。例如x1 x2 10和x1 x2 20显然矛盾。边界约束与不等式/等式约束矛盾lb和ub指定的范围可能与A*xb或Aeq*xbeq定义的范围没有交集。调试方法逐一注释法这是最实用的方法。暂时注释掉一部分约束比如先注释所有等式约束Aeq, beq运行看是否有解。然后逐步加入约束直到加入某一条约束后问题变得不可行那么这条约束很可能就是矛盾的源头。检查数据输入错误非常常见仔细核对A,b,Aeq,beq,lb,ub中的每一个数值。一个数字输错比如把100打成1000就可能导致可行域为空。检查单位一致性确保所有约束方程中的量纲是一致的。例如一边是“公斤”另一边不能是“吨”。使用linprog的‘diagnostics’选项设置options optimoptions(linprog, Display, iter, Diagnostics, on)更详细的输出有时能提供线索。5.3 其他常见错误与技巧结果不直观或不符合预期即使exitflag1解也可能很奇怪。比如投资比例加起来超过了1或者产量是小数但现实中必须是整数。这时你需要检查模型是否完整是否漏掉了“所有投资比例之和为1”这样的等式约束考虑整数规划如果变量必须取整如生产设备台数线性规划允许非整数解的特性会导致结果失真。这时应使用混合整数线性规划函数是intlinprog。检查目标函数系数确认系数的正负号是否正确最大化 vs 最小化。大规模问题的性能对于变量和约束成千上万的大规模问题默认的内点法可能内存消耗较大。可以尝试切换算法为‘dual-simplex’它对某些稀疏结构问题更有效。同时确保你的A和Aeq矩阵以稀疏矩阵格式存储使用sparse函数可以极大提升求解速度和降低内存使用。灵敏度分析与影子价格linprog的输出不直接提供影子价格对偶变量。但你可以通过求解对偶问题或者使用linprog的完整输出格式[x, fval, exitflag, output, lambda] linprog(...)来获取lambda结构体。lambda.ineqlin和lambda.eqlin分别给出了不等式和等式约束的影子价格这在经济解释和“如果资源增加一单位利润能增加多少”的边际分析中非常有用。6. 进阶从线性规划到整数规划与实战建议线性规划假设变量可以取任意实数但现实中的很多问题要求变量是整数比如“购买几台机器”、“派遣几个员工”、“一个项目是否启动0或1”。这时就需要用到整数规划或混合整数线性规划。6.1 使用intlinprog处理整数约束Matlab中求解混合整数线性规划的函数是intlinprog它的语法和linprog非常相似多了一个指定哪些变量需要取整数的参数intcon。语法[x, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub)参数intcon: 一个整数向量指定哪些决策变量必须取整数值。例如如果x(1)和x(3)是整数变量则intcon [1, 3]。例子背包问题假设有一个容量为10的背包有4件物品其价值和体积如下表。每件物品只能选择带或不带0-1变量如何选择物品使得总价值最大物品价值体积152263341474% 0-1背包问题 f -[5; 6; 4; 7]; % 最大化价值转为最小化负价值 intcon 1:4; % 所有变量都是0-1整数 % 体积约束 A [2, 3, 1, 4]; b 10; % 无等式约束 Aeq []; beq []; % 0-1变量边界 lb zeros(4,1); ub ones(4,1); [x_opt, fval_opt] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub); disp(最优选择1表示携带0表示不携带); disp(x_opt); disp([最大总价值, num2str(-fval_opt)]);intlinprog的求解比linprog复杂得多计算时间随问题规模指数级增长。对于大规模整数规划问题可能需要设置更长的求解时间限制options.MaxTime或调整分支定界法的相关选项。6.2 数学建模竞赛中的实战建议结合我参加和辅导数学建模竞赛的经验在比赛中使用线性/整数规划时有几点特别重要模型第一编程第二花足够的时间把问题理清建立正确的数学模型。一个错误的模型即使用最高级的算法也得不到正确的答案。在写代码前最好先在纸上把决策变量、目标函数、所有约束清晰地列出来。从简单开始逐步复杂化不要试图一口气建立一个包含所有细节的完美模型。先建立一个最简化的核心模型比如忽略一些次要约束用linprog快速求解验证基本逻辑是否正确。然后逐步添加更复杂的约束如整数约束、非线性约束的线性近似、多阶段决策等。重视结果的分析与解释数学建模不是解出答案就结束了。你需要分析结果敏感性分析如果某个参数如资源量、价格变化了10%最优解会变化多少这可以通过改变b或f中的值重新求解来实现。影子价格的经济意义对于资源约束其影子价格代表了该资源的边际价值。在论文中解释这个值能为你的方案提供深刻的洞察。结果的合理性得到的解是否在物理意义或常识上是合理的如果产量是负数或者投资比例加起来不等于1一定要回头检查模型。注意性能对于变量较多比如成千上万的问题在构建A和Aeq矩阵时尽量使用稀疏矩阵sparse(i, j, v, m, n)来存储可以节省大量内存和计算时间。做好错误处理在你的脚本中一定要检查exitflag。如果求解失败可以尝试输出中间变量或者给出一个备用的启发式方案并在论文中说明求解的局限性。一个健壮的代码能让你在紧张的比赛时间里更从容。最后线性规划是数学建模中极其强大的工具但它只是一个工具。真正考验你的是如何将一个复杂的现实问题抽象成一个简洁而准确的数学模型。Matlab的优化工具箱为你解决了计算上的难题让你能更专注于建模本身。多练习多思考下次当你遇到“在有限条件下寻求最优”的问题时不妨先想想“这能不能用一个线性规划来描述”
Matlab线性规划建模实战：从生产计划到投资组合优化
1. 从一道简单的生产计划题说起如果你正在准备数学建模竞赛或者你的专业课程里涉及到运筹优化那么“线性规划”这个词你肯定不陌生。它听起来有点学术但说白了就是一种在给定限制条件下寻找最优方案的方法。比如一个工厂生产两种产品每种产品需要不同的原料和工时原料和工时是有限的每种产品利润不同。厂长想知道怎么安排生产计划才能让总利润最高这就是一个典型的线性规划问题。我第一次在数学建模国赛里用上线性规划是处理一个资源调度问题。当时手忙脚乱对着理论公式推导了半天代码写得磕磕绊绊。后来才发现有了Matlab这个工具很多复杂的计算和求解过程都能被极大地简化。Matlab里的linprog函数就像是一个封装好的“规划求解器”你只需要把问题“翻译”成它认识的语言——也就是目标函数和约束条件的系数矩阵它就能给你算出最优解。这不仅仅是省去了手算的麻烦更重要的是它让你能把精力从繁琐的计算中解放出来专注于问题本身的分析和建模。所以这篇内容我想从一个建模者的角度而不是纯数学的角度来聊聊怎么用Matlab搞定线性规划。我们会从最标准的形式入手理解每个参数的意义然后通过几个从易到难的例子看看它如何解决实际问题比如资源分配、投资组合甚至一些简单的机器学习问题。最后我还会分享一些我踩过的坑比如解不出来、解无界或者结果不符合常识时该怎么去排查和调整你的模型。无论你是数学建模的新手还是想更高效地使用Matlab进行优化计算希望这些经验能帮到你。2. 线性规划的标准形式与Matlab的“输入语言”在让Matlab帮你求解之前你们俩得先对一下“暗号”。Matlab的linprog函数只认一种标准形式的线性规划问题。如果你手里的问题长得不一样那就必须先给它“化化妆”变成标准样子。2.1 线性规划的标准形式Matlab的linprog函数要求问题必须是如下形式最小化一个线性目标函数服从于线性等式和不等式约束。 具体写出来是这样的求x最小化f^T * x满足A * x bAeq * x beqlb x ub这里每一个符号都代表一个矩阵或向量x这是我们要求解的决策变量向量。比如生产计划里x1和x2就代表两种产品的产量。f目标函数的系数向量。f^T * x展开就是f1*x1 f2*x2 ...。注意linprog默认是求最小值。如果你的原始问题是求最大利润那很简单把利润系数都乘以-1求这个新目标的最小值等价于求原目标的最大值。A和b线性不等式约束的系数矩阵和右端向量。A * x b可以表示很多限制比如“原料消耗总量不能超过库存”、“工作时间不能超过8小时”。Aeq和beq线性等式约束的系数矩阵和右端向量。比如“两种产品的产出比例必须为1:1”或者“所有投资资金必须全部用完”。lb和ub决策变量的下界和上界向量。这通常代表了物理意义比如产量不能为负lb [0, 0]或者某个设备的产能有上限。注意这是最基本、最完整的形式。在实际使用时A,Aeq,lb,ub这些参数如果不存可以用空矩阵[]来代替。但f是必须的。2.2 一个简单的翻译例子假设我们有一个经典的生产计划问题 一家工厂生产桌子和椅子。生产一张桌子需要4单位木材和2单位工时利润为50元。生产一把椅子需要2单位木材和1单位工时利润为30元。现有木材100单位工时60单位。 问如何安排生产生产多少桌子和椅子使得总利润最大第一步定义决策变量设x1为桌子产量x2为椅子产量。第二步建立数学模型目标最大化总利润Z 50*x1 30*x2约束木材约束4*x1 2*x2 100工时约束2*x1 1*x2 60非负约束x1 0,x2 0第三步翻译成Matlab标准形式因为linprog求最小我们将最大化问题转化为最小化minimize -Z -50*x1 -30*x2。所以f [-50; -30]。不等式约束A [4, 2; 2, 1]每一行是一个约束的系数b [100; 60]没有等式约束所以Aeq [], beq []。变量非负所以下界lb [0; 0]上界ub不设限用[]表示。现在这个生产计划问题就变成了Matlab能听懂的“语言”linprog(f, A, b, Aeq, beq, lb, ub)。接下来我们就可以在Matlab里求解了。3. 核心函数linprog详解与基础求解掌握了“翻译”规则我们就可以召唤Matlab的求解器了。linprog函数是优化工具箱里的核心它的基本调用语法就是我们上面提到的。但要想用好它还得了解它返回什么以及有哪些选项可以调整。3.1 linprog的基本调用与输出我们接着用上面的生产计划例子在Matlab中实现求解。% 定义参数 f [-50; -30]; % 目标函数系数求最小所以利润取负 A [4, 2; 2, 1]; % 不等式约束系数矩阵 b [100; 60]; % 不等式约束右端项 Aeq []; % 无等式约束 beq []; lb [0; 0]; % 变量下界 ub []; % 变量上界无限制 % 调用linprog求解 [x, fval, exitflag, output] linprog(f, A, b, Aeq, beq, lb, ub); % 输出结果 disp(最优生产计划); disp([桌子数量 x1 , num2str(x(1))]); disp([椅子数量 x2 , num2str(x(2))]); disp([最大利润为, num2str(-fval)]); % 注意fval是最小化目标函数的值取负得到最大利润 disp([求解器退出状态 exitflag , num2str(exitflag)]); disp(output.message);运行这段代码你会得到类似下面的结果最优生产计划 桌子数量 x1 20 椅子数量 x2 10 最大利润为1300 求解器退出状态 exitflag 1 Optimization terminated.exitflag 1表示求解器成功找到了最优解。输出参数的含义如下x最优解向量即决策变量的值。这里x [20; 10]表示生产20张桌子和10把椅子利润最大。fval在最优解x处目标函数的值。因为我们求的是min -50x1-30x2所以fval -1300。原问题的最大利润就是-fval 1300。exitflag求解器的退出状态。这是非常重要的诊断信息。1函数收敛到最优解x。0迭代次数超过options.MaxIterations或函数计算次数超过options.MaxFunctionEvaluations。-2无可行解。这意味着你给出的约束条件互相矛盾找不到任何一个点能同时满足所有约束。比如你要求x1 x2 5同时又要求x1 10这显然不可能。-3问题无界。这意味着在满足约束的条件下目标函数值可以无限地减小对于最小化问题。比如你的约束只有x1 0目标函数是min -x1那么x1越大目标函数值就越小没有下限。其他负值代表求解过程中遇到了其他错误。output一个结构体包含关于优化过程的详细信息比如迭代次数、算法、收敛信息等。output.message通常包含了可读的退出原因。3.2 求解选项options设置默认情况下linprog会使用内点法求解。但对于不同规模、不同特性的问题你可能需要调整求解选项。这通过optimoptions函数来设置。% 创建优化选项 options optimoptions(linprog, ... Display, iter, ... % 显示每次迭代信息 Algorithm, dual-simplex, ... % 使用对偶单纯形法 MaxIterations, 1000, ... % 最大迭代次数 OptimalityTolerance, 1e-8); % 最优性容差 % 使用选项重新求解 [x, fval, exitflag, output] linprog(f, A, b, Aeq, beq, lb, ub, [], options);常用的选项包括Display: 控制输出信息的详细程度。off不显示final默认只显示最终结果iter显示每次迭代信息用于调试。Algorithm: 选择求解算法。对于线性规划主要有dual-simplex默认对偶单纯形法。通常对中等规模问题非常高效尤其在重新求解一个类似问题时。interior-point-legacy内点法旧版。对于大规模稀疏问题可能更快。interior-point内点法新版R2015b以后。MaxIterations/MaxFunctionEvaluations: 设置迭代或函数计算的上限防止程序陷入无限循环。OptimalityTolerance/ConstraintTolerance: 优化终止的容差。如果结果对精度要求极高可以适当调小这些值如1e-10但可能会增加计算时间。对于初学者大部分问题使用默认选项即可。当你遇到求解失败、速度慢或者对结果精度有疑问时才需要考虑调整这些选项。4. 从理论到实战三类典型建模案例拆解懂了基本操作我们来看几个更贴近实际数学建模赛题的例子。这些例子会涉及到如何将文字描述转化为数学模型以及处理一些常见的变体。4.1 案例一资源分配与成本最小化不等式约束为主问题描述一个饲料公司需要生产一种混合饲料要求每公斤饲料中营养成分A至少含18单位B至少含12单位C至少含16单位。现有三种原料可供选择其每公斤成本及营养成分含量如下表原料成本元/公斤营养A单位营养B单位营养C单位甲4321乙2213丙3132问如何配比这三种原料才能在满足营养要求的前提下使每公斤饲料的成本最低建模与求解决策变量设每公斤饲料中使用甲、乙、丙三种原料的量分别为x1, x2, x3公斤。目标函数最小化总成本min Z 4*x1 2*x2 3*x3。约束条件营养A要求3*x1 2*x2 1*x3 18营养B要求2*x1 1*x2 3*x3 12营养C要求1*x1 3*x2 2*x3 16非负约束x1, x2, x3 0隐含约束这里没有要求总重量必须等于1公斤因为x1, x2, x3本身就是重量目标是最小化单位重量的成本。如果要求“配制恰好1公斤饲料”则需要增加等式约束x1 x2 x3 1此时x1, x2, x3代表的是比例。本题没有这个要求所以不添加。Matlab标准化注意约束是“大于等于”而linprog标准形式是“小于等于”。我们需要两边乘以-1来转换。例如3x12x2x3 18等价于-3x1 -2x2 -x3 -18。% 饲料配比问题 - 成本最小化 f [4; 2; 3]; % 成本系数直接求最小 % 不等式约束 ( 转化为 ) A [-3, -2, -1; % 对应营养A: -3x1-2x2-x3 -18 -2, -1, -3; % 对应营养B: -2x1-x2-3x3 -12 -1, -3, -2]; % 对应营养C: -x1-3x2-2x3 -16 b [-18; -12; -16]; % 无等式约束 Aeq []; beq []; % 变量非负 lb [0; 0; 0]; ub []; [x_opt, fval_opt, exitflag] linprog(f, A, b, Aeq, beq, lb, ub); if exitflag 0 disp(最优原料配比公斤/每公斤饲料); disp([甲: , num2str(x_opt(1))]); disp([乙: , num2str(x_opt(2))]); disp([丙: , num2str(x_opt(3))]); disp([最低成本元/每公斤饲料: , num2str(fval_opt)]); else disp(求解失败请检查模型和约束。); end运行后你可能会得到类似x_opt [2; 0; 6]的结果表示每公斤饲料使用2公斤甲和6公斤丙总成本为4*23*626元。注意这里x1x2x38意味着为了满足高营养要求需要8公斤原料混合出“1公斤”符合要求的饲料这在逻辑上是“浓缩”饲料模型本身是合理的。如果要求总重为1则需要添加等式约束结果和意义会完全不同。4.2 案例二投资组合优化含等式约束与边界问题描述你有100万资金打算投资到三种资产股票A、债券B、基金C。已知期望年收益率股票A为10%债券B为5%基金C为8%。风险系数假设股票A为8债券B为2基金C为5。你希望总投资风险不超过500风险系数加权和并且为了分散风险规定对基金C的投资比例不低于20%。同时你与银行有一个协议必须将总资金的10%用于购买该银行发行的一款特殊债券D收益率3%这款债券必须单独考虑。问如何分配资金包括债券D在满足风险、比例和全额投资的要求下使期望年收益最大建模与求解决策变量设投资于股票A、债券B、基金C、特殊债券D的资金分别为x1, x2, x3, x4万元。目标函数最大化总期望收益max Z 0.10*x1 0.05*x2 0.08*x3 0.03*x4。约束条件总投资额x1 x2 x3 x4 100全部资金用完等式约束风险约束8*x1 2*x2 5*x3 500特殊债券D风险假设为0基金C比例约束x3 0.2 * 100即x3 20特殊债券D约束x4 0.1 * 100即x4 10这也是一个等式约束非负约束x1, x2, x3, x4 0Matlab标准化目标函数取负求最小。注意现在我们有两个等式约束总资金约束和债券D约束。% 投资组合问题 f -[0.10; 0.05; 0.08; 0.03]; % 收益系数求最大转为求最小 % 不等式约束风险约束和基金C下限 A [8, 2, 5, 0; % 风险约束系数 0, 0, -1, 0]; % x3 20 等价于 -x3 -20 b [500; -20]; % 等式约束总资金和债券D Aeq [1, 1, 1, 1; % x1x2x3x4 100 0, 0, 0, 1]; % x4 10 beq [100; 10]; % 变量非负 lb [0; 0; 0; 0]; ub []; % 无上限 [x_opt, fval_opt] linprog(f, A, b, Aeq, beq, lb, ub); disp(最优投资方案万元); disp([股票A: , num2str(x_opt(1))]); disp([债券B: , num2str(x_opt(2))]); disp([基金C: , num2str(x_opt(3))]); disp([特殊债券D: , num2str(x_opt(4))]); disp([最大年期望收益万元: , num2str(-fval_opt)]);这个例子展示了如何同时处理等式和不等式约束以及如何用矩阵Aeq和beq来表达多个等式关系。4.3 案例三数据拟合与简单机器学习线性回归作为线性规划线性规划不仅能解决资源分配问题还可以用于某些类型的数据拟合。考虑一个简单的中位数回归L1回归问题。与普通最小二乘法L2最小化误差平方和不同L1回归最小化误差的绝对值之和对异常值更不敏感。问题描述给定一组数据点(x_i, y_i)我们想用一条直线y a*x b来拟合。最小二乘法的目标是min Σ|y_i - (a*x_i b)|^2而L1回归的目标是min Σ|y_i - (a*x_i b)|。这个绝对值目标函数不是线性的但可以通过一个经典技巧转化为线性规划。转化技巧对于每个数据点i我们引入两个非负辅助变量u_i和v_i令y_i - (a*x_i b) u_i - v_i 且u_i 0, v_i 0。 那么绝对值|y_i - (a*x_i b)|就可以用u_i v_i来等价表示。因为对于任意实数它都可以表示成两个非负数的差而其绝对值就是这两个非负数的和。这样原问题就转化为决策变量a, b, u_1, v_1, u_2, v_2, ..., u_n, v_n目标函数min (u_1v_1 u_2v_2 ... u_nv_n)约束条件对于每个i有a*x_i b u_i - v_i y_i且u_i, v_i 0。a, b无约束可为负。% 生成带异常值的示例数据 rng(0); % 固定随机种子使结果可重复 x (1:20); y_true 2*x 5; y y_true randn(20,1)*3; % 加入正态噪声 y(10) y(10) 30; % 在第10个点加入一个异常值 % 使用线性规划进行L1回归 n length(x); % 决策变量顺序: [a; b; u1; v1; u2; v2; ...; un; vn] % 目标函数系数: a和b的系数为0u_i和v_i的系数为1 f [0; 0; ones(2*n, 1)]; % 等式约束: a*x_i b u_i - v_i y_i Aeq zeros(n, 2 2*n); for i 1:n Aeq(i, 1) x(i); % a 的系数 Aeq(i, 2) 1; % b 的系数 Aeq(i, 22*(i-1)1) 1; % u_i 的系数 Aeq(i, 22*(i-1)2) -1; % v_i 的系数 end beq y; % 不等式约束: 无 A []; b []; % 边界约束: a, b 无界u_i, v_i 0 lb [-inf; -inf; zeros(2*n, 1)]; ub []; % 求解线性规划 options optimoptions(linprog, Display, off); [z, fval] linprog(f, A, b, Aeq, beq, lb, ub, [], options); a_l1 z(1); b_l1 z(2); % 对比普通最小二乘法L2回归 X [x, ones(size(x))]; coeff_l2 X \ y; a_l2 coeff_l2(1); b_l2 coeff_l2(2); % 绘图对比 figure; scatter(x, y, bo, DisplayName, 数据点含异常值); hold on; plot(x, a_l1*x b_l1, r-, LineWidth, 2, DisplayName, L1回归线性规划); plot(x, a_l2*x b_l2, g--, LineWidth, 2, DisplayName, L2回归最小二乘); plot(x, y_true, k:, LineWidth, 1.5, DisplayName, 真实关系); legend(Location, best); xlabel(x); ylabel(y); title(L1回归 vs L2回归对异常值的鲁棒性); grid on; hold off; disp([L1回归结果: y , num2str(a_l1), * x , num2str(b_l1)]); disp([L2回归结果: y , num2str(a_l2), * x , num2str(b_l2)]);运行这段代码你会直观地看到由于异常值的存在绿色的L2回归线被明显拉偏而红色的L1回归线则更接近真实的黑色虚线。这个例子展示了线性规划在稳健统计中的应用也体现了将非线性的绝对值问题转化为线性规划的巧妙思路。在数学建模中这种转化能力非常重要。5. 求解失败问题无界或无解常见错误排查指南在实际使用中你可能会遇到linprog返回负的exitflag比如-2或-3这意味着求解失败。别慌这通常是你的模型出了问题而不是Matlab的错。下面是一些排查思路。5.1 问题无界exitflag -3表现求解器告诉你目标函数值可以无限小对于最小化问题。可能原因与排查缺少关键约束这是最常见的原因。检查你是否漏掉了对决策变量的必要限制。例如在最大化利润的生产问题中如果你忘记了原料或工时的约束产量就可以无限大利润也就无限大对应最小化问题的负利润无限小。约束方向错误检查你的不等式约束符号。如果你本意是消耗 资源但不小心写成了消耗 资源并且资源量很小这可能导致变量为了满足这个“宽松”的约束而变得非常大。变量无下界对于求最小化的问题如果目标函数系数f中存在负数且对应的变量没有下界lb为-inf则该变量趋向负无穷时目标函数值也会趋向负无穷。调试方法首先简化问题。尝试去掉一些约束或者固定一些变量的值看问题是否变得有解。这能帮你定位是哪个约束缺失或错误。其次检查模型假设。回到问题描述确保每一个现实中的限制条件都在模型中有对应的数学表达式。最后可视化对于二维或三维问题。你可以手动绘制约束条件围成的区域可行域看看它是否在一个方向上是不封闭的。例如只有x0和y0两个约束可行域就是第一象限对于目标函数min -x-y解就是无界的。5.2 无可行解exitflag -2表现求解器找不到任何一个点能同时满足所有约束条件。可能原因与排查约束条件互相矛盾这是根本原因。比如你要求x1 x2 5同时又要求x1 10且x2 0。这两个条件不可能同时成立。等式约束过于严格多个等式约束可能共同定义了一个空集。例如x1 x2 10和x1 x2 20显然矛盾。边界约束与不等式/等式约束矛盾lb和ub指定的范围可能与A*xb或Aeq*xbeq定义的范围没有交集。调试方法逐一注释法这是最实用的方法。暂时注释掉一部分约束比如先注释所有等式约束Aeq, beq运行看是否有解。然后逐步加入约束直到加入某一条约束后问题变得不可行那么这条约束很可能就是矛盾的源头。检查数据输入错误非常常见仔细核对A,b,Aeq,beq,lb,ub中的每一个数值。一个数字输错比如把100打成1000就可能导致可行域为空。检查单位一致性确保所有约束方程中的量纲是一致的。例如一边是“公斤”另一边不能是“吨”。使用linprog的‘diagnostics’选项设置options optimoptions(linprog, Display, iter, Diagnostics, on)更详细的输出有时能提供线索。5.3 其他常见错误与技巧结果不直观或不符合预期即使exitflag1解也可能很奇怪。比如投资比例加起来超过了1或者产量是小数但现实中必须是整数。这时你需要检查模型是否完整是否漏掉了“所有投资比例之和为1”这样的等式约束考虑整数规划如果变量必须取整如生产设备台数线性规划允许非整数解的特性会导致结果失真。这时应使用混合整数线性规划函数是intlinprog。检查目标函数系数确认系数的正负号是否正确最大化 vs 最小化。大规模问题的性能对于变量和约束成千上万的大规模问题默认的内点法可能内存消耗较大。可以尝试切换算法为‘dual-simplex’它对某些稀疏结构问题更有效。同时确保你的A和Aeq矩阵以稀疏矩阵格式存储使用sparse函数可以极大提升求解速度和降低内存使用。灵敏度分析与影子价格linprog的输出不直接提供影子价格对偶变量。但你可以通过求解对偶问题或者使用linprog的完整输出格式[x, fval, exitflag, output, lambda] linprog(...)来获取lambda结构体。lambda.ineqlin和lambda.eqlin分别给出了不等式和等式约束的影子价格这在经济解释和“如果资源增加一单位利润能增加多少”的边际分析中非常有用。6. 进阶从线性规划到整数规划与实战建议线性规划假设变量可以取任意实数但现实中的很多问题要求变量是整数比如“购买几台机器”、“派遣几个员工”、“一个项目是否启动0或1”。这时就需要用到整数规划或混合整数线性规划。6.1 使用intlinprog处理整数约束Matlab中求解混合整数线性规划的函数是intlinprog它的语法和linprog非常相似多了一个指定哪些变量需要取整数的参数intcon。语法[x, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub)参数intcon: 一个整数向量指定哪些决策变量必须取整数值。例如如果x(1)和x(3)是整数变量则intcon [1, 3]。例子背包问题假设有一个容量为10的背包有4件物品其价值和体积如下表。每件物品只能选择带或不带0-1变量如何选择物品使得总价值最大物品价值体积152263341474% 0-1背包问题 f -[5; 6; 4; 7]; % 最大化价值转为最小化负价值 intcon 1:4; % 所有变量都是0-1整数 % 体积约束 A [2, 3, 1, 4]; b 10; % 无等式约束 Aeq []; beq []; % 0-1变量边界 lb zeros(4,1); ub ones(4,1); [x_opt, fval_opt] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub); disp(最优选择1表示携带0表示不携带); disp(x_opt); disp([最大总价值, num2str(-fval_opt)]);intlinprog的求解比linprog复杂得多计算时间随问题规模指数级增长。对于大规模整数规划问题可能需要设置更长的求解时间限制options.MaxTime或调整分支定界法的相关选项。6.2 数学建模竞赛中的实战建议结合我参加和辅导数学建模竞赛的经验在比赛中使用线性/整数规划时有几点特别重要模型第一编程第二花足够的时间把问题理清建立正确的数学模型。一个错误的模型即使用最高级的算法也得不到正确的答案。在写代码前最好先在纸上把决策变量、目标函数、所有约束清晰地列出来。从简单开始逐步复杂化不要试图一口气建立一个包含所有细节的完美模型。先建立一个最简化的核心模型比如忽略一些次要约束用linprog快速求解验证基本逻辑是否正确。然后逐步添加更复杂的约束如整数约束、非线性约束的线性近似、多阶段决策等。重视结果的分析与解释数学建模不是解出答案就结束了。你需要分析结果敏感性分析如果某个参数如资源量、价格变化了10%最优解会变化多少这可以通过改变b或f中的值重新求解来实现。影子价格的经济意义对于资源约束其影子价格代表了该资源的边际价值。在论文中解释这个值能为你的方案提供深刻的洞察。结果的合理性得到的解是否在物理意义或常识上是合理的如果产量是负数或者投资比例加起来不等于1一定要回头检查模型。注意性能对于变量较多比如成千上万的问题在构建A和Aeq矩阵时尽量使用稀疏矩阵sparse(i, j, v, m, n)来存储可以节省大量内存和计算时间。做好错误处理在你的脚本中一定要检查exitflag。如果求解失败可以尝试输出中间变量或者给出一个备用的启发式方案并在论文中说明求解的局限性。一个健壮的代码能让你在紧张的比赛时间里更从容。最后线性规划是数学建模中极其强大的工具但它只是一个工具。真正考验你的是如何将一个复杂的现实问题抽象成一个简洁而准确的数学模型。Matlab的优化工具箱为你解决了计算上的难题让你能更专注于建模本身。多练习多思考下次当你遇到“在有限条件下寻求最优”的问题时不妨先想想“这能不能用一个线性规划来描述”