微网多电源容量配置：两阶段鲁棒优化与CCG求解实践

简介面向微电网与电力系统优化研究者的MATLAB源代码包聚焦基于两阶段鲁棒优化算法的微网多电源容量配置问题适合具备一定优化理论基础的学者、工程师用于算法复现与改进。压缩包共426个文件约91.42MB主要包含276个Excel数据文件用于负荷与新能源历史数据整理110个MAT数据文件存储仿真输入输出MATLAB源码.m实现两阶段鲁棒优化模型CSV文件记录实测环境数据另有rar压缩包、Word文档、PDF及CAJ文献等辅助资料。资源在设计上覆盖数据输入、不确定性建模、优化求解、结果可视化与仿真验证等模块便于对照论文理解算法流程。目前已有150人学习下载内容完整度较高可作为微网容量配置研究的参考实现也可用于教学演示与二次开发。1. 两阶段鲁棒优化做微网多电源容量配置先想清楚“最坏场景”再装机做微网多电源容量配置的人大概率都经历过这种场景按典型日出力曲线算出来的容量方案安装完第一年运行就翻车——风电连续几天低出力负荷又顶着尖峰最后只能切负荷经济性账完全对不上。这套代码走的是另一条路用两阶段鲁棒优化算法把“明天风多大、负荷多高”这类不确定性显式建模第一阶段定各电源装多少容量第二阶段在最坏不确定性场景下求最小运行成本。压缩包里是完整的MATLAB代码、六个CSV数据文件、晏鸣宇关于含风电鲁棒机组组合的CAJ文献和一份docx格式的代码说明整套正好覆盖数据—建模—求解—验证闭环适合做微电网规划、电力系统优化方向的工程师和研究生对照复现。2. 从两阶段决策到不确定集合把容量配置拆成可求解的结构2.1 为什么容量配置必须分两层投资决策与运行模拟容量配置问题本质上是一个双层决策外层决定风机、光伏、储能、燃气轮机的装机容量内层模拟这些机组在不同场景下的运行方式。如果只做单层规划把风电出力取一个典型值、负荷取均值得到的配置方案在平均工况下很好看但在极端天气和负荷同时叠加时毫无招架之力。这就是为什么这份代码采用两阶段结构——第一阶段变量在不确定性实现之前就必须定下来第二阶段变量等不确定性参数暴露后再做调整这种“先决策、后调整”的结构非常适合描述容量规划的物理过程。代码里对应的变量划分大概是这样的阶段变量含义决策时机第一阶段x各电源装机容量、机组类型选择不确定性实现前第二阶段y机组出力、启停、储能充放电、切负荷量看到风电/负荷实际值后不确定性u风电出力、负荷需求、温度自然产生无法控制第一阶段目标通常是投资成本加运行成本的期望或最坏值第二阶段目标是在给定x和u下找到最经济的运行方案。两阶段之间通过耦合约束连接比如容量上限约束机组出力不能超过装机容量这一条把x和y绑在一起也是容量配置问题的核心物理约束。2.2 盒式不确定集合与预算参数Γ保守度的旋钮鲁棒优化处理不确定性的方式不是给风电出力一个概率分布而是定义一个不确定集合保证集合内所有场景都能被满足。这套代码最常用的集合是盒式集合加预算约束u_i u_i_nominal ± δ_i * u_i_nominal sum(|u_i - u_i_nominal| / (δ_i * u_i_nominal)) ≤ Γ第一行给出每个不确定参数的单点波动范围δ是波动幅度比例比如风电场实际出力偏离预测值最多30%第二行是整个系统层面的预算约束Γ限制同一时刻最多有几个参数同时取到极端值。Γ0时退化为确定性模型Γ1时允许所有参数同时恶化是最保守的配置。我在复现时习惯先把Γ设成不确定参数数量的三分之一跑通后再做敏感性分析。代码里有专门的参数定义区直接改Gamma_value这个变量就行。需要说明的是Γ并不是越大越好——保守度越高、装机容量越大、投资成本越高实际运行时大量容量闲置经济性反而下降。这个旋钮决定了鲁棒解和名义解之间的平衡点。2.3 机组组合约束与禁止运行区间容量配置的隐藏边界容量配置如果只算“装多少”不进去“怎么开”很容易高估系统的可用容量。常规模型认为机组出力可以在0到额定容量之间连续调节但真实机组存在最小出力限制、爬坡约束还有禁止运行区间。晏鸣宇那篇CAJ文献的核心就是含风电的鲁棒机组组合其中机组禁止运行区间指机组在某些出力区间内无法稳定运行比如一台燃气轮机可能无法在20%到40%额定出力之间连续运行。这类约束在模型里是用二进制变量实现的给状态变量z加一条逻辑约束当机组处于开机状态且出力落在禁止区间时要么增大出力越过区间要么降回最小出力以下。加入禁止运行区间后第二阶段问题的求解难度明显上升因为纯粹的线性规划变成了混合整数规划。复现代码时会看到主循环里第一阶段变量只有容量x但第二阶段子问题里带了整数变量这就是为什么求解时间比普通两阶段鲁棒问题要长。这个约束虽然让模型复杂了但很有必要忽略禁止运行区间的容量配置可能在某个场景下计算出机组能输出35%额定功率而实际机组在这个区间无法稳定运行运行成本和安全边界全部失真。所以拿到代码后第一步建议先确认这个约束是否激活再看它如何与第二阶段出力变量耦合。3. 数据文件怎么用六个CSV、一份CAJ、一份docx的对齐顺序3.1 CSV文件的角色分工与读取顺序压缩包里数据文件数量不少但角色分工很清晰。zhenjiang_power.csv和Tianchi_power.csv是两组电力系统公开数据前者是镇江地区的负荷或功率序列后者是天池电力比赛的数据集一般包含多个机组的实际出力记录。Load_history.csv是微网本地负荷历史temperature_history.csv是温度序列Benchmark.csv提供基准场景参数weights.csv存成本权重系数submission_template.csv和test.csv是预留给验证或测试阶段的格式化数据。我的读取顺序建议是先读Load_history.csv和temperature_history.csv把微网自身的负荷与温度基线建立起来再看zhenjiang_power.csv和Tianchi_power.csv用外部数据集补充风电出力或区域负荷的波动特征。最后打开weights.csv看成本项怎么定义——是燃料成本权重还是碳排放权重直接决定目标函数里的系数矩阵怎么填。在设计数据读取模块时一般是按列名索引而不是按列位置这样即使CSV列顺序调整也不会出错。一个稳妥的读取模板长这样% 读取负荷与温度历史按时间对齐 load_data readtable(Load_history.csv); temp_data readtable(temperature_history.csv); power_data readtable(zhenjiang_power.csv); % 统一时间戳格式便于后续按小时聚合 load_data.Time datetime(load_data.Time, InputFormat, yyyy-MM-dd HH:mm:ss); temp_data.Time datetime(temp_data.Time, InputFormat, yyyy-MM-dd HH:mm:ss); % 检查采样间隔常见坑两个文件一个15分钟一个1小时 fprintf(负荷数据采样间隔: %d 分钟\n, minutes(median(diff(load_data.Time)))); fprintf(温度数据采样间隔: %d 分钟\n, minutes(median(diff(temp_data.Time))));这段代码做的事情很基础但很关键先把三个文件的时间戳解析成统一的datetime类型再检查采样间隔。因为zhenjiang_power.csv和Tianchi_power.csv来自不同的采集系统采样频率很可能不一致直接横向合并会制造大量NaN。代码里Data_Align.m那个函数就是干这个的内部用retime或interp1做重采样建议不要跳过。3.2 特征工程温度如何进入风电出力模型temperature_history.csv不是摆设它至少有两层用途。第一层是负荷预测——空调和采暖负荷对温度高度敏感温度序列可以直接作为负荷模型的回归特征。第二层是风电出力的间接影响因素部分代码实现里用温度替代风速数据因为公开的风电场数据往往缺少测风塔记录而温度在多数地区与风电出力存在季节性相关。weights.csv这个文件容易被忽略但它决定了优化模型的量纲。比如权重系数配置成[1, 0.5, 0.2]含义可能是一单位投资成本、半单位燃料成本、五分之一单位碳排放成本之间的相对重要性。因此读weights.csv后最好先归一化确认权重和代码中目标函数矩阵的维度一致。代码中objective函数构造成本向量C_inv和C_op时会自动把weights.csv的系数乘上去如果维数对不上MATLAB会在矩阵运算时报维度错误——这是最常见的启动报错之一。还有一个实用技巧把Load_history和temperature_history合并后做一次相关性分析确认负荷和温度之间的滞后关系再决定是否给温度特征加滞后项。代码里的Feature_Engineer.m算是比较朴素的做法滞后阶数取3到4个小时用过去几小时温度预测当前负荷整体上比直接用同期温度效果稳定一些。3.3 参考论文CAJ的读法禁止运行区间如何变成约束CAJ文件是知网格式的文献需要用CAJViewer打开。这篇晏鸣宇的《考虑机组禁止运行区间的含风电鲁棒机组组合》建议先读第三节的模型构建部分重点看禁止运行区间的数学表达。论文里的表达通常是P_min * z ≤ p ≤ P_low * z P_up * (1 - z)或者用多个二进制变量把可行域拆成互斥区间。这种约束在MATLAB代码里对应的是机组组合约束部分注释里一般会写“禁止运行区间”或“禁区约束”。读懂这篇论文再回头看代码能理解为什么子问题里多了那么多二进制变量——它们不是冗余的而是把连续可行域切成多个子区间的必要手段。docx格式的MATLAB代码说明文档则是对照着看的教程里面应该有变量表、约束编号和参数含义。建议把docx里的参数命名和实际代码里的变量名做一个对照表因为压缩包里的代码有可能在搬运过程中改过变量名和文档不完全一致。4. CCG求解流程与MATLAB代码结构主问题、子问题、割平面4.1 CCG主循环骨架上界、下界与终止gap两阶段鲁棒优化最经典的求解方法是列与约束生成核心思路是主问题给出一个容量方案子问题找出这个方案最害怕的场景然后把该场景对应的约束加回主问题循环直到上下界收敛。代码里的主循环一般在Main_CCG.m或Run_CCG.m里结构是下面这样% CCG主循环初始化上界为无穷大下界为负无穷 UB inf; LB -inf; Gamma 4; % 预算参数控制保守度 MaxIter 20; % 最大迭代次数 GapTol 0.01; % 终止gap阈值1% iter 1; while iter MaxIter (UB - LB) / abs(LB) GapTol % 求解主问题给定已有割平面得到容量方案x和theta [x_sol, theta_sol, MP_obj] solve_MP(cut_set); LB MP_obj; % 主问题是松弛问题给出下界 % 固定x_sol求解子问题找最坏场景u和运行成本 [u_worst, SP_obj] solve_SP(x_sol); UB min(UB, c_inv * x_sol SP_obj); % 子问题给出上界 % 把最坏场景对应约束作为割平面加入主问题 cut_set{end1} generate_cut(u_worst); iter iter 1; end这段代码的逻辑很清楚主问题因为放宽了部分场景约束目标值偏小所以是下界子问题在固定容量后求最坏场景下的运行成本接受容量配置的最优性所以目标值偏大构成上界。终止条件是上下界相对gap小于1%或者迭代达到20次。实际运行时Gamma值越大需要迭代的次数越多因为更保守的配置意味着更多场景约束会被激活。如果发现20次迭代还没收敛优先检查是不是Gamma设得太大或者子问题里的整数变量导致割平面生成不稳定。另外注意UB和LB的更新顺序不能反否则会出现gap为负的诡异情况。4.2 子问题max-min的强对偶转换与大M化子问题是整个求解流程里最麻烦的环节——它内部又是一个两阶段先找最坏场景u再在该场景下求最小运行成本。直接求解这个max-min问题很困难代码里通常用强对偶定理把内层min问题转成max问题合并成一个单层max问题。转换如果不到位最常见的翻车方式是出现双线性项也就是对偶变量乘以不确定变量u两个都是变量乘积无法直接用线性求解器处理。标准解法是引入大M法把双线性项线性化代码里对应的片段是% 双线性项 lambda * u 的线性化引入辅助变量 w % 其中 u 是连续不确定变量lambda 是对偶变量连续 % 约束条件改为 % w M * z; % w -M * z; % w u M * (1 - z); % w u - M * (1 - z); % z 是引入的二进制辅助变量M 取合理上界 Constraints_SP [Constraints_SP, w M * z]; Constraints_SP [Constraints_SP, w -M * z]; Constraints_SP [Constraints_SP, w u M * (1 - z)]; Constraints_SP [Constraints_SP, w u - M * (1 - z)];这段代码本身不是完整的子问题而是线性化模板核心逻辑是用二进制变量z把lambda和u的乘积拆成四段线性约束。M的取值是个技术活太小会切掉可行解太大导致数值病态CPLEX求解时容易出现精度警告。我一般把M设成u上界的10倍再根据求解日志微调。如果不想做对偶转换另一条路是直接枚举极值场景。盒式不确定集合的极值点有限当不确定参数数量少时可以暴力枚举但参数一多就指数爆炸只适合小规模验证。压缩包里的代码走的是对偶路线这也是大规模问题最常用的处理方式。4.3 求解器选择与参数配置CPLEX、MIP gap、场景轮数这套代码底层用的是YALMIP建模加外部求解器推荐用CPLEX或Gurobi。用YALMIP的sdpsettings配置求解器参数ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 0.001; % MIP gap 0.1% ops.cplex.mip.tolerances.integrality 1e-6; % 整型变量容差 ops.cplex.parallel 4; % 并行核数按机器调整 ops.cplex.timelimit 3600; % 单次求解时限单位秒mipgap设0.001可以在求解速度和精度之间取得较好平衡如果设太严第二阶段问题的二进制变量会让单次求解动辄数小时。parallel参数按机器核心数设置4核是起步值。timelimit建议设3600秒防止单次子问题把总循环拖死——我在复现时遇到过子问题求解超过两小时的情况最后发现是某个割平面导致数值病态求解器在里面空转。运行整份代码的方式有两种。一种是在MATLAB图形界面里直接打开主脚本F5运行另一种是命令行批量跑方便后续调参matlab -batch run_CCG;参数修改集中在脚本顶部的Config块系统容量上下限、储能容量离散步长、风电渗透率、不确定度delta等。每次修改参数后建议清空工作区再运行防止旧变量污染新计算。5. 复现避坑这份代码最容易翻车的五个位置5.1 现象一第一次跑就报“Infeasible problem”第一次运行主脚本求解器直接返回不可行连下界都算不出来。原因通常是第一阶段容量变量和第二阶段运行约束之间存在不可满足组合。最常见的是风电装机容量上限设置太低而第二阶段强制要求满足固定负荷风不够、其他机组容量又补不上的时候问题直接不可行。另一个高频原因是数据文件里存在NaN读取进来后约束变成无效不等式。解决的办法是先给模型加一个切负荷变量让系统在最坏情况下允许少量甩负荷并给切负荷设置高昂罚成本。这样问题从“必须满足负荷”变成“尽量满足负荷”可行性大幅改善。同时在读取CSV后统一用fillmissing填充或剔除含NaN的行。5.2 现象二子问题对偶转换后出现无界解加了割平面进入下一轮迭代后子问题的目标值突然变成负无穷求解器提示“Unbounded objective”。这几乎都是大M参数设置不当引起的。M太大导致CPLEX的数值容差无法精确区分0和1对偶变量在边界处振荡M太小又会让原本可行的双线性组合被错误剪枝。另一个原因是引入了冗余的自由变量比如对偶变量漏掉了非负约束。解决方法是先做一次量纲分析u的取值区间是0到1之间还是0到100之间对偶变量的量级是多少M取u上界的5到10倍最稳妥。排查时可以暂时固定u为名义值如果子问题仍然无界说明问题不在双线性项而是对偶约束本身有缺漏——检查互补松弛条件是否漏写。5.3 现象三禁止运行区间约束不起作用优化结果中机组出力正好落在禁止区间内甚至完全没有避开的意思。原因大概率是代码里禁止运行区间约束和机组最小开机时间约束被分开写了但两者共用同一个状态变量时逻辑出现了漏洞。比如机组在禁止区间内本来应该强制切换状态但代码里只写了区间约束忘记把区间状态和启停状态做绑定。更隐蔽的是约束里有小数误差机组出力区间边界是[30, 50]兆瓦但求解器算出的出力是29.9999恰好漏过约束。解决方法是把禁止运行区间统一写成互斥形式声明两个互斥的区间状态保证任何情况下最多一个状态生效。同时把边界值处理成绝对容差——约束写成p ≤ 30 - 1e-4或p ≥ 50 1e-4不给数值误差留缝隙。5.4 现象四结果对Γ不敏感怎么调都一样把Gamma从0调到10容量配置结果几乎不变看起来鲁棒优化没有起作用。这说明不确定集合的边界设置有问题。最常见的原因是delta设得太小比如风电出力波动只有5%那么Γ再大集合覆盖的场景和名义场景差别也有限。另一个原因是不确定参数只进了目标函数没进约束——第二阶段约束里根本没有风电出力的影子那么最坏场景再坏也影响不了结果。解决方法是先检查不确定参数是否真实耦合到约束里比如发电平衡约束P_wind P_gas P_storage P_load确认P_wind确实来自不确定变量。然后做一个单参数敏感性测试固定其他参数把delta从0.05扫到0.3画出装机容量随delta变化的曲线曲线没有单调变化就说明模型哪里断了。5.5 现象五两个CSV时间轴错位配置结果偏到离谱配置出来的储能容量是常理推断的两倍以上怎么看都不对。这几乎都是时间戳问题。zhenjiang_power.csv的时区可能和Load_history.csv不一致或者一个用UTC一个用本地时间错位8小时后风电出力和负荷峰值的相关性被彻底打乱。模型为了“稳妥”只能靠堆容量来兜底结果自然离谱。解决方法是先画时间序列对比图把两个数据源的前100小时曲线叠加人工确认峰值是否对齐。再做一个跨相关分析如果最佳滞后时间不是0小时说明时间坐标系需要校准。最后统一转换为UTC时间戳或统一到北京时间再进模型这一步做完再去配置容量结果通常会回归正常区间。6. 验证鲁棒解蒙特卡洛回测才是最后一道关代码跑通、结果收敛不代表这份容量配置就真的可靠。我拿到这套代码后做的最重要一件事是把鲁棒解和名义解同时丢进蒙特卡洛回测里比谁更抗造。方法很简单用历史数据的统计特征生成几百组随机场景每组场景里让所有机组的实际出力和负荷都偏离预测值然后分别用两套容量配置去计算运行成本、切负荷量和弃风量。% 蒙特卡洛回测对比鲁棒解与名义解 N_sim 500; cost_robust zeros(N_sim, 1); cost_nominal zeros(N_sim, 1); shortfall_robust zeros(N_sim, 1); for k 1:N_sim % 从历史数据中抽样生成一组随机场景 u_sim generate_scenario(load_data, temp_data); % 固定容量配置运行仿真模块求实际成本 [cost_robust(k), shortfall_robust(k)] simulate(u_sim, x_robust); cost_nominal(k) simulate(u_sim, x_nominal); end fprintf(鲁棒解平均成本: %.2f, 切负荷概率: %.2f%%\n, ... mean(cost_robust), 100 * mean(shortfall_robust 0));回测观察两个指标平均运行成本和切负荷概率。鲁棒解通常平均成本略高但切负荷概率会显著低于名义解。如果两者成本差超过20%说明Gamma设大了容量配置过度保守如果切负荷概率没有实际下降说明不确定集合本身构建有问题。从那以后我每次拿到微网规划的代码都不急着跑主程序而是先花半小时把数据对齐、名义解校准、蒙特卡洛基线这三件事做完再上两阶段鲁棒循环。这个习惯帮我筛掉过不少表面光鲜、实际脆弱的模型也省下了大量调参时间。希望帮到你。本文还有配套的精品资源点击获取