基于马尔可夫链与NSGA-II的水资源公平分层供需建模

简介这份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, alpha1.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(axis1, keepdimsTrue)返回矩阵的第 i 行第 j 列含义是“今天处于状态 i明天转到状态 j”的条件概率。平滑系数一般取 1历史序列上万天时可以取 0.5 让数据说话。这里有个容易被忽略的细节一定要按季节或月份分段估计。把夏季和冬季的数据混在一起估出来的矩阵两边都不像后面蒙特卡洛的结果自然偏。估完先做一次一致性检查解平稳分布 π满足 πP π 且各分量和为 1看 π 里的干旱占比和样本里实际的干旱天数占比差多少。差得远说明状态划分或分段方式有问题别急着往下走。这一步花不了十分钟但能省掉后面返工的两天。2.3 蒙特卡洛跑出连续干旱天数分布有了转移矩阵就可以模拟。每次从平稳分布抽一个初始状态然后按矩阵逐日抽样走完一个月的窗口记录这个窗口里最长的一段连续干旱天数。重复十万次得到游程长度的经验分布取分布里最有代表性的那个值作为更新周期。def simulate_dry_runs(trans, days30, trials100000, seed42): 模拟一个月窗口内最长连续干旱天数。 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, dtypeint) for t in range(trials): state rng.choice(2, ppi) best cur 0 for _ in range(days): state rng.choice(2, ptrans[state]) cur cur 1 if state 0 else 0 # 0 表示干旱 best max(best, cur) runs[t] best return runs参数怎么调days 跟着决策窗口走题目按“一个月内”考虑就取 30trials 十万次在普通笔记本上几十秒能跑完想更快可以把两层循环改成批量抽样。种子必须固定否则答辩现场复现不出来这是血泪经验。从游程数组里取代表值最直觉的是取众数但如果分布很平比如 9 天占 15%、10 天也占 15%众数就不稳。我一般取“累计概率首次超过 50%”的那个天数比单纯取众数更接近业务上“多数情况下够用”的语义也更好解释。2.4 结果落地分时段更新时长表按季节分段估计转移矩阵后两库在不同时段的更新周期差别很明显时段P 库干旱概率P 库降雨概率P 库更新天数M 库干旱概率M 库降雨概率M 库更新天数1—5 月0.89260.1074100.91950.0805106—8 月0.87910.1209110.89010.1099109—12 月0.91840.0816110.95900.041017这张表的用法不是把数字抄进论文就完事而是要能解释清楚 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: 水位序列单位 ftabove MSL storages: 对应蓄量序列单位 maf 返回斜率 b、截距 a、R2以及由水位反解库容的函数 x np.log(np.asarray(storages, dtypefloat)) y np.log(np.asarray(levels, dtypefloat)) 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, theta0.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, eta0.85): q: 发电流量m3/shead: 净水头meta: 机组综合效率 返回功率单位 kW return eta * RHO * G * q * head / 1000.0q 受最小发电流量和最大过机流量两头夹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_varn, n_obj2, n_ieq_constr2, xlnp.zeros(n), xudemand.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 hrefhttps://download.csdn.net/download/nvlnlf37/86729634 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p