做数学建模题的时候,只要模型里冒出“绝对偏差之和最小”“到各个需求点的距离之和最短”“总调整量绝对值最小”这类字眼,很多人第一反应就是打开 MATLAB 写个abs()丢进linprog,然后被一句报错打回原形——linprog只认线性表达式,abs在它眼里就是非线性。于是要么改用手写遗传算法碰运气,要么直接换成fmincon求个局部解,最后结果忽好忽坏,论文里也不敢写“全局最优”。
我这些年帮人看过的建模代码里,含绝对值的线性规划(L1 范数最小化、绝对值约束)几乎是出错率最高的一类。问题不在于 MATLAB 不会算,而在于建模那一层没把绝对值“拆干净”,矩阵拼装时行数列数对不上、符号写反、辅助变量忘了给下界,这些坑一个比一个隐蔽,而且报错信息往往只告诉你“维度不一致”,根本指不到病灶。下面我就按“为什么要拆——怎么拆——拆完怎么喂进去——跑通之后还有哪些坑”这条线,把含绝对值线性规划的 MATLAB 求解方法完整讲一遍,中间给两个可以直接复现的完整案例:一个是带离群点的直线拟合,一个是曼哈顿距离意义下的选址问题。基础一般的读者能照着抄作业,有经验的读者可以重点看关于解不唯一和数值稳定性的那两节。
1. 绝对值让模型“写不下去”的根源到底在哪
1.1 linprog 接受的输入长什么样
先把linprog的标准形式摆出来,这是后面所有推导的基准:
min f' * x s.t. A * x <= b Aeq* x == beq lb <= x <= ub它接受的只有四种东西:一个线性目标、一组线性不等式、一组线性等式、一组变量的上下界。判断标准非常机械——目标函数和约束里,每个决策变量只能以“常数乘一次”的形式出现,不能有平方、不能有乘积、不能有两个变量相乘,当然也不能有绝对值。
绝对值的麻烦之处在于它是分段线性的:|x|在x >= 0时等于x,在x < 0时等于-x。它不是一条直线,而是一条折线。你没法用一个线性函数在整条实数轴上表示它——这是数学上的硬约束,不是 MATLAB 的能力问题。所以正确的思路不是“让 MATLAB 认识 abs”,而是“把这条折线用额外的变量和约束铺平成平面上的一个凸区域”。
这里有个很重要的前提需要说清楚:绝对值函数是凸函数。这个性质决定了两件事。第一,如果绝对值出现在最小化目标里,整个问题仍然是凸的,可以用线性规划精确求解,改写是等价改写,不是近似。第二,如果绝对值出现在需要“大于等于”的约束里(比如|x| >= 3),可行域会变成两条射线,是非凸的,纯线性规划就搞不定了,只能上整数变量。很多人一开始不区分这两种情况,写出来的模型看着像线性规划,其实可行域已经碎成两块,linprog给出的解要么是其中一块的顶点,要么直接告诉你不可行。
1.2 目标里带绝对值和约束里带绝对值,处理手法完全不同
我把常见的位置分成三类,处理难度递增:
| 出现位置 | 典型场景 | 是否保持凸性 | 处理手法 |
|---|---|---|---|
| 最小化目标中的 `sum | r_i | ` | L1 回归、最小绝对偏差、总距离最小 |
| 约束中的 ` | a'x + b | <= c` | 偏差限幅、容差带、区间限制 |
| 约束中的 ` | a'x + b | >= c` | 最小偏离量、必须“跳出”某个区间 |
第一类是最常见的,做法是“上境图法”(epigraph form):把|r_i|用一个新变量t_i顶替,然后在约束里逼着t_i不小于r_i和-r_i。因为目标是最小化sum t_i,最优化过程会自动把t_i压到刚好等于|r_i|的位置,所以最终sum t_i就等于sum |r_i|,等价性成立。
第二类最省事,|a'x + b| <= c直接等价于-c <= a'x + b <= c,拆成两条不等式塞进A*x <= b里就行,一个辅助变量都不用加。很多人习惯性地一看到绝对值就引入u, v分解,结果变量数平白多了一倍,规模大一点的问题立刻就慢了。
第三类必须用整数变量,|a'x + b| >= c等价于a'x + b >= c或a'x + b <= -c,两个区域中间断开,需要一个 0-1 变量z做选择器:a'x + b >= c - M*z和a'x + b <= -c + M*(1-z)。这时候求解器要从linprog换成intlinprog。建模比赛里这类约束出现频率不高,但一旦出现,如果还硬用linprog,得到的解几乎肯定不满足原约束。
2. 目标函数里的绝对值:辅助变量怎么加才不漏不重
2.1|x| = u - v分解法的适用边界
网上流传最广的做法是x = u - v,u, v >= 0,然后|x|写成u + v。这个写法数学上没错,因为对任意实数x,总能找到非负的u, v使u - v = x且u + v >= |x|(等号在u和v至少有一个为 0 时取到)。
但它有个非常容易踩的陷阱:只有当你确实在最小化u + v的系数为正时,这个分解才等价。如果u + v只是被塞进某条约束里,目标函数里根本没有它,那么求解器完全可以让u和v同时加上一个大常数,u - v的差不变,但u + v被撑大了。这时候你拿到的解在数学上满足约束,但对应的绝对值根本不是你想要的。
我见过一个队友在约束里写u + v <= 10想表示|x| <= 10,结果模型给出的x值是 3,而u - v算出来是 3 的同时u + v被撑到了 9.99——看起来没问题,实际上是碰巧。如果那条约束右边是 100,解就完全不靠谱了。
结论很清楚:u - v分解法适合“决策变量本身可正可负、且绝对值要进入最小化目标”的场合;而如果绝对值是出现在约束的“上界”一侧,直接用双向不等式-c <= x <= c,别绕这个弯子。
2.2 上境图法的完整推导:为什么t >= ±r就够了
设我们要最小化sum |r_i|,其中r_i是决策变量的线性组合,比如r_i = y_i - (a*x_i + b)。引入t_i,构造:
min sum t_i s.t. t_i >= r_i t_i >= -r_i t_i >= 0 (其实由上面两条可推出,但显式加上更稳)为什么这两条约束加上最小化目标就能保证t_i = |r_i|?反过来想:假设最优解里某个t_k > |r_k|,那么我把t_k下调到|r_k|,约束t_k >= r_k和t_k >= -r_k依然满足(因为|r_k|就是这两个下界的最大值),而目标函数值严格变小了。这与“最优解”矛盾。所以最优解里必然t_k = |r_k|,等价性得证。
这个证明虽然短,但它揭示了一个关键条件:目标函数里t_i的系数必须是正数。如果你写的是min sum t_i - 5*t_1,那t_1就会被无限制地推大,模型直接无界。我在审别人代码时,看到目标向量f里出现负系数配辅助变量,基本可以直接判定这模型有问题。
顺带说一句,t_i >= 0这条其实可以从t_i >= r_i和t_i >= -r_i推出(因为max(r, -r) = |r| >= 0),但显式写进lb里是零成本的好习惯——它能帮求解器更快地做预处理,也能在你拼错了符号时让问题表现为“不可行”而不是“无界”,报错信息更容易读懂。
2.3 一次完整拼装:从数学式子到linprog入参
拿最典型的 L1 直线拟合举例。给n个数据点(x_i, y_i),要拟合y = a*x + b,最小化残差绝对值之和。
决策变量排成列向量z = [a; b; t_1; t_2; ...; t_n],一共n + 2个。目标向量f = [0; 0; 1; 1; ...; 1]。
约束t_i >= r_i = y_i - a*x_i - b,移项得a*x_i + b + t_i >= y_i,转成<=形式:-x_i*a - b - t_i <= -y_i。
约束t_i >= -r_i = -(y_i - a*x_i - b) = a*x_i + b - y_i,移项得t_i - a*x_i - b >= -y_i,转成<=:x_i*a + b - t_i <= y_i。
所以A是2n行、n + 2列。第i行(对应第一条)在a位置是-x_i,在b位置是-1,在t_i位置是-1,其余为 0;第n+i行在a位置是x_i,在b位置是1,在t_i位置是-1,其余为 0。右端项b_vec前n个是-y_i,后n个是y_i。
下界lb = [-inf; -inf; zeros(n,1)],a和b是自由变量,必须给-inf,不能给 0——给了 0 就相当于强行要求截距非负,很多拟合结果会因此偏得离谱。
对应的 MATLAB 代码:
function [a, b, res] = l1fit(x, y) x = x(:); y = y(:); n = numel(x); f = [0; 0; ones(n,1)]; % 上半部分: -x_i*a - b - t_i <= -y_i I1 = [ -x, -ones(n,1), -eye(n) ]; % 下半部分: x_i*a + b - t_i <= y_i I2 = [ x, ones(n,1), -eye(n) ]; A = [I1; I2]; bvec = [-y; y]; lb = [-inf; -inf; zeros(n,1)]; ub = []; opts = optimoptions('linprog', 'Algorithm', 'dual-simplex', ... 'Display', 'off'); [z, fval, exitflag] = linprog(f, A, bvec, [], [], lb, ub, opts); if exitflag ~= 1 warning('linprog exitflag = %d,结果需要复核', exitflag); end a = z(1); b = z(2); res = fval; end这段代码里的eye(n)构造对角块是个小技巧,-eye(n)直接把-t_i的对角结构铺出来,比循环填矩阵快得多,也避免了索引写错。数据量上到几千个点时,这种向量化写法能把拼装时间从秒级压到毫秒级。
3. 约束里的绝对值:展开方式决定问题规模
3.1|a'x + b| <= c直接双向展开,别加变量
|a'x + b| <= c与-c <= a'x + b <= c完全等价,这是绝对值最基本的定义推论。拆开就是两条线性不等式:
a'x + b <= c -a'x - b <= c两条都塞进A*x <= b里就完事,不需要任何辅助变量。这个做法在建模中特别有用,例如:
- 生产计划中要求某产品的产量与基准值的偏差不超过 5%,直接写成
产量 - 基准 <= 0.05*基准和基准 - 产量 <= 0.05*基准。 - 投资组合中要求某资产权重偏离目标权重的幅度受限,同款展开。
- 曲线拟合中要求拟合值落在观测值的某个容差带内,每个观测点贡献两条约束。
这里唯一的坑是右端常数c必须是非负的。如果c < 0,那么|某东西| <= 负数永远不成立,问题必然不可行。我遇到过一次,队员把c写成了一个含变量的表达式,展开后忘了检查符号,模型跑出来exitflag = -2,查了半天才发现是容差写成了负号。
3.2|a'x + b| >= c为什么不能照搬
|a'x + b| >= c(c > 0)展开是a'x + b >= c或a'x + b <= -c,注意是“或”不是“且”。线性规划的约束天然是“且”的关系,所有约束同时成立。你没法在纯线性规划里表达“或”。
硬要写成“且”,比如同时加上a'x + b >= c和a'x + b <= -c,那这两条本身互相矛盾,可行域直接空掉,linprog返回exitflag = -2。
正确的做法是引入 0-1 变量z ∈ {0, 1}和大 M 常数:
a'x + b >= c - M*z a'x + b <= -c + M*(1 - z)z = 0时第二条退化成a'x + b <= -c + M(M 足够大时无约束力),第一条生效;z = 1时反过来。两条约束里恰好有一条在起作用。M 的取值要覆盖a'x + b可能达到的最大绝对值,取太小会误砍掉可行解,取太大会让松弛变差、数值不稳定。经验值是把|a'x + b|在上界和下界处的取值都算一遍,取两者绝对值的较大者再乘 1.5 左右。
求解器要换成intlinprog,调用格式和linprog几乎一样,只是多一个intcon参数指定哪些变量是整数:
intcon = n + 1; % 假设 z 排在最后 [z_sol, fval] = intlinprog(f, intcon, A, b, Aeq, beq, lb, ub);intlinprog内部走的是分支定界,规模稍大(比如 0-1 变量超过几百个)就会明显变慢。建模比赛里如果时间紧,要提前评估:这类约束是不是真的必要,能不能用其他建模方式绕开。
3.3 等式约束带绝对值:几乎都是建模错误
偶尔会看到有人写|a'x + b| = c这样的约束。这个等式意味着a'x + b只能取c或-c两个值,可行域是两个超平面的并集,同样非凸。处理方式和上面一样,也要 0-1 变量。
但我更想说的是:九成情况下出现这种约束,是因为建模时想表达的东西没说清楚。比如“产量必须正好等于目标值”写成绝对值等于零就行了(a'x + b = 0),“偏差必须恰好是某个数”这种需求本身就很少见。碰到这种约束,先回头看看原始问题描述,是不是把“不超过”或者“不少于”误写成了“等于”。
4. 两个能直接跑通的实战案例
4.1 案例一:带离群点的直线拟合
构造一组数据:真实关系是y = 2.5*x + 1,在前 20 个点上加小幅噪声,然后在最后 3 个点上故意加一个很大的偏移(模拟记录错误或者异常样本)。这种数据下,普通最小二乘会被离群点拽偏,而 L1 拟合对离群点不敏感,这正是绝对值建模的价值所在。
rng(7); x = linspace(0, 10, 23)'; y = 2.5*x + 1 + 0.3*randn(size(x)); y(end-2:end) = y(end-2:end) + 12; % 三个离群点 % L1 拟合 [a1, b1, res1] = l1fit(x, y); % 对照:最小二乘 p = polyfit(x, y, 1); fprintf('L1 : a = %.4f, b = %.4f, 绝对残差和 = %.4f\n', a1, b1, res1); fprintf('LS : a = %.4f, b = %.4f\n', p(1), p(2));在类似的数据设定下,最小二乘的斜率会被离群点明显拉低,截距被抬高;而 L1 拟合的a和b会非常接近真实值2.5和1。这个对比本身就是建模论文里很有说服力的一段。
这里有个容易被忽略的细节:res1是目标函数的最优值,也就是sum |r_i|,它和sum(r_i.^2)不是一个量纲上的东西,写论文时别把两者放在一张表里比大小,没有意义。要比较就比“估计参数与真实参数的偏差”,或者用留出集上的绝对误差。
4.2 案例二:曼哈顿距离意义下的单设施选址
某园区要新建一个配送站,服务m个需求点,坐标(p_i, q_i),需求量w_i。园区道路是网格状的,车辆只能沿横纵方向走,所以实际行驶距离是曼哈顿距离|x - p_i| + |y - q_i|。目标是让加权总距离最小。
模型:
min sum w_i * (u_i + v_i) s.t. u_i >= x - p_i u_i >= -(x - p_i) v_i >= y - q_i v_i >= -(y - q_i) u_i, v_i >= 0其中u_i承担|x - p_i|,v_i承担|y - q_i|。注意这里x和y的绝对值是分开处理的,因为曼哈顿距离本身就是两个一维绝对值之和。这一点很关键:如果误用欧氏距离的写法,模型就变成二阶锥规划,linprog直接失效。
决策变量z = [x; y; u_1..u_m; v_1..v_m],共2m + 2个。约束转成<=形式:
x - u_i <= p_i (来自 u_i >= x - p_i) -x - u_i <= -p_i (来自 u_i >= -(x - p_i)) y - v_i <= q_i -y - v_i <= -q_iMATLAB 实现:
function [xs, ys] = manhattan_loc(p, q, w) p = p(:); q = q(:); w = w(:); m = numel(p); f = [0; 0; w; w]; % 目标: sum w_i*(u_i + v_i) % 构造 A, 每 2 个需求点贡献 4 行 rowX1 = [ 1, 0, -eye(m), zeros(m)]; rowX2 = [-1, 0, -eye(m), zeros(m)]; rowY1 = [ 0, 1, zeros(m), -eye(m)]; rowY2 = [ 0,-1, zeros(m), -eye(m)]; A = [rowX1; rowX2; rowY1; rowY2]; bvec = [p; -p; q; -q]; lb = [-inf; -inf; zeros(2*m,1)]; opts = optimoptions('linprog','Algorithm','dual-simplex','Display','off'); [z, fval, exitflag] = linprog(f, A, bvec, [], [], lb, [], opts); xs = z(1); ys = z(2); end拿一组数据试一下:p = [2; 8; 5; 9; 3],q = [7; 3; 9; 1; 6],w = [1; 2; 1; 3; 2]。跑出来的最优点会落在某个需求点附近,或者落在两个需求点之间的一段区域上——曼哈顿距离下的中位数性质决定了最优解可能是一个区域而不是一个点,这一点后面第 5 节会详细说。
4.3 用第二种方法做交叉验证,别只信一次求解
上面两个案例跑完之后,一定要做交叉验证。最简单的办法是用fmincon直接对原始的非线性模型(带abs)求一次,看两者结果是否一致。fmincon求的是局部解,但在这种凸问题上,局部最优就是全局最优,两者应该对得上。
% 案例二的对照求解 obj = @(z) sum(w .* (abs(z(1) - p) + abs(z(2) - q))); z0 = [mean(p), mean(q)]; [z_nl, fval_nl] = fmincon(obj, z0, [], [], [], [], [-inf;-inf], [], []);如果fval和fval_nl相差超过 1e-6,那一定有一边写错了。这个对比步骤看起来麻烦,但在比赛里能救命的——线性化矩阵拼错是最常见也最难肉眼发现的错误,靠单元测试比对是最快的定位方式。
5. 跑通之后才开始的坑:规模、数值与解不唯一
5.1Algorithm选dual-simplex还是interior-point
linprog提供三个算法选项:'dual-simplex'、'interior-point'、'interior-point-legacy'(新版 MATLAB 里还多了'highs'系列)。实际用下来:
dual-simplex适合中小规模(变量数几百到几千),收敛稳、对退化问题处理得好,缺点是内存占用随约束数增长较快。绝对值线性化后约束数量往往翻倍甚至三倍,规模上去以后这个算法会变慢。interior-point适合大规模稀疏问题,内存友好、速度快,但它求出的解不一定是顶点解,对“需要整数解或者需要明确顶点”的场合不友好。另外内点法在问题接近退化时容易报exitflag = -3或者给出“可行但目标值略差”的结果。highs系列(新版 MATLAB 默认)综合了前两者的优点,一般直接用它就行,除非你在老版本 MATLAB 上跑。
我个人的习惯是:建模比赛先用默认选项跑通,再用dual-simplex复核一遍目标值。两者差得多就说明模型有数值问题,通常是因为系数矩阵里存在数量级差异极大的行(比如一行的系数是1e-6量级,另一行是1e5量级)。
5.2 变量翻倍之后,稀疏矩阵能救回多少性能
含绝对值的线性化会把变量数和约束数都放大两到三倍。案例一里n = 23就只有 25 个变量,问题不大;但如果换成图像处理里n = 250000的像素级问题,变量数直接上到 25 万,A矩阵如果是全稠密存储,2n × (n+2)就是 50 万乘 25 万,内存直接爆掉。
解决办法是用sparse构造:
I1 = sparse([1:n]', [ones(n,1), 2*ones(n,1), (3:n+2)'], ... [-x, -ones(n,1), -ones(n,1)], 2*n, n+2); I2 = sparse([n+1:2*n]', [ones(n,1), 2*ones(n,1), (3:n+2)'], ... [ x, ones(n,1), -ones(n,1)], 2*n, n+2); A = [I1; I2];稀疏矩阵在linprog里是原生支持的,A和Aeq传sparse类型以后,求解器会自动走稀疏路径。实测下来,对于 10 万个变量级别的 L1 拟合,稀疏存储比稠密存储快一到两个数量级,内存占用也能从几十 GB 降到几百 MB。
还有一个减少变量的小技巧:如果能把问题拆成若干个独立子问题(比如按时间分段、按区域分块),拆开来分别求解,规模能降好几个数量级,比硬上一个巨型问题划算得多。
5.3exitflag不是 1 的时候,按这个顺序排查
linprog的exitflag是排查问题的第一手信息:
| exitflag | 含义 | 首先该检查什么 |
|---|---|---|
| 1 | 找到最优解 | 用fmincon交叉验证目标值 |
| 0 | 迭代次数上限,未收敛 | MaxIterations、OptimalityTolerance是否过严 |
| -1 | 被输出函数终止 | 一般不会出现,检查是否误设OutputFcn |
| -2 | 问题不可行 | 约束符号、右端项符号、lb/ub是否冲突 |
| -3 | 问题无界 | 目标向量f是否含负系数、辅助变量lb是否漏给 |
-2和-3是最常见的两个。-2的原因九成是符号写反:在转A*x <= b形式时,把>=直接照抄成<=是最经典的错误,记住“乘 -1 再变号”这条铁律。另一个常见原因是lb和ub交叉,比如某个变量lb = 5而ub = 3,这种问题必然不可行,但如果你是从别的代码改过来的,很可能残留了旧的边界值。
-3的原因基本都是目标向量的符号问题。回到第 2.2 节那个结论:最小化问题里,如果辅助变量的系数是负的,求解器会把它推到无穷大,然后告诉你无界。检查一遍f里辅助变量对应的位置是不是全为正。
5.4 解不唯一:为什么你的结果和队友的对不上
这个问题在曼哈顿选址里特别典型。曼哈顿距离意义下的中位数,最优解可能是一整段区域——只要x落在需求点横坐标的加权中位数区间内,目标值都一样。这意味着linprog给出的x只是这个区间里的某一个点,换一台机器、换一个算法版本、甚至只是改一下Algorithm选项,得到的x都可能不一样。
这不是程序错了,是问题本身的最优解不唯一。处理办法有几种:
一种是接受它,在论文里明确写“最优解位于[x_lo, x_hi] × [y_lo, y_hi]区域内,模型给出的是其中一个代表点”,然后补充说明该区域的计算方法——把问题再解两次,一次加约束x <= x_opt - ε求左边界,一次加x >= x_opt + ε求右边界。
另一种是加一个次级目标打破平局。比如在两个目标值相同的最优解中,选一个“离现有设施最近”或者“坐标总和最小”的,把次级目标加一个小系数到主目标上:min sum w_i*(u_i + v_i) + 1e-6*(x + y)。系数要足够小,不能改变主目标的最优值,但足以在平局解里做出确定性选择。这个技巧我在实际项目里用得很多,好处是结果可复现,坏处是论文里得解释清楚这个小系数是干什么的,别让评阅的人以为是拍脑袋加的。
6. 建模论文里怎么写才不会被追问
6.1 线性化的推导必须写出来,不能只写结果
比赛论文里写“将含绝对值的模型转化为线性规划模型”然后直接给最终矩阵,是很常见的写法,但也最容易在答辩环节被追问。评委真正想看的是:你怎么证明改写是等价的。
我的习惯是三步走:第一,写出辅助变量法的一般形式min sum t_i, s.t. t_i >= ±r_i;第二,用一句话说明等价性成立的理由(即第 2.2 节那个“t 调小则目标变小”的反证);第三,给出矩阵的维度说明和每个块的物理含义(哪几行对应上界、哪几行对应下界、t_i放在第几列)。
这三步加起来不超过半页纸,但能省掉后面一堆麻烦。尤其是“每个块的物理含义”这一条,在变量数上百的时候,你自己回头看代码全靠这个说明。
6.2 敏感性分析和模型检验不能省
含绝对值的模型有一个很有意思的敏感性特征:L1 拟合的解对数据点的“抖动”不敏感(某点动一点点,只要不跨过中位数位置,解不变),但对离群点的位置很敏感(离群点移出集合,解会跳变)。这个性质和最小二乘完全相反,是可以在论文里展开写的一段。
具体做法:把每个数据点逐个剔除,重新求解,记录参数变化幅度。变化幅度大的那些点就是“高影响点”,论文里列个表,说明模型对这些点的依赖性。这个过程写个循环几十行代码就能跑完,说服力却很强。
另外,曼哈顿选址这类问题,可以拿最优解和“重心法”(欧氏距离下的近似解)做对比,说明两种距离度量下选址结果的差异,并解释为什么网格道路条件下曼哈顿距离更合适。这种对比分析往往比多算几个算例更能体现建模思维。
6.3 代码交付时的几个小习惯
最后分享几个我在交付建模代码时坚持的习惯,都是踩坑踩出来的。
第一个,所有linprog调用都包一层exitflag检查,~= 1时直接warning或者error,绝不静默返回。比赛现场赶时间,很容易把一个不可行的解当成最优解写进论文,后面发现时已经来不及重来。
第二个,把“建模函数”和“求解调用”分开写。l1fit只负责拼矩阵和求解,数据生成、结果可视化放在脚本里。这样换组数据时只改脚本,逻辑不变,出错的概率小很多。
第三个,凡是涉及大 M 的地方,把 M 的值打出来看看。M = 1e9这种写法在数值上很危险,会让系数矩阵的条件数变得极大,求解器可能给出“看起来收敛但误差很大”的结果。能用1e4就用1e4,实在不行先估一下变量的实际取值范围再定。
第四个,辅助变量的物理含义在注释里写清楚。t(1:n)是残差绝对值、u(1:m)是横向距离,这种注释在半个月后回来看代码时价值连城。我自己有过一次经历,一个含绝对值的模型放了两个月再改,因为没写注释,对着A矩阵的列索引数了四十分钟才想起来第 37 列是哪个变量。
说到底,含绝对值的线性规划本身不难,难的是从“看到绝对值”到“条件反射地判断它该用双向展开还是辅助变量还是整数变量”这个判断过程。这个判断做对了,剩下的都是矩阵拼装的机械活;判断做错了,代码写得再漂亮也是在解一个错的问题。我个人的经验是,每次动手写代码之前,先在纸上把“哪些是决策变量、哪些是辅助变量、每个辅助变量被哪两条约束夹住”列一遍,这五分钟的纸上功夫能省掉后面几个小时的调试。