简介《2022年美赛B题M奖优秀论文》面向备战美国大学生数学建模竞赛、关注水资源优化配置与多目标建模的读者完整呈现一道美赛B题M奖论文的写作范式与建模思路。资源为1个PDF文件压缩包约1.16MB篇幅紧凑便于在电脑或平板上直接阅读与批注。论文围绕美国西部干旱危机构建FHSD公平分层供需模型结合马尔科夫链与蒙特卡洛模拟确定模型更新周期并以多目标非线性规划刻画公平性与经济效益采用NSGA2算法求解帕累托解集读者可从中学习赛题拆解、供应与需求函数设计、约束设置、敏感性分析及论文摘要与目录组织方式。目前已有3065人学习下载适合需要参考获奖论文结构、积累优化建模表达与复盘美赛B题解法的参赛者和指导教师。1. 从一份沿用百年的分水协议说起两座梯级水库给下游五个州同时供农业、工业、居民用水还要靠落差发电这是 2022 年美赛 B 题的题面也是这份 M 奖论文要回答的事。多数人的第一版做法是把年可用水量按各州需求比例一分了事数字看着漂亮可只要遇上连续枯水年水位跌破最低发电水位下游连居民用水都得限供模型可能直接失效。这套 FHSD 方案Fair Hierarchical Supply and Demand把问题拆成三层先用马尔可夫链加蒙特卡洛确定模型多久更新一次再用水位—库容回归和分层需求函数搭出供需两侧最后用多目标非线性规划配合 NSGA-II 求一组帕累托解。它适合正在做水资源调度、多目标优化类毕设或课设的人也适合想把“公平”这个偏软的指标写进数学模型的从业者下面按这条主线把每一步的实现和验收都摊开。2. 更新周期建模为什么用连续干旱天数而不是年均降雨水库调度成本高、灵活性低模型多久更新一次直接决定方案的稳健性。用年均降雨概率去推更新周期看着省事但年均值把丰枯波动抹平了方案在枯水年几乎必然出问题。这篇论文换了个思路从五州历史日尺度气象数据里统计“连续干旱天数”取出现概率最大的那个 T 作为更新间隔相当于按最坏情形设计模型宁可保守一点。2.1 状态定义与转移矩阵的统计方法把天气压成两个状态就够用0 表示无雨1 表示有雨。无后效性是这套方法的全部前提今天下不下雨只跟昨天有关跟更早的历史无关。实际数据往往不满足严格的无后效性但日尺度天气的马尔可夫性质在工程上已经被验证过足够好用作为更新周期的粗估工具完全够。import numpy as np def estimate_transition(seq: np.ndarray) - np.ndarray: seq: 日尺度天气序列0 无雨(dry)1 有雨(rain) 返回 2x2 转移矩阵 PP[a][b] P(明天b | 今天a) P np.zeros((2, 2), dtypefloat) for a, b in zip(seq[:-1], seq[1:]): P[a, b] 1.0 # 按行归一化某状态样本过少时用拉普拉斯平滑避免出现 0/0 row_sum P.sum(axis1, keepdimsTrue) row_sum[row_sum 0] 1.0 return P / row_sum这段代码只做统计不引入任何主观假设。seq是清洗后的日序列缺失值先剔除再统计P[a][b]的含义是“今天状态 a、明天状态 b”的频率估计。行归一化这一步不能省否则后续蒙特卡洛里pP[state]传进去的不是概率分布采样会直接抛错。样本量小于 200 天时建议做拉普拉斯平滑分子加 1、分母加状态数代价是概率被轻微拉平换来的是数值稳定。平稳分布可以顺手算出来它对应长期的干旱占比是后面敏感性分析的基准输入def stationary(P, tol1e-10, max_iter10000): 幂迭代求平稳分布返回 [P(dry), P(rain)] pi np.array([0.5, 0.5]) for _ in range(max_iter): nxt pi P if np.abs(nxt - pi).sum() tol: return nxt pi nxt return pitol控制收敛精度日尺度天气矩阵一般几十次迭代就收敛max_iter只是兜底真跑到上限说明矩阵可能有多解先回去检查是不是有状态根本没出现过。2.2 蒙特卡洛把一个月的最长连续干旱天数分布跑出来有了转移矩阵接下来模拟一个月 30 天的天气路径统计每条路径上“最长连续干旱天数”重复两万次得到一张分布表。取众数作为更新周期 T取分位数作为保守档两种用法都合理看你要多稳。def longest_dry_run(seq): 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, n_days30, n_paths20000, seed20220201): 返回长度 n_days1 的数组dist[t] 最长连续干旱恰好为 t 天的概率 rng np.random.default_rng(seed) counter np.zeros(n_days 1, dtypenp.int64) for _ in range(n_paths): state 0 if rng.random() 0.9 else 1 # 初始状态严谨做法按平稳分布采样 seq [] for _ in range(n_days): state rng.choice(2, pP[state]) seq.append(state) counter[longest_dry_run(seq)] 1 return counter / n_pathsn_paths决定蒙特卡洛噪声大小两万次时众数的抖动基本在 1 天以内seed固定住是为了结果可复现答辩时能当场重跑。初始状态那一行是很多人的失分点直接用固定值会让前几天的路径分布偏离平稳态正确做法是按stationary(P)采样路径短于 10 天时这个偏差尤其明显。论文把全年切成三段分别统计得到的更新周期如下表时段P 库无雨概率P 库有雨概率P 库更新周期 TM 库无雨概率M 库有雨概率M 库更新周期 T1–5 月0.89260.107410 天0.91950.080510 天6–8 月0.87910.120911 天0.89010.109910 天9–12 月0.91840.081611 天0.95900.041017 天2.3 分段结果的读法与异常值处理这张表里最值得警惕的是 M 库 9–12 月的 17 天。同一座库在其他两段都是 10 天上下突然跳到 17 天通常有两种解释一是这段历史样本里真的出现过长时间连续无雨二是样本量偏小导致分布右偏。我一般的处理是两条腿走路先看这段历史序列长度是否与其他段持平样本量不足就合并相邻时段重新统计样本量够就保留 17 天但要在模型里注明“该时段按保守档设计”别让一个极端的 T 把整个调度频率拖成季度级。分段还有一个现实好处农业用水集中在 4–5 月的融雪期6–8 月是发电高峰9–12 月进入枯水期三段的用水结构差异很大。按段设置更新周期等于让模型在丰水期更勤快地重算在枯水期更谨慎地少动这个节奏比全年一个 T 要贴近真实调度。3. 供需两侧分层建模从水位—库容曲线到四系数需求函数3.1 水位与库容为什么用双对数线性回归题给数据是水位英尺海拔和库容百万英亩-英尺两者的关系不是直线直接线性拟合的残差会随水位升高越来越大外推到高水位段基本没法用。双对数能把它压回一条直线形式是 ln S a·ln H b物理上对应幂律关系拟合出来再做指数还原。import numpy as np from scipy import stats def fit_level_storage(level, storage): level: 水位(ft above MSL)storage: 库容(maf) 返回斜率 a、截距 b、R²、p 值 ln_h np.log(np.asarray(level, dtypefloat)) ln_s np.log(np.asarray(storage, dtypefloat)) res stats.linregress(ln_h, ln_s) return {a: res.slope, b: res.intercept, r2: res.rvalue ** 2, p: res.pvalue} def storage_from_level(level, a, b): 反解给定水位返回库容 return np.exp(b) * np.power(level, a)a是弹性系数含义是水位上涨 1% 对应的库容上涨百分比两座库的a都应落在 1 到 3 之间偏离太多说明数据里有异常点先查原始表。r2是拟合优度p是回归系数的显著性检验论文里两座库的 p 值都远小于 0.01意味着幂律关系在 99% 置信度下显著不为零这个结论是后面所有库容推算的地基。storage_from_level里的np.exp(b)是幂律的前置常数别把指数写反。3.2 水位分档与水电出力的接线方式水位不是一个连续变量就能打发的它有三条硬线最大库容对应水位、最低发电水位、死水位。三条线把水位切成四档每档的供水和发电策略完全不同水位区间供水策略发电策略高于最大库容对应水位弃水或跨库调蓄满发最低发电水位 ~ 最大库容正常供水按落差与下泄流量发电死水位 ~ 最低发电水位限量供水优先居民不发电低于死水位停止供水不发电分档之后水电出力才有意义否则会出现“水位已经跌到死水位以下还在算发电收入”这种明显不合理的解。出力的物理公式是三段式流量、落差、效率。def hydro_power(flow_m3s, head_m, eta0.88, rho1000.0, g9.81): flow_m3s: 下泄流量 m³/s head_m: 上下游水位差 m eta: 综合效率0~1大型水电站取 0.85~0.90 返回功率 kW return rho * g * flow_m3s * head_m * eta / 1000.0rho是水密度g是重力加速度eta把水轮机效率、发电机效率、管路损失一次性打包取值范围 0.85 到 0.90 之间比较稳妥。head_m要用上游库水位减去下游尾水位不是减死水位这个减错的人不少结果就是发电量偏高 20% 以上。接线顺序是水位 → 库容 → 年可调蓄量 → 下泄流量 → 功率 → 发电收入中间任何一环用了错的水位后面全错。3.3 需求函数的四个系数怎么落到参数上需求侧不是把各州历史用水量直接抄进去那样等于假设需求刚性价格、成本、公平性全都进不了模型。论文的做法是把需求写成弹性形式用四个系数调公平-效益系数alpha、价格系数price、成本系数cost、满意度s。import numpy as np def demand_matrix(D0, alpha, price, price_ref, eps): D0: (n_state, n_sector) 基准需求量 alpha: (n_state,) 公平-效益系数反映各州在公平原则下的权重 price: (n_sector,) 各行业水价系数 eps: (n_state, n_sector) 需求价格弹性 返回调整后的需求矩阵 ratio (price / price_ref)[None, :] ** (-eps) return alpha[:, None] * D0 * ratioalpha用各州人口占比和历史用水量占比的加权平均来定取值在 0.5 到 1.5 之间超过这个范围说明权重设计失衡要么有州被过度照顾要么被过度压制。price_ref是基准水价用全流域平均水价即可作用是让ratio无量纲化。eps是需求价格弹性农业通常低于工业居民用水最低量级上农业 0.2–0.4、工业 0.5–0.8、居民 0.1–0.3这个顺序不能反反了模型会得出“涨价后农业用水反而增加”这种反直觉结论。成本系数cost来自“水库到州的距离 × 单位输水成本”距离用地理直线距离做近似单位输水成本按行业分档农业低、工业高。满意度s_i min(1, x_i / D_i)是供水量与需求量之比上限截断在 1否则超额供水会被算成“满意度 1.5”把公平性指标带偏。3.4 分层结构怎么体现在约束里分层指的是“年总量 → 更新周期内总量 → 单州单行业配额”三层往下拆。年层确定总可供水量是硬约束周期层把年量按 2.3 节算出的 T 切成若干段每段受当年水文情势约束州层再按需求函数分配。三层之间的传递关系是上层出总量、下层出结构不能反过来反过来就变成“各州先报需求、总量再兜底”遇到枯水年直接无解。4. 多目标非线性规划两个目标、三类约束与 NSGA-II 的完整程序4.1 经济效益与公平性怎么同时写进目标函数经济效益是收入减成本对每个州、每个行业求和公平性用“最差州的满意度”来衡量最大化最小值也就是 maximin。两个目标天然冲突把钱全砸给高水价行业经济效益上去了农业州的满意度会掉到地板平均主义地分公平性好看经济效益又难看。import numpy as np from pymoo.core.problem import ElementwiseProblem class Alloc(ElementwiseProblem): def __init__(self, D, price, cost, W_year): self.D, self.price, self.cost D, price, cost self.W_year W_year self.ns, self.nk D.shape super().__init__( n_varself.ns * self.nk, n_obj2, n_ieq_constrself.ns 1, xl0.0, xu1.0, ) def _evaluate(self, x, out, *args, **kwargs): X x.reshape(self.ns, self.nk) * self.D # 决策变量缩放到实际水量 profit np.sum((self.price[None, :] - self.cost) * X) s X.sum(axis1) / self.D.sum(axis1) # 各州满意度 out[F] [-profit, -np.min(s)] # pymoo 默认最小化 g_total X.sum() - self.W_year # 年总量约束 g_state X.sum(axis1) - self.D.sum(axis1) # 单州不超需求 out[G] np.concatenate(([g_total], g_state))决策变量x归一化到 [0,1]乘上需求矩阵D之后才是实际水量这样做的好处是搜索空间规整交叉变异不会把某个州的水量顶到天上去。out[F]里两个目标都取负号因为 pymoo 统一按最小化处理取负等价于最大化收益和最大化最差满意度。out[G]是三个约束的合并第一个是年总量不超过可供水量后面ns个是单州供水量不超过自身需求全部写成g(x) 0的形式。4.2 约束怎么写才不会让解全落在不可行域约束写法有两个常见问题。一是把库容下限写成等式水位必须正好等于最低发电水位这种约束在连续优化里几乎不可能满足解会全部被判为不可行正确做法是写成不等式水位 - 最低发电水位 0允许有裕量。二是把“不发电”这种离散决策塞进连续变量里NSGA-II 处理不了得用分档函数近似或者在评估阶段根据水位直接截断出力。年总量约束W_year的来源是供水函数两座库的年初库容加上年内来水减去蒸发和生态基流剩下的才是可分配量。这个值在敏感性分析里会上下浮动W_year一变帕累托前沿整体平移所以它不能写死在代码里要做成可配置参数。4.3 NSGA-II 的参数选择与帕累托解挑选算法参数没有普适最优但有几条经验线种群规模不低于决策变量数的 10 倍这里 15 个变量取 120 到 150 之间代数看前沿是否稳定跑到 400 代左右前沿形状基本不再变化就可以停。from pymoo.algorithms.moo.nsga2 import NSGA2 from pymoo.optimize import minimize algorithm NSGA2(pop_size120, eliminate_duplicatesTrue) res minimize(Alloc(D, price, cost, W_year), algorithm, (n_gen, 400), seed1, verboseFalse) F, X res.F, res.Xpop_size太小前沿会稀疏太大收敛慢120 是 15 维问题比较合适的折中eliminate_duplicatesTrue去掉重复个体避免种群早熟seed固定住答辩现场重跑能复现同一张前沿图。res.F是目标值矩阵res.X是对应的决策变量两者行对齐。从帕累托前沿里挑最终解常用两种办法。拐点法找前沿上曲率最大的点对应“再让一点公平要付出很多经济代价”的转折TOPSIS 按与理想解的相对距离排序适合有明确偏好的场景。我一般两个都跑一遍如果两种方法选出的解落在同一条解附近说明这个解比较稳健如果分歧很大就把两个解都拿出来做敏感性对比看哪个在参数扰动下更稳。论文最终给出的结果量级是经济效益约 287.01 亿美元按一年期供给P 库需要 10 个月满足既定需求M 库需要 12 个月剩余约 20 亿立方米的余量最终流入下游湾口。当水资源紧缺、水电需求无法满足时通过收紧供水函数的取值区间来实现节流这一步本质上是把“不发电”档的水位区间扩大让模型自己把水留住。5. 避坑与常见问题五个最容易让方案失真的地方5.1 单位不统一目标函数量级离谱现象经济效益算出来是 10¹² 或者 10³ 这种量级跟论文里的百亿级对不上帕累托前沿两个目标的数值差十几个数量级图根本画不出来。原因水位用英尺、库容用百万英亩-英尺、功率用千瓦、电价用美分每千瓦时四套单位混在一个目标函数里量纲没有抵消。解决建模前先定一套单位制全部换成 SI转换只在数据入口做一次目标函数内部不再出现任何转换系数。具体就是英尺乘 0.3048 变米百万英亩-英尺乘 1.2335×10⁶ 变立方米美分每千瓦时除 100 再乘汇率口径变元每千瓦时。转换表贴在代码注释里谁接手都能核对。5.2 更新周期用了年均值现象模型在验证集上表现不错一放到枯水年情景就断供水位直接击穿死水位方案失效。原因用全年平均降雨概率推 T丰枯波动被抹平算出来的更新周期在枯水段偏长模型反应迟钝。解决按 2.3 节那样分段统计至少切成丰水、平水、枯水三段每段单独算转移矩阵和 T。如果数据量只够分两段就按汛期和非汛期切别硬凑三段。5.3 满意度直接比绝对供水量现象优化结果里大州拿走了七成水量小州满意度接近零公平性指标却没报警。原因公平性直接用供水量x_i比较没做归一化需求基数大的州天然占优。解决一律用比例s_i min(1, x_i / D_i)先除需求再比。这个改动看起来小但它把“公平”从绝对量拉回到相对量是整套模型里最容易被忽略、也最影响结论的一行代码。5.4 提前加权把多目标压成单目标现象帕累托前沿退化成一个点或者是一条斜率固定的直线NSGA-II 跑出来的解全挤在一起。原因在目标函数里直接写w1*profit w2*fairness权重提前定死支配关系被破坏。解决目标函数保持两个独立分量权重留到选解阶段再用。真要加权也得跑一组权重扫描看前沿形状随权重怎么变而不是拍一个 0.6/0.4 就完事。5.5 罚函数权重量级拍脑袋现象所有解都落在不可行域或者全部贴在约束边界上可行解一个都找不到。原因把硬约束写成了罚项罚项系数跟目标函数量级不匹配要么罚得太轻等于没约束要么罚得太重把搜索空间压扁。解决能写硬约束就别用罚函数pymoo 的n_ieq_constr支持不等式约束直接写进去最干净。非要用罚函数罚项系数取目标函数典型值的 1 到 10 倍跑之前先用随机解试一遍看违约解的适应度是否明显劣于可行解。6. 验收与答辩怎么证明这组解不是凑出来的验收的核心不是“结果好看”而是“结果经得起问”。我一般把验收拆成三组检查每组都有明确的通过标准任何一组不过就回去改模型而不是改结论。检查项通过标准不通过时先看哪里回归拟合R² 0.95p 0.01原始数据是否有异常点、是否该分段拟合帕累托前沿解数量 ≥ 50分布均匀无断档种群规模是否过小、代数是否不足敏感性三档情景下最优解排序不变参数扰动幅度是否超出合理区间敏感性分析至少跑三档农业用水量 ±10%、可再生能源占比 ±20%、总需求 -10%。论文里做的正是这三组看的是需求侧变化怎么传导到供给侧。判断标准不是“结果完全不变”而是“解的相对排序不变”排序变了说明模型对某个参数过度敏感得回去检查那个参数是不是被当成了硬编码。答辩现场五分钟起手别讲模型多复杂先讲最坏情形两座库在连续干旱 17 天的情况下会怎样把评委的注意力拉到问题的难点上。然后过渡到两层时间尺度年层定总量、周期层定节奏一句话说清为什么不用日尺度因为日尺度会让模型维度爆掉、求解时间从分钟级涨到小时级。最后落到帕累托前沿那张图上指着拐点说“这个解是公平性和经济效益的平衡点”比念一堆公式有用得多。从那以后我每次做多目标优化题都强制自己先跑一遍单位检查再跑模型因为见过太多方案栽在英尺和米的换算上明明思路是对的结果一公布就被问住。希望这套结构能帮到你把方案讲成一个能被追问、也能被复现的完整东西。本文还有配套的精品资源点击获取
