简介：这份2022年美赛B题M奖优秀论文围绕美国西部千年大旱下的水资源分配问题，给出了可复用的建模范式，适合备战美赛、学习多目标优化的本科生与研究生参考。论文提出FHSD公平分层供需模型：先用马尔科夫链分析历史气候数据，结合蒙特卡洛模拟求出连续干旱最大概率，以确定模型更新周期；再以一年为尺度构建供需函数，引入公平效益、价格、成本与满意度系数，并用NSGA2求解多目标非线性规划，在公平与经济效益间取得帕累托折中。压缩包仅含1个PDF文件，约1.16MB，即完整论文全文，含摘要、引言、模型假设、求解与灵敏度分析等章节，可据此复盘M奖论文的写作结构与建模链条。目前已有3065人学习下载。

1. 从一份拿奖的水资源分配建模方案说起

做水资源分配这类题，最容易翻车的地方不是算法写不出来，而是把年水量按比例一分，答辩时被问一句“凭什么这么分”就卡住了。公平分层供需模型（FHSD）的解法是把问题拆成两层：先算清楚模型该多久更新一次，再在年尺度和更新周期尺度上做多目标非线性规划，用 NSGA-II 求出一族帕累托解，而不是一个拍脑袋的比例。这套结构适合正在做水资源、能源调度、公共资源分配方向毕设或课设的人，也适合手里有一份数学建模论文、想把它改成可交付方案的从业者。接下来的顺序是周期确定、供需函数构建、多目标求解、避坑排查，最后落到验收和答辩。

2. 更新周期怎么定：马尔可夫链加蒙特卡洛把“多久重算一次”变成可算的

2.1 为什么更新周期要先定

水库调度的现实约束很硬：闸门动作成本高，下游用水户的生产计划不能天天改，一份分配方案下发之后总要维持一段时间。如果按天甚至按小时建模，变量数会随受水区、行业、时段的组合成倍膨胀，求解时间和可解释性一起崩。反过来，周期定得太长，遇到连续干旱就会失真，方案还没执行完就已经不适用。

通常的做法是把“更新周期”当成一个待定参数，用历史上最可能出现的连续干旱天数来标定它。逻辑是：如果两次决策之间真的出现了比这更长的连续无雨期，那属于极端情景，本来就该触发人工干预；而在绝大多数月份里，模型不会因为低估干旱而被推翻。这是一种偏保守的鲁棒性设计，比直接拍一个“每月一次”要经得起追问。上游那段主库记作 P 库，下游那段记作 M 库，两库的入流和蓄量状态不同，更新周期也要分开算。

2.2 马尔可夫链估计降雨转移概率

把每天的天气压成两个状态：无雨（干旱）和有雨。假设明天状态只和今天有关，用历史日序列统计转移频次，得到 2×2 的转移矩阵。这个假设在气象学上当然不严格，天气有持续性也有季节循环，但作为“更新周期”的输入够用，因为我们只用它估连续干旱长度的分布，不拿它做天气预报。

import numpy as np def estimate_transition(series, alpha=1.0): """用历史日天气序列估计一阶马尔可夫转移矩阵。 series: 元素为 0(干旱) / 1(有雨) 的序列，建议按季节切分后分别调用 alpha: 拉普拉斯平滑系数，防止某个转移对样本为 0 时概率被估成 0 """ counts = np.full((2, 2), alpha) for prev, cur in zip(series[:-1], series[1:]): counts[prev, cur] += 1 return counts / counts.sum(axis=1, keepdims=True)

返回矩阵的第 i 行第 j 列，含义是“今天处于状态 i，明天转到状态 j”的条件概率。平滑系数一般取 1，历史序列上万天时可以取 0.5 让数据说话。这里有个容易被忽略的细节：一定要按季节或月份分段估计。把夏季和冬季的数据混在一起，估出来的矩阵两边都不像，后面蒙特卡洛的结果自然偏。

估完先做一次一致性检查：解平稳分布 π（满足 πP = π 且各分量和为 1），看 π 里的干旱占比和样本里实际的干旱天数占比差多少。差得远说明状态划分或分段方式有问题，别急着往下走。这一步花不了十分钟，但能省掉后面返工的两天。

2.3 蒙特卡洛跑出连续干旱天数分布

有了转移矩阵就可以模拟。每次从平稳分布抽一个初始状态，然后按矩阵逐日抽样，走完一个月的窗口，记录这个窗口里最长的一段连续干旱天数。重复十万次，得到游程长度的经验分布，取分布里最有代表性的那个值作为更新周期。

def simulate_dry_runs(trans, days=30, trials=100000, seed=42): """模拟一个月窗口内最长连续干旱天数。 trans: 2x2 转移矩阵，行是今天、列是明天 days: 仿真窗口，取一个月的天数 trials: 重复次数，决定经验分布的平滑程度 """ rng = np.random.default_rng(seed) # 初始状态按平稳分布抽样，避免固定初值带来的偏差 eigvals, eigvecs = np.linalg.eig(trans.T) pi = np.real(eigvecs[:, np.argmax(np.real(eigvals))]) pi = pi / pi.sum() runs = np.empty(trials, dtype=int) for t in range(trials): state = rng.choice(2, p=pi) best = cur = 0 for _ in range(days): state = rng.choice(2, p=trans[state]) cur = cur + 1 if state == 0 else 0 # 0 表示干旱 best = max(best, cur) runs[t] = best return runs

参数怎么调：days 跟着决策窗口走，题目按“一个月内”考虑就取 30；trials 十万次在普通笔记本上几十秒能跑完，想更快可以把两层循环改成批量抽样。种子必须固定，否则答辩现场复现不出来，这是血泪经验。

从游程数组里取代表值，最直觉的是取众数，但如果分布很平（比如 9 天占 15%、10 天也占 15%），众数就不稳。我一般取“累计概率首次超过 50%”的那个天数，比单纯取众数更接近业务上“多数情况下够用”的语义，也更好解释。

2.4 结果落地：分时段更新时长表

按季节分段估计转移矩阵后，两库在不同时段的更新周期差别很明显：

这张表的用法不是把数字抄进论文就完事，而是要能解释清楚 M 库在 9—12 月为什么跳到 17 天。该时段降雨概率只有 4.1%，连续干旱几乎成为常态，最长游程的分布被压到窗口上限附近，代表值自然偏大。换句话说，这段时间模型更新得“懒”一点可以接受，因为天气本来就不太给意外；反过来 1—5 月两库都在 10 天，说明这段时间天气切换频繁，更新要勤。

工程上还有一个用法：把更新周期和库容状态挂钩。水位接近死水位时，不管周期到没到都强制重算一次。周期是常态下的节奏，不是免检牌照，这一点在写模型假设时要单独说明，否则容易被理解成“干旱期反而不更新”的逻辑漏洞。

3. 供需函数怎么建：水位库容回归与四系数分层需求

3.1 供给端：水位—库容对数线性回归

水库的“可供水量”不是固定数，它随水位变化。要把它写进优化模型，先得有一条水位到库容的连续关系。常见做法是查两库的历史水位—蓄量对应表，然后拟合。

直接把水位对库容做线性回归是不行的：低水位段库容变一点水位就掉很多，高水位段又变得很平，线性拟合在两端都系统性偏。取双对数之后，两段的曲率被拉平，残差分布更均匀，外推到接近死水位时也不容易出现负库容这种荒谬结果。

import numpy as np def fit_level_storage(levels, storages): """水位-库容对数线性回归：ln(H) = b * ln(V) + a levels: 水位序列，单位 ft（above MSL） storages: 对应蓄量序列，单位 maf 返回斜率 b、截距 a、R2，以及由水位反解库容的函数 """ x = np.log(np.asarray(storages, dtype=float)) y = np.log(np.asarray(levels, dtype=float)) b, a = np.polyfit(x, y, 1) # numpy 返回顺序是 [斜率, 截距] y_hat = b * x + a r2 = 1 - np.sum((y - y_hat) ** 2) / np.sum((y - y.mean()) ** 2) def storage_from_level(h): return np.exp((np.log(h) - a) / b) return b, a, r2, storage_from_level

参数说明：b 是弹性系数，物理含义是“水位每变化 1%，库容变化 1/b%”，两库的 b 值不同，不能互相套用。R² 用来判断拟合质量，两条回归的 R² 都在 0.99 以上、回归系数 p 值远小于 0.01，说明关系显著。反解函数 storage_from_level 是后面所有约束的入口——优化模型里的约束都是对库容写的，而调度员看的是水位，两边靠这条函数互转。

注意截断：回归只在有观测数据的区间可信，外推到死水位以下必须夹住，写成 max(死库容, 反解值) 再和最大蓄量取 min。三个特征水位要分清楚：最大蓄水位决定上限，最小发电水位决定还能不能开机，死水位决定还能不能放水。很多模型把这三个混成一个“最低水位”，结果发电约束和供水约束打架，解出来的方案根本没法执行。

3.2 需求端：四个系数把公平和效益写进函数

需求侧不能只写一个“总需水量”，那样优化出来一定把水全给产值最高的行业。要刻画公平，得把需求拆到“受水区 × 行业”的粒度，再给每个格点配上四个系数。

公平效益系数：每个受水区一个权重，综合历史用水量、人口、耕地面积得到，体现基本生活用水不能砍。价格系数：不同行业水价不同，工业高于居民，居民高于农业，这是效益目标的来源。成本系数：由输水距离和扬程决定，距离越远、提水越高，单位成本越大。满意度：实际分配量除以需求量，截断到 0 到 1 之间。

把满意度做成带指数的方式进入目标函数，是这套模型比较关键的一步：

import numpy as np def social_welfare(x, d, w, theta=0.6): """分层需求的社会效益目标。 x: 分配量矩阵，形状 (受水区数, 行业数) d: 需求量矩阵，同形状 w: 各受水区公平效益系数，长度等于受水区数 theta: 公平偏好参数，越小越倾向均等分配 """ u = np.clip(x / np.maximum(d, 1e-9), 0.0, 1.0) # 满意度 return float(np.sum(w[:, None] * (u ** theta) / theta))

theta 的作用要讲清楚：小于 1 时满意度函数的边际收益递减，把一单位水从满意度 0.9 的用户挪给满意度 0.3 的用户，总福利上升，算法自然把水往缺水方倾斜；theta 等于 1 退化成线性加权，公平性最弱；theta 继续往 0 走，分配越来越平均，但经济效益掉得也快。答辩时最好给出 0.3、0.6、1.0 三档对比，说明为什么选中间那档，而不是只报一个数。

“分层”体现在这里：先保居民基本用水，再保农业，最后才是工业增量。实现方式不是把顺序硬编码进约束，而是给不同层级设不同的 theta 和权重下限，让优化器在满足硬约束的前提下自己找平衡。这样模型既可解释，又不会因为写死而失去调整空间。

3.3 水电出力与“水—电”耦合

水电收益是这套模型的第二个经济来源，出力公式本身不复杂：

RHO, G = 1000.0, 9.81 def hydro_power(q, head, eta=0.85): """q: 发电流量，m3/s；head: 净水头，m；eta: 机组综合效率 返回功率，单位 kW """ return eta * RHO * G * q * head / 1000.0

q 受最小发电流量和最大过机流量两头夹，head 是上下游水位差减去水头损失。电价按题给假设取 0.1041 美元每千瓦时，这个数字必须写进假设清单，因为它直接决定“发电”和“供水”谁更值钱，答辩一定会追问来源和时效性。

真正容易出问题的是耦合关系：发电用的水和供给农业、工业、居民的水是同一股水。取水口在上游还是下游，决定了这股水能不能被算两次。合理的处理是先发电、尾水回归河道，下游取水只计一次，同时保证河道内还剩生态基流。如果模型里两边都按满额计收益，水量平衡一定对不上，效益还会虚高一大截，这种错误在初稿里非常常见。

4. 多目标求解：NSGA-II 在年尺度定总量、在更新周期尺度分水量

4.1 决策变量与约束怎么定

变量设计上，年尺度和更新周期尺度要分开。年尺度上，决策变量是每个受水区每个行业的年分配量，加上两库的年放水量，5 个受水区 × 3 个行业 = 15 维，再加 2 维放水量，一共 17 维。这个维度对 NSGA-II 来说很轻，种群 100 到 200 就够用。

约束有这几类：水量平衡（年初蓄量 + 区间入流 − 蒸发渗漏 − 放水 = 年末蓄量）；库容上下界（由 3.1 的反解函数把水位约束翻译过来）；发电约束（放水量不低于最小发电流量折算的时段水量）；需求上限（每个格点的分配量不超过它的需求量）；非负与总量（所有分配量之和不超过可用水量）。

单位时间尺度上，把年分配量按更新周期切成若干段，决策变量变成“每段每个格点的分配量”，约束换成当段来水。两层之间靠年结果做边界：年总量是上界，周期内分配之和不能超过它。这样拆的好处是既避免了按日建模的维度爆炸，又保留了季节来水差异带来的细节。

4.2 目标函数：经济效益与社会效益

两个目标是有冲突的，这是模型的核心张力。经济效益 = 售水收入 + 发电收入 − 输水成本，售水收入按行业水价乘分配量，发电收入按出力乘电价再乘时间，输水成本按成本系数乘距离乘水量。社会效益就是 3.2 里的满意度加权函数，它不看钱，只看各格点的缺水程度，天然偏向把水给需求大、当前满意度低的一方。

如果只做一个目标，用加权求和也能解，但要试很多组权重，而且权重怎么定说不清。多目标的价值在于一次给出一族帕累托解，让决策者在前沿上看“多花一亿成本能换来多少公平”，这个权衡过程本身就是答辩的加分项。反过来说，如果答辩时你只能拿出一个解，评委很难判断你是真的权衡过，还是随手调了个权重。

4.3 用 NSGA-II 求帕累托前沿

NSGA-II 的用法不复杂，关键是目标和约束写对。

import numpy as np from pymoo.core.problem import Problem from pymoo.algorithms.moo.nsga2 import NSGA2 from pymoo.optimize import minimize class WaterAllocation(Problem): def __init__(self, demand, price, cost, w, theta, avail, eco_flow): self.d, self.p, self.c = demand, price, cost self.w, self.theta = w, theta self.avail, self.eco_flow = avail, eco_flow n = demand.size super().__init__(n_var=n, n_obj=2, n_ieq_constr=2, xl=np.zeros(n), xu=demand.flatten()) def _evaluate(self, X, out, *args, **kwargs): X3 = X.reshape(-1, *self.d.shape) # 目标一：经济效益，取负值以适配 pymoo 的最小化约定 revenue = np.sum(self.p[None] * X3) - np.sum(self.c[None] * X3 <p> <a href="https://download.csdn.net/download/nvlnlf37/86729634" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>