简介面向微电网规划与电力系统优化领域的研究者这套资料围绕含风电、光伏、储能及燃气轮机的微网电源容量配置问题给出基于两阶段鲁棒优化算法的完整MATLAB实现。压缩包共19个文件大小3.76MB包含5个.m源码文件、6个docx建模与推导文档、2个pdf及1个caj参考论文、1个典型日数据xlsx以及3张矩阵推导截图文件类型覆盖代码、文档、数据与图示。目前已有1328人学习。资源从两阶段鲁棒构建过程、CCG列与约束生成算法求解到MP/SP主子问题代码均有细致呈现并配有实现效果截图与readme说明可帮助读者厘清一阶段容量决策与二阶段运行出力之间的迭代关系可直接参考用于学术复现或实际微网规划项目。1. 两阶段鲁棒优化算法在微网电源容量配置中的切入点微网电源容量配置最直接的做法是给定风、光、负荷典型曲线通过单层优化求解风电、光伏、储能和燃气轮机容量。但这样得到的方案无法回答一个关键问题当风电出力低于预测、光伏被云层遮挡、负荷又同时攀升时容量是否仍然够用。两阶段鲁棒优化算法把容量决策放在第一阶段把最坏不确定场景下的运行调度放在第二阶段用min-max-min结构的CCG算法迭代求解。这个项目提供的正是这样一套完整代码包含MP、SP、addC等核心文件适合微网规划、鲁棒优化和综合能源方向的研究者做基线参考。2. 不确定集与两阶段决策框架容量配置模型的矩阵化2.1 一阶段容量决策与二阶段运行调度的变量边界两阶段鲁棒优化必须先明确哪些变量在一阶段确定哪些在二阶段调整。对于微网电源容量配置一阶段变量是风电、光伏、储能和燃气轮机的安装容量这类变量在不确定性实现前就要敲定二阶段变量是各个时段的实际出力、充放电功率、切负荷量以及储能SOC这类变量可以在看到实际风、光、负荷后进行实时调度。这个程序里一阶段变量x可以写成4×1向量风电机组装机容量、光伏装机容量、储能额定功率、燃气轮机装机容量。二阶段变量y的维度要按时段展开例如24时段下至少包含风电出力、光伏出力、储能充电功率、储能放电功率、燃气轮机出力、切负荷量以及储能SOC递推中的每个时段能量状态。把变量分层的好处是后续用列约束生成时主问题只需要管容量和已发现的最坏场景子问题再对不同容量的运行成本做评估。在模型文档中目标函数通常写作min_x { c_inv * x max_{u∈U} min_{y∈F(x,u)} c_ope * y }其中c_inv是单位容量投资成本c_ope是单位运行成本。去掉max-min部分就退化为确定性规划加上后就变成一个三层结构不能直接用求解器一步解出必须分解为主问题MP和子问题SP。下面的符号表是后续建模的基础。符号含义维度x一阶段容量决策4×1y二阶段出力决策(4T1)×1u归一化不确定变量T×3Γ不确定预算1×1U不确定集由u的约束定义F(x,u)给定x和u的可行域多面体2.2 风电光伏出力与负荷的不确定集构建不确定集决定了鲁棒优化在多大范围内“防患于未然”。微网中风电、光伏出力可由预测曲线加偏差区间描述负荷则围绕预测值上下波动。项目采用盒式不确定集加预算约束表达式如下U { u_w(t), u_pv(t), u_load(t) ∈ [0,1] : Σ_t (u_w(t) u_pv(t) u_load(t)) ≤ Γ }其中u_w(t)1时对应风电出力取区间下界u_load(t)1时对应负荷取区间上界。Γ用于限制最坏情况同时出现的总次数避免每个时刻都对抗不确定性而过度保守。一个典型参数组合是风电偏差δ_wt取15%光伏偏差δ_pv取15%负荷偏差δ_load取10%Γ取6到8。四个典型日数据分别代表冬季、夏季、过渡季和极端天气针对不同典型日可以单独调整偏差比例。下面是构建不确定集的MATLAB函数function U build_uncertainty_set(P_forecast, delta, Gamma, T) % P_forecast: T*3矩阵[风电预测; 光伏预测; 负荷预测] % delta: [风电偏差, 光伏偏差, 负荷偏差] % Gamma: 不确定性预算0表示退化为确定性 U.w_up P_forecast(:,1) .* (1 delta(1)); U.w_lo P_forecast(:,1) .* (1 - delta(1)); U.pv_up P_forecast(:,2) .* (1 delta(2)); U.pv_lo P_forecast(:,2) .* (1 - delta(2)); U.load_up P_forecast(:,3) .* (1 delta(3)); U.load_lo P_forecast(:,3) .* (1 - delta(3)); U.Gamma Gamma; end这段代码把每个时段的预测值按比例展开成上下界Gamma作为全局参数存到结构体U中。在SP子问题里不确定变量u会被限制在这六个上下界之间并通过预算约束关联。一个常见误用是直接让u取0或1这会丢失中间值信息实际上第二阶段目标函数通常线性最优解天然在顶点处因此最终取到的还是0或1但连续化处理有助于求解器保持稳定。2.3 目标函数与约束的矩阵形式两阶段鲁棒规划中约束矩阵的组织方式决定了CCG能否顺利迭代。功率平衡是最主要的等式约束各电源出力之和减去切负荷量等于负荷。储能部分有SOC递推约束和充放电功率约束燃气轮机有出力上下限和爬坡约束。这些约束写成紧凑形式为A_y * y ≤ b_y B_x * x B_u * u其中A_y对应二阶段变量之间的耦合系数B_x把上限容量乘到出力变量上B_u把不确定量放入右端项。项目内矩阵推导图就是在把上面的等式和不等式逐条整理成该标准形。整理时最容易出问题的是储能SOC递推SOC(t) SOC(t-1) η_c * P_ess_c(t) - P_ess_d(t) / η_d如果漏了SOC变量就会导致子问题无界。建议写代码前先确认每个约束的右端项是否含x或u。功率平衡右端项只含负荷不含x但风电和光伏出力上限由容量x约束因此要把这类约束的右端项整理成由x控制的项。例如P_wt(t) ≤ P_wt_cap * C_wt(t)这里的C_wt(t)是风电归一化出力系数。经过这种整理SP就是一个标准LP可以直接用强对偶求解。3. CCG求解MP与SP交替迭代的列约束生成实现3.1 为什么选择CCG而不是Benders分解两阶段鲁棒优化的经典求解算法是Benders分解它通过对偶变换把内层min问题转化为割平面加入主问题。但Benders在处理min-max-min结构时每次迭代只增加一个约束外部max部分依赖极点枚举收敛速度慢且上下界可能震荡。CCG则直接把最坏场景对应的全部二阶段变量和约束加入主问题相当于把原问题不断“外延扩大”理论上迭代次数等于需要辨识的极端场景数工程上通常10次以内即可收敛。MP和SP的分工很清晰主问题负责决策容量x并给出下界子问题固定x后寻找使运行成本最大的不确定场景并给出上界。下面表格列出两者在实际编程中的差异便于对照项目代码。对比项CCG主问题CCG子问题求解目标容量x和辅助变量η最坏场景u和运行成本数值作用提供下界提供上界变量规模随迭代增长固定约数百个核心难点动态加列max-min线性化3.2 子问题SP.m的max-min线性化过程SP.m接收主问题传来的x求解max_{u∈U} min_{y∈F(x,u)} c_ope y。内层min是一个以y为变量的线性规划强对偶成立时可以将其转化为对偶变量的max问题再与外层max合并。最终SP变成单层max问题决策变量包含u、y的对偶变量π目标函数变为π (b_y B_x * x B_u * u)约束为对偶可行域A_y π ≤ c_ope以及π ≥ 0。对偶转化的前提是内层min可行且目标函数与约束没有不可控的非线性。对于本项目中的功率平衡、储能SOC和出力上下限全部满足LP条件。SP.m核心实现如下function [worst_u, obj_sp] solve_SP(x, U) % 输入x为当前容量配置 % 输出worst_u为最坏场景obj_sp为子问题目标值 % 使用强对偶将max-min转为max问题 cvx_begin variable u(T,3) variable pi(n_cons) maximize( pi * (b_y B_x * x B_u * u(:)) ) subject to A_y * pi c_ope; pi 0; sum(u(:)) U.Gamma; 0 u 1; cvx_end worst_u reshape(u, T, 3); obj_sp cvx_optval; end代码中A_y、B_x、B_u与第2.3小节的矩阵对应。参数n_cons是第二阶段不等式约束的数量需要与A_y的行数一致。如果子问题无解常规做法是在SP中加入人工切负荷变量并设置很高的惩罚系数这样既能保证可行性又能通过惩罚成本暴露容量不足的时段。调试时还会把u的初值设成全0观察目标值是否等于确定性场景以此检查对偶方向是否正确。3.3 主问题MP.m与addC.m的动态加列实现主问题的初始模型只包含一个初始场景通常取预测场景u0。迭代一次后SP返回一个最坏场景u*main.m调用addC.m把对应的一组y变量和约束加入主问题。主问题形式如下min_{x, y_i, η} c_inv x η s.t. η ≥ c_ope y_i, i 1...k A_y * y_i ≤ b_y B_x * x B_u * u_i*, i 1...k x ∈ X, y_i ∈ Y在MATLAB中addC.m通过YALMIP的sdpvar动态拼接变量。代码片段如下function MP addC(MP, u_new, params) y_new sdpvar(params.n_y, 1); MP [MP, params.eta params.c_ope * y_new]; MP [MP, params.A_y * y_new params.b_y ... params.B_x * params.x params.B_u * u_new]; endmain.m里的循环控制上下界收敛。下界是主问题目标值上界是当前x对应的投资成本加上SP目标值。当gap小于1e-3时停止迭代。如果迭代次数超过15次仍不收敛需要检查SP返回的最坏场景是否与上一轮重复重复说明主问题已经包含足够场景但上下界仍未闭合通常是因为主问题中的辅助变量η没有和y_i正确耦合或SP目标值漏加了投资成本。4. 四个典型日数据如何驱动容量配置结果4.1 从Excel到MATLAB的典型日读取四个典型日数据.xlsx是模型输入的关键。工程上通常把表结构设计为每个Sheet代表一种设备每列是一个典型日的24小时数据这样读取时不需要反复切换维度。下表给出一种常见的表格布局具体数值随项目参数变化时段h负荷典型日1负荷典型日2风电典型日1光伏典型日11120013500.3502115013000.400...............24122013800.300读取代码如下data readtable(四个典型日数据.xlsx); T 24; P_load data{1:T, 1:4}; % 四列负荷 P_wt data{1:T, 5:8}; % 四列风电归一化出力 P_pv data{1:T, 9:12}; % 四列光伏归一化出力这里的P_wt和P_pv是归一化到额定容量的出力系数不是绝对功率。在目标函数和约束中实际出力等于容量乘以该系数。如果表格中给出的是绝对功率则需要先除以规划基准容量否则不确定集和容量约束会冲突。4.2 典型日与不确定性预算的交叉测试四个典型日覆盖了微网全年负荷和新能源出力的不同形态。选择典型日的依据往往是最冷日、最热日、过渡季和极端新能源日。规划时通常对每个典型日单独做一组Gamma测试形成下面的表格典型日负荷峰值(kW)风电容量(MW)光伏容量(MW)储能容量(MW)燃气轮机容量(MW)日1145018.29.54.210.8日2160016.011.05.011.5日3120014.58.53.09.0日4180020.06.06.513.0从趋势看日4这种极端场景会显著抬高储能和燃气轮机容量而日3过渡季则允许系统依靠更多的光伏和较少的储能。这个结果符合鲁棒优化的直觉最差场景决定了容量下限常态场景只影响经济性排序。增加Gamma后不确定集变大SP会找到更严格的最坏场景主问题被迫扩大容量。一般规则是Gamma每增加2总成本上升约3%~6%。如果成本增幅过大说明偏差δ取值偏大需要根据历史预测误差重新统计。4.3 从容量配置结果反推运行瓶颈拿到容量配置结果后直接入库容易错过潜在问题。我习惯把x固定回二阶段调度中检查每小时功率和SOC曲线。这个步骤类似半仿真用于发现两个常见问题一是储能容量绰绰有余但功率上限限制放电导致晚间峰值要靠燃气轮机硬扛二是光伏容量大但中午无法消纳弃光率超过20%说明储能容量应该增大或者燃气轮机最小出力限制导致无法调峰。具体做法可以复用SP.m中的约束把x固定为最优解令u取预测场景然后求解确定性调度LP。如果此时出现切负荷说明容量配置实际不足如果燃气轮机在大部分时段贴近上限说明需要增加储能或对峰值负荷进行转移。这个验证结果配合SOC曲线图比单纯贴CCG迭代次数更有说服力。5. 扩展机组禁止运行区间与addC.m的鲁棒机组组合改造5.1 禁止运行区间的0-1建模燃气轮机的运行区间并非完全连续部分机组在某个出力范围存在振动过大、效率骤降等禁止运行区。参考文件提到的含风电鲁棒机组组合就是在两阶段模型中加入0-1变量来避开这些区间。假设燃气轮机有m个禁止区间每个区间上下界为[L_j, U_j]则需引入二进制变量z_j(t)Σ_j z_j(t) 1, z_j(t) ∈ {0,1} P_gt(t) ≤ L_j M * (1 - z_j(t)) P_gt(t) ≥ U_j - M * (1 - z_j(t))这里M取一个足够大的正数保证z_j0时对应约束自动松弛。实际工程中禁止区间通常不是绝对不能用而是避免长时间停留所以也可以把约束写成“最多连续运行k小时”但0-1区间拆分更直观。5.2 addC.m如何扩展包含二进制变量的场景加入0-1变量后主问题中的每个场景都要额外维护一组z变量addC.m不再只添加连续变量。同时子问题内层min从LP变为MILP强对偶失效。处理这种扩展有两条路线第一种是继续使用SP对偶但把0-1变量作为参数在SP外层枚举第二种是直接使用KKT条件把内层min用KKT互补松弛条件等价转化形成的模型是混合整数非线性规划一般要用Gurobi的bilinear或big-M线性化。实际扩展时我偏向于第一种路线把0-1变量决策从SP剥离在每个不确定场景u下先枚举可行的开停机组合再对每个组合求解LP取其中最优值作为该场景的子问题目标。修改后的addC.m示意如下function MP addC(MP, u_new, z_init, params) y_new sdpvar(params.n_y, 1); z_new sdpvar(params.n_gt, T, full); MP [MP, params.eta params.c_ope * y_new]; MP [MP, params.A_y * y_new params.b_y ... params.B_x * params.x params.B_u * u_new]; for t 1:params.T MP [MP, sum(z_new(:,t)) 1]; MP [MP, P_gt(t) L params.M * (1 - z_new(1,t))]; MP [MP, P_gt(t) U - params.M * (1 - z_new(2,t))]; end end通过把z_new加入主问题CCG主问题从LP变成MILP但规模仍然可控。每一轮新增的z变量数量是机组数×时段数比如2台机组×24小时增加48个0-1变量对Gurobi来说压力不大。5.3 求解规模控制与参数收敛经验加入禁止运行区间后求解时间会成倍上升。下表给出不同规模下的经验配置模型规模求解器间隙设置参考时间24时段纯LPlinprog/intlinprog1e-4秒级24时段含禁止区间Gurobi/Cplex0.5%分钟级8760时段CCGGurobi/Cplex1%小时级时间只是参考与机器配置和矩阵稀疏度相关。一般不建议在容量配置阶段就加入所有禁止区间约束而是先用纯LP鲁棒模型确定容量大致范围再加入禁止区间对燃气轮机容量进行微调。这样两步走的方案比一次性求解MILP更快且容量结果差异通常在5%以内。6. 用场景回带法验证CCG容量配置可靠性CCG收敛后最容易被忽视的一步是验证容量配置在外层不确定集中是否真正可行。场景回带法通过固定容量x遍历不确定集代表性的极端场景逐一求解确定性调度问题以此检验模型有没有漏约束。我的做法是构造M2T1个极端场景。前T个场景分别让每个时段的风电取下界、光伏取下界、负荷取上界其余时段取预测值第T1个场景让所有时段同时达到最坏组合最后T个场景让每个时段的偏差组合达到预算上限。对这M个场景调用一个纯调度函数solve_dispatch(x_curr, u_test)返回最优运行成本或不可行标志。核心验证代码如下for m 1:M u_test generate_extreme_scenario(U, T, m); [~, exitflag] solve_dispatch(x_curr, u_test); if exitflag 1 fprintf(场景%d不可行\n, m); end end如果所有场景都可行再对比CCG的UB与这M个场景调度成本的算术关系。正常情况下UB应该等于M个场景中最大的调度成本加投资成本。如果UB明显小于该最大值说明SP中定义的U与验证时用的generate_extreme_scenario不一致常见原因是SP中的u被限制为0到1连续而验证时直接取了端点值并让负荷同时取上界超出了预算约束。修正的方法是让生成场景的函数也遵守同样的预算约束。这个方法还能顺带发现储能约束问题如果某个极端场景下SOC越界导致不可行需要检查二阶段中的SOC上下界是否乘以了容量x。如果x固定后SOC上限初始化成了固定值而不是x的倍数就会在容量增大时出现“容量利用不上”的奇怪结果。把回带验证加入常规流程后CCG迭代结果才真正具备可交付性。本文还有配套的精品资源点击获取