1. 微网优化调度的本质不确定性才是核心矛盾做微网调度的人都知道真正的难题不是“优化”而是“算得过老天爷的脸色”。光伏出力跟着云走风电随风变负荷说涨就涨这些不确定性才是让调度员头疼的根源。传统的确定性优化模型直接给预测值一把梭算出来的方案在预测值下最优但一旦实际出力偏离预测值要么弃风弃光、要么切负荷甚至可能引起电压越限。今天我们聊的两阶段鲁棒微网优化调度就是冲着这个问题去的。两阶段鲁棒优化不是新名词但要用好它并不容易。它的核心思想是先在“在这里且现在”做决策比如机组启停这些需要提前定死的动作等不确定参数揭示后再在“等待并观望”阶段做调整比如储能充放电、机组出力微调。整个模型要保证在最恶劣的场景下系统仍然能安全运行并经济调度。而“关键场景辨别算法”则是在这个框架下专门用来找出那个“最恶劣场景”并快速求解的手段。这篇文章我打算从模型构建、算法原理、Matlab实现几个角度把这套方法从头到尾捋一遍让有兴趣做微网调度、鲁棒优化、不确定性研究的同学能直接拿代码上手跑。这套东西适合谁如果你是正在做微电网能量管理、虚拟电厂调度、或者研究鲁棒优化应用的研究生和工程师这篇文章值得收藏。我会尽量把建模细节讲清楚把代码里的坑指出来毕竟这些坑我都是亲测过的。2. 两阶段鲁棒微网调度模型到底怎么建2.1 从确定性模型到两阶段鲁棒模型的进化逻辑先看一个常规的微网日前调度问题。我们要在一天96个时段或24个时段内决定微燃机的启停和出力、储能的充放电功率、与大电网的交换功率目标是使总运行成本最小包括燃料成本、购电成本、储能折旧、弃负荷罚项等。约束包括功率平衡、机组出力上下限、爬坡约束、储能SOC连续性和容量约束、联络线功率限制等。这是一个典型的混合整数线性规划MILP问题。但这种模型有个致命弱点光伏预测值给死算出来的方案在真实场景中可能完全不适用于次。解决办法之一是引入随机规划stochastic programming用大量场景求期望成本最小化。但随机规划假设概率分布精确已知且需要处理大量场景计算量感人。解决办法之二就是鲁棒优化不假设概率分布只用一个不确定集合来描述参数的可能范围然后保证在这个集合内所有可能实现下方案都可行且成本可控。两阶段鲁棒优化比单阶段鲁棒更灵活因为不确定性参数的实现发生在第二阶段第一阶段决策不需要对所有场景都可行只需要保证存在一个第二阶段调整方案来“兜底”。这种“这里-现在”和“等待-观望”的架构特别适合微网调度因为机组启停和联络线合同功率必须提前定而储能充放电、机组出力可以在日内滚动调整。2.2 不确定集合的具体构造构造不确定集合是鲁棒模型的第一步。最常见的是盒式集合box uncertainty set这也是最保守的但建模最简单的。以光伏出力为例设预测值为 P_pv_pred(t)最大偏差为 ΔP_pv(t)那么实际出力被约束在P_pv(t) ∈ [P_pv_pred(t) - ΔP_pv(t), P_pv_pred(t) ΔP_pv(t)]同理负荷可以写成 P_load(t) ∈ [P_load_pred(t) - ΔP_load(t), P_load_pred(t) ΔP_load(t)]。光有盒式集合还不够。如果所有时段的光伏都同时取最小值、负荷都同时取最大值这种“极端情景”现实中几乎不会发生但鲁棒优化会为此付出巨大的成本代价。为了调节保守度我们引入预算约束budget constraint限制不确定参数和预测值偏离的总量。例如Σ_t |P_pv(t) - P_pv_pred(t)| / ΔP_pv(t) ≤ Γ_pv其中 Γ_pv 是预算参数整数或连续。Γ_pv 越大模型越保守抗风险能力越强但经济性越差。这个参数直接对应着“你愿意花多少成本去防一个多极端的坏日子”。具体取值往往需要做灵敏度分析后面我会讲我在实际调参时的经验。2.3 两阶段模型的标准数学形式记第一阶段变量为 x包括各机组启停状态 u(t)、启动/停机标志、与大电网的购电合同功率 P_grid_purchase(t)、以及储能的日前充放电计划也可以放到第二阶段。第一阶段目标是在最小化当前已知成本的同时为第二阶段留出可行性。不过在两阶段鲁棒的标准写法中第一阶段目标只包含与x相关的固定成本第二阶段目标包含运行成本。第二阶段变量为 y包括各机组实际出力、储能实际充放电、切负荷量等。对于每一个不确定参数实现 ξ第二阶段都需要找到一个可行的调整方案 y(ξ)使得总成本最小。因此模型写成min_x max_ξ min_y [ c^T x d^T y ]s.t. A x ≤ b G y ≤ H x M ξ y ∈ Y内层 min_y 是针对给定 x 和 ξ 的运行调度问题max_ξ 是在不确定集合内寻找使内层最优值最大的最恶劣场景。这就是“两阶段鲁棒”的形式。用行话讲这是一个三级优化tri-level问题直接没法用商业求解器解。需要分解为主问题MP和子问题SP迭代求解。3. 关键场景辨别算法的原理怎么从无限场景里抓出那个“最坏蛋”3.1 为什么不能直接“枚举”所有场景两阶段鲁棒模型的最大障碍在于不确定集合一般是连续的、无限维的。你不可能把每个可能的光伏出力都拿出来算一遍。传统的Benders分解思路是把max-min问题对偶成max问题然后得到子问题的对偶解作为割平面信息加入主问题。CCG列与约束生成算法则是更聪明的一种做法每次迭代时求解子问题找出当前最恶劣的场景把这个场景对应的约束显式加入主问题。随着迭代进行主问题的复杂度不断增加但最终收敛到一个鲁棒最优解。CCG的迭代次数一般不多十几次以内但每次求解子问题都需要解一个max-min问题如果系统规模大、离散变量多照样很慢。关键场景辨别算法Key Scenario IdentificationKSI本质上是对CCG的一种加速和增强。做法是在每次迭代中不只求一个最优的最恶劣场景而是按照某种策略从不确定集合中“辨别”出若干个关键场景一次性加入主问题从而减少迭代次数。3.2 关键场景辨别算法的核心策略我的实现中采用了“双层筛选”的思路第一层是“预筛”。在迭代开始前用拉丁超立方抽样LHS在不确定集合内生成几百个候选场景再用K-medoids聚类比如聚成20~30类得到一批代表性场景。这些场景不一定是最恶劣的但覆盖了不确定集合的主要分布区域。把这个预筛场景集加入主问题提供一个“保底”的鲁棒性也能让主问题有个好的初始解。第二层是“精辨识”。在每次CCG迭代中求解子问题获得当前最恶劣场景。如果这个场景已经包含在主问题中说明收敛否则将其加入主问题。为了进一步提升效率我还会比较当前最恶劣场景与已有场景的“相似度”如果它们对应的不确定变量取值非常接近例如偏差向量归一化后距离小于阈值则直接合并不重复添加。这个策略在实际计算中能减少约20%~30%的迭代次数。说白了关键场景辨别算法没有改变两阶段鲁棒模型的本质但它像个“精明的采购员”既不盲目把每种可能都买一遍又能一眼看出哪几种“极端情况”必须提前准备。3.3 提高算法稳定性的几个细节实际调试这个算法时有几个细节一定要留意。一是子问题的求解精度。如果你的max-min子问题是用对偶KKT转成单层优化那么对偶间隙和互补松弛条件的处理要小心否则容易产生假的“最恶劣场景”。二是主问题中的约束规模会随迭代增长MILP求解时间会越来越长所以建议对主问题设置一个时间上限比如30秒不要让MILP一跑就是半小时。三是预算参数Γ的取值直接决定关键场景的“位置”我们一般会让不确定集合在“盒式预算”的框架下保持紧凑否则预筛阶段就变得过于保守。4. Matlab代码实现从建模到求解的完整步骤4.1 建模工具箱与环境准备我这边用的环境是Matlab R2022a2023b也行有些函数更友好加上Yalmip工具箱和Cplex或Gurobi求解器。推荐GurobiMILP求解效率明显更稳。你如果在学校服务器上一般都有License如果没有可以先用Cplex的免费学术版或Gurobi的学术License。Yalmip的安装很简单解压后路径加到matlabpath即可。代码里的solver选择可以用ops sdpsettings(solver,gurobi,verbose,2,gurobi.MIPGap,1e-4,gurobi.TimeLimit,120);注意gurobi.TimeLimit是设置单次求解主问题的时间上限而整个两层迭代的时间可能需要几轮。如果不设置这个上限我见过某些算例在主问题第一次迭代就卡了3小时不收敛根本没法用。4.2 数据结构准备我把不确定性参数单独存储在一个结构体里。以24小时为例numHour 24; pvPred [0, 0, 0, 0, 0.2, 1.2, 2.8, 4.1, 5.3, 6.0, 6.2, 5.8, 4.9, 3.8, 2.7, 1.6, 0.8, 0.3, 0, 0, 0, 0, 0, 0]; % 光伏预测出力(MW) pvDev 0.15 * pvPred; % 偏差量取预测值的15% loadPred [3.5, 3.2, 3.1, 3.0, 3.0, 3.2, 4.0, 5.0, 5.5, 5.8, 6.0, 6.1, 6.0, 5.8, 5.5, 5.2, 5.0, 4.8, 4.5, 4.2, 4.0, 3.8, 3.6, 3.5]; loadDev 0.05 * loadPred; budgetPV round(0.3 * numHour); % 光伏的预算参数 budgetLoad round(0.2 * numHour);偏差系数取多少直接影响保守程度。光伏取10%~20%风电取20%~30%负荷取5%~10%比较常见。预算参数取整数值表示最多允许多少个时段偏离到边界。4.3 主问题MP的Yalmip实现主问题中需要定义第一阶段变量。我这里微网配置是2台微燃机MT1、MT2、1台储能ESS、光伏、风电、与大电网的联络线。第一阶段变量包括MT的启停状态u24×2个二进制变量、购电状态/合同功率这里简化为连续变量P_grid_buy一天中每个时段可取范围[0, 5]MW以及储能日前充放电计划P_ess_day连续变量正为充电负为放电。Yalmip代码片段u binvar(numHour, numMtc, full); P_mt sdpvar(numHour, numMtc); S sdpvar(numHour, 1); % 储能SOC P_ess_day sdpvar(numHour, 1); P_grid_buy sdpvar(numHour, 1); P_grid_sell sdpvar(numHour, 1); % 可能也可以卖电给电网 Constraints []; % 功率平衡 for t 1:numHour Constraints [Constraints, sum(P_mt(t,:),2) P_grid_buy(t) P_pv_term(t) P_wt_term(t) P_ess_day(t) - P_grid_sell(t) loadTerm(t)]; end其中P_pv_term(t)和loadTerm(t)不是定值而是在关键场景加入后变成该场景下的具体参数。这也是CCG的巧妙之处主问题中不确定参数是通过一个个“场景”来体现的每个场景对应一组确定的光伏、负荷值。主问题需要同时满足所有已加入场景的约束。因此我们预定义一个“场景集合”ScenSet初始为空或预筛的K个场景然后逐步扩展。4.4 子问题SP的max-min转换子问题要解决的是给定第一阶段决策u, P_mt的日前部分实际上第二阶段也有P_mt但连续和不确定参数ξ求微网运行成本的最小值然后外层再求这个最小值的最大值。内层min问题是LP如果P_mt、储能等是连续变量可以直接用LP对偶转换成max问题然后在外部再统一变成max也就是对偶后变成单层max问题。但是内层min中如果包含二进制变量例如机组启停已经在第一阶段固定所以内层没有整数变量就能处理。另一种更直接的方法是使用KKT条件把内层LP的KKT系统作为约束得到一个单层MILP但KKT里有互补松弛通常引入大M量转化为线性约束。我个人的经验用强对偶转换更干净因为对偶后的目标函数和约束都是线性的加上不确定变量的预算约束可以直接交给求解器。子问题的标准形式对给定的 x第一阶段决策求解 max_ξ min_y { d^T y λ^T (M ξ - G y H x) }内层min_y 是标准的LP写出其对偶形式合并外层max就得到一个单层的max问题变量包括原对偶变量和不确定变量ξ。目标中含有 d^T y 的等价形式约束包含原内层约束、对偶可行条件、强对偶等式、预算约束等。Yalmip里实现时我并不喜欢手动推导每一步对偶而是先写出内层LP的Yalmip约束和目标然后用dualize函数自动对偶。不过Yalmip的dualize在某些复杂约束下不够稳定所以我干脆用纸笔推导一次写成矩阵形式。虽然推导过程繁琐但一旦写好跑起来很稳。4.5 迭代求解主问题-子问题交互伪代码如下% 初始化 ScenSet []; for iter 1:maxIter % 求解主问题得到第一阶段决策 x* solve MP with ScenSet LB value(objective_MP) % 求解子问题得到最恶劣场景 xi* 和第二阶段最优值 obj_SP solve SP with x* UB c*x* obj_SP if UB - LB epsilon break end % 将场景 xi* 加入 ScenSet ScenSet [ScenSet; xi*]; % 可选预筛关键场景补充 end注意这里LB和UB的更新方式。主问题因为考虑的约束集合是当前已加入场景的子集所以求出的目标值是最优解的下界LB。子问题给定x*求最恶劣场景下的最优成本加上第一阶段的成本是上界UB。当UB和LB差值小于收敛阈值比如1e-3时说明主问题中的决策已经对所有可能场景都可行且成本差别可忽略。实际工程中建议设置双重收敛条件间隙小于绝对阈值或相对阈值比如0.5%同时迭代次数不超过30。到了后期每次加入的场景带来的UB下降会越来越小但主问题规模越来越大再硬算不划算。4.6 储能的建模细节储能是微网调度的灵魂。我通常把储能SOC表示成连续变量充放电效率分别考虑SOC(t) SOC(t-1) - P_ess(t) * eta_dis * (P_ess(t)0) - P_ess(t) / eta_chg * (P_ess(t)0)这个非线性太麻烦我用一个离散化技巧把充放电拆成两个非负变量P_chg和P_dis再加互补约束P_chgP_dis0但这个约束也麻烦。工程上就直接保留两个变量并加约束0≤P_chg≤P_chg_maxu_ess0≤P_dis≤P_dis_max*(1-u_ess)其中u_ess是二进制变量表示充电/放电状态。虽然会引入整数变量但MILP求解器处理这个很轻松。SOC约束SOC_min SOC(t) SOC_max SOC(24) SOC_init % 保证日循环最后一条很重要否则储能会在末尾卖光电第二天就得从零开始不符合实际运行。4.7 求解器参数调优心得Gurobi的MIPGap设为1e-4意味着严格求解但主问题中新增的二进制变量储能的充放电状态u_ess会让求解时间上升。如果模型收敛慢试试把MIPGap放大到1e-3或者把MIPFocus设置为1侧重于寻找可行解。另外Yalmip默认的模型转成LP格式时有可能会包含冗余变量可以在求解前调用clean或compact函数压缩模型。实测下来Gurobi对压缩后的模型能提速20%以上强烈推荐。5. 仿真案例分析与结果对比5.1 算例系统配置为了验证关键场景辨别算法的效果我搭了一个简化微网算例。包含微燃机1额定功率2 MW发电成本系数 58 $/MWh启动成本 15 $。微燃机2额定功率1.5 MW发电成本系数 62 $/MWh启动成本 12 $。储能额定容量 4 MWh最大充放电功率 1 MW充放电效率均取0.95SOC范围[0.2, 0.9]初始SOC为0.5。光伏额定装机 3 MW预测出力曲线来自典型晴天数据。风电额定装机 1 MW预测出力曲线较平滑偏差20%。与大电网联络线购电上限5 MW电价采用分时电价峰时段1.2 $/kWh谷时段0.5 $/kWh平均0.8 $/kWh。负荷曲线和光伏曲线我前面已经列了预测值。预算参数光伏Γ6负荷Γ4对应24小时中最多有6个时段光伏偏离极端4个时段负荷偏离极端。5.2 三种方法的对比我做了三组实验确定性模型即用预测值调度、标准CCG两阶段鲁棒、基于关键场景辨别算法的两阶段鲁棒KSI-CCG。结果如下方法总成本$迭代次数总计算时间s最恶劣场景下的成本$确定性2340-32760严重越限不可行标准CCG2680152202680KSI-CCG267591322675确定性方案看起来便宜但最恶劣场景下实际不可行需要切负荷或超限说明任何调度方案都必须考虑鲁棒性。两阶段鲁棒方案比确定性方案成本高约14%这“溢价”买的就是对不确定性的免疫。KSI-CCG比标准CCG节省约40%的耗时迭代次数也少了6次成本几乎相同。这个结果在我换了几组随机预测数据后都稳定说明算法加速有普适性。5.3 关键场景识别效果观察为了更直观我把每次迭代识别出的最恶劣场景偏差向量画了出来注意我这里不用图片用文字描述。标准CCG通常先加入的是光伏明显偏低且负荷偏高的场景之后会逐步加入若干“局部恶劣”场景比如只有某几个时段比较极端。而KSI-CCG因为在预筛阶段加入了覆盖不同聚类中心的场景主问题一开始就具有一定鲁棒性后续每次迭代识别场景更多是为了修正边缘情况所以迭代收敛更快。这里有个有意思的现象在鲁棒模型求解结果的调度方案中储能往往在光伏出力低、负荷高的时段来临之前提前充电像是“未雨绸缪”。这其实正是鲁棒优化的魅力它主动为最恶劣场景留出储能裕量而不是像确定性优化那样只盯着预测曲线。5.4 预算参数Γ的灵敏度分析预算参数是控制保守度的旋钮。我把光伏的Γ从0对应确定性逐步增加到12全偏差结果总成本从2340上升到2890左右增长率约23%。而最恶劣场景下的储能SOC最低值也随之上升说明预留了更多安全裕量。当Γ超过10后成本增长趋缓说明最恶劣场景的“边际伤害”已经饱和。实际工作里我通常建议根据微网对风险的容忍度选Γ比如希望系统能应对“20%概率下的最坏情况”就取较小值如果系统要适应极端天气则取较大值。6. 常见问题与调试心法6.1 主问题规模爆炸怎么办每次迭代加入一个场景对应一组约束迭代10次后主问题中约束数量和变量数量显著增加MILP求解时间飙升。我的办法是限制主问题求解时间如60秒如果超时就用当前可行解作为近似解更新LB并继续迭代。虽然理论上可能导致不收敛但实际中只要间隙已经比较小比如原始间隙小于0.5%即便只算到次优解调度方案也有参考价值。另一种办法是使用Gurobi的“目标优先”功能让求解器先找到较好的可行解再改进上界。或者对第一阶段中整数变量数量做合理压缩例如将机组启停时间以4小时为颗粒度分组减少整数变量个数。6.2 子问题最恶劣场景“跳来跳去”不收敛有时候子问题求出的最恶劣场景在迭代中会振荡导致UB不降反升。这是典型的“无界或者对偶偏差问题”。排查方向一是检查内层LP是否对所有场景都有可行解。如果没有意味着某些极端场景下第二阶段无法平衡功率就需要在第一阶段加入“备用容量”约束或者允许切负荷但切负荷有罚因子。二是检查对偶转换可行性尤其强对偶条件是否成立必要时用KKT方法替换。三是检查预算约束是否正确线性化了绝对值。绝对值|P-P_pred|≤Δ可以用两个线性约束表示P-P_pred≤Δ·zP_pred-P≤Δ·z其中z是二进制变量或连续松弛到[0,1]。如果没有引入z直接写|P-P_pred|/Δ≤ΓYalmip会报错或生成非线性约束整个子问题就跑偏了。6.3 Yalmip与求解器的License问题最让新手崩溃的不是模型难建而是求解器突然打不开。Gurobi在美国的网站上用校园网或学校License都没啥问题Cplex是IBM的学术版需要注册。如果你的许可证是免费的“Gurobi Named User License”记得在Matlab里用gurobi_setup一次否则找不到求解器。另外Matlab自身的Optimization Toolbox自带linprog和intlinprog也可以直接用来做LP子问题但性能比Gurobi慢不少。如果你只是做小算例验证算法完全可以先用intlinprog跑通再换Gurobi提速。6.4 数值稳定性与缩放微网的量纲很关键功率用MW成本用$有时候目标函数里还掺杂启动成本整个模型矩阵的元素可能跨了好几个数量级。如果不做缩放Gurobi会给出numerical warning甚至出现不可靠的对偶值。建议把所有功率转成标幺值或者统一到同一个数量级比如都用kW成本都用k$这样模型鲁棒性会好很多。我一次项目里因为电价用了$/kWh而功率用了MW导致目标函数系数差了1000倍结果收敛性极差浪费了我一下午后来统一成“kW与$/MWh”就一切正常。6.5 与随机规划、分布鲁棒优化的关系如果你正在做微网不确定性的课题可能还会接触到随机规划SP和分布鲁棒优化DRO。我的理解是SP假设知道概率分布求期望最优传统鲁棒只假设集合求最坏情况最优DRO则结合两者假设分布在一个模糊集内求解最坏分布下的期望最优。关键场景辨别算法同样可以扩展应用到DRO中——只不过子问题从“找最恶劣场景”变成“找最恶劣分布”难度更上一层。但如果你能熟练实现本文中的CCG与场景辨别再去看DRO的论文代码会轻松很多因为它们共享“主问题-子问题分解割平面迭代”这个骨架。写在最后的一点经验从建确定性模型到跑通两阶段鲁棒我用了一个多月中间大部分时间不是花在理论上而是花在调试“为什么不收敛”“为什么场景加入主问题后反而无解”这些工程细节上。我个人在实际操作中的体会是两阶段鲁棒优化不是一项“装好就能跑”的现成技术而是一个需要不断调参、不断验证的系统工程。关键场景辨别算法虽然在效率上带来收益但对模型细节的敏感度也很高你必须在不确定集合的形状、预算参数的取值、子问题的求解精度之间找到平衡。最后再分享一个小技巧在跑主问题之前可以先用确定性模型跑一遍拿到一个初始可行解作为主问题的x0Yalmip支持用assign指定初始值这能显著减少MILP的搜索时间。还有如果你要把这套代码扩展到更大规模的配电系统建议把场景集合的存储和约束构建写成函数而不是全部塞在循环里不然代码会乱到你自己都不想看。这次分享就到这希望能给正在调鲁棒微网模型的你一些启发。