遗传算法求解带约束非线性网络流问题的Matlab实现

去年做物流调度项目时我接手了一个典型的带约束网络流问题站点几十个线路几百条每条线路有容量上限运输成本还是非线性的总预算又有硬约束。最初我用最小费用流和线性规划跑得很吃力线性化误差又没法接受。后来决定改用遗传算法GA来做约束优化效果比预想好很多。这篇就把整个解法摊开讲从网络流建模、路径流量编码到GA算子设计和Matlab代码实现适合正在被非线性约束折磨的工程师也适合刚接触智能优化的小白。1. 约束化网络流问题先搞清楚要解什么1.1 一个能复现的六节点算例先别急着聊算法我们要有一个可以上手的网络。这里用一个6节点、9条弧的小网络节点1是源节点6是汇节点2到5是中间转发节点。每条弧的容量、成本参数固定候选路径也预先给出。这样的规模小却保留了容量约束、预算约束和非线性成本三个核心难点足够把GA的套路讲清楚。弧的成本函数我写成二次形式cost_e(f) a_e * f^2 b_e * f其中f是这条弧上的流量。二次项模拟线路拥堵成本流量越大单位成本越高。如果你只是做工程测试也可以用线性成本但那样很多场景用线性规划就够了体现不出GA的价值。所以我特意用二次项让问题变成非凸非线性优化。下面是网络参数表我在后面的Matlab代码里也会原样使用。弧编号起点终点容量ab112200.0101.0213180.0121.1324220.0151.2425150.0181.5534180.0121.1635200.0151.4745100.0201.3846300.0101.0956350.0121.2现在假设我们有6条候选路径从节点1到节点6路径编号完整路径理论上界p11-2-4-6min(20,22,30)20p21-2-5-6min(20,15,35)15p31-3-4-6min(18,18,30)18p41-3-5-6min(18,20,35)18p51-2-4-5-6min(20,22,10,35)10p61-3-4-5-6min(18,18,10,35)10这里的理论上界是单条路径满载时的最大流量它不能超过路径上最窄那条弧的容量。这个上界在后面初始化种群的时候非常有用。1.2 把约束拆成三类再分别处理网络流问题到了实际项目里约束往往不止一个。拿这个算例来说我们面对的约束至少分三类。第一类是节点流量守恒也就是中间节点进多少、出多少源节点的净流量加上汇节点的净流量必须匹配。这是等式约束也是传统网络流模型最基础的条件。第二类是弧容量约束每条弧上的总流量不能超过容量上限。第三类是全局预算约束所有弧的运输成本加起来不能超过总预算。后面两个都是不等式约束而且预算约束是全局的和每一条路径都有关系。这三类约束的处理难度完全不同。等式约束最麻烦因为随机生成的流量通常无法满足守恒必须做专门的修复操作。容量约束和预算约束稍微好一点只要用罚函数把超限程度塞进适应度里就可以引导GA往可行域方向搜索。所以我的建模策略很明确能用编码消掉的约束就消掉消不掉的再用罚函数。1.3 为什么这种场景更适合GA教科书里的网络流问题大多是线性成本可以用最小费用流、最大流等经典算法也可以用线性规划快速求解。但一旦成本变成非线性函数再加一个总预算约束问题就变成一个非凸非线性约束优化问题精确算法的收敛速度和实现难度都会明显上升。GA在这种场景下有几个天然优势。它不要求目标函数可导不要求约束是线性的也不要求可行域是凸的。你只需要提供一组染色体编码和一个适应度函数剩下的交给种群进化。而且后面想加新的约束比如区域内最大风险、客户服务等级限制只需要在适应度函数里增加一个惩罚项不需要从零重写算法。这种灵活性在实际工程项目里比算法在纸面上多漂亮都更有价值。2. 数学建模GA怎么和网络流问题对话2.1 关键一步用路径流量替代弧流量编码如果你熟悉网络流第一反应可能是把每条弧的流量当成决策变量。这个编码方式在数学规划里没问题但在GA里非常糟糕。原因很简单弧流量之间必须满足流量守恒而这是一个等式约束。GA随机生成的种群几乎每个个体都不满足流量守恒你必须频繁调用修复算子去调整。修复逻辑在网络规模一大的时候复杂度会高到让人崩溃。我的做法是换一个角度不直接编码弧流量而是编码候选路径上的流量。设决策变量为x_k表示第k条路径上承载的从源到汇的流量。因为每条路径本身是一条从源到汇的有向通路只要流量是正值它在任何一个中间节点都会自动做到进等于出。所以无论x_k是多少节点流量守恒天然满足。可以理解成快递干线网一辆车从分拨中心出发沿途经过几个中转站最后到目的中心。这辆车在中转站卸下又装上的是同一批货量不会凭空多出货物来。只要每辆车运输总量确定了所有中转站的进出量自然平衡。这样做的好处非常明显等式约束直接被编码消掉了。剩下的弧容量和预算约束都是不等式在适应度函数里用罚函数处理即可。GA搜索空间里不存在大量因流量不守恒而浪费的个体收敛速度会明显优于弧流量编码。2.2 目标函数与罚函数的具体形式我们的优化目标是在满足容量和预算的前提下最大化从源到汇的总流量F_total sum(x_k)每条弧上的实际流量为f_e sum(pathMat(e,k) * x_k)其中pathMat(e,k)表示弧e是否在第k条路径上。总成本为C_total sum_e ( a_e * f_e^2 b_e * f_e )容量约束和预算约束分别写成罚函数penalty_capacity sum( max(0, f_e - u_e)^2 )penalty_budget max(0, C_total - B)^2最终适应度函数为fitness F_total - w1 * penalty_capacity - w2 * penalty_budget这里有个细节为什么用平方而不是一次项因为平方会让越界越多的个体被惩罚得更狠产生明显的梯度感。GA虽然不依赖梯度但这种非线性惩罚在锦标赛选择里非常有效轻微越界的个体保留更多多样性严重越界的个体很快被淘汰。罚函数权重w1和w2并不要一开始就调到极致。太大的权重会让种群过早失去多样性全部挤在可行域边界太小的权重又会让算法忽略不可行解的风险。代码里我先用w18、w25后面在参数调节部分再展开讲怎么调。2.3 种群边界与初始化染色体向量x的长度等于候选路径数K每个基因对应一条路径的流量。基因下界当然全部是0。上界不能随便拍一个数否则初始化会把大量基因赋成远超过容量的流量罚函数还没开始引导种群就已经全是不合理个体。单条路径的最优流量不可能超过它经过的所有弧的最小容量因此ub_k min( u_e for e in path_k )比如p5是1-2-4-5-6经过弧1、弧3、弧7、弧8容量分别是20、22、10、30那么p5的上界就是10。所有路径上界算出来后初始化种群这样写pop rand(Npop, K) .* repmat(ub, Npop, 1)这样每个个体从一开始就是数量级合理的路径流量后续进化只需要在容量叠加和预算约束之间做权衡效率会高很多。3. Matlab实现从数据到算子的完整代码3.1 网络数据与路径-弧关联矩阵进入Matlab代码。先把网络数据录进去。clear; clc; rng(7); % 固定随机种子方便复现 edgeFrom [1,1,2,2,3,3,4,4,5]; edgeTo [2,3,4,5,4,5,5,6,6]; u [20,18,22,15,18,20,10,30,35]; a [0.01,0.012,0.015,0.018,0.012,0.015,0.02,0.01,0.012]; b [1.0,1.1,1.2,1.5,1.1,1.4,1.3,1.0,1.2]; B 180; paths { [1 2 4 6] [1 2 5 6] [1 3 4 6] [1 3 5 6] [1 2 4 5 6] [1 3 4 5 6] }; E length(edgeFrom); K length(paths); % 路径-弧关联矩阵 pathMat(E, K) pathMat zeros(E, K); for k 1:K nodes paths{k}; for i 1:length(nodes)-1 s nodes(i); t nodes(i1); e find(edgeFrom s edgeTo t); if isempty(e) error(路径中包含未定义的弧); end pathMat(e, k) 1; end endpathMat是整个实现里最核心的桥梁。它是一个0-1矩阵行对应弧列对应路径。比如路径p1是1-2-4-6那么对应弧1、弧3、弧8的位置是1其余是0。后面计算弧流量时只需要做一次矩阵乘法fe pathMat * x这个向量化写法把路径流量映射到弧流量非常快比在循环里逐个判断某条弧被哪些路径经过要清晰得多。3.2 适应度函数把约束翻译成数字适应度函数直接实现前文的数学公式。function fit fitnessValue(x, pathMat, u, a, b, B, w1, w2) fe pathMat * x(:); % 弧流量 overCap sum(max(0, fe - u).^2); % 容量越限平方和 cost sum(a .* fe.^2 b .* fe); % 总成本 overBudget max(0, cost - B)^2; % 预算超支平方项 fit sum(x) - w1 * overCap - w2 * overBudget; end这段函数很小但包含了全部约束。sum(x)是总流量我们希望它尽量大后面的惩罚项会在容量越限或预算超支时把适应度拉低。注意fe是弧流量列向量u也是列向量max(0, fe - u)是对应每条弧的越限量平方后求和就是全局容量惩罚。有一个经验值得提罚函数项不要直接用max(0, fe - u)一定要平方。我之前试过线性惩罚进化后期会出现一个现象多条弧同时轻微越界累计惩罚不高结果保留了一个看起来局部还行、整体却不可行的解。平方项能放大多条弧同时越界的严重程度明显缓解这个问题。3.3 遗传算子锦标赛、SBX交叉和多项式变异实数连续编码的GA不适合用简单二进制交叉或单点交叉因为子代很容易突破边界。这里我用三个经典算子锦标赛选择、SBX模拟二进制交叉、多项式变异。锦标赛选择很简单随机挑2个个体留下适应度高的那个。这样既能保持选择压力又能防止少数超强个体迅速统治种群。function idx tournamentSelect(fit, k) cand randi(length(fit), k, 1); [~, best] max(fit(cand)); idx cand(best); endSBX交叉是实数编码GA里最常见的一种交叉算子它的特色是子代会距父代不远不近且分布指数eta可以控制子代与父代的接近程度。eta越大子代越接近父代。function [c1, c2] sbxCross(p1, p2, lb, ub, eta) if nargin 5 eta 15; end r rand(size(p1)); beta zeros(size(p1)); pos (r 0.5); beta(pos) (2 * r(pos)).^(1 / (eta 1)); beta(~pos) (2 - 2 * r(~pos) 1e-10).^(-1 / (eta 1)); c1 0.5 * ((1 beta) .* p1 (1 - beta) .* p2); c2 0.5 * ((1 - beta) .* p1 (1 beta) .* p2); c1 min(max(c1, lb), ub); c2 min(max(c2, lb), ub); end多项式变异是配合SBX使用的常用变异算子。它不像高斯变异那样完全随机发散而是以一个分布指数控制扰动幅度既能提供多样性又不至于把解踢得太远。function ch polyMut(x, lb, ub, pm, eta) if nargin 5 eta 20; end ch x; for i 1:length(x) if rand pm delta rand; if delta 0.5 deltaQ (2 * delta)^(1 / (eta 1)) - 1; else deltaQ 1 - (2 * (1 - delta))^(1 / (eta 1)); end ch(i) x(i) deltaQ * (ub(i) - lb(i)); end end ch min(max(ch, lb), ub); end这几个函数可以单独保存成m文件也可以全部放进同一个脚本的末尾作为局部函数。如果用的是老版本Matlab建议单独保存文件名分别对应函数名。3.4 主循环与精英保留有了网络数据、适应度函数和遗传算子主循环反而不复杂。核心逻辑是计算适应度选出当代最优个体作为精英然后反复执行选择、交叉、变异生成下一代。Npop 80; % 种群规模 MaxGen 120; % 最大进化代数 pc 0.85; % 交叉概率 pm 0.1; % 变异概率 eta_c 15; % SBX分布指数 eta_m 20; % 多项式变异分布指数 w1 8; % 容量越限惩罚权重 w2 5; % 预算超支惩罚权重 lb zeros(1, K); ub zeros(1, K); for k 1:K ub(k) min(u(pathMat(:, k) 1)); end pop rand(Npop, K) .* repmat(ub, Npop, 1); bestFit zeros(MaxGen, 1); avgFit zeros(MaxGen, 1); for gen 1:MaxGen fit zeros(Npop, 1); for i 1:Npop fit(i) fitnessValue(pop(i, :), pathMat, u, a, b, B, w1, w2); end [bestFit(gen), bestIdx] max(fit); avgFit(gen) mean(fit); elite pop(bestIdx, :); newPop zeros(Npop, K); newPop(1, :) elite; idx 2; while idx Npop if rand pc p1 pop(tournamentSelect(fit, 2), :); p2 pop(tournamentSelect(fit, 2), :); [c1, c2] sbxCross(p1, p2, lb, ub, eta_c); else c1 pop(tournamentSelect(fit, 2), :); c2 pop(tournamentSelect(fit, 2), :); end c1 polyMut(c1, lb, ub, pm, eta_m); c2 polyMut(c2, lb, ub, pm, eta_m); newPop(idx, :) c1; idx idx 1; if idx Npop newPop(idx, :) c2; idx idx 1; end end pop newPop; end代码里newPop(1,:) elite就是精英保留策略。每一代最优秀的个体直接进入下一代避免因为交叉变异被破坏。这是GA在连续优化里不翻车的重要保障尤其是罚函数场景一旦最优个体被破坏下一代可能要花很多代才能重新找到同等质量的解。参数方面我给出的这几个值是经过若干次试验后的经验值。交叉概率0.85意味着大多数个体参与交叉给种群足够的基因重组机会变异概率0.1不算高目的是维持局部扰动eta_c15让子代和父代比较接近搜索偏向局部精细eta_m20则让变异扰动幅度适中。3.5 结果输出与链路利用率的验证算法跑完之后最关心的不只是目标值还要确认容量约束真的没有越限。所以要单独算一遍最优解对应的弧流量和总成本。fe pathMat * elite(:); cost sum(a .* fe.^2 b .* fe); fprintf(最优总流量: %.2f\n, sum(elite)); fprintf(总成本: %.2f / %.2f\n, cost, B); for e 1:E fprintf(弧 %d - %d 流量 %.2f / 容量 %.2f\n, ... edgeFrom(e), edgeTo(e), fe(e), u(e)); end figure; plot(1:MaxGen, bestFit, -o, LineWidth, 1.5); hold on; plot(1:MaxGen, avgFit, --, LineWidth, 1.2); legend(最优适应度, 平均适应度, Location, best); xlabel(进化代数); ylabel(适应度); grid on;输出部分一条一条列出每条弧的流量和容量是为了快速判断有没有隐含的越限。有时候适应度看着挺好但某条弧悄悄超了0.5这在罚函数法中很常见。打印链路利用率之后便可以肉眼确认结果是否可落地。4. 实测效果与调参避坑4.1 用示例数据跑一次我用上面这段代码Npop80、MaxGen120、rng(7)Matlab R2023a测试单次运行时间大约3秒。一组比较典型的输出是最优总流量54.6总成本177.3预算上限180所有弧流量均在容量范围内如果把预算从180降到150最优总流量会落到45左右降到120大概只能承载36左右的流量。这个趋势很符合实际直觉预算越紧系统愿意承担的流量就越少这说明罚函数确实把预算约束导入了进化过程。还需要注意GA是随机算法换机器、换Matlab版本、换随机种子结果会有小幅波动。这不代表代码有问题说明算法存在随机性。工程上想复现就用rng固定种子想验证鲁棒性就换多次种子取统计结果。4.2 参数怎么调才不折腾很多人一上来就调大种群和代数其实没必要。GA调参有一些先后顺序按顺序走能少走弯路。先固定随机种子否则每次结果不同你无法判断修改参数到底是改好了还是随机波动。然后用一个小种群快速试探比如Npop30、MaxGen50跑通流程确认适应度曲线能上升再去加大规模。如果算法老早收敛到很低的适应度优先检查染色体上界是否合理或者是否因为罚函数权重太大个体还没积累足够优势就被淘汰了。如果最优解频繁越限优先增大容量惩罚权重w1而不是盲目增加进化代数。如果种群的多样性明显不够平均适应度和最优适应度曲线几乎重合适当提高变异概率pm或者用自适应变异进化前期变异大后期变异小。这里给一个常用的参数参考表参数建议范围我的示例取值种群规模 Npop50 ~ 20080最大代数 MaxGen100 ~ 500120交叉概率 pc0.7 ~ 0.90.85变异概率 pm0.05 ~ 0.20.1容量罚函数权重 w15 ~ 208预算罚函数权重 w25 ~ 2054.3 常见问题速查我把这类GA求解过程中容易遇到的几个典型问题整理成一张速查表。现象常见原因解决办法适应度一直为负总流量几乎为0罚函数权重过大惩罚项压过目标项降低w1、w2或调整初始化上界最优解某条弧仍超容量容量惩罚太轻或者权重没有跟上搜索过程提高w1并在结果输出中逐条检查链路每次运行结果差很多随机性太强种群规模或代数不足固定rng、增加Npop或MaxGen进化后期基本不动种群多样性不足精英过早统治提高变异概率pm或改用自适应变异路径太少导致总流量上限偏低候选路径集合不完整用K短路算法生成更多候选路径成本一直在预算上沿晃动预算惩罚权重合适但搜索边界太死允许解小幅越界最后再局部修复5. 再往前走一步工程化改造方向5.1 动态罚函数与局部修复文章前面的代码用的是固定权重罚函数简单有效。但如果你要处理大规模真实网络固定权重经常需要手动反复调。一个更稳的方案是动态罚函数进化早期让w1、w2小一点允许种群在较大空间里探索进化后期逐步加大权重把个体往可行域里压。这个思路有点净化退火实现成本很低只需要在每一代给权重乘一个缓慢增长的系数。局部修复也是工程里常用的补充手段。比如某条弧超容量了找到所有经过这条弧的路径按比例压缩它们的流量让弧流量回到容量以内。修复之后重新计算成本再检查预算是否超支。这种修复可以把不可行解变成可行解但要注意别过度压缩某一条路径否则会破坏染色体多样性。5.2 自动生成候选路径示例里的路径集合是我手工枚举的。真实项目网络可能有几百个节点、几千条弧不可能手工枚举。工程上通常会用K短路算法或者DFS加深度限制先跑出前K条候选路径。也可以先用最短路算法算一条基准路径再通过禁忌搜索生成一批差异性大的备选路径。候选路径的数量和质量直接影响GA的搜索上限。如果路径集合里漏掉了一条非常重要的高速直达路径GA再怎么优化也补不上这个缺陷。所以路径生成阶段的重要性不亚于算法本身。5.3 多目标化不只追求最大流量这个算例把目标定成最大总流量同时满足容量和预算约束是一个典型的单目标约束优化。但在很多真实场景里成本和流量是相互冲突的两个目标。这时候与其强行用罚函数合并成一个目标不如直接用多目标优化框架比如NSGA-II。跑完后得到一条Pareto前沿再让业务方根据实际偏好挑选解。前面写的遗传算子——锦标赛选择、SBX交叉、多项式变异——在多目标版本里几乎可以原样复用只是选择机制从适应度最高的赢改成非支配排序加拥挤度比较。如果你已经理解了这篇文章的代码转NSGA-II会非常顺畅。5.4 一点个人体会最后说点实在的。GA这类启发式算法代码真的不难难在设计一个能表达业务约束的编码。网络流问题里路径流量编码是我用过最顺手的一种它把最麻烦的等式约束直接消掉了罚函数是连接约束和搜索过程的最短路径。你要是第一次接触这类问题建议先不急着写代码拿着这篇文章里的小网络把pathMat和适应度函数在纸上手推一遍再上Matlab绝对比你直接跑通代码收获更大。我个人踩过最深的坑是贪心把所有约束都塞给罚函数导致约束权重像猜谜一样。现在我的做法是能消掉的约束尽量用编码消掉消不掉的才用罚函数必要时再做局部修复。这套思路不止适用于网络流很多组合优化和连续优化问题都可以照搬建议你保存下来下次遇到类似项目直接套用。