简介本资源面向微电网容量配置方向的毕业设计与算法研究者围绕《基于两阶段鲁棒优化算法的微网多电源容量配置》一文提供完整程序实现。针对可再生能源与负荷不确定性构建min-max-min结构的两阶段鲁棒优化模型在储能、需求侧负荷及可控分布式电源约束下求解最恶劣场景的最低运行成本并引入不确定性调节参数灵活控制调度保守性采用列约束生成算法与强对偶理论将原问题分解为主、子问题交替求解。压缩包共423个文件约89.47MB以276个xls数据表、110个mat数据文件为主辅以m脚本、csv与xlsx算例、docx说明及caj文献覆盖建模、算例与结果分析全流程。已有148人学习适合需要复现鲁棒优化调度、验证分时电价下储能调度边界条件的读者参考。1. 两阶段鲁棒优化做微网容量配置这份源程序把论文里的公式变成了能跑的代码微网容量配置这件事最怕的不是算不出来而是算出来的结果一遇到光伏波动、负荷突增就崩了。确定性优化给出的那套配置方案在仿真里跑得漂漂亮亮真到工程评审时被问一句「阴雨天连续三天怎么办」就哑火。这份源程序对应的论文《基于两阶段鲁棒优化算法的微网多电源容量配置》走的是另一条路先把最恶劣的场景找出来再在这个场景下求最优配置用「最坏情况下的最优解」替代「平均情况下的最优解」。源程序把两阶段鲁棒优化的建模、对偶转化、列与约束生成CCG求解流程全部落地成了可运行的代码适合做毕业设计、课程设计或者想从确定性优化跨到鲁棒优化的同学。论文在知网可下载配套博客有逐段解读源程序本身是复现论文结果的完整工程。2. 两阶段鲁棒优化的数学骨架从 min-max-min 到可求解的主子问题2.1 为什么微网容量配置需要两阶段结构微网容量配置的核心矛盾是投资决策要在不确定性揭晓之前做运行调度要在不确定性揭晓之后调。风机、光伏的出力是波动的负荷也不是一条平线如果用一个确定性的典型日曲线去优化得到的容量方案本质上是在赌天气。两阶段鲁棒优化把这个问题拆成两层第一阶段决定装多少块光伏板、配多大储能、买多大容量的柴油机这些是「here-and-now」变量一旦定了就不改第二阶段是在给定容量下面对实际出现的不确定场景去调度各电源出力这些是「wait-and-now」变量可以随场景调整。数学形式写出来就是 min-max-min 三层结构。最外层 min 是对投资变量求最小化中间 max 是让不确定变量取到使运行成本最大的那个场景最内层 min 是在这个最坏场景下求运行成本最小。这个结构直接求解是不可能的因为中间那层 max 让问题变成了半无限规划。常见做法是用对偶理论把内层 min 转成 max然后和中间的 max 合并得到一个单层的 max 问题再用 CCG 算法交替求解主问题MP和子问题SP。注意两阶段鲁棒优化的「两阶段」指的是决策时序不是指算法有两个步骤。很多初学者会把两阶段和 CCG 的两步混淆实际上 CCG 是求解两阶段模型的一种算法不是模型本身。2.2 不确定集怎么选盒式、多面体还是数据驱动不确定集的选取直接决定鲁棒优化的保守程度。最常用的是盒式不确定集每个不确定参数独立地在区间内波动形式简单但过于保守因为所有参数同时取到最坏值的概率极低。稍微精细一点的是多面体不确定集引入一个预算参数 Γ 来控制同时波动的参数个数Γ 越大越保守Γ 等于参数总数时退化为盒式。源程序里用的是盒式不确定集光伏出力和负荷各自有上下界这是论文的设定。如果你要改成多面体需要在子问题里增加对偶变量和预算约束代码改动量不小。我一般会建议先跑通盒式版本理解 CCG 的迭代逻辑之后再尝试把不确定集换成多面体或者基于历史数据构建的数据驱动不确定集。数据驱动的方法需要你有光伏和负荷的历史时序数据用聚类或者分位数回归提取不确定边界这一步在源程序里没有属于进阶改造。2.3 CCG 主问题与子问题的代码实现主问题是在已知的有限个最坏场景集合下求投资变量和运行变量的最优解。随着迭代进行场景集合不断扩充主问题的目标值单调递增给出下界。子问题是在给定投资变量下寻找使运行成本最大的不确定场景给出上界。上下界收敛时停止。# 主问题给定场景集优化投资变量和运行变量 def solve_master(scenarios, params): scenarios: 当前已识别的最坏场景列表 params: 包含投资成本系数、运行成本系数、设备参数等 返回投资变量值、运行变量值、主问题目标值下界 m gp.Model(master) # 投资变量光伏容量、储能容量、柴油机容量 x_pv m.addVar(lb0, namex_pv) x_ess m.addVar(lb0, namex_ess) x_dg m.addVar(lb0, namex_dg) # 运行变量按场景索引 p_pv {} p_ess {} p_dg {} for s in range(len(scenarios)): p_pv[s] m.addVar(lb0, namefp_pv_{s}) p_ess[s] m.addVar(lb-params[ess_max], ubparams[ess_max], namefp_ess_{s}) p_dg[s] m.addVar(lb0, namefp_dg_{s}) # 目标投资成本 最坏场景下的运行成本 invest_cost params[c_pv] * x_pv params[c_ess] * x_ess params[c_dg] * x_dg run_cost 0 for s in range(len(scenarios)): run_cost params[c_dg_run] * p_dg[s] params[c_ess_run] * gp.abs_(p_ess[s]) m.setObjective(invest_cost run_cost, gp.MINIMIZE) # 约束功率平衡、设备出力上限、储能SOC等 for s, scenario in enumerate(scenarios): m.addConstr(p_pv[s] p_ess[s] p_dg[s] scenario[load] - scenario[pv_actual]) m.addConstr(p_pv[s] x_pv * scenario[pv_available]) m.addConstr(p_dg[s] x_dg) m.optimize() return x_pv.X, x_ess.X, x_dg.X, m.ObjVal上面这段是主问题的骨架实际代码里还要加储能 SOC 连续性约束、柴油机爬坡约束、光伏出力不能超过可用容量等。参数说明几个关键的c_pv、c_ess、c_dg是单位容量投资成本c_dg_run是柴油机单位发电运行成本ess_max是储能最大充放电功率。scenario[pv_actual]是不确定变量在当前场景下的取值scenario[pv_available]是光伏可用出力系数。子问题部分需要把内层 min 对偶化得到一个 max 问题。源程序里用的是对偶转化加线性化因为内层有绝对值项和双线性项。对偶之后子问题变成一个混合整数线性规划可以用 Gurobi 或 CPLEX 直接求解。# 子问题给定投资变量寻找最坏场景 def solve_subproblem(x_pv_val, x_ess_val, x_dg_val, params): 输入主问题给出的投资变量值 返回最坏场景下的不确定变量取值、子问题目标值上界 m gp.Model(subproblem) # 不确定变量光伏实际出力系数、负荷波动系数 u_pv m.addVar(lbparams[pv_min], ubparams[pv_max], nameu_pv) u_load m.addVar(lbparams[load_min], ubparams[load_max], nameu_load) # 对偶变量 lambda1 m.addVar(lbNone, namelambda1) lambda2 m.addVar(lb0, namelambda2) # 目标最大化运行成本对偶后的形式 m.setObjective(lambda1 * (u_load * params[load_base] - u_pv * params[pv_base]) lambda2 * x_pv_val, gp.MAXIMIZE) # 对偶约束 m.addConstr(lambda1 - lambda2 params[c_dg_run]) m.addConstr(lambda1 params[c_ess_run]) m.addConstr(lambda1 -params[c_ess_run]) m.optimize() return u_pv.X, u_load.X, m.ObjVal子问题的对偶变量lambda1对应功率平衡约束lambda2对应光伏出力上限约束。pv_min、pv_max是光伏出力的不确定区间load_min、load_max是负荷的不确定区间。子问题求出的u_pv和u_load就是当前投资方案下最恶劣的光伏和负荷组合把这个场景加回主问题进入下一轮迭代。2.4 迭代收敛判据与上下界更新CCG 的收敛判据是上下界间隙小于阈值。主问题给出下界 LB子问题给出上界 UB当 (UB - LB) / LB ε 时停止。源程序里 ε 默认取 0.01可以调小到 0.001 获得更精确的解但迭代次数会增加。# CCG 主循环 LB -float(inf) UB float(inf) scenarios [initial_scenario] # 初始场景通常取不确定区间的中点 max_iter 50 gap_tol 0.01 for k in range(max_iter): # 求解主问题 x_pv_val, x_ess_val, x_dg_val, mp_obj solve_master(scenarios, params) LB max(LB, mp_obj) # 求解子问题 u_pv_val, u_load_val, sp_obj solve_subproblem(x_pv_val, x_ess_val, x_dg_val, params) UB min(UB, mp_obj sp_obj) # 注意这里要加上投资成本 # 收敛判断 if (UB - LB) / abs(LB) gap_tol: print(f收敛于第 {k1} 次迭代间隙 {(UB-LB)/abs(LB):.4f}) break # 添加新场景到主问题 new_scenario {pv_actual: u_pv_val, load: u_load_val * params[load_base], pv_available: u_pv_val} scenarios.append(new_scenario)这段循环里有个容易翻车的地方UB的更新不能直接用子问题目标值因为子问题只算了运行成本上界应该是投资成本加最坏运行成本。源程序里把投资成本单独存下来在更新 UB 时加上。另外初始场景的选择会影响收敛速度取中点是比较稳妥的做法取边界值可能导致主问题初始可行域过窄。3. 把源程序跑起来环境配置、数据替换与结果验证3.1 依赖安装与求解器配置源程序是 Python 写的依赖 Gurobi 求解器。Gurobi 需要许可证学校邮箱可以申请学术版否则只能用试用版变量规模受限。如果不想折腾 Gurobi可以把求解器换成 CBC 或 GLOP但混合整数问题的求解速度会慢很多CCG 迭代几次就可能卡住。# 创建虚拟环境 python -m venv venv source venv/bin/activate # Windows 用 venv\Scripts\activate # 安装依赖 pip install gurobipy numpy pandas matplotlib # 验证 Gurobi 许可证 python -c import gurobipy; m gurobipy.Model(); print(Gurobi OK)如果gurobipy导入报错大概率是许可证没配好。学术版许可证需要把gurobi.lic文件放到用户目录下Windows 是C:\Users\你的用户名\Linux 是/home/你的用户名/。环境变量GRB_LICENSE_FILE也可以指定许可证路径。3.2 输入数据格式与替换方法源程序的数据放在data/目录下光伏出力、负荷、电价各一个 CSV 文件。光伏和负荷是 24 小时时序数据单位是标幺值或者 kW具体看论文里的基准值。替换数据时注意列名要和代码里的读取逻辑一致源程序用的是pv、load、price三列。# 数据读取与预处理 import pandas as pd import numpy as np def load_data(data_dir): 读取光伏、负荷、电价数据 返回包含各时序数组的字典 pv pd.read_csv(f{data_dir}/pv.csv)[pv].values load pd.read_csv(f{data_dir}/load.csv)[load].values price pd.read_csv(f{data_dir}/price.csv)[price].values # 归一化处理论文里通常用标幺值 pv_norm pv / pv.max() load_norm load / load.max() return {pv: pv_norm, load: load_norm, price: price, pv_base: pv.max(), load_base: load.max()}替换数据时最常见的坑是量纲不统一。光伏数据如果是 W 而负荷是 kW功率平衡约束直接崩掉。我一般会在读取之后打印各序列的最大值和均值确认量级一致再往下跑。另外时间分辨率也要注意源程序默认是 1 小时一个点如果你的数据是 15 分钟一个点需要先做聚合或者修改代码里的时段数。3.3 结果解读容量配置方案与运行成本构成跑完之后源程序会输出三个投资变量的值光伏容量、储能容量、柴油机容量以及最坏场景下的运行成本。结果解读时重点看两个地方一是储能容量和光伏容量的比例鲁棒优化通常会配更多储能来应对光伏波动二是柴油机容量如果柴油机容量接近零说明在最坏场景下可再生能源加储能已经能覆盖负荷。# 结果可视化 import matplotlib.pyplot as plt def plot_results(result): result: 包含投资变量、运行成本、迭代历史的字典 fig, axes plt.subplots(1, 2, figsize(12, 4)) # 容量配置柱状图 axes[0].bar([PV, ESS, DG], [result[x_pv], result[x_ess], result[x_dg]]) axes[0].set_ylabel(Capacity (kW)) axes[0].set_title(Optimal Capacity Configuration) # 迭代收敛曲线 axes[1].plot(result[LB_history], labelLB) axes[1].plot(result[UB_history], labelUB) axes[1].set_xlabel(Iteration) axes[1].set_ylabel(Cost) axes[1].legend() axes[1].set_title(CCG Convergence) plt.tight_layout() plt.savefig(results.png, dpi150)收敛曲线如果出现上下界震荡不收敛通常是子问题对偶转化有误或者不确定集边界设得太宽导致主问题可行域变化剧烈。可以先检查子问题目标值是否始终大于等于主问题目标值如果不是说明对偶方向搞反了。4. 避坑与排查CCG 迭代不收敛、对偶符号错误、数据量纲混乱4.1 迭代震荡不收敛现象上下界在迭代十几次后仍然来回跳动间隙不缩小。原因通常是子问题求出的最坏场景没有正确加回主问题或者主问题里场景对应的运行变量没有正确索引。解决方法是检查scenarios列表是否在每次迭代后都追加了新场景以及主问题里p_pv[s]、p_ess[s]、p_dg[s]的索引是否和场景列表对齐。另一个可能原因是子问题的对偶变量符号搞反了导致求出的不是最坏场景而是最好场景。4.2 对偶转化后目标函数符号错误现象子问题目标值为负或者主问题下界大于上界。原因是对偶理论里 min 转 max 时目标函数要取反约束的符号也要跟着变。源程序里内层是 min 运行成本对偶之后变成 max 对偶目标如果忘记取反整个 CCG 的上下界逻辑就反了。排查方法是手动算一个简单场景的对偶值和原问题对比。4.3 数据量纲不统一导致功率平衡崩溃现象求解器报 infeasible或者结果里储能充放电功率异常大。原因是光伏、负荷、设备容量的单位不一致。源程序里所有功率单位是 kW如果你替换的数据是 MW需要乘以 1000。另外标幺值的基准值也要统一光伏和负荷不能用不同的基准。4.4 Gurobi 许可证失效或变量规模超限现象跑小规模数据正常换成全年 8760 小时数据后报内存不足或许可证错误。原因是试用版 Gurobi 限制 2000 个变量和约束全年数据轻松超过。解决方法是先用典型日或者聚类后的几个场景跑通逻辑确认无误后再用学术版许可证跑全年数据。如果只能用试用版可以把时间分辨率降到 4 小时或者 6 小时。4.5 储能 SOC 约束遗漏导致充放电无限制现象储能容量配置结果为零但运行成本却很低。原因是储能 SOC 约束没加储能可以无限充放电而不受容量限制优化器自然不愿意投资储能。源程序里 SOC 约束在constraints.py里检查soc_min、soc_max、soc_init三个参数是否都设置了以及 SOC 连续性约束是否写对。5. 从跑通到改对不确定集预算参数 Γ 的调参与结果验证盒式不确定集跑通之后最值得动手改的是把不确定集换成多面体引入预算参数 Γ 来控制保守程度。Γ 的物理含义是在所有不确定参数中最多允许 Γ 个同时取到最坏值。Γ 等于 0 时退化为确定性优化Γ 等于参数总数时退化为盒式鲁棒优化。调参时可以从 Γ1 开始逐步增加到参数总数观察投资成本和运行成本的变化。# 多面体不确定集下的子问题改造关键部分 def solve_subproblem_polyhedral(x_pv_val, params, Gamma): Gamma: 预算参数控制同时波动的不确定参数个数 m gp.Model(subproblem_poly) T params[T] # 时段数 # 不确定变量每个时段的光伏和负荷波动 u_pv m.addVars(T, lb-1, ub1, nameu_pv) u_load m.addVars(T, lb-1, ub1, nameu_load) # 预算约束所有不确定变量绝对值之和不超过 Gamma m.addConstr(gp.quicksum(u_pv[t] u_load[t] for t in range(T)) Gamma) m.addConstr(gp.quicksum(u_pv[t] u_load[t] for t in range(T)) -Gamma) # 后续对偶转化和盒式类似但需要对每个时段分别处理 # ...省略对偶变量定义和约束 m.optimize() return m.ObjValΓ 的调参结果可以画成一条曲线横轴 Γ纵轴总成本。曲线通常先快速上升后趋于平缓拐点对应的 Γ 就是性价比最高的保守程度。如果 Γ 从 1 增加到 5 成本只涨了 2%但从 5 增加到 10 成本涨了 15%那 Γ5 附近就是合理选择。验证鲁棒优化结果是否真的鲁棒可以用蒙特卡洛方法在不确定集内随机生成 1000 个场景把优化出的容量配置代入每个场景做运行调度统计运行成本超过子问题最坏成本的场景比例。如果超过 5%说明不确定集边界设窄了或者 Γ 太小。我一般会把这个验证步骤作为最后一道关卡从那以后每次改完不确定集参数都强制跑一遍蒙特卡洛不然心里没底。希望帮到你。本文还有配套的精品资源点击获取