做虚拟电厂调度项目有一阵子了最让我觉得有意思的其实是处理不确定性这件事本身。光伏出力和负荷波动就像两个不听话的变量单拎出来任何一个做调度都还有章可循一旦同时把它们丢进优化模型问题就从“确定性的线性规划”直接升级成“既要抗扰动、又要保经济性”的鲁棒优化。这篇就把我在实现过程中用到的核心思路、数学模型和代码整理出来给正在搞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, nonnegTrue) p_sell cp.Variable(T, nonnegTrue) p_gen cp.Variable(T) p_ch cp.Variable(T, nonnegTrue) p_dis cp.Variable(T, nonnegTrue) 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(solvercp.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, nonnegTrue) p_sell cp.Variable(T, nonnegTrue) p_gen cp.Variable(T, nonnegTrue) p_ch cp.Variable(T, nonnegTrue) p_dis cp.Variable(T, nonnegTrue) soc cp.Variable(T) u_buy cp.Variable(T, booleanTrue) # 参数 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(solvercp.GUROBI, TimeLimit60) print(最优成本:, prob.value) print(储能SOC曲线:, soc.value)有几个细节必须提一下。首先是充放互斥约束我这里直接写了 p_ch_t * p_dis_t 0这在CVXPY里是非线性非凸约束GUROBI处理起来会报错或转成MIP求解速度慢很多。更好的做法是引入布尔变量表示充放电状态然后加互斥约束。实际调试时我先用线性大M方式u_ch cp.Variable(T, booleanTrue) u_dis cp.Variable(T, booleanTrue) 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.412/246.28312.715/247.116338.918/248.024365.220/248.448完全保守412.624/249.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%但极端场景下几乎没出现过功率失衡备用机组响应次数也在可控范围内。对虚拟电厂这种对可靠性和经济性都敏感的系统来说这个代价换来的安全裕度是值得的。希望这份实现笔记能帮你少走几步弯路。