做微电网调度的人大概率都经历过这种场景日前把机组启停和储能充放电计划排得明明白白结果第二天一阵云遮住光伏或者风突然掉下去整个计划就得推倒重来更麻烦的是如果备用留得不够还得在高峰时段切负荷这事放在实际园区里没人能接受。既然预测做不到百分百准确那不如换个思路不追求预测准而是让调度方案在预测偏差的“最坏情况”下依然能扛得住。这正是含可再生能源和储能的区域微电网最优运行转向鲁棒优化的核心动机。我最近完整复现了一套考虑鲁棒性和不确定性的区域微电网多阶段鲁棒调度模型用Matlab把数学建模、CCG列与约束生成算法、求解器配置一路跑通代码从零搭起没有依赖现成工具箱里的黑盒函数。整个过程比想象中坑多但也确实把鲁棒优化的来龙去脉摸透了。这篇博文把模型推导、算法逻辑、代码实现和踩过的坑全部写出来给正在做微电网调度、储能优化、鲁棒优化或者准备复现类似论文的朋友做个参考尤其是对“论文公式能看懂、动手写代码就卡壳”的同学应该会有帮助。1. 调度困境的根源风光不确定性如何扰动微电网运行1.1 确定性日前调度的局限在哪里传统微电网优化调度最常见的形式是确定性模型也就是把光伏、风电的预测出力当成固定参数在满足各种运行约束的前提下求解总成本最小化问题。这种模型实现简单求解也快但它有一个天然缺陷——预测值本身是错的。以我搭建的某园区微电网测试系统为例风电额定装机800 kW光伏额定装机600 kW。日前预测某时刻风电出力为500 kW实际运行中可能只有380 kW偏差接近24%。如果调度计划按500 kW安排负荷供应那120 kW的缺口就必须靠其他机组或储能顶上。可如果在同一时刻联络线功率已经顶到上限储能SOC又偏低切负荷就成了唯一选择。很多工程人员应对这个问题的方法是留旋转备用所有可调机组按最大出力减去当前出力来预留空间储能也强制设置SOC下限以上的一段“安全区”。这种方法逻辑上没错问题在于“备多少”完全靠经验。备用容量留大了机组的运行点被迫抬高燃料成本和购电成本同步上升留小了遇到大风大光或大波动依然出问题。尤其对于区域微电网这种可调资源有限、联络线功率受限的系统盲目留备用的代价非常可观。这也是为什么这两年“用显式建模方式处理不确定性”成为主流方向。不再把预测值当作确定参数而是把预测误差建模成一个不确定集合让调度方案在这个集合里的任意取值下都满足运行约束。这就是鲁棒优化的基本逻辑也是标题里“考虑鲁棒性和不确定性”这句话的实质含义。1.2 不确定性的数学描述从预测误差到盒式集合把不确定性写进模型的第一步是定义“误差长什么样”。使用最广泛的是盒式不确定集形式非常直白[ W_t^w \in [\bar{W}_t^w - \hat{W}_t^w,\ \bar{W}_t^w \hat{W}_t^w] ]其中 (\bar{W}_t^w) 是预测值(\hat{W}_t^w) 是允许的偏差半径可以取预测值的某个百分比也可以取历史预测误差的90%分位数。光伏同理。这个区间就是“最坏情况”的边界。但只给区间还不行。如果所有时段都同时取边界值这个场景物理上几乎不可能出现求解出来的调度方案会保守到离谱运行成本暴涨几十个百分点。所以标准做法是再加一条预算约束限制整个调度周期内总偏差的大小引入保守度参数Gamma。它的含义可以理解为“允许多少个时段出现接近最坏偏差”也可以写成累计偏差不超过某个总量[ \sum_t \frac{|W_t^w - \bar{W}_t^w|}{\hat{W}_t^w} \le \Gamma ]Gamma越大不确定集合越“胖”调度方案越保守Gamma等于0模型就退化成确定性模型。这里有一个容易忽略的细节偏差半径 (\hat{W}_t^w) 必须和预测时长的尺度匹配。15分钟级滚动调度和24小时日前调度的预测误差分布完全不同不能照搬同一个百分比。我复现时的典型配置是日前风电偏差取15%光伏偏差取20%Gamma取调度点数的25%到40%。这个尺度把握直接决定结果有没有工程参考价值。2. 多阶段鲁棒调度模型目标函数、分层决策与约束2.1 目标函数怎么建模型的目标是最小化调度周期内的总运行成本但要区分不同决策阶段。第一阶段是“预调度”在不确定性尚未揭晓前决定机组启停、储能基准充放电计划和联络线基准功率第二阶段是“调整调度”在最坏场景出现后允许各机组和储能做出相应的再调度响应。写成两阶段鲁棒优化的标准形式是[ \min_{x} \left( c^T x \max_{u \in U} \min_{y} d^T y \right) ]其中 (x) 是第一阶段决策变量(u) 是归结于不确定集 (U) 的风光出力真实值(y) 是第二阶段针对特定场景的调整变量。外层min处理预调度成本内层max-min处理“最坏场景下的最小调整成本”。实际建模时目标函数包含的项主要有四类微型燃气轮机的燃料成本按出力二次曲线近似或分段线性化从上级电网购电的成本按分时电价计算弃风弃光的惩罚成本鼓励尽量消纳清洁能源切负荷的惩罚成本这是模型中最重的惩罚项理论上优先级最高。为什么把切负荷惩罚设得非常高因为在工程语义里弃风弃光是可以接受的是“可控的浪费”而切负荷意味着用户侧断电是运行事故。所以目标函数中切负荷惩罚系数至少要比购电成本高一个数量级否则优化器会选择切负荷而不是高价买电结果不符合运行逻辑。2.2 决策变量的分层结构“多阶段”三个字在具体模型里常常被简化处理但决策变量的阶段归属必须清晰。第一阶段决策在不确定性揭晓前就要确定包括燃气轮机的启停状态0-1变量和基准出力储能系统的基准充放电状态与功率与上级电网联络线的基准交换功率。第二阶段决策是针对每个可能场景的“响应动作”包括燃气轮机的实际出力调整量爬坡和上下限约束下储能实际充放电功率修正切负荷量弃风弃光量。这种分层结构的意义在于日前计划要提前下发但具体怎么执行要等实际风光出力出来之后再说。CCG算法迭代过程中每识别出一个新的最坏场景就会给主问题增加一组与该场景对应的第二阶段约束这一步正好体现了“逐步逼近真实不确定场景”的多阶段决策特征。2.3 关键约束功率平衡、储能SOC与网络边界模型的约束体系按功能划分成四组。功率平衡约束是每一时刻都必须满足的物理守恒机组出力、光伏出力、风电出力、储能放电功率、购电功率加总等于负荷加上储能充电功率和售电功率。注意在第二阶段这个等式里所有的功率项都要带上“场景索引”因为不同场景下风光出力不同平衡关系自然不同。储能约束是模型最容易写错的地方。SOC递推方程、充放电功率上下限、同一时刻不能同时充放电等约束必须严格写全。SOC的初值一般取0.2或0.5容量调度周期末再约束回落到接近初值否则储能会在最后一个时段把电全放光结果没有工程意义。我用的典型参数是储能容量1000 kWh最大充放电功率200 kWSOC上下限0.1到0.9末端SOC约束0.5初始0.5。机组约束包括出力上下限、爬坡速率约束和最小启停时间约束。如果只做日前经济调度最小启停时间可以简化处理但如果要考虑机组的物理运行特性这一段必需完整建模否则调度计划在执行层面会遇到机组根本无法按计划启停的情况。联络线约束限制与上级电网交换功率不超过给定上限这一点在区域微电网里非常关键。联络线功率上限设得越小对系统本地调节能力的要求越高储能的角色就越重要。我在测试系统里把联络线功率上限设为300 kW相当于负荷峰值的不到三分之一这样能充分暴露鲁棒调度的价值。3. 列与约束生成算法为什么它能破解这类模型3.1 为什么随机优化不够鲁棒优化需要专门算法有人会问既然有随机优化为什么要用鲁棒甚至有人直接把场景削减后的随机规划当成鲁棒模型来用。二者的本质区别在于随机优化要求的是“期望成本最低”它对极端场景的容忍度较高鲁棒优化要求的是“最坏情况下成本最低”它对极端场景更敏感。如果你的微电网服务的是医院、数据中心这类对供电可靠性要求极高的负荷显然更希望用鲁棒尺度来评估方案。模型建好之后直接求解是不现实的。min-max-min三层结构本身就是NP难问题不可能扔给商用MIP求解器一把梭。工程上常用的破解手段有两类Benders对偶分解和列与约束生成Column-and-Constraint GenerationCCG。Benders方法通过对偶把子问题的极端场景信息反馈给主问题迭代较慢CCG方法则是把极端场景对应的完整变量和约束直接加入主问题主问题规模增长更快但收敛性远远好于Benders尤其是第二阶段有大量连续变量时CCG的优势非常明显。3.2 主问题与子问题的迭代机制CCG算法的骨架并不复杂。主问题MP是最小化第一阶段成本加上一个辅助变量theta并满足已经识别出来的所有最坏场景对应的第二阶段约束。每轮迭代给MP增加一个场景实际上是给MP“塞”了一组新的连续变量和约束[ \min_{x, y, \theta} c^T x \theta ]s.t. [ \theta \ge d^T y^{(k)}, \quad A x B y^{(k)} \le C u^{(k)} ]其中 (u^{(k)}) 是第k轮从子问题中解出的最坏场景。子问题SP在给定第一阶段解 (x^*) 的情况下寻找使调整成本最大的风光出力场景[ \max_{u \in U} \min_{y} d^T y \quad \text{s.t.} \quad A x^* B y \le C u ]整个迭代流程可以写成四步求解MP得到当前最优的第一阶段决策和LB固定第一阶段决策求解SP得到最坏场景和最优调整成本若SP求出的成本加上第一阶段成本构成新的UB比较UB和LB的间距收敛则停止否则把SP找到的场景加入MP重新求解。为什么这一步叫“多阶段”因为每一轮CCG迭代都会识别一个新的极端场景并且把该场景下的完整调度决策纳入主问题。随着迭代进行主问题实际上是在对一个不断扩大的“场景决策树”做联合优化这正是多阶段决策思想的数值实现。3.3 SP内部的双层问题强对偶的处理逻辑SP是一个max-min双层问题不能直接求解。标准处理方式是把内层的min问题看作一个线性规划利用强对偶定理将其转化为max问题于是外层max和内层max合并SP变成一个单层max问题。这里需要特别留意内层LP的对偶变量维度必须与原问题约束一一对应。对偶后的目标函数会出现“第一阶段决策变量乘对偶变量”的双线性项。虽然它们都是连续变量但乘积项在标准MIP框架下不可直接处理。常见的做法是用大M法做线性化[ x \cdot \lambda z,\quad 0 \le z \le M \cdot x,\quad \lambda - M(1-x) \le z \le \lambda ]这一步是整个算法实现中最容易出错的地方也是后面提到的坑。如果对偶方向写错或大M取值不当SP求解出的场景根本代表不了“最坏情况”CCG迭代自然会收敛到错误的结果。4. Matlab代码复现从数学公式到可运行程序的完整过程4.1 代码文件结构与数据准备我在复现时把代码拆成几个文件不搞一个几千行的main脚本这样调试起来效率高很多。核心文件分工如下main_ccg.m主入口负责初始化参数、循环调用MP和SPcase_data.m存放所有系统参数包括机组参数、负荷曲线、风光预测出力与偏差半径build_uncertainty.m生成不确定集矩阵包括Box区间和预算约束mp_model.m构建并求解主问题MPsp_model.m构建并求解子问题SPplot_result.m结果汇总与绘图。数据准备的关键是基准场景。我生成了24个时段的负荷预测曲线、风电预测曲线和光伏预测曲线。风电曲线设定为夜间稍有出力、白天相对平稳的模式光伏曲线设定为中午高、早晚低的“馒头形”。这样的典型场景能较好展示鲁棒模型在风光波动最剧烈时段的应对能力。4.2 主问题MP的核心实现主问题用YALMIP建模求解器调用Gurobi。核心代码的关键片段如下% 主问题变量 x binvar(ng, T); % 机组启停 p_g sdpvar(ng, T); % 机组基准出力 p_ess_c sdpvar(1, T); % 储能基准充电功率 p_ess_d sdpvar(1, T); % 储能基准放电功率 theta sdpvar(1, 1); % 辅助变量 soc sdpvar(1, T1); % 储能SOC % 目标函数第一阶段成本 theta objective sum(sum(fuel_cost .* p_g)) ... sum(grid_price .* p_buy) theta; % 基础约束 constraints []; constraints [constraints, p_g_min .* x p_g p_g_max .* x]; constraints [constraints, soc(:, t1) soc(:, t) eta_c * p_ess_c(:, t) - p_ess_d(:, t) / eta_d]; % 把识别出的最坏场景约束加入主问题 for k 1:num_scenarios u_k worst_scenarios{k}; % 为第k个场景定义调整变量 y_k % 加入约束: theta cost_k, Ax By_k C u_k end注意theta cost_k这一组约束是CCG迭代的核心每轮识别出一个新场景就往主问题里加一组带场景索引的变量和约束。随着k增加主问题的规模会变大这是正常的。4.3 子问题SP的核心实现与对偶转换SP的建模比MP复杂得多。我把SP写成对偶后的单层max问题这样Gurobi可以直接求MIP。简化版代码如下% 内层LP的原始变量 y sdpvar(size_y, T); % 第二阶段调整变量 % 对偶变量 lambda sdpvar(size_lambda, T); % 对应等式约束 mu sdpvar(size_mu, T); % 对应不等式约束 u_w sdpvar(1, T); % 风电实际出力 u_pv sdpvar(1, T); % 光伏实际出力 % SP目标max 对偶目标函数包含双线性项 x* * lambda sp_obj sum(b * lambda) sum(c * mu); % 用大M线性化双线性项 for t 1:T % z x* * lambda 的线性化 constraints [constraints, z(:, t) M * x_fixed(:, t)]; constraints [constraints, z(:, t) lambda(:, t)]; constraints [constraints, z(:, t) lambda(:, t) - M * (1 - x_fixed(:, t))]; constraints [constraints, z(:, t) 0]; end % 不确定集约束 constraints [constraints, u_w wind_pred - wind_dev]; constraints [constraints, u_w wind_pred wind_dev]; constraints [constraints, sum(abs(u_w - wind_pred) ./ wind_dev) Gamma];这里有两处需要注意。第一x_fixed是MP求解得到的第一阶段决策值传给SP后要当作常数不能再作为变量第二大M的取值不能随便设太大会导致求解器数值病态我测试下来取机组最大出力的10倍比较稳妥。4.4 求解器配置与运行耗时模型用YALMIP作为建模层求解器采用GurobiMIPGap设为0.01%时间上限设为300秒。测试系统规模是24时段、3台燃气轮机、1台储能、1条联络线MP加SP一轮迭代的耗时大约10到20秒整体CCG迭代8到10轮即可收敛总耗时在3分钟以内。这个规模在普通笔记本上完全可以接受。如果模型规模更大比如时段细化到96点迭代轮数和每轮耗时都会显著上升。这时候建议先关闭一些不关键的场景约束或者放宽MIPGap优先保证CCG迭代能跑完再逐步收紧。5. 复现过程中的三个典型坑每一处都值得记录5.1 对偶变量方向混乱导致的SP无界第一次写SP时我把内层min问题转max的时候对偶变量方向写反了直接后果是SP的目标函数可能在可行域的某个方向上趋向无穷Gurobi报出“模型无界”的错误。排查后发现问题出在约束含义与对偶变量正负号的匹配上内层min问题里“≤”型约束对应的对偶变量应限制为非负而等式约束对应的对偶变量是自由变量。如果搞反了强对偶后的上界性质彻底失效SP求出的“最坏场景”数值完全失真。排查这类问题有一个好习惯先固定一组已知的第一阶段解手算一个小规模SP比如只用2个时段、1台机组把对偶后的结果和暴力枚举精确解对比符号错乱和方向反了立刻就能暴露。5.2 大M线性化的数值病态双线性项线性化时我最初把大M设为10000结果Gurobi一直在报数值警告MIP求解时间从十几秒暴涨到几分钟而且迭代出来的场景约束精度很差。后来把大M调整到与对应变量量级匹配的范围问题立刻缓解。经验是大M只要能保证约束覆盖变量可行域即可越小越好。对这台测试系统机组出力、联络线功率、切负荷量都在几百到一千千瓦的量级大M取2000就能覆盖所有情况。数值稳定之后每轮耗时降到原来的三分之一。5.3 收敛判据的“假收敛”与场景重复CCG迭代中如果UB和LB的间距除以UB小于某个阈值就停止看起来很快收敛但这可能是“假收敛”。我遇到过一种情况SP找到的最坏场景和前一轮几乎相同加入MP后对解没有实际改善UB和LB虽然贴得近但其实迭代陷入了一个局部场景的死胡同。原因出在SP求解时使用了松弛处理导致SP返回的并不是严格的最坏场景。解决方法是限定SP求解的严格性MIPGap不放宽同时在每次迭代后检查新场景与历史场景的差异度。如果连续两轮场景完全相同就说明SP需要重新审视约束是否写全而不是直接宣布收敛。6. 调度结果解读鲁棒性背后到底付出了什么代价6.1 典型日调度结果与最坏场景识别CCG迭代收敛后主问题得到一组“能应对所有已识别最坏场景”的调度方案SP则给出了系统面临的最严苛风光场景。在我的测试案例中最坏场景并不是“全天都没风没光”而是出现在晚间负荷高峰时段——光伏已经归零风电同时跌到预测下限而此时负荷正处在一天中的最高点。这个场景虽然概率不高但后果很严重鲁棒模型正是为此预留了裕度。对应的调度结果中储能表现很有意思确定性模型倾向于在光伏大发的中午时段把SOC充满然后在傍晚放电到接近下限鲁棒模型则在中午保留了一部分充电空间傍晚放电也更加温和最终的SOC曲线整体抬高相当于提前为“晚间最坏场景”预留了额外的放电能力。这个行为特征在图像上非常直观也是鲁棒调度“未雨绸缪”属性的最好体现。6.2 保守度参数Gamma的成本敏感性我把Gamma从0逐步调到16观察总成本的变化Gamma值总运行成本元最大切负荷量kW储能最大放电深度0确定性1286022082%4134509570%8140802263%1614890058%可以看到从Gamma0到Gamma8成本只上升了9.5%最大切负荷量却从220 kW降到了22 kW换来了可靠性的大幅提升。而从Gamma8到Gamma16切负荷已经被压到0成本却还要再增加近6%。这说明Gamma的取值在中间区间收益最高过了某个点之后就是纯粹的“为不存在的最坏情况买单”。实际工程中Gamma不必取太大结合历史误差分布取中间值往往性价比最高。6.3 鲁棒调度与确定性调度的对照实验对照实验是复现时最值得做的一件事它能把鲁棒调度的价值量化。我的做法是用同一套系统参数分别跑确定性模型和鲁棒模型然后把两组调度方案放到同一批随机抽样场景里做后验仿真。仿真结果显示确定性方案在约35%的抽样场景中会出现不同程度的切负荷鲁棒方案只有不到5%的场景需要少量切负荷且切负荷量非常小。代价是鲁棒方案的日前发电计划更“保守”购电成本平均高出约8%。但这笔钱买来的是用户侧供电可靠性和运行风险的大幅下降。对供电可靠性要求高的场景这个溢价完全值得如果负荷允许一定程度的中断那么确定性方案配合更灵活的日内滚动调度可能也是合理选择。最后分享一个调试经验把不确定性预算参数从0开始逐步调大分别跑一遍完整模型然后把总成本和储能SOC曲线放在同一张图上对比。Gamma增大后成本曲线应当平滑上升储能曲线应当越来越“谨慎”。如果哪一段出现成本骤降或SOC曲线异常跳变那大概率不是模型变聪明了而是某个约束在特定场景下被写松了。这个方法帮我排查过至少两个隐蔽的建模错误远比翻公式调试效率高建议你也试试。