1. 项目背景与核心挑战从“碳中和”到数学建模的跨越去年五一杯数学建模竞赛的C题把“碳中和”这个宏大的国家战略直接摆在了我们这些参赛者面前。说实话当时拿到赛题团队里先是兴奋紧接着就是一阵头大。兴奋的是这个题目太“潮”了紧扣时代脉搏谁都能聊上两句头大的是它太“虚”了——碳中和涉及能源、工业、交通、碳汇等无数系统边界模糊数据庞杂怎么把它变成一个能用数学语言清晰描述、能用算法求解的模型这中间的鸿沟就是建模竞赛最核心的挑战也是最能体现我们价值的地方。这道题的精髓不在于让你去设计一个国家的碳中和路线图那太不切实际。它的核心是考察你如何将一个复杂的现实问题通过合理的假设、抽象和简化转化为一个结构化的数学模型并设计有效的算法进行求解或分析。简单说就是“大题小做”在有限的篇幅和时间内展现你定义问题、拆解问题和解决问题的能力。很多新手队伍容易犯的错误就是试图面面俱到构建一个包罗万象的超级模型结果往往是参数过多、数据难寻、求解失败论文变成了一锅“理论浆糊”。我们的策略恰恰相反抓住一个具体场景做深做透。我们当时选择的核心场景是“区域电网的低碳化调度与规划”。为什么选这个第一电力系统是碳排放的大户也是实现碳中和的主战场之一重要性毋庸置疑。第二电力系统的运行有相对成熟的物理模型比如直流潮流、机组组合数据相对规范负荷曲线、机组参数、可再生能源出力预测这为建模提供了坚实的基础。第三这个问题天然地融合了优化、预测和多目标决策非常适合数学建模发挥。接下来我就把我们当时的思考过程、模型构建的细节、代码实现的坑以及最后论文的亮点毫无保留地拆解一遍。这不是标准答案而是我们趟出来的一条路希望能给你带来启发。2. 问题定义与模型框架设计如何将“碳中和”装进数学公式面对“碳中和”这样的大命题第一步也是最关键的一步就是划定你的战场。我们将其具体化为在满足一个区域未来若干年比如5年电力需求增长的前提下如何规划新增发电机组风电、光伏、储能、气电等的容量和布局并优化其年度运行方式使得在规划期内的总成本投资成本运行成本最低同时碳排放总量不超过一个逐年递减的硬性约束。这个定义一下子就把问题收紧了。它包含了几个关键要素时间尺度中长期规划年为单位与短期运行小时为单位的结合。决策变量一是容量规划变量整数单位MW二是机组出力变量连续单位MWh。核心约束电力供需平衡、机组技术约束爬坡率、最小技术出力、网络约束简化、以及最重要的——碳排放总量约束。目标函数总成本最小化这是一个典型的经济性目标。基于此我们构建了一个两阶段随机规划模型。为什么用随机规划因为风电、光伏的出力具有强不确定性直接用历史平均值或典型日曲线会严重失真必须考虑其随机性对系统运行和规划的影响。2.1 模型核心数学表达我们的模型主体是一个混合整数线性规划MILP问题下面是其核心部分的白话解读和数学表达第一阶段投资决策Here-and-Now在规划期初我们需要决定建什么电厂、建多大。这些决策一旦做出就不可更改且必须在所有可能的风光出力场景下都可行。决策变量是各种类型机组i的新增容量x_iMW。第二阶段运行模拟Wait-and-See在每一个具体的风光出力场景s我们通过历史数据生成了大量场景下我们需要做出最优的运行决策每台机组在每个时段t的出力P_{i,t,s}储能设备的充放电功率P_{ch/dis,t,s}和荷电状态E_{t,s}等。这些决策依赖于第一阶段的投资和当前场景的具体情况。目标函数最小化总期望成本Minimize: 总投资成本 期望运行成本总投资成本 Σ (单位容量投资成本_i * x_i) 期望运行成本 Σ (场景概率_s * Σ (机组运行成本 启停成本 碳排放成本)) 这里的关键是引入了“碳排放成本”或者说“碳约束的影子价格”。我们并没有直接把它作为一个成本项加进去而是通过下面的约束来实现。核心约束一碳排放总量上限这是体现“碳中和”目标的核心约束。对于每一年yΣ_{t, i, s} [概率_s * (机组i的碳排放强度 * P_{i,t,s})] 年度碳排放限额_y年度碳排放限额_y是一个逐年递减的数值比如每年比上一年减少5%。这个约束将环境目标直接转化为对系统运行的硬性限制。在求解时这个约束会对应一个拉格朗日乘子影子价格其经济意义就是该年度的“碳价”。碳价越高说明减排压力越大模型会更倾向于使用可再生能源。核心约束二电力平衡与机组运行对于每一个时间点t和每一个场景sΣ (所有机组出力 储能放电 - 储能充电) 该时段负荷需求同时每个机组的出力必须在其最小技术出力和最大容量现有容量新增容量x_i之间并且要满足爬坡率约束。核心约束三可再生能源消纳我们强制要求风电和光伏的发电量必须全部上网除非技术性弃电这体现了优先消纳清洁能源的政策导向。这个模型框架将长期的容量规划、短期的随机运行和碳排放的刚性约束有机地耦合在了一起。它的输出不仅仅是“该建多少风电和光伏”还包括了在不同减排力度下系统的最优电源结构演变、系统的运行成本变化、以及隐含的碳价信号。注意这里有一个重要的建模技巧。直接求解这样一个包含大量场景可能成千上万的两阶段随机规划MILP问题计算量是灾难性的。我们采用了“场景削减”技术通过K-means聚类或前向选择算法从大量场景中选取几十个具有代表性的典型场景来近似代替完整的概率分布从而在计算精度和求解时间之间取得平衡。3. 数据准备、处理与关键假设模型建得再漂亮没有数据支撑就是空中楼阁。五一杯这类竞赛通常不提供现成数据需要自己搜集、加工这部分的工作量和技术含量绝不亚于建模本身。3.1 数据需求清单电力负荷数据目标区域的历史逐时负荷数据用于预测未来负荷。我们使用了该区域过去3-5年的数据通过时间序列分析如SARIMA模型预测未来5年的负荷曲线并考虑了年均增长率。可再生能源数据风电场和光伏电站的历史逐时出力数据。同样用于生成未来出力的随机场景。我们从气象数据库获取了风速和辐照度数据通过功率转换模型得到出力。技术经济参数各类电源单位容量投资成本元/kW、固定运维成本元/kW/年、可变运维成本元/MWh、发电效率、碳排放强度gCO2/kWh、最小技术出力、爬坡率、寿命等。储能系统投资成本元/kWh和元/kW、循环效率、充放电速率、寿命。碳排放限额根据国家或区域的碳中和目标自己设定一个合理的逐年递减路径。例如以某基准年为起点设定每年减排百分比。3.2 数据处理中的坑与技巧风光出力场景生成这是随机规划的灵魂。我们采用的方法是首先对历史风光数据进行聚类得到几种典型的天气类型如“大风晴天”、“小风阴天”等。然后针对每一种天气类型用Copula函数描述风速和辐照度的联合分布再通过蒙特卡洛模拟生成大量相关性的出力场景。最后用场景削减算法得到代表性场景。# 示例使用K-means进行场景削减伪代码 from sklearn.cluster import KMeans # historical_profiles 形状为 (n_samples, n_timesteps) kmeans KMeans(n_clusters20, random_state42) cluster_labels kmeans.fit_predict(historical_profiles) representative_scenarios kmeans.cluster_centers_ # 得到20个典型场景 scenario_probabilities np.bincount(cluster_labels) / len(cluster_labels) # 计算每个场景的概率负荷预测的波动性不仅要预测平均负荷还要预测负荷的波动范围。我们在确定性预测的基础上叠加了一个基于历史误差分布的随机扰动以体现负荷的不确定性。成本参数的归一化与贴现投资成本是期初一次性投入运行成本是每年发生。为了在目标函数中相加必须将未来每年的运行成本贴现到规划期初。我们使用了标准的净现值计算方法设定了一个社会折现率如5%。关键假设必须明确在论文中我们专门用一小节罗列所有重要假设例如“忽略输电网络阻塞”、“风电机组和光伏组件的成本在未来五年内按年均下降率X%考虑”、“碳排放限额为外生给定且必须严格执行”。清晰的假设是模型合理性的护城河。4. 模型求解算法选择与代码实现核心模型是MILP求解器自然首选Gurobi、CPLEX或COPT这类商业求解器它们对大规模MILP的求解效率远超开源工具。竞赛环境通常允许使用这些求解器的免费学术版。4.1 编程框架与代码结构我们选择Python作为编程语言搭配gurobipy库。代码结构清晰是团队协作和调试的基础。项目目录/ ├── data/ # 存放所有原始和处理后的数据 │ ├── load.csv │ ├── wind_scenarios.csv │ └── tech_economic_params.json ├── src/ │ ├── data_preprocessing.py # 数据清洗、场景生成 │ ├── model_builder.py # 构建Gurobi模型对象定义变量、约束、目标 │ ├── solver_config.py # 求解器参数设置如时间限制、MIPGap │ └── post_processing.py # 结果解析、可视化 ├── main.py # 主程序串联整个流程 └── results/ # 输出结果图表、报告4.2model_builder.py核心代码解析这里展示最关键的建模部分代码片段并附上详细注释。import gurobipy as gp from gurobipy import GRB import numpy as np def build_two_stage_stochastic_model(data, scenarios, probabilities): 构建两阶段随机规划模型。 data: 包含所有参数成本、技术参数、负荷等的字典 scenarios: 列表每个元素是一个字典代表一个风光出力场景 probabilities: 列表每个场景对应的概率 model gp.Model(Carbon_Neutrality_Power_Planning) # 第一阶段变量投资决策 # x[tech]: 技术类型tech的新增容量 (MW) x model.addVars(data[technologies], lb0, vtypeGRB.CONTINUOUS, namex_capacity) # 第二阶段变量每个场景下的运行决策 # 这里以火电机组为例其他类似 # p[scenario_idx, tech, time]: 机组出力 # u[scenario_idx, tech, time]: 机组启停状态 (0/1) p {} u {} for s_idx, sc in enumerate(scenarios): for tech in data[dispatchable_techs]: # 可调度机组如火电、气电 for t in range(data[num_time_periods]): p[s_idx, tech, t] model.addVar(lb0, namefp_s{s_idx}_{tech}_t{t}) u[s_idx, tech, t] model.addVar(vtypeGRB.BINARY, namefu_s{s_idx}_{tech}_t{t}) # 目标函数总投资成本 期望运行成本 # 1. 投资成本 investment_cost gp.quicksum(data[inv_cost][tech] * x[tech] for tech in data[technologies]) # 2. 期望运行成本 expected_op_cost 0 for s_idx, (sc, prob) in enumerate(zip(scenarios, probabilities)): scenario_cost 0 # 燃料成本、运维成本 for tech in data[dispatchable_techs]: for t in range(data[num_time_periods]): scenario_cost data[var_om_cost][tech] * p[s_idx, tech, t] scenario_cost data[startup_cost][tech] * model.addVar(...) # 启停成本变量需额外定义 # 碳排放“成本”通过约束体现此处不直接加入目标 # 将场景成本乘以概率加到总期望成本中 expected_op_cost prob * scenario_cost model.setObjective(investment_cost expected_op_cost, GRB.MINIMIZE) # 核心约束添加 # 约束1: 每个场景每个时刻的电力平衡 for s_idx, sc in enumerate(scenarios): for t in range(data[num_time_periods]): # 总发电火电气电风电出力光伏出力储能放电 total_generation gp.quicksum(p[s_idx, tech, t] for tech in data[dispatchable_techs]) \ sc[wind][t] sc[pv][t] \ (discharge_power[s_idx, t] - charge_power[s_idx, t]) # 储能 model.addConstr(total_generation data[load][t], namefbalance_s{s_idx}_t{t}) # 约束2: 机组出力上下限约束与投资变量x关联 for s_idx, sc in enumerate(scenarios): for tech in data[dispatchable_techs]: existing_cap data[existing_capacity][tech] max_cap existing_cap x[tech] # 最大出力不能超过现有新增容量 for t in range(data[num_time_periods]): model.addConstr(p[s_idx, tech, t] max_cap * u[s_idx, tech, t], namefmax_cap_s{s_idx}_{tech}_t{t}) model.addConstr(p[s_idx, tech, t] data[min_output][tech] * u[s_idx, tech, t], namefmin_cap_s{s_idx}_{tech}_t{t}) # 约束3: 年度碳排放总量约束最关键 yearly_emissions {} for year in data[planning_years]: yearly_emissions[year] 0 # 计算该年份在所有场景下的期望碳排放量 # 简化处理假设每个场景代表一年中的一种可能运行情况 for s_idx, (sc, prob) in enumerate(zip(scenarios, probabilities)): for tech in data[dispatchable_techs]: tech_emission gp.quicksum(data[carbon_intensity][tech] * p[s_idx, tech, t] for t in range(data[num_time_periods])) yearly_emissions[year] prob * tech_emission # 添加约束期望碳排放 该年限额 model.addConstr(yearly_emissions[year] data[carbon_cap][year], namefcarbon_cap_{year}) # ... 其他约束爬坡、储能动态、可再生能源消纳等 return model, x, p, u # 返回模型和关键变量便于后续处理4.3 求解策略与调参经验设置合理的求解时间与MIP Gap对于复杂模型追求最优解可能耗时过长。我们通常设置一个时间限制如2小时和一个可接受的MIP Gap如0.5%或1%。这样能在有限时间内得到一个高质量的可行解。model.setParam(TimeLimit, 7200) # 2小时 model.setParam(MIPGap, 0.005) # 0.5%的间隙利用回调函数记录进度在长时间求解过程中使用回调函数输出当前最优解、间隙等信息既能监控进度也能在意外中断后有一个可用的解。分步求解策略如果模型太大可以尝试先求解一个简化版比如减少时间分辨率从逐时变为4小时一点或者减少场景数得到一个初始解然后将其作为初始解提供给完整模型能大大加速求解过程。并行计算Gurobi支持多线程并行求解MIP问题。确保你的代码环境能利用多核CPU设置model.setParam(Threads, 0)让Gurobi自动使用所有可用线程。5. 结果分析与可视化让模型“说话”模型求解完成后输出是一堆数字。如何将这些数字转化为有说服力的结论和直观的图表是论文拿高分的关键。5.1 核心结果输出最优投资方案每种技术类型的新增容量x_i是多少这是最直接的规划建议。系统总成本与成本构成总投资是多少运行成本是多少碳排放约束导致了多少额外的“系统成本”隐含碳价通过查询碳排放总量约束的对偶变量影子价格我们可以得到每一年每吨CO2的隐含价格。这是一个非常有力的经济信号指标。典型日运行模拟选取一个代表性场景如夏季高峰日绘制该日的电力供需平衡图。图中应清晰展示负荷曲线、各类电源的出力堆叠、储能充放电状态。这张图能直观验证模型的合理性和系统的运行特性。敏感性分析这是体现模型稳健性和你思考深度的部分。改变关键参数看结果如何变化。碳排放限额的敏感性如果减排目标更激进限额下降更快最优电源结构如何变化系统成本增加多少可再生能源成本的敏感性如果风电/光伏/储能成本下降速度超预期规划结果会怎样负荷增长的敏感性如果经济增长超预期电力需求更高如何调整规划5.2 可视化代码示例使用Matplotlibimport matplotlib.pyplot as plt import pandas as pd def plot_dispatch_for_typical_day(scenario_idx, results, data): 绘制某个典型场景下一天的调度结果。 fig, ax plt.subplots(figsize(14, 7)) times range(24) # 获取结果数据假设results是一个包含所有变量值的字典 p_coal [results[p][scenario_idx, coal, t] for t in times] p_gas [results[p][scenario_idx, gas, t] for t in times] p_wind [data[scenarios][scenario_idx][wind][t] for t in times] p_pv [data[scenarios][scenario_idx][pv][t] for t in times] load [data[load][t] for t in times] # 创建堆叠面积图 ax.stackplot(times, p_coal, p_gas, p_wind, p_pv, labels[燃煤, 燃气, 风电, 光伏], colors[#8B4513, #FFA500, #87CEEB, #FFD700]) # 绘制负荷曲线 ax.plot(times, load, colorblack, linewidth2, label负荷需求, markero) ax.set_xlabel(小时, fontsize12) ax.set_ylabel(功率 (MW), fontsize12) ax.set_title(典型日电力系统调度模拟, fontsize14, fontweightbold) ax.legend(locupper left) ax.grid(True, linestyle--, alpha0.6) ax.set_xlim(0, 23) plt.xticks(times) plt.tight_layout() plt.savefig(typical_day_dispatch.png, dpi300) plt.show() def plot_capacity_expansion(results, data): 绘制规划期内电源结构演变图。 techs [coal, gas, wind, pv, storage] existing_cap data[existing_capacity] new_cap results[x_optimal] # 最优新增容量 total_cap {} for tech in techs: total_cap[tech] existing_cap.get(tech, 0) new_cap.get(tech, 0) df pd.DataFrame(list(total_cap.items()), columns[Technology, Capacity]) df df.sort_values(Capacity, ascendingFalse) fig, ax plt.subplots(figsize(10, 6)) bars ax.bar(df[Technology], df[Capacity], color[#8B4513, #FFA500, #87CEEB, #FFD700, #9370DB]) ax.set_ylabel(装机容量 (MW), fontsize12) ax.set_title(规划期末最优电源结构, fontsize14, fontweightbold) # 在柱子上标注数值 for bar in bars: height bar.get_height() ax.text(bar.get_x() bar.get_width()/2., height 10, f{int(height)}, hacenter, vabottom) plt.tight_layout() plt.savefig(optimal_capacity_mix.png, dpi300) plt.show()5.3 敏感性分析示例我们以碳排放限额为例展示如何编程实现批量计算和可视化。def sensitivity_analysis_carbon_cap(base_data, cap_reduction_rates): 分析不同碳排放下降速率对结果的影响。 cap_reduction_rates: 列表如 [0.03, 0.05, 0.07, 0.10] 表示每年减排3%,5%,7%,10% results_summary [] for rate in cap_reduction_rates: # 1. 修改数据中的碳排放限额 modified_data base_data.copy() base_year_cap modified_data[carbon_cap][2023] for year in modified_data[carbon_cap]: years_from_base year - 2023 modified_data[carbon_cap][year] base_year_cap * ((1 - rate) ** years_from_base) # 2. 重新构建并求解模型 model, _, _, _ build_two_stage_stochastic_model(modified_data, scenarios, probs) model.optimize() if model.status GRB.OPTIMAL or model.status GRB.TIME_LIMIT: total_cost model.ObjVal new_cap_wind model.getVarByName(x_capacity[wind]).X # ... 获取其他关键结果 results_summary.append({ reduction_rate: rate, total_cost: total_cost, wind_capacity: new_cap_wind, # ... 其他指标 }) else: print(f求解失败减排率: {rate}) results_summary.append({reduction_rate: rate, total_cost: None}) # 3. 可视化 df_sens pd.DataFrame(results_summary) fig, ax1 plt.subplots(figsize(10, 6)) ax1.plot(df_sens[reduction_rate], df_sens[total_cost], b-o, linewidth2, label系统总成本) ax1.set_xlabel(碳排放年下降速率, fontsize12) ax1.set_ylabel(系统总成本 (亿元), colorb, fontsize12) ax1.tick_params(axisy, labelcolorb) ax2 ax1.twinx() ax2.plot(df_sens[reduction_rate], df_sens[wind_capacity], r-s, linewidth2, label风电新增容量) ax2.set_ylabel(风电新增容量 (MW), colorr, fontsize12) ax2.tick_params(axisy, labelcolorr) fig.suptitle(碳排放约束强度对系统规划的影响, fontsize14, fontweightbold) fig.legend(locupper left, bbox_to_anchor(0.1, 0.9)) plt.tight_layout() plt.savefig(sensitivity_carbon_cap.png, dpi300) plt.show() return df_sens通过这样的分析我们可以得出清晰的结论随着减排力度加大系统总成本会上升同时电源结构会向风电、光伏等清洁能源显著倾斜。这种定量的结论比空泛的论述有力得多。6. 论文撰写要点与竞赛心得模型和代码是骨架论文才是呈现给评委的血肉。一篇好的数模论文逻辑清晰、图表美观、表达准确三者缺一不可。1. 摘要要像“电梯演讲”摘要必须在半页纸内讲清楚针对什么问题、建立了什么模型、用了什么方法、得到了什么结论、有什么特色亮点。避免细节突出整体逻辑和创新点。我们当时的摘要第一句就是“本文针对‘碳中和’目标下的区域电力系统规划问题构建了一个考虑风光不确定性的两阶段随机规划模型……”2. 模型假设部分不能少明确列出所有主要假设这是模型成立的前提也能体现你思考的严谨性。例如“假设输电网络无阻塞”、“假设未来五年燃料价格保持不变”、“假设所有机组均能可靠运行”。3. 模型部分重逻辑轻公式不要堆砌公式。先文字描述清楚模型的整体框架、决策变量、目标、约束的核心思想再用清晰的公式表达。对于复杂的约束如储能动态最好配以简单的示意图或流程图说明。4. 灵敏度分析是加分项它展示了模型的稳健性和你对问题理解的深度。不要只做一个参数的敏感性最好能做2-3个关键参数如碳限额、风光成本、折现率的分析并讨论其现实意义。5. 图表质量决定第一印象多用组合图比如将电源结构演变和系统成本变化放在一张图的两个Y轴上。颜色要专业使用区分度高的颜色并保持一致性如煤电用棕色、气电用橙色、风电用蓝色、光伏用黄色。标注要清晰坐标轴标签、单位、图例必须一目了然。图中重要的数据点可以标注具体数值。避免截图尽量导出矢量图如PDF、SVG格式在论文中插入以保证清晰度。6. 代码附录与可重复性在附录中提供核心算法的伪代码或流程图并说明主要函数的功能。虽然不要求提交全部代码但清晰的说明能让评委相信你的工作是扎实、可重复的。7. 团队协作是效率关键三个人一定要有明确分工一人主攻模型与算法负责model_builder.py和求解一人主攻数据处理与可视化负责data_preprocessing.py和post_processing.py一人主攻论文撰写与整合。每天至少同步两次用Git管理代码和论文版本避免最后时刻合并冲突。回过头看这道“碳中和”赛题考察的远不止数学和编程。它要求我们从庞杂的现实问题中提炼出科学问题用严谨的数学工具进行刻画再用可靠的工程方法求解最后用清晰的逻辑和专业的表达呈现出来。这个过程本身就是一次微缩的科研训练。最大的收获不是那个奖项而是这套从“问题”到“解决方案”的完整思维框架和实战能力。如果你也在准备类似的竞赛我的建议是尽早选定一个具体场景把模型做“厚”把故事讲“薄”。深度永远比广度更有力量。