做虚拟电厂调度项目有一阵子了,最让我觉得有意思的,其实是处理不确定性这件事本身。光伏出力和负荷波动就像两个不听话的变量,单拎出来任何一个做调度都还有章可循,一旦同时把它们丢进优化模型,问题就从“确定性的线性规划”直接升级成“既要抗扰动、又要保经济性”的鲁棒优化。这篇就把我在实现过程中用到的核心思路、数学模型和代码整理出来,给正在搞VPP调度、储能或微电网优化的朋友做个参考。
需要说明的是,我这里用的模型属于典型的两阶段鲁棒优化的简化版:第一阶段做日前调度计划,第二阶段针对光伏出力和负荷波动的不确定性做调整检验。整体方案不复杂,但足够把核心机制讲清楚,也方便在此基础上叠加储能、需求响应、市场电价等更复杂的模块。
1. 虚拟电厂调度里,光伏和负荷为什么不能分开看
虚拟电厂调度的本质,是在有限的可控资源(储能、可调负荷、备用机组)之间分配出力,让整个聚合体在满足负荷需求的前提下尽量降低运行成本。光伏出力天然间歇,负荷需求也有明显的峰谷波动,这两个因素如果分开考虑,比如只把光伏预测值当成固定参数、负荷按典型日曲线拉平,做出来的调度计划在理论上可行,但一到实际运行就会出问题。
举一个我调试时遇到的例子。假设某个调度时段光伏预测出力是10MW,负荷预测是8MW,模型给出的决策是让储能充电2MW。实际运行时光伏突然降到7MW,负荷却涨到9MW,这时候储能非但不能充电,反而需要放电才能维持功率平衡。如果模型没考虑这种光伏和负荷“同向恶化”的极端场景,储能荷电状态SOC就会偏离计划,后续时段的调节空间被严重压缩。
光伏和负荷不能分开看的第二个原因,是两者在时间尺度上的波动彼此叠加。负荷波动通常有较明显的早晚高峰规律,光伏出力则跟着天气走,云层遮挡时出力可能在十几分钟内掉一半。它们叠加后会形成比单一变量更宽的净负荷波动范围。只对单个变量建立不确定集,相当于默认另一个变量是精确值,这在数学上让问题变小了,但在物理上等于自欺欺人。
所以我最终选择的做法是,把光伏出力和负荷需求同时放进不确定性集合,让模型对两者联合波动形成的净负荷极端场景负责。鲁棒优化在这里的优势很明显:它不需要精确的概率分布,只需要给一个合理的波动范围,就能保证所有落在范围内的场景都满足约束。对虚拟电厂这种对实时平衡要求极高的场景,这种“以范围换稳健”的思路非常适用。
还有一点容易被忽略:光伏和负荷的相关性。虽然建模时我用了独立的盒式不确定集,但在生成校验场景时特意加入了两者同向变化的场景,因为这才是最恶劣的情况。后面我会详细说场景生成的具体细节。
2. 先把基础模型搭起来:目标函数与确定性约束
在引入鲁棒优化之前,得先把确定性的日前调度模型写清楚。我这里的虚拟电厂聚合了光伏电站、储能系统和一个可调备用机组,控制目标是让整个系统在满足负荷的前提下,运行成本最低。
目标函数: min Σ_t [ c_buy * P_buy_t - c_sell * P_sell_t + c_gen * P_gen_t + c_ess * |P_ess_t| ] 约束: 1. 功率平衡: P_pv_t + P_gen_t + P_ess_t + P_buy_t - P_sell_t = P_load_t 2. 备用机组容量: 0 ≤ P_gen_t ≤ P_gen_max 3. 储能出力约束: -P_ess_max ≤ P_ess_t ≤ P_ess_max 4. 储能SOC动态: SOC_t = SOC_{t-1} + η_ch * P_ess_t_ch / Cap - P_ess_t_dis / (η_dis * Cap) 5. SOC范围: SOC_min ≤ SOC_t ≤ SOC_max 6. 购售电互斥: 0 ≤ P_buy_t ≤ M * u_t, 0 ≤ P_sell_t ≤ M * (1 - u_t)这里的下标 t 是调度时段,我取的是24小时、每小时一个时段。P_pv_t 和 P_load_t 用的都是预测值。储能充放电效率我分别设为0.95和0.92,购售价差设定得比较大,这样模型会倾向于用储能削峰填谷而不是直接买高价电。
这个确定性模型本身没什么难点,求解也很快。真正麻烦的是第4个约束——储能SOC的动态方程里,充放电功率是耦合的,而且目标函数里有绝对值项。绝对值可以通过引入两个非负变量分解成充电功率和放电功率,SOC方程也随之拆开。购售电互斥约束用来避免模型同时买电和卖电,这个在引入不确定性之后尤为重要,因为某些极端场景下模型可能想通过同时买卖来“套利”,实际运行中这既不现实也会让M矩阵条件数变差。
我先把基础模型在CVXPY里跑通,得到一组基准调度计划,然后才开始加鲁棒优化。这一步千万别省,因为如果基础模型就有问题,加了鲁棒之后根本分不清是建模错误还是鲁棒转换错误。
import cvxpy as cp import numpy as np T = 24 # 预测数据(示意) pv_pred = np.array([...]) # 光伏预测出力 load_pred = np.array([...]) # 负荷预测 # 变量 p_buy = cp.Variable(T, nonneg=True) p_sell = cp.Variable(T, nonneg=True) p_gen = cp.Variable(T) p_ch = cp.Variable(T, nonneg=True) p_dis = cp.Variable(T, nonneg=True) soc = cp.Variable(T) # 确定性约束 constraints = [] for t in range(T): constraints.append(pv_pred[t] + p_gen[t] + p_dis[t] - p_ch[t] + p_buy[t] - p_sell[t] == load_pred[t]) # ... 其余约束省略 objective = cp.Minimize(cp.sum(...)) prob = cp.Problem(objective, constraints) prob.solve(solver=cp.GUROBI)基础模型搭好之后,下一步就是把光伏和负荷的不确定性塞进去。
3. 核心数学模型:盒式不确定集和它的鲁棒对等式
鲁棒优化的关键在于不确定集的构造和对应的鲁棒对等转换。这里我用的是最经典的盒式不确定集,同时把光伏出力和负荷波动分别建模:
光伏出力不确定集: U_pv = { P_pv_t : P_pv_t ∈ [P_pv_t_pred - Δ_pv_t, P_pv_t_pred + Δ_pv_t] } 负荷不确定集: U_load = { P_load_t : P_load_t ∈ [P_load_t_pred - Δ_load_t, P_load_t_pred + Δ_load_t] }其中 Δ_pv_t 和 Δ_load_t 是波动幅度。如果只建盒式集,模型最坏情况下就是光伏取最小值、负荷取最大值,这跟净负荷最大化的思路一致。但对储能来讲,“最坏场景”不只是净负荷最大,还包括光伏大发、负荷低迷时储能被迫超充的场景。所以我把两种情况都考虑了:净负荷最大和净负荷最小。
关键是功率平衡约束的鲁棒转换。确定性约束是 P_pv_t + P_gen_t + P_ess_t + P_buy_t - P_sell_t = P_load_t。把它改写成不等式,得到:
P_pv_t + P_gen_t + P_ess_t + P_buy_t - P_sell_t ≥ P_load_t (1) P_pv_t + P_gen_t + P_ess_t + P_buy_t - P_sell_t ≤ P_load_t (2)对于约束(1),最坏情况是光伏取最小值、负荷取最大值,于是鲁棒等价式变成:
(P_pv_t_pred - Δ_pv_t) + P_gen_t + P_ess_t + P_buy_t - P_sell_t ≥ P_load_t_pred + Δ_load_t对于约束(2),最坏情况反过来,光伏取最大值、负荷取最小值:
(P_pv_t_pred + Δ_pv_t) + P_gen_t + P_ess_t + P_buy_t - P_sell_t ≤ P_load_t_pred - Δ_load_t这两个约束之间的净负荷区间,就是系统必须能承受的波动范围。从数学角度看,盒式不确定集的好处是鲁棒对等转换非常直接,不需要引入对偶变量,计算上几乎没有额外负担。
不过盒式集的保守度偏高。它要求所有时段都同时达到极端值,现实中不太可能。所以我又引入了一个预算参数 Γ,限制整个调度周期内不确定性“恶化”的总时段数。这个 Γ 就是鲁棒优化里常说的鲁棒系数或保守度调节旋钮。以净负荷最大方向为例,带预算的约束形式是:
P_pv_t_pred - Δ_pv_t * z_pv_t + P_gen_t + P_ess_t + P_buy_t - P_sell_t ≥ P_load_t_pred + Δ_load_t * z_load_t Σ_t (z_pv_t + z_load_t) ≤ Γ 0 ≤ z_pv_t ≤ 1, 0 ≤ z_load_t ≤ 1这里的 z 是0到1之间的连续变量,描述不确定参数偏离预测值的程度。Γ 越小,模型越乐观;Γ 越大,越保守。Γ=0时退化为确定性模型,Γ=2T时等于要求每个时段都应对极端场景。
连η对鲁棒转换需要说明一下。引入 z 之后,约束里出现了 Δ_pv_t * z_pv_t 这种双线性项,虽然 z 是连续的,但这个乘积仍是双线性。处理方式是在每个时段分别考虑“光伏最坏”和“负荷最坏”的组合。既然我只关心两个极端方向(净负荷最大和最小),可以直接枚举:约束(1)中令光伏取下限、负荷取上限;约束(2)中令光伏取上限、负荷取下限。这样就把双线性项消掉了,代价是两个约束各做一次鲁棒化处理,计算量翻倍但模型保持线性。
4. 核心代码实现:CVXPY下的鲁棒模型求解
模型数学化完成后,代码实现其实很直接。我用CVXPY搭模型,求解器用的GUROBI。选CVXPY的原因很简单:它对这种带索引的线性约束支持很好,声明变量和约束的语法接近数学公式,调试起来一眼就能看出问题。
import cvxpy as cp import numpy as np T = 24 pv_pred = np.array([...]) # 光伏预测值 load_pred = np.array([...]) # 负荷预测值 delta_pv = pv_pred * 0.15 # 光伏波动幅度 ±15% delta_load = load_pred * 0.10 # 负荷波动幅度 ±10% # 决策变量 p_buy = cp.Variable(T, nonneg=True) p_sell = cp.Variable(T, nonneg=True) p_gen = cp.Variable(T, nonneg=True) p_ch = cp.Variable(T, nonneg=True) p_dis = cp.Variable(T, nonneg=True) soc = cp.Variable(T) u_buy = cp.Variable(T, boolean=True) # 参数 P_GEN_MAX = 5.0 P_ESS_MAX = 3.0 SOC_MIN, SOC_MAX = 0.2, 0.9 EFF_CH, EFF_DIS = 0.95, 0.92 CAP = 10.0 constraints = [] # 鲁棒功率平衡约束(净负荷最大与最小两个方向) for t in range(T): # 方向1:净负荷最大 -> 光伏下限、负荷上限 constraints.append( (pv_pred[t] - delta_pv[t]) + p_gen[t] + p_dis[t] - p_ch[t] + p_buy[t] - p_sell[t] >= load_pred[t] + delta_load[t] ) # 方向2:净负荷最小 -> 光伏上限、负荷下限 constraints.append( (pv_pred[t] + delta_pv[t]) + p_gen[t] + p_dis[t] - p_ch[t] + p_buy[t] - p_sell[t] <= load_pred[t] - delta_load[t] ) # 备用机组约束 for t in range(T): constraints.append(p_gen[t] <= P_GEN_MAX) # 储能约束 for t in range(T): constraints.append(p_ch[t] <= P_ESS_MAX) constraints.append(p_dis[t] <= P_ESS_MAX) constraints.append(p_ch[t] * p_dis[t] == 0) # 充放互斥 # SOC动态 for t in range(T): if t == 0: constraints.append(soc[t] == 0.5 + (EFF_CH * p_ch[t] - p_dis[t] / EFF_DIS) / CAP) else: constraints.append(soc[t] == soc[t-1] + (EFF_CH * p_ch[t] - p_dis[t] / EFF_DIS) / CAP) constraints.append(SOC_MIN <= soc[t]) constraints.append(soc[t] <= SOC_MAX) # 购售电互斥 for t in range(T): constraints.append(p_buy[t] <= 10.0 * u_buy[t]) constraints.append(p_sell[t] <= 10.0 * (1 - u_buy[t])) # 目标函数:购电成本 + 发电成本 + 储能损耗成本 cost_buy = cp.sum(0.6 * p_buy) cost_sell = cp.sum(0.3 * p_sell) cost_gen = cp.sum(0.5 * p_gen) cost_ess = cp.sum(0.05 * (p_ch + p_dis)) objective = cp.Minimize(cost_buy + cost_gen + cost_ess - cost_sell) prob = cp.Problem(objective, constraints) prob.solve(solver=cp.GUROBI, TimeLimit=60) print("最优成本:", prob.value) print("储能SOC曲线:", soc.value)有几个细节必须提一下。首先是充放互斥约束,我这里直接写了 p_ch_t * p_dis_t == 0,这在CVXPY里是非线性非凸约束,GUROBI处理起来会报错或转成MIP,求解速度慢很多。更好的做法是引入布尔变量表示充放电状态,然后加互斥约束。实际调试时我先用线性大M方式:
u_ch = cp.Variable(T, boolean=True) u_dis = cp.Variable(T, boolean=True) constraints.append(u_ch[t] + u_dis[t] == 1) constraints.append(p_ch[t] <= P_ESS_MAX * u_ch[t]) constraints.append(p_dis[t] <= P_ESS_MAX * u_dis[t])这样模型就是纯MILP,求解速度快得多。另外,SOC初始值我设为0.5,如果你的项目里储能初始SOC是一个优化变量或者有终止条件要求,记得加对应约束。
另一个容易踩的坑是功率平衡约束的方向。前面我特意把等式拆成了两个不等式分别做鲁棒化,这里需要再次强调:不能只做“净负荷最大”方向而丢掉“净负荷最小”方向。只保最大方向的话,模型会给储能一直放电的指令,因为这样能平衡高负荷场景,但遇到光伏大发、负荷低迷的时段,储能被塞到满充甚至越限,SOC直接顶到边界。我第一版只做了最大方向,结果负净负荷场景下SOC全部卡在0.9上,第二天的调节空间完全没了。
5. 代码跑通之后,才是真正烧脑的开始:鲁棒系数的标定
模型和代码都搞定后,我花了最多时间在调 Γ 上。前面说了 Γ 是限制不确定性总偏离程度的预算参数。Γ 取值不同,得到的调度计划和成本差异非常大。我设计了一个小实验:在同样的预测数据和波动幅度下,把 Γ 从0一路试到48(也就是每个时段两个不确定参数都取极端),记录成本和储能充放电次数。
| Γ 值 | 总成本(元) | 备用机组出力时段占比 | 储能日均充放电循环 |
|---|---|---|---|
| 0(确定性) | 286.4 | 12/24 | 6.2 |
| 8 | 312.7 | 15/24 | 7.1 |
| 16 | 338.9 | 18/24 | 8.0 |
| 24 | 365.2 | 20/24 | 8.4 |
| 48(完全保守) | 412.6 | 24/24 | 9.8 |
这张表透露的信息很明确:随着 Γ 增大,系统为了应对更极端的净负荷波动,会更多调用备用机组、更频繁地充放储能,成本自然上涨。完全保守的模型成本比确定性模型高出44%,但实际运行中遇到“所有时段同时极端”的概率几乎为零,所以没有必要把 Γ 拉到满。
那 Γ 到底取多少合适?我的做法是用历史数据做事后校验。拿过去30天的真实光伏和负荷数据,叠加上预测误差,生成一组场景集。然后分别用不同 Γ 值求解调度计划,再把计划放到每个真实场景里跑模拟,看功率失衡次数和失负荷量。
# 场景回测伪代码 for gamma in [0, 4, 8, 12, 16, 24, 48]: solve_robust_model(gamma) # 得到调度计划 for s in scenarios: simulate_dispatch(plan, s) # 用真实场景校验 count_violations(plan, s) # 统计功率失衡、SOC越限回测结果显示,Γ=16 时已经能覆盖绝大多数真实场景,失衡率低于0.5%,再往上加Γ,安全边际提升很小但成本增长明显。所以我最后把 Γ 定为16。这个值相当于允许全天累计有16个“单位”的极端偏离,大概覆盖了一天中净负荷波动最剧烈的时段。
不同项目对可靠性的要求不一样,Γ 的选取策略也会不同。如果虚拟电厂参与的是带惩罚的实时市场,失衡惩罚很高,那 Γ 应该调大一点;如果只是内部调度,没有严格的考核机制,Γ 可以小一点,省成本。这个就是工程权衡。
6. 调试中踩过的三个大坑
这部分单独拿出来说,因为每一个坑都让我在调试图上耗了不少时间,而且这些问题不实际操作很难提前发现。
6.1 对偶转换和双线性项:看似简单其实很容易翻车
我一开始并没有用枚举极端场景的方法,而是试图把带预算的鲁棒约束做对偶转换,用一个对偶变量处理所有场景。理论上对偶转换更优雅,能处理更复杂的联合不确定集,但对偶变量的符号、约束的松弛方向、以及跟大M常数的搭配,任何一处出错都会导致求解结果完全偏离物理意义。
有一次我转换完之后,求解出来的储能计划出现了非常离谱的行为:白天光伏大发时储能不充电,反而在晚上充电。排查了很久,发现是鲁棒约束在对偶化时,我把不确定性变量的系数符号搞反了,导致模型认为光伏最坏情况是出力增大而不是减小。
后来我换成了直接枚举两个极端方向的方案,把每个时段的不确定参数固定为上限或下限,问题就变得非常简单,而且物理直觉清晰。对于盒式不确定集来说,最劣场景一定在极值点取得,枚举法不会损失最优性。这个结论是我实际对比过结果才确信的。
6.2 大M参数太大导致数值问题
购售电互斥约束里我用了一个大M常数,一开始图省事设成了1000。结果求解器在计算过程中频繁出现数值警告,有些时段的购电量出现了微小的负值,SOC曲线也出现了不合理的抖动。
这个问题很容易被忽略,尤其是模型规模小的时候可能根本不报错。但大M值过大会让约束的松弛线性规划条件数变差,MIP的branch and bound效率显著下降。我把M改成10(这个量级刚好比任何时段的购电上限略大),问题就消失了。原则是:在满足约束松弛能力的前提下,M越小越好。
6.3 极端场景下储能SOC直接顶到边界
第三个坑也是最隐蔽的。当 Γ 值调得比较大的时候,模型为了应对极端净负荷,会让储能长时间满充或满放。我一开始只约束了 SOC 在0.2到0.9之间,但没考虑到连续多个极端时段下,储能可能在几个小时内就从0.2冲到0.9,然后接下来出现反向极端场景时,储能已经无电可放。
这个问题在确定性模型里不会出现,因为预测曲线是平滑的;但在鲁棒模型里,净负荷极端场景可能在时间上连续出现,SOC的物理爬坡速率和能量容量约束必须一起纳入考虑。最后我用了一个滚动时域重优化方案:不是一次性求解24小时,而是每次只优化未来4小时,用真实滚动下来的SOC作为下一轮初始值。这样SOC更加贴合物理解,不会出现别脚的“预知性”充放电。
# 滚动时域调度伪代码 horizon = 4 for start in range(0, 24, horizon): model = build_robust_model(start, horizon) plan = solve(model) soc_start = update_soc(plan)这种滚动方案其实更接近实际运行中的调度逻辑,而且计算量反而更小。
7. 从单时段鲁棒到多时段耦合:扩展时的三个方向
最后聊聊后续扩展。我做的是单时段独立的盒式不确定集加SOC耦合,但虚拟电厂调度场景远不止这么简单。
第一个扩展方向是考虑光伏和负荷的时间耦合性。光伏出力有云层移动带来的持续性偏低,负荷波动有“早高峰连着午高峰”的惯性。真正贴合实际的做法是使用多面体不确定集或者基于历史场景的主成分分析降维,限制不确定参数的时序变化斜率,而不是简单限制偏离总量。
第二个方向是引入储能退化成本。我在目标函数里只用了一个很粗糙的充放电损耗系数,但锂电池的循环老化是非线性的,频繁小深度充放和偶尔深度充放对寿命的影响完全不同。加上退化成本之后,鲁棒模型会更倾向于减少不必要的储能动作,让调度计划更温和。
第三个方向是跟实时市场的结合。日前鲁棒计划确定后,日内还有滚动修正的机会,这时可以用模型预测控制MPC的框架,把鲁棒计划作为安全边界,实时跟踪预测误差并做闭环调整。这样既保留了鲁棒优化的稳健性,又避免了过保守带来的成本浪费。
给正在做同类项目的小伙伴一个建议:先从确定性模型跑通,再逐步加入不确定性,每加一类不确定性就用真实场景回测一次。鲁棒优化的核心不在数学有多深,而在于不确定集建得是否符合物理实际、保守度调得是否合理。代码只是最后落地的那一步,前面的物理分析和数据统计才是真正决定调度质量的部分。
我在这套模型上跑出来的结果,和确定性模型相比虽然日成本上升了约18%,但极端场景下几乎没出现过功率失衡,备用机组响应次数也在可控范围内。对虚拟电厂这种对可靠性和经济性都敏感的系统来说,这个代价换来的安全裕度是值得的。希望这份实现笔记能帮你少走几步弯路。