做微电网调度的人,大概率都经历过这种场景:日前把机组启停和储能充放电计划排得明明白白,结果第二天一阵云遮住光伏,或者风突然掉下去,整个计划就得推倒重来;更麻烦的是,如果备用留得不够,还得在高峰时段切负荷,这事放在实际园区里没人能接受。既然预测做不到百分百准确,那不如换个思路:不追求预测准,而是让调度方案在预测偏差的“最坏情况”下依然能扛得住。这正是含可再生能源和储能的区域微电网最优运行转向鲁棒优化的核心动机。
我最近完整复现了一套考虑鲁棒性和不确定性的区域微电网多阶段鲁棒调度模型,用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 kW,SOC上下限0.1到0.9,末端SOC约束0.5,初始0.5。
机组约束包括出力上下限、爬坡速率约束和最小启停时间约束。如果只做日前经济调度,最小启停时间可以简化处理;但如果要考虑机组的物理运行特性,这一段必需完整建模,否则调度计划在执行层面会遇到机组根本无法按计划启停的情况。
联络线约束限制与上级电网交换功率不超过给定上限,这一点在区域微电网里非常关键。联络线功率上限设得越小,对系统本地调节能力的要求越高,储能的角色就越重要。我在测试系统里把联络线功率上限设为300 kW,相当于负荷峰值的不到三分之一,这样能充分暴露鲁棒调度的价值。
3. 列与约束生成算法:为什么它能破解这类模型
3.1 为什么随机优化不够,鲁棒优化需要专门算法
有人会问:既然有随机优化,为什么要用鲁棒?甚至有人直接把场景削减后的随机规划当成鲁棒模型来用。二者的本质区别在于:随机优化要求的是“期望成本最低”,它对极端场景的容忍度较高;鲁棒优化要求的是“最坏情况下成本最低”,它对极端场景更敏感。如果你的微电网服务的是医院、数据中心这类对供电可靠性要求极高的负荷,显然更希望用鲁棒尺度来评估方案。
模型建好之后,直接求解是不现实的。min-max-min三层结构本身就是NP难问题,不可能扔给商用MIP求解器一把梭。工程上常用的破解手段有两类:Benders对偶分解和列与约束生成(Column-and-Constraint Generation,CCG)。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和SP;
- case_data.m:存放所有系统参数,包括机组参数、负荷曲线、风光预测出力与偏差半径;
- build_uncertainty.m:生成不确定集矩阵,包括Box区间和预算约束;
- mp_model.m:构建并求解主问题MP;
- sp_model.m:构建并求解子问题SP;
- plot_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, T+1); % 储能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(:, t+1) == 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作为建模层,求解器采用Gurobi,MIPGap设为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(确定性) | 12860 | 220 | 82% |
| 4 | 13450 | 95 | 70% |
| 8 | 14080 | 22 | 63% |
| 16 | 14890 | 0 | 58% |
可以看到,从Gamma=0到Gamma=8,成本只上升了9.5%,最大切负荷量却从220 kW降到了22 kW,换来了可靠性的大幅提升。而从Gamma=8到Gamma=16,切负荷已经被压到0,成本却还要再增加近6%。这说明Gamma的取值在中间区间收益最高,过了某个点之后就是纯粹的“为不存在的最坏情况买单”。实际工程中,Gamma不必取太大,结合历史误差分布取中间值往往性价比最高。
6.3 鲁棒调度与确定性调度的对照实验
对照实验是复现时最值得做的一件事,它能把鲁棒调度的价值量化。我的做法是:用同一套系统参数,分别跑确定性模型和鲁棒模型,然后把两组调度方案放到同一批随机抽样场景里做后验仿真。仿真结果显示,确定性方案在约35%的抽样场景中会出现不同程度的切负荷,鲁棒方案只有不到5%的场景需要少量切负荷,且切负荷量非常小。
代价是鲁棒方案的日前发电计划更“保守”,购电成本平均高出约8%。但这笔钱买来的是用户侧供电可靠性和运行风险的大幅下降。对供电可靠性要求高的场景,这个溢价完全值得;如果负荷允许一定程度的中断,那么确定性方案配合更灵活的日内滚动调度,可能也是合理选择。
最后分享一个调试经验:把不确定性预算参数从0开始逐步调大,分别跑一遍完整模型,然后把总成本和储能SOC曲线放在同一张图上对比。Gamma增大后成本曲线应当平滑上升,储能曲线应当越来越“谨慎”。如果哪一段出现成本骤降或SOC曲线异常跳变,那大概率不是模型变聪明了,而是某个约束在特定场景下被写松了。这个方法帮我排查过至少两个隐蔽的建模错误,远比翻公式调试效率高,建议你也试试。