简介面向数学建模竞赛选手与水资源优化研究者这份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_state2, alpha1.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(axis1, keepdimsTrue) 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, tol1e-13, max_iter100000): 幂迭代求平稳分布避免直接做特征分解带来的数值问题 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, days30, n200000): 模拟 n 个月返回最长连续干旱天数的经验分布 hist np.zeros(days 1, dtypenp.int64) for _ in range(n): seq np.empty(days, dtypenp.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 万次在普通笔记本上是几秒的量级代价可以接受。这里还必须说明取众数而非期望是因为更新周期要的是一个“最可能出现的整段干旱长度”期望会被长尾拉高导致更新周期定得过长反而在最常见的中等干旱段里失去响应能力。分季节建模之后两个水库、三个时段各得一组结果形态大致如下时段水库 P 干旱概率水库 P 更新周期天水库 M 干旱概率水库 M 更新周期天1–5 月0.8926100.9195106–8 月0.8791110.8901109–12 月0.9184110.959017表格里最值得琢磨的是下游水库在 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, k9.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 分层需求四个系数如何把“公平”变成可计算的量需求侧不能简单地写成“每个区需求量是常数”。真实的需求随价格、输水成本、缺水程度变化而且不同行业对缺水的敏感度差异巨大。分层结构就是在这里起作用的农业、工业、居民三层每层单独建需求函数层与层之间的权重由政策偏好决定。四个系数各有分工系数符号作用取值逻辑公平-效益系数α决定该区在该行业上的分配优先级人口、产值、历史用水权综合价格系数p反映需求对水价的弹性行业平均水价工业最高、农业最低成本系数c输水距离决定的单位成本与水库到该区的距离成正比满意度sat供给/需求的凹化函数越缺水边际痛感越强满意度这一步是整套需求建模里最容易被做坏的地方。如果直接用线性的供给需求比那么总社会效益就是各区的加权和优化器的必然选择是把水集中给“权重高、人口多”的区其余区被压到零这显然不是公平。用凹函数处理就能改变这个结论def satisfaction(supply, demand, eps1e-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 是目标数这里是 2N 是种群规模。种群取 200 时一次排序大约几万次比较完全不是瓶颈。拥挤度距离负责在同一层内部排序。做法是按每个目标分别排序边界解的距离设为无穷大保证极端解一定被保留中间解的距离取相邻两解在该目标上的差值再归一化对所有目标求和。def crowding_distance(front, objs): 同一非支配层内的拥挤度值越大说明周围越稀疏越该保留 dist {i: 0.0 for i in front} for k in range(objs.shape[1]): order sorted(front, keylambda 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%就说明某一层的假设过强需要把第二级的可行性约束回灌到第一级或者干脆把两级的粒度调成一致。另一个值得说的取舍是更新周期取众数而不是均值。用均值会让周期偏长因为连续干旱段长度的分布是右偏的长尾会把均值拉高结果是模型对中等长度干旱段的响应变迟钝。取众数保证了模型在“最常发生的那种干旱”里保持最佳响应代价是对极端长旱的响应不足——这个缺口靠单独设置应急触发条件来补而不是靠拉长周期。这套模型最值钱的地方其实是那套“先定最坏情况、再做优化”的次序。绝大多数优化模型假定参数是确定的把不确定性留给敏感性分析去事后处理而这里把不确定性前置成了模型结构的一部分更新周期本身就是不确定性的函数。这个思路可以迁移到很多资源分配问题上——库存补货周期、应急物资轮换周期、设备检修窗口都可以用同样的逻辑从“最坏情况持续多久”倒推出来。从那以后我每次做带时间维度的调度模型都会先花半天时间把“更新周期从哪来”这个问题单独回答一遍而不是默认按月、按周切开。哪怕最后结论还是按月至少知道这个“月”是怎么来的出了问题该往哪里查。希望帮到你。本文还有配套的精品资源点击获取
