Matlab动态规划求解沙漠越野：状态定义、剪枝与回溯

简介：这份2020年数学建模B题的Matlab源代码面向参加数学建模竞赛、需要复盘优化类赛题的学生与培训者，围绕矿区路线规划、天气不确定下的资金最大化等问题给出可运行的求解实现。第二关针对双地图与双村庄场景，先以距离矩阵整合数据，再用单目标优化模型比较矿山开采与村庄补给方案；后续关卡引入天气分布概率预测、动态规划、数学期望和随机分布等方法，并讨论收益参数、沙暴条件对挖矿决策的影响。压缩包共114个文件，约66KB，以93个m脚本为核心，辅以15个mat数据、5个asv备份和1个xlsx表格，便于按关卡查看建模过程、运行结果与方案对比。已有1240人学习，适合用于理解从数据处理、模型构建到Matlab编程验证的完整思路，也可作为培训后的代码参考与排错对照。

1. 沙漠、水与天数：2020 年数学建模 B 题的决策骨架

2020 年数学建模 B 题给出一张沙漠地图，玩家要背着有限重量出发，在限定天数内抵达终点，沿途靠补给点补充水和食物，同时尽量少花钱。很多队伍拿到题后第一反应是找一份 Matlab 源代码跑通，结果只得到一串数字：换个地图或天气文件就不工作，询问“为什么这么走”也答不上来。这道题真正的骨架是“状态定义 + 按天推进 + 剪枝”，代码只是把这个骨架写下来。下面按我复现时的顺序走一遍：先把三张表读进来，再处理状态数量，最后用一份能自检的 Matlab 代码把策略输出。适合已经读过赛题、卡在“知道是动态规划但写不出来”这一步的人。

2. 把赛题摊成三张表：地图、天气、消耗在 Matlab 里怎么存

这道题最容易低估的是输入部分的杂乱。规则分散在几张附件表里：一张说哪些区域相邻，一张说每个区域的类型和地形，一张给未来若干天的天气。很多人一上来就写搜索，把if node == 3这类字面量散落在函数里，第三天发现补给点价格要分绿洲和村庄两档，只能全文件搜替换。我一般先把输入处理干净：把所有会变的数字关进一个结构体，其余函数只读不改，后面调参就变成改一个文件。

2.1 邻接关系：27 个节点用邻接矩阵还是边表

赛题地图的规模不大，节点数是两位数级别，邻接矩阵和边表两种都放得下。但两者在后续写法上差别明显。邻接矩阵的好处是“判断能不能走”是一行索引，P.adj(pos, to)直接出逻辑值，状态转移里不用写循环查找；边表的好处是省内存、方便按边遍历，但每次转移都要扫一遍边列表。

我选邻接矩阵，理由很实际：这道题的瓶颈不在内存，在状态数量。邻接矩阵占的那点空间和状态矩阵比可以忽略，而它能让我在转移函数里把“可达节点”写成一行：

% 取当前节点所有相邻节点，返回列向量 nextNodes = find(P.adj(pos, :))';

注意两个容易出错的地方。第一，地图是无向图，读边表时必须补对称边，只填一半会导致回程走不回来，而这类错误在前期测试中往往因为起点离终点很近而被掩盖。第二，不要把“相邻”和“可通行”混为一谈，沙暴天气下相邻也可能不能移动，这是转移函数的事，不是邻接矩阵的事，两者分开写，便于修改。

还有一种处理是给边加权重，比如把地形代价预先算到边上。我不建议这么做，因为这道题的代价不只取决于地形，还同时取决于当天天气和玩家携带量，预先算好等于把两个维度压成一个，后面想加“高温天沙地额外消耗”就得重来。邻接矩阵只存连通性，代价留给转移函数算，边界更清晰。

2.2 天气与地形：把规则从 if 里赶进查表

赛题的消耗规则可以写成一张二维表：行是天气，列是地形，格子里是当天该扣多少水、多少食物。一旦写成查表，转移函数里就不该再出现任何天气名称的字面量，只留下下标。

% consume(weatherIdx, terrainIdx, :) -> [waterBox, foodBox] % weatherIdx: 1 晴朗 / 2 高温 / 3 沙暴 % terrainIdx: 1 普通 / 2 沙地 P.consume = zeros(3, 2, 2); P.consume(1,1,:) = [1, 1]; P.consume(1,2,:) = [1, 1]; P.consume(2,1,:) = [2, 1]; P.consume(2,2,:) = [2, 2]; P.consume(3,1,:) = [2, 2]; P.consume(3,2,:) = [3, 2];

上面这些数值是占位的示意值，不是赛题原值，落地的时候按附件逐格替换就行，代码结构不用动。真正要守住的是“三维数组的下标语义”：第一维天气、第二维地形、第三维是水和食物。语义定死之后，转移函数里就只有squeeze(P.consume(w, t, :))这一种写法，谁也不会再写错顺序。

天气本身要按天索引，P.weather(day)返回当天的天气编号。如果赛题给的是多套天气序列，我一般把它们整理成矩阵，一行一套，求解时外层套一层循环，内层逻辑完全复用，这样问题二、问题三的扩展就是加一层循环，而不是复制一份代码改。

提示：查表的下标一定要在初始化函数里做一次合法性检查，比如assert(all(P.weather >= 1 & P.weather <= 3))。越界下标在 Matlab 里会直接抛错，但如果数组刚好够大，它不会报错，只会读到错误的值，这类问题查起来非常费时间。

2.3 一份只改数字不改逻辑的初始化函数

把三张表归一到一个函数里，是整个 Matlab 源代码里最值得先写好的部分。我的习惯是三份 csv 放在data/下，初始化函数只负责读表和拼结构体，不做任何业务判断。

function P = loadParams(dataDir) % 只做读表与结构体拼装，不做任何业务判断 nodeT = readtable(fullfile(dataDir, 'node.csv')); edgeT = readtable(fullfile(dataDir, 'edge.csv')); wthT = readtable(fullfile(dataDir, 'weather.csv')); P.N = max(nodeT.id); P.adj = false(P.N, P.N); for i = 1:height(edgeT) P.adj(edgeT.from(i), edgeT.to(i)) = true; P.adj(edgeT.to(i), edgeT.from(i)) = true; % 无向图必须补对称边 end P.type = zeros(P.N, 1); % 1 普通 2 起点 3 终点 4 补给 P.terrain = zeros(P.N, 1); % 1 普通 2 沙地 P.type(nodeT.id) = nodeT.type; P.terrain(nodeT.id) = nodeT.terrain; P.D = height(wthT); P.weather = wthT.weather; P.start = find(P.type == 2, 1); P.goal = find(P.type == 3, 1); P = loadRules(P); % 负重上限、单价、消耗表都在这里塞进来 assert(P.N > 0 && P.D > 0, '输入表为空，检查路径'); end

逻辑上有两处值得说。第一，P.type和P.terrain都是按节点编号做下标赋值的，所以node.csv里的id必须是连续正整数，中间断号会让数组错位；读表后加一句assert(isequal(sort(nodeT.id), (1:P.N)'))能挡住这个问题。第二，loadRules单独拆出来，是为了把“结构”和“数值”分开：改负重上限、改单价只动loadRules，读表逻辑一行不用碰。参数说明也简单，P.Wcap是负重上限（单位跟赛题一致，别中途换成克），P.price用长度为三的向量存三档价格，下标顺序写进注释里，后面查错时比翻文档快。

3. 状态怎么定义才不炸：从五元组到按天推进的层式 DP

写不出这道题的人，绝大多数卡在同一个地方：状态里到底放什么。放少了，转移不合法；放多了，状态数指数级增长，跑一晚上出不来结果。这一章讲的就是这条线的画法。

3.1 钱不进键：状态去重时最容易犯的错

最直观的状态是五元组：(位置, 第几天, 剩余水, 剩余食物, 剩余钱)。如果把它整个当去重键，问题立刻出现：两个状态位置、天数、水、食物都一样，只有钱差一点，它们会被当成两个不同的状态分别扩展；而实际上钱少的那个已经被钱多的完全支配了，后面的路一模一样，再展开纯属浪费。

正确的做法是把钱从键里拿出来，当值用。键是(位置, 天数, 水, 食物)，值是“在这个键下能拿到的最大剩余钱”。这样同一组资源状态只保留最优的那个，状态数量直接掉一个数量级。这条规则听起来简单，但我在不止一份流传的 Matlab 源代码里见过反例：作者把所有维度拼成一个字符串当containers.Map的键，跑小规模数据没问题，一上完整地图就卡死。

这里还有一个连带的好处：目标函数变成“最后一层里位置等于终点、天数不超过上限的那些状态中，钱最大的那个”，取最大值是一行代码，不需要额外设计评价函数。

注意：去重不能只看键相同就随便保留一个。必须显式地“按键分组，取组内最大资金”，如果只是用unique去重，Matlab 保留的是排序后的第一行，结果依赖于行顺序，跑两次可能不一样。

3.2 三条硬剪枝：负重、资金、返程下界

状态压下来之后，还需要剪枝，否则大地图上依然吃不住。我固定用三条，都是“活着但已经被支配”的状态。

第一条是负重剪枝。水加食物的总重量不能超过负重上限，购买动作产生新状态时要立刻判定，超了就不生成。这条剪枝的正确性来自物理约束，没有商量余地，但要注意单位换算：如果水的单位是“箱”、每箱固定重量，那么重量是箱数乘单箱重量，不能拿箱数直接和上限比。

第二条是资金剪枝。购买之后钱不能为负，为负的状态直接丢弃。这条看似显然，但在“先扣消耗再补给”的顺序下容易漏判，因为扣消耗本身不涉及钱，很多人就忘了在购买分支里再查一次。

第三条是返程下界剪枝，也是收益最大的一条：预先算出每个节点到终点的最短步数（忽略天气，按最少消耗估算），如果当前剩余天数小于这个步数，或者当前剩余水粮连“最乐观情况下走完这段路”都不够，直接丢掉。最乐观情况可以用“每天都按最低消耗”来估，虽然偏松，但足以砍掉大量在远处徘徊的状态。

function keep = pruneByBound(P, S, d) % S: [pos, day, water, food, money] % 预计算 P.dist2goal(pos) 为最少步数，P.minCost(pos) 为乐观最小消耗 remainDay = P.D - d; needDay = P.dist2goal(S(:,1)); needWater = P.minCostWater(S(:,1)); needFood = P.minCostFood(S(:,1)); keep = (S(:,3) >= needWater) & (S(:,4) >= needFood) & (remainDay >= needDay); end

参数说明：dist2goal用 BFS 在邻接矩阵上求一遍就行，几十个节点是毫秒级；minCost我一般按“剩余路段全是最低消耗地形加最低消耗天气”估，偏保守但计算快，如果发现剪枝效果不够，可以把它换成“按实际天气序列做一次反向 DP”，代价是初始化时间变长。

第三条剪枝要小心一个坑：它必须是“必要条件”而不是“充分条件”。也就是说，不满足就一定能丢，满足了不代表一定能到。如果把估算写得太乐观，会误杀可行解，结果比不剪还糟——你会得到一个看起来像最优解的错误答案，而且没有任何报错提示。

3.3 unique 加 accumarray：一晚上跑完所有天气组合

状态推进按天分层，每层是一个矩阵，每行一个状态。层与层之间做三件事：扩展、合并、剪枝。合并这一步用unique拿唯一键，用accumarray取组内最大值，两行就能写完，比手写哈希表快得多，也不用担心containers.Map的字符串键开销。

function [S, pick] = mergeStates(cand) % cand: 候选状态 [pos, day, water, food, money] [U, ~, g] = unique(cand(:, 1:4), 'rows'); % 前四列做键 bestM = accumarray(g, cand(:, 5), [], @max); % 组内最大资金 [~, ord] = sortrows([g, -cand(:, 5)]); % 按组号升序、资金降序 [~, first] = unique(g(ord)); % 每组第一行即最大资金 pick = ord(first); S = cand(pick, :); S(:, 5) = bestM(g(pick)); end

逐句解释一下。unique的第三个输出g是“每个候选行属于哪个唯一键”的分组号，长度和候选行数一致。accumarray用g做分组、用第五列做值、用@max做聚合，得到每个唯一键对应的最大资金。sortrows把分组号和负资金拼在一起排序，效果是“同组的排在一起，组内资金从大到小”，再取每组第一行，就拿到了那一行原始数据。最后把这一行的资金改写成组内最大值，保证矩阵里的值和自己算出来的一致。

这一步最容易出的错是sortrows的列顺序写反。如果写成[g, cand(:,5)]，组内是资金升序，取第一行拿到的是最小值，结果会莫名其妙偏小，而且小得不多，肉眼很难发现。跑完先拿一组小数据手算对比，能很快定位。

4. Matlab 源代码骨架：四个函数把求解流程串成一条线

前面两章是数据层和算法层，这一章落到文件结构。我不喜欢把所有东西塞进一个.m文件，出了问题只能靠dbstop一行行看。四个函数分工，每个都能单独测。

4.1 主脚本的调用顺序与目录约定

目录就三层：data/放三份 csv，src/放函数，根目录放main.m和validate.m。main.m里只做四件事，顺序固定，中间任何一步失败都能一眼看出问题在哪。

% main.m clear; clc; P = loadParams('data'); % 1 读表与参数 P = precompute(P); % 2 预计算 dist2goal / minCost [histS, histP] = runDP(P); % 3 层式 DP，保留每层状态与父指针 [plan, best] = backtrack(P, histS, histP); % 4 回溯出策略表 disp(best); writetable(plan, 'plan.csv');

precompute干两件事：在邻接矩阵上 BFS 求每个节点到终点的最少步数，以及按最乐观假设求最少消耗。这两样东西和天气无关，一次算完反复用。runDP是消耗资源的主要部分，返回每层状态矩阵和父指针数组，内存占用与状态数成正比。backtrack只做查表，很快。

把runDP和backtrack拆开的好处是：调参的时候只需要重跑后半段，不用整个重算；想换一种回溯输出格式（比如按天出一张明细表），也不影响求解。

4.2 状态转移函数：停留、移动、补给三类动作怎么枚举

转移函数是整个代码里唯一需要理解赛题规则的地方，逻辑写错基本都写在这里。我的做法是按动作分类枚举，每一类产出候选状态，最后统一合并。

function C = expandState(P, s, day) % s: [pos, day, water, food, money]，返回候选状态矩阵 pos = s(1); w = s(3); f = s(4); money = s(5); C = zeros(0, 5); acts = [0, find(P.adj(pos, :))]; % 0 表示停留，其余为相邻节点 for k = 1:numel(acts) to = acts(k); if to == 0 tIdx = P.terrain(pos); % 停留按当前所在地形消耗 nPos = pos; else tIdx = P.terrain(to); % 移动按目标地形消耗 nPos = to; end c = squeeze(P.consume(P.weather(day), tIdx, :))'; nw = w - c(1); nf = f - c(2); if nw < 0 || nf < 0, continue; end if P.type(nPos) == 4 % 到达补给点，枚举购买 C = [C; buyBranch(P, nPos, nw, nf, money, day)]; %#ok<AGROW> else C = [C; [nPos, day, nw, nf, money]]; %#ok<AGROW> end end end

buyBranch里枚举购买量，量程用负重上限卡住：ceil((P.Wcap - nw*P.wWater - nf*P.wFood) / P.wWater)是水的理论上限，食物同理。这个上限在大地图上依然偏大，我一般再加一刀：“单次购买超过若干天正常消耗的部分不枚举”，因为它几乎不可能出现在最优解里。

参数上要注意两点。第一，P.weather(day)用的是当天的天气，转移的目标层是day之后的状态，所以写循环时day的取值要和状态里的天数对齐，别一边从 0 开始一边从 1 开始。第二，消耗表的列顺序必须和c(1)/c(2)的语义一致，水和食物反过来会让结果看起来“能跑但总是差一点”，这类错误最难查。

提示：expandState我强烈建议单独写测试，固定一个小的三节点地图，手算三种天气、三个动作的候选状态，和函数输出逐一比对。这一步花二十分钟，后面能省掉几个小时的排查。

4.3 用 parent 索引回溯出每天走哪条路

DP 只给出最后一层的最优值，要输出策略表还需要回溯。回溯靠的是每层保存的父指针：新状态的每一行记录它由上一层哪一行扩展而来。

function [histS, histP] = runDP(P) histS = cell(P.D + 1, 1); histP = cell(P.D + 1, 1); histS{1} = [P.start, 0, P.water0, P.food0, P.money0]; for d = 1:P.D C = []; CP = []; for i = 1:size(histS{d}, 1) newRows = expandState(P, histS{d}(i, :), d); if isempty(newRows), continue; end C = [C; newRows]; %#ok<AGROW> CP = [CP; repmat(i, size(newRows, 1), 1)]; %#ok<AGROW> end [histS{d+1}, pick] = mergeStates(C); histP{d+1} = CP(pick); % 与合并后的状态一一对应的父行号 end end

这里CP(pick)是关键：合并时把一批候选行压成一行，父指针必须跟着一起压，用同一个pick索引。如果只对状态矩阵做合并、忘了同步父指针，回溯出来的路径会指向毫不相干的状态，表现是“策略表能生成，但逐天重放时水量对不上”。

回溯本身从最后一层的最优行开始，一行一行往上找父节点：

function [plan, best] = backtrack(P, histS, histP) last = histS{end}; cand = find(last(:,1) == P.goal); [best, k] = max(last(cand, 5)); row = cand(k); plan = zeros(P.D, 5); for d = P.D + 1:-1:2 plan(d-1, :) = [d-1, histS{d}(row, 1), histS{d}(row, 3), histS{d}(row, 4), histS{d}(row, 5)]; row = histP{d}(row); end end

plan每行是“第几天、在哪个节点、剩余水、剩余食物、剩余钱”。有了这张表，人工核对就变得很轻：随便挑一天，看看水量变化是否等于当天规则该扣的量加减补给量。

4.4 三条自检断言：跑完先别急着交

代码跑出结果不等于结果对。我会在main.m末尾加三条断言，任何一条不通过就说明逻辑有问题，而不是参数有问题。

第一条，终点可达：assert(~isempty(find(histS{end}(:,1) == P.goal, 1)), '末层没有终点状态')。这条防的是“跑完了但根本没走到终点”却因为取最大值时没加条件而误判成功。

第二条，逐天重放一致：把plan从第一天重新算一遍，每一天的水、食物、钱必须和表里完全相等。这条能抓住所有“转移函数和回溯不同步”的问题。

第三条，约束全通过：逐天检查负重不超限、资金不为负、每天所在地和相邻关系一致。这条最啰嗦但最值得写，边界条件基本都在这露头。

5. 五类高频翻车点：从状态爆炸到回溯错位的排查清单

下面这五条，是我在复现这道题和帮别人看代码时反复遇到的。每一条按“现象—原因—解决”写，照着对号入座比看文档快。

5.1 现象：跑十分钟内存爆掉，任务管理器一片红

原因几乎都是状态没有真正合并。常见写法是把状态拼成字符串当哈希键，或者把资金也放进了键，导致同一资源组合被当成多个状态。还有一种隐蔽情况：合并做了，但unique之后没有重新赋返回值，状态矩阵还是老的那一份。

解决分三步查。先打印每层状态数，正常应该先涨后稳，如果一路线性涨而不回落，说明合并没生效。再看键的组成，确认资金不在键里。最后在合并函数里加一句assert(size(U,1) < size(cand,1) || size(cand,1) == 0)，候选行数没被压缩就说明键选得太细了。

5.2 现象：结果比手算的还差，明显有一条更省的路

原因通常是剪枝写成了充分条件。返程下界那条最容易出问题：估算剩余消耗时用了“最省”的假设，但同时又把“剩余天数”用“最短步数”卡死，如果玩家计划在补给点停留一天，最短步数够但实际天数不够，可行解就被误杀了。

解决是给剪枝留冗余。天数下界乘一个小于 1 的松弛系数再比较，或者干脆关掉这条剪枝跑一遍小数据，打开再跑一遍，两次结果必须一致；不一致就说明剪枝有问题。

注意：剪枝的验证不能只看最终答案的数字，因为误杀可能刚好不改变最优值。要在小规模数据上把“剪枝前状态集合”和“剪枝后状态集合”都打出来，确认被丢的确实被支配了。

5.3 现象：策略表里某天买了负数的水

原因在购买枚举的边界没有卡严。负重上限算剩余可承载量时，如果用了“当前已装箱数”而不是“当前重量”，或者单位换算时把箱数当重量，就可能算出负的量程，循环里0:负数直接不执行，或者更糟，产生了负数下标。

解决是把量程单独算出来并断言行非负：ub = floor((P.Wcap - nw*P.wWater - nf*P.wFood) / P.wWater); assert(ub >= 0);。这类问题查起来很快，但如果没有断言，它会悄悄地让某些状态永远不被生成，结果偏保守。

5.4 现象：回溯出来的路径和状态表对不上，第一天就错

原因是父指针没有跟着合并一起压。mergeStates返回的pick是“合并后第 i 行来自候选矩阵的哪一行”，父指针数组必须用同一个pick索引，写成histP{d+1} = CP(pick)。写成histP{d+1} = CP就会长度对不上，Matlab 有时不报错，只取前若干个，回溯到后面全乱。

解决是把mergeStates的接口改成同时返回状态和索引，中间不拆开传。再配合 4.4 里的逐天重放断言，任何错位都会在第一天暴露。

5.5 现象：换一组天气文件，结果全变但看不出哪里错了

原因是天气索引和天数索引错位。赛题的天气是按天给的，状态里的day从 0 开始还是从 1 开始，两条线必须统一。常见写法是状态里day从 0 开始、循环变量从 1 开始，然后P.weather(d)取到了前一天的天气，整体偏移一天。

解决是在读表后立刻把天气序列转成从 1 开始索引的向量，并在转移函数开头加assert(day >= 1 && day <= P.D)。另一个做法是把“第几天”直接写进状态矩阵并在合并时保留，虽然多一列，但调试时肉眼可见，值得。

6. 从能跑到可信：上界对拍与策略重放校验

代码能出数字，和数字可信，中间隔着一整套验证。我习惯做两件事：给结果配一个上界，配一个独立重放器。

上界怎么来？把问题松弛掉一部分约束，让它变简单但一定不会更差。比如忽略天气最劣的影响、允许水粮按小数购买、允许在任意节点补给，这样求出来的最优值一定是原问题的上界。把自己的解和这个上界比，差距在合理范围内就说明没算错方向；如果自己的解反而超过了上界，那一定是约束写漏了，这个信号比任何日志都直接。

重放器是另一个函数，故意不复用expandState，按plan.csv逐天重算，用最简单直白的写法：

function ok = replay(P, plan) % plan 每行: 天数 节点 水 食物 钱；独立实现，不复用求解期函数 w = plan(1,3); f = plan(1,4); m = plan(1,5); for i = 2:height(plan) dw = w - plan(i,3); df = f - plan(i,4); dm = plan(i,5) - m; assert(dw >= 0 && df >= 0, '第 %d 天消耗为负', i); assert(w*P.wWater + f*P.wFood <= P.Wcap, '第 %d 天超重', i); if plan(i,2) ~= plan(i-1,2) % 发生了移动 assert(P.adj(plan(i-1,2), plan(i,2)), '第 %d 天走了不存在的边', i); end w = plan(i,3); f = plan(i,4); m = plan(i,5); end ok = plan(end,2) == P.goal && m >= 0; end

这份重放器的价值在于“独立实现”四个字。求解期的转移函数有一堆剪枝和合并，逻辑绕；重放器只用最笨的方式顺着走一遍，两边都算出同一个结果，才敢说数据没问题。我一般会拿它跑三组数据：一组手算的小地图、一组完整地图、一组随机生成的天气序列。前两组用来验正确性，第三组用来验稳定性。

最后留一个习惯。每次改完参数，先把plan.csv导出来，挑中间三天人工看一遍水量变化，比看任何汇总指标都可靠。血泪经验是：越野类题目的错误很少出现在第一天和最后一天，几乎都藏在中段那几个补给点前后的收支平衡上。希望帮到你。

本文还有配套的精品资源，点击获取