简介：面向数学建模竞赛选手与水资源优化研究者，这份2022年美赛B题M奖论文围绕美国西部千年大旱背景下的水资源公平分配与调度展开。论文提出FHSD公平分层供需模型，以多目标非线性规划为核心，先用马尔科夫链与蒙特卡洛模拟分析连续干旱天数并确定模型更新频率，再构建以年为周期的供需函数，引入公平效益、价格、成本与满意度等参数，并采用NSGA2求取帕累托非劣解，在公平性与经济效益之间权衡；文中还给出年度供水量分配、泄流控制及敏感性分析。压缩包共1个PDF文件，约1.16MB，保留完整论文的摘要、模型推导、求解结果与结论，适合参考美赛论文结构、多目标建模表达及算法落地写法。已有3065人学习下载，可供备赛者拆解赛题思路、对照图表与结论组织自己的论文框架。

1. 从连续干旱天数反推更新周期：一个反直觉的鲁棒性起点

先给一个容易被忽略的结论：在这套水资源分配模型里，最关键的一个参数——模型多久更新一次——不是调度经验拍出来的，也不是按自然月固定下来的，而是从“连续干旱天数出现概率最大是多少”倒推出来的。换句话说，先假设最坏情况会持续多久，再决定模型多久重算一次。这份围绕多目标非线性规划构建的公平分层供需模型，把鲁棒性放在了精度前面：水库调度成本高、调整灵活性差，更新过于频繁会让调度方案反复变动，更新太慢又可能在水位已经逼到最小发电水位时才发现约束不可行。这个次序上的取舍，恰恰是很多同类模型容易忽略的地方。

2. 马尔可夫链与蒙特卡洛：干旱概率如何被压缩成一个更新周期

更新周期这件事，本质上是在回答一个概率问题：未来一段时间里，连续多少天不下雨的可能性最大？这个问题看起来需要长期天气预报，但长期预报本身不可靠，于是建模思路转向了统计结构——不预测具体哪天下雨，只描述“今天的状态如何影响明天的状态”。

2.1 状态空间、转移矩阵与无后效性假设

把每天的天气压成两个状态：0 表示无雨，1 表示有雨。状态空间只有两个元素，转移概率矩阵是一个 2×2 的方阵，第 i 行第 j 列表示“今天处于状态 i 时，明天转到状态 j 的概率”。矩阵的每一行加起来必须等于 1，这是概率归一化的硬要求，也是后面所有推导成立的前提。

用历史日序列估计这个矩阵，做法非常直接：数频次，再按行归一化。

import numpy as np def build_transition(series, n_state=2, alpha=1.0): """由历史日状态序列估计一阶转移矩阵。 series: 0/1 序列，0=无雨，1=有雨 alpha : 拉普拉斯平滑系数，防止某行样本为零导致归一化失败 """ counts = np.full((n_state, n_state), alpha) # 先铺平滑底数 for prev, cur in zip(series[:-1], series[1:]): counts[prev][cur] += 1.0 # 统计相邻日共现频次 row_sum = counts.sum(axis=1, keepdims=True) return counts / row_sum # 行归一化，每行和为 1 P = build_transition(history) print(P) # P[i][j] = P(明天=j | 今天=i)

这里有个必须说清的假设：无后效性。它意味着给定今天的状态，前天、大前天的状态对明天没有任何额外信息量。对日尺度天气来说这个假设是粗糙的——真实的降水系统有天气尺度的持续性，也有季节内的低频振荡，一阶链会把“连续第 5 天干旱”和“连续第 1 天干旱”的次日转移概率算成同一个值。但建模要的不是精确复现气候系统，而是一个能稳定输出“连续干旱段长度分布”的生成器。一阶链的参数只有 2 个自由度，估计方差小，在只有几年日观测数据的条件下，比硬上高阶链更不容易过拟合。这是一个典型的偏差-方差取舍：宁可系统性低估持续性，也不要在样本外抖动。

平滑系数 alpha 不是可有可无的装饰。某些季节分段里可能整月没有有效降水，“有雨→有雨”这个格子频次为 0，归一化之后该行出现 0/0，整个矩阵变成 NaN，后面所有计算全部污染。加一个极小的先验（拉普拉斯平滑或杰弗里斯先验）就避开了这个数值陷阱。

2.2 收敛性与稳态分布：为什么可以跳过逐日模拟

一阶链有个很好的性质：只要矩阵非周期且任意两个状态互通（从任何状态出发都能以正概率到达任何状态），那么 P 的 n 次幂当 n 趋于无穷时会收敛到一个所有行都相同的矩阵，每一行都是同一个平稳分布 π。π 满足 π = πP，且与初始状态无关。

这条性质的实际价值在于：如果只关心“长期来看干旱日占比多少”，可以完全不跑模拟，直接对矩阵做幂迭代，几毫秒就能得到答案。

def stationary(P, tol=1e-13, max_iter=100000): """幂迭代求平稳分布，避免直接做特征分解带来的数值问题""" pi = np.full(P.shape[0], 1.0 / P.shape[0]) for _ in range(max_iter): nxt = pi @ P # 左乘转移矩阵，向前推一步 if np.abs(nxt - pi).sum() < tol: # L1 残差足够小即认为收敛 return nxt pi = nxt raise RuntimeError("未收敛，检查矩阵是否可约或周期")

收敛速度由矩阵的第二大特征值模长决定，两个状态的链一般几步就收敛，所以这个迭代几乎不花代价。真正的用途，是给后续蒙特卡洛模拟提供一个交叉校验：模拟出来的长期干旱频率应该和 π 中的干旱分量接近，如果差得远，说明模拟代码里的状态转移写反了，或者随机数种子被复用导致样本相关。这种“解析解校验模拟器”的习惯，在写任何随机仿真时都值得保留。

2.3 从日状态序列到“连续干旱段长度”的分布

知道干旱日的长期占比还不够，需要的是“连续干旱段有多长”的分布。这一步交给蒙特卡洛：反复模拟长度为一个月的日天气序列，统计每条序列里最长连续无雨段的天数，得到一张直方图。

rng = np.random.default_rng(20220213) def longest_dry_run(seq): """返回序列中最长连续 0 的长度""" best = cur = 0 for s in seq: cur = cur + 1 if s == 0 else 0 best = max(best, cur) return best def simulate_month(P, days=30, n=200000): """模拟 n 个月，返回最长连续干旱天数的经验分布""" hist = np.zeros(days + 1, dtype=np.int64) for _ in range(n): seq = np.empty(days, dtype=np.int8) seq[0] = 0 if rng.random() < 0.9 else 1 # 初始状态按季节先验抽样 for t in range(1, days): seq[t] = 0 if rng.random() < P[seq[t - 1]][0] else 1 hist[longest_dry_run(seq)] += 1 return hist / n dist = simulate_month(P) T_star = int(np.argmax(dist)) # 概率最大的连续干旱天数，即更新周期

样本量取到 20 万次，是因为关注的是分布的众数而不是均值。众数对抽样噪声很敏感，样本少的时候直方图顶端会随机跳来跳去，跑两次得到不同的 T。20 万次在普通笔记本上是几秒的量级，代价可以接受。这里还必须说明：取众数而非期望，是因为更新周期要的是一个“最可能出现的整段干旱长度”，期望会被长尾拉高，导致更新周期定得过长，反而在最常见的中等干旱段里失去响应能力。

分季节建模之后，两个水库、三个时段各得一组结果，形态大致如下：

表格里最值得琢磨的是下游水库在 9–12 月那个 17 天。干旱概率高达 0.959，更新周期反而最长。常见的理解是：这个时段下游水库几乎完全依赖上游泄流补给，本地降水对入库流量的边际贡献趋近于零，连续干旱段的长度主要由上游调度节奏决定而不是由本地天气的随机性决定，因此统计出来的“最可能连续干旱段”被拉长。这也意味着，下游水库在这个季节的模型更新，本质上是在跟随上游的决策周期，而不是跟随天气。把这件事在模型里显式化，比把它藏在参数里要诚实得多。

3. 供需函数：水位-库容回归怎么变成一条可求导的供给曲线

有了更新周期，下一个问题是每一期到底能供多少水、发多少电。供给侧的物理约束全部围绕一个量展开：水库水位。而需求侧要的是水量，优化器要的是可导的目标函数。中间这座桥，就是把“水位”翻译成“可用水量”的函数关系。

3.1 双对数线性回归与它的反解

水库的水位-库容关系本质上是地形积分，理论上没有解析式，实测数据画出来是一条单调上凸的曲线。直接做多项式拟合会在端点外剧烈震荡，工程上的常规做法是先取对数再做线性回归：ln h = a + b·ln S，其中 h 是水位，S 是库容。

import numpy as np from scipy import stats def fit_level_storage(S, h): """双对数线性回归：ln(h) = intercept + slope * ln(S)""" x, y = np.log(S), np.log(h) res = stats.linregress(x, y) # 返回斜率、截距、R²、p 值、标准误 return {"slope": res.slope, "intercept": res.intercept, "r2": res.rvalue ** 2, "p": res.pvalue} def storage_from_level(h, fit): """反解：给定水位，求对应库容。优化器里真正要用的是这个方向""" return np.exp((np.log(h) - fit["intercept"]) / fit["slope"])

回归的方向和水位、库容的因果方向相反：拟合时以库容为自变量、水位为因变量，因为地形决定了库容是“累积量”，水位是它的函数；但优化时要以水位为约束（水位不能低于死水位、不能低于最小发电水位），库容是被算出来的量，所以必须反解。这两个方向的残差结构完全不同——正向回归最小化的是水位残差，反解之后库容的相对误差会被斜率放大或压缩，斜率绝对值小于 1 时误差被压缩，大于 1 时被放大。这一点在设定容差时必须心里有数。

双对数模型的系数有明确的物理含义。斜率 b 反映库容随水位增长的快慢：b 越小，说明单位水位上升带来的库容增量越大，库区越“扁平”；b 越大，说明库区越“陡峭”，水位稍微掉一点库容就掉很多。R² 高、p 值远小于 0.01，只能说明在观测区间内拟合得好，完全不能保证外推到死水位附近仍然成立——这是下一章要专门讲的坑。

3.2 水电站出力公式与分段供给结构

水电站的实际出力由一段很经典的式子给出：P = ρ·g·Q·H·η，其中 ρ 是水密度，g 是重力加速度，Q 是过机流量，H 是上下游水位差（水头），η 是综合效率。把常数打包成一个系数 k，就得到 P = k·Q·H。这里面 H 是上游水位减去下游水位，而上游水位又通过上一步的反解和库容挂钩，于是出力最终变成了库容和泄流量的二元函数。

关键在于这个函数是分段的，断点就在最小发电水位：

• 上游水位高于最小发电水位时，机组满发，出力与 Q、H 的乘积成正比；
• 水位落在这条线和死水位之间时，可开机机组减少，出力随水位近似线性衰减；
• 水位低于死水位时，出力为零，水库只剩生态基流和生活供水的功能。

def hydropower(Q, h_up, h_down, h_gen, h_dead, k=9.81 * 0.85): """分段出力：h_gen 为最小发电水位，h_dead 为死水位""" if h_up <= h_dead: return 0.0 # 死水位以下，无电可发 if h_up >= h_gen: head = h_up - h_down # 满发段，水头取全落差 else: head = (h_up - h_dead) / (h_gen - h_dead) * (h_gen - h_down) return k * Q * head # kW 量级，视单位而定

这个分段结构直接决定了供给函数的形状：供水量和发电量不是两个独立的目标，它们共享同一个决策变量 Q。多放水给农业，就少了一份可用于发电的水头；为了保证发电而压低下泄流量，农业和工业就得让步。这个耦合关系是多目标优化能成立的根本原因，也是后面所有冲突的源头。

供给函数最终被写成“给定水位 h，在满足最小生态流量和最大泄流能力的前提下，可用于各行业分配的水量上限”。取水价按每千瓦时约 10 美分的量级折算成单位水量的发电机会成本，让“水”和“电”在同一个量纲下可比。

3.3 分层需求：四个系数如何把“公平”变成可计算的量

需求侧不能简单地写成“每个区需求量是常数”。真实的需求随价格、输水成本、缺水程度变化，而且不同行业对缺水的敏感度差异巨大。分层结构就是在这里起作用的：农业、工业、居民三层，每层单独建需求函数，层与层之间的权重由政策偏好决定。

四个系数各有分工：

满意度这一步是整套需求建模里最容易被做坏的地方。如果直接用线性的供给需求比，那么总社会效益就是各区的加权和，优化器的必然选择是把水集中给“权重高、人口多”的区，其余区被压到零，这显然不是公平。用凹函数处理就能改变这个结论：

def satisfaction(supply, demand, eps=1e-9): """凹化满意度：缺水初期边际痛感强，后期钝化。 用平方根而非线性，避免优化器把水全部堆给高权重区。""" ratio = np.clip(supply / (demand + eps), 0.0, 1.0) return np.sqrt(ratio) def transport_cost(distance, unit_cost, volume): """输水成本：距离线性计入，忽略地形起伏与泵站扬程""" return unit_cost * distance * volume

平方根在这里不是随手选的。它的导数在 ratio 趋近 0 时趋于无穷，意味着“从完全没水到有一点水”带来的社会效益增量极大，优化器会优先保证每个区都不至于断水；而在 ratio 接近 1 时导数趋于零，多给的水几乎不增加社会效益，于是余水会被让给发电或其他区。这个形状和决策者的真实偏好是对得上的，也是把“公平”这个模糊词变成可优化目标的最直接手段。

4. 多目标非线性规划与 NSGA-II：帕累托前沿是怎么被“养”出来的

到这里，供给函数和需求函数都有了，剩下的事是把它们塞进一个优化框架里，并且承认一件事：不存在唯一最优解。给农业多分一点水，经济效益可能上升而社会公平下降；保发电，农业就得休耕。目标之间是真正冲突的，所以输出必须是一组解，而不是一个解。

4.1 两个目标、一组约束与决策变量规模

经济效益目标可以写成：售水收入减去输水成本，再减去因下泄流量减少而损失的发电收益。后一项是关键，它把水权和发电权放到了同一个账本上。社会效益目标则是各区各行业满意度的加权和，权重由公平-效益系数给出。两个目标都取最大化，统一转成最小化形式后写成向量目标 F = (−经济效益, −社会效益)。

约束包括：总供水量不超过可用水量、各水库库容在死库容和最大库容之间、水位在死水位和最大蓄水位之间、下泄流量不超过泄流能力且不低于最小生态流量。

决策变量的规模值得算一下：5 个受水区 × 3 个行业 × 2 座水库 = 30 个连续变量。如果按日或按月分配，变量数要乘上时间步数，一年按日算就是 10950 维，NSGA-II 的种群规模和迭代代数都要跟着涨，计算量会失控。这正是分层模型存在的直接动因——先在年度层面把总量分掉，再在更新周期层面细化，两级各自只有几十维，都落在遗传算法能高效搜索的区间里。

4.2 非支配排序、拥挤度与精英保留

NSGA-II 的核心是两件事：用非支配排序把种群分层，用拥挤度距离在同一层内保持多样性。

支配关系的定义是：解 p 支配解 q，当且仅当 p 在所有目标上都不劣于 q，且至少在一个目标上严格更优。第一层是不被任何解支配的解集合，也就是帕累托前沿；把它们剥掉之后剩下的解里再找不被支配的，构成第二层，依次类推。

def fast_nondominated_sort(objs): """objs: (N, M) 全部按最小化处理。返回各层索引，fronts[0] 即帕累托前沿""" N = len(objs) dominated_by = [[] for _ in range(N)] # 被 p 支配的解 dom_count = [0] * N # 支配 p 的解的个数 fronts = [[]] for p in range(N): for q in range(N): if p == q: continue pq = all(objs[p] <= objs[q]) and any(objs[p] < objs[q]) qp = all(objs[q] <= objs[p]) and any(objs[q] < objs[p]) if pq: dominated_by[p].append(q) elif qp: dom_count[p] += 1 if dom_count[p] == 0: fronts[0].append(p) # 无人支配，进入第一层 i = 0 while fronts[i]: nxt = [] for p in fronts[i]: for q in dominated_by[p]: dom_count[q] -= 1 # 剥掉一层后重新计数 if dom_count[q] == 0: nxt.append(q) i += 1 fronts.append(nxt) return fronts[:-1]

复杂度是 O(MN²)，M 是目标数（这里是 2），N 是种群规模。种群取 200 时，一次排序大约几万次比较，完全不是瓶颈。

拥挤度距离负责在同一层内部排序。做法是按每个目标分别排序，边界解的距离设为无穷大（保证极端解一定被保留），中间解的距离取相邻两解在该目标上的差值再归一化，对所有目标求和。

def crowding_distance(front, objs): """同一非支配层内的拥挤度，值越大说明周围越稀疏，越该保留""" dist = {i: 0.0 for i in front} for k in range(objs.shape[1]): order = sorted(front, key=lambda i: objs[i, k]) dist[order[0]] = dist[order[-1]] = float("inf") # 边界解必留 span = objs[order[-1], k] - objs[order[0], k] + 1e-12 for a, b, c in zip(order[:-2], order[1:-1], order[2:]): dist[b] += (objs[c, k] - objs[a, k]) / span # 归一化后累加 return dist

拥挤度对水资源问题的重要性，比在一般测试函数上要大得多。决策者拿到前沿之后要回答的问题是“多给农业 10% 的水，会牺牲多少发电收益”，这需要在整条前沿上有均匀分布的点可供比较。如果种群全部塌缩到前沿的一端，那张表就没法用。

选择阶段用二元锦标赛：先比非支配层级，层级相同再比拥挤度。交叉用模拟二进制交叉（SBX），变异用多项式变异，两者都是为连续变量设计的、能控制子代靠近父代的邻域算子。精英保留策略把父代和子代合并后再做一次非支配排序和截断，保证已经找到的好解不会因为随机性丢失。

4.3 两级分配：年度总量与单元时间量的衔接

整个求解被拆成两级。第一级在年度尺度上跑 NSGA-II，决定每个区每个行业拿到的年度总量，得到一个年度帕累托前沿。第二级在更新周期的尺度上，把年度量拆到各个时段，用“泄流—径流—供水”的链条做细化：上游水库按调度曲线泄流，径流沿河道演进到下游水库，再按第一级给定的份额分配。

这一级拆分的代价是层级之间可能不一致：第一级给农业 40% 的年度份额，第二级在某个干旱时段因为来水不足只能兑现 30%，总效益会打折扣。工程上的处理是把第二级的可行性反馈成第一级的约束，或者接受一定比例的衰减。

论文报告的量化结果是：经济效益约 287.01 亿量级（按原文货币单位），上游水库需要 10 个月才能满足既定需求，下游水库需要 12 个月，期间有约 20 亿立方米的多余水量下泄。这组数字里最值得关注的不是效益值本身，而是“10 个月”和“12 个月”的差异——它说明在既定需求下，下游水库的补给链路比上游长两个月，这两个月就是整条链路的风险窗口。任何一个月的来水低于预期，压力都会集中在这个窗口里释放。缺水且水电需求无法满足时，通过压缩供给函数的取值区间来实现节流，这个操作在模型里就是把 Q 的上界往下调，重新求解，而不是修改目标函数——改目标会改变帕累托前沿的形状，改约束只改变可行域，后者更容易解释给决策者听。

5. 常见问题与排查：这类分层优化模型最容易翻车的五个地方

5.1 数据与建模层面的三个坑

转移矩阵出现全零行，稳态分布算出 NaN。现象是幂迭代跑了几步之后 π 变成 nan，或者某一行的和不为 1。原因几乎总是某个季节分段里“雨天→雨天”的样本数为零，归一化时除零。解决方式是加平滑先验，或者在分段之前先检查每个格子的样本量，样本量低于阈值就把相邻时段合并。血泪经验是：不要在归一化之后再检查，那时候 NaN 已经扩散到整条链路，回溯起来很费劲。

双对数回归 R² 很高，反解却在端点外崩掉。现象是把水位设到死水位附近，反解出来的库容变成负数或者暴涨几个数量级。原因是回归只在观测区间内有效，对数模型在区间外的行为完全不受数据约束。解决方式是对反解函数做 clip，把输出限制在 [死库容, 最大库容] 之间，或者干脆用单调分段线性插值替代回归——插值在端点外可以按端点斜率外推，行为可控得多。

单位混乱导致经济效益差三个数量级。现象是同一套代码换一批输入数据，经济效益从几十亿变成几万。原因是库容的单位在不同数据源里可能是立方米、亿立方米、英亩·英尺、百万英亩·英尺，水头可能是米或英尺，电价可能是每千瓦时或每兆瓦时。这类论文里也常见单位笔误，比如把水量的单位写成面积单位。解决方式是在读入层就统一转成 SI，模型内部只认一套单位，输出时再转回业务单位。

5.2 算法与结果层面的两个坑

帕累托前沿不收敛，或者全部挤在一端。现象是跑到 500 代，前沿上还是零星几个点，而且集中在一侧。常见原因是约束用罚函数处理，罚系数比目标函数量级大三个数量级，搜索被彻底压平；也可能是两个目标量级差太多，拥挤度计算被大目标主导。解决方式是改用可行性优先的约束支配原则（可行解无条件优于不可行解，可行解之间再比目标），并且在计算拥挤度之前对每个目标做归一化。

公平性目标被平均主义绑架。现象是前沿上社会效益最高的解，是把水按人口比例平均分，工业区缺水到停产。原因是满意度用了线性求和，边际效用恒定，最优解必然向“人多”的方向倾斜。解决方式是把满意度凹化（平方根、对数），或者把社会效益改成 max-min 形式，即最大化最差那个区的满意度，这在福利经济学里对应罗尔斯准则，和“保基本民生”的政策取向是一致的。

6. 设计取舍：为什么“分层”比“精细化”更值得押注

回头看整套模型，最核心的设计决策不是用了 NSGA-II，也不是用了马尔可夫链，而是选择了分层而不是扁平化。如果把这 30 个变量放到日尺度上做联合优化，理论上能得到更精确的解，实际上会得到三个后果：计算量涨到遗传算法难以承受的量级；帕累托前沿维度爆炸，决策者无从选择；任何一处数据误差都会通过上万维的耦合关系放大到全局。

分层的代价是层级之间可能不一致。第一级在年度尺度上承诺的分配方案，第二级在某个干旱时段兑现不了，这个缺口不会自动消失，必须靠回代校验来发现。一个可操作的验证方法是：把第一级输出的年度解作为第二级的输入，跑完之后统计总效益相对第一级估计值的衰减率。衰减率在 5% 以内，说明年度尺度的聚合没有丢失关键约束；超过 10%，就说明某一层的假设过强，需要把第二级的可行性约束回灌到第一级，或者干脆把两级的粒度调成一致。

另一个值得说的取舍是更新周期取众数而不是均值。用均值会让周期偏长，因为连续干旱段长度的分布是右偏的，长尾会把均值拉高，结果是模型对中等长度干旱段的响应变迟钝。取众数保证了模型在“最常发生的那种干旱”里保持最佳响应，代价是对极端长旱的响应不足——这个缺口靠单独设置应急触发条件来补，而不是靠拉长周期。

这套模型最值钱的地方，其实是那套“先定最坏情况、再做优化”的次序。绝大多数优化模型假定参数是确定的，把不确定性留给敏感性分析去事后处理；而这里把不确定性前置成了模型结构的一部分，更新周期本身就是不确定性的函数。这个思路可以迁移到很多资源分配问题上——库存补货周期、应急物资轮换周期、设备检修窗口，都可以用同样的逻辑从“最坏情况持续多久”倒推出来。

从那以后，我每次做带时间维度的调度模型，都会先花半天时间把“更新周期从哪来”这个问题单独回答一遍，而不是默认按月、按周切开。哪怕最后结论还是按月，至少知道这个“月”是怎么来的，出了问题该往哪里查。希望帮到你。

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