先交代一下项目背景。前阵子我在做一套虚拟电厂微网系统的日前调度方案聚合了分布式光伏、储能、柴油机、工商业可调负荷和上级电网的交互接口。各个环节都有成熟模型但真正让人头疼的只有一件事明天光伏发多少、负荷用多少根本没办法提前精确知道。光伏出力忽高忽低负荷曲线受商户开工率和天气影响飘忽不定按预测值做出来的“最优”计划在实际执行时经常变成次优甚至不可行方案。后来我把思路从确定性优化转成随机优化用蒙特卡洛模拟生成源荷不确定性的海量场景再用快速概率距离削减法压缩场景规模最终建立了一套可落地的日前随机优化调度模型。这篇文章就把整个项目的思路、数学原理、建模细节和踩坑过程完整复盘一遍希望能帮到正在做虚拟电厂、微网或园区级能源调度的同行。1. 确定性日前调度在虚拟电厂场景下的三个致命短板1.1 点预测掩盖了“明天可能有很多种明天”传统日前调度最常用的做法是拿一组预测曲线——比如光伏出力预测、负荷预测——丢进优化模型算出一套功率分配和机组启停计划。这个套路在源荷波动小的时候问题不大但到了虚拟电厂和微网场景至少三个地方会失效。第一个短板是光伏出力预测误差的分布往往不是对称的。上午一次云层过境光伏出力可能从预测值的80%直接掉到20%傍晚天气转晴又可能冲到预测值的120%。你用均值或者中位数预测参与调度等于默认“未来只有一条路”但实际上未来是带状分布甚至双峰分布。更麻烦的是预测误差存在明显的时间相关性——连续阴天时段误差不是独立的而是成片出现。这种相关性用单点预测根本看不到。第二个短板是负荷侧的不确定性跟电源侧叠加不是简单地加一个“备用系数”就能应付。微网里面的工商业负荷有随机性空调启停、产线调整、临时订单都会让负荷曲线偏离基线。哪怕是同一时刻光伏出力偏低和负荷偏高这两个坏消息也可能同时发生。如果调度模型只按预测值安排功率一旦出现这种“双低”或者“双高”组合要么电压越限要么需要临时切负荷甚至会被上级调度按偏差电量考核。第三个短板是备用容量的配置没有概率依据。很多项目做日前计划时习惯性地按照最大负荷的10%或者光伏装机容量的20%留旋转备用但这个百分比从哪来没人说得清。备用留多了经济性变差留少了风险暴露。确定性模型给不出“风险有多大”这个信息因为它内部根本没有概率结构。随机优化模型则可以直接回答问题在某个置信水平下预留多少备用、切负荷概率控制在多少需要付出多少成本。1.2 随机优化不是“多算几次”而是把不确定性写进模型我最早也走过弯路以为随机优化就是拿预测的乐观、中性、悲观三组曲线各算一次然后加权平均。后来发现这个做法属于“多场景确定性优化”它虽然比单场景好一点但本质上没有解决“决策要适应场景变化”的问题。真正的随机优化尤其是两阶段随机规划强调的是“现在做决策时就要考虑未来所有可能的实现并且不能被某个单一场景绑架”。在这个模型里第一阶段决策比如机组启停、日前申报的联络线计划需要在看到光伏和负荷实际值之前就定下来因为电网调度在日前就要收到申报第二阶段决策比如储能实际充放电、切负荷量、柴油机出力调整则可以根据场景实际发生情况进行调整。这种结构天然适合描述虚拟电厂和微网的日前调度你先向系统报一个基准计划但实际运行时可以在一定范围内修正内部资源修正的主体就是第二阶段的灵活性变量。处理源荷双重不确定性时最核心的问题就是如何把连续的概率分布变成有限个离散场景。我的做法就是用蒙特卡洛模拟把不确定性“场景化”再用快速概率距离削减法在保持概率特征的条件下压缩场景数量。下面两节分别展开。2. 蒙特卡洛模拟把不确定的光伏和负荷变成一千个平行世界2.1 先从概率模型说起光伏和负荷应该用什么样的分布要做蒙特卡洛模拟第一步是确定随机变量的概率分布。光伏出力我用的比较多的是贝塔分布。光伏出力的标幺值实际出力/额定容量范围在0到1之间贝塔分布形状灵活可以通过两个形状参数拟合“晴天接近1多云偏中低雨天贴近0”这类偏态分布。负荷预测误差一般用正态分布描述均值取预测值标准差取预测时段对应的历史误差标准差。不过要注意负荷误差在早晚高峰时段往往变异明显最好按小时分别统计标准差而不是全时段用一个固定值。另外还要考虑源荷相关性。光伏偏小的时候负荷未必偏小但如果园区里面有大量空调负荷晴天午后光伏出力高、负荷也高制冷需求阴雨天光伏低、负荷也可能比较低。这种相关性如果无视生成出来的场景会失真。简单处理可以用Cholesky分解先把两个不相关的标准正态随机数做线性组合让组合后的相关系数接近历史样本的相关系数再把组合后的数分别转成光伏和负荷的随机样本。工程上如果历史数据有限可以先算皮尔逊相关系数如果相关性不显著就按独立抽样处理别强行引入相关性。2.2 抽样生成与场景树的形成为什么抽1000次而不是抽10次蒙特卡洛模拟的思路很直白从每个随机变量的分布里抽取一组样本每个样本代表一种“未来”。对于24小时日前调度一个场景就是一条完整的光伏出力和负荷曲线。假设光伏预测为 (P_{pv,forecast,t})实际光伏场景可以写成 实际值 预测值 × 贝塔分布采样值负荷场景可以写成 实际值 预测值 正态误差。把每个时刻的抽样值拼接起来就得到一个场景。抽样次数为什么选1000这是计算精度与代价的折中。按照大数定律样本均值收敛于期望值的速度是 (\frac{\sigma}{\sqrt{N}})也就是说要提升一位精度样本量要扩大100倍。但优化模型不是无限放大生成本事越大越好因为场景越多模型里变量和约束成倍增加。我一般分两步先抽1000个原始场景然后用削减法压到10到20个既保留概率信息又让MILP可解。这比直接抽20个更稳妥因为直接抽少量样本很容易漏掉小概率但又带来高成本的极端场景比如傍晚负荷尖峰叠加光伏出力骤降。抽样时还有几个操作细节一是随机数种子要固定否则每次运行模型结果都变没法对比方案二是对光伏和负荷的时间序列误差建议做上下限截断比如光伏标幺值限定在0.05到1.05之间负荷偏差限定在±3σ避免把明显不合理的场景灌进模型三是每个场景都要记录发生概率初始抽样时所有原始场景的概率默认相等即 (1/N)但削减之后概率就不一样了这个后面会讲到。3. 快速概率距离削减法怎么从一千个平行世界挑出最有代表性的几个3.1 为什么不能直接保留随机抽出的20个场景我一开始图省事想直接保留抽样的前20个或随机挑20个结果发现模型对抽样种子非常敏感换一下随机种子两个方案的调度成本能差出十几个百分点。原因很简单随机抽样的20个点不保证落在原始分布的“高概率区域”里更不保证它们之间的概率距离最小。所谓“概率距离”就是两个概率分布之间的一种差异度量场景削减的目标就是找到一个更小的场景集合使得新旧集合之间的概率距离尽可能小同时保留原有场景的概率权重。业界用得比较成熟的是基于Kantorovich距离简称KD的场景削减方法。Kantorovich距离可以理解为“把一个概率分布的物资搬到另一个概率分布所需的最小运输成本”。在离散场景集合里KD定义为所有场景对之间距离乘概率权重的最优传输结果。削减算法就是在保持KD尽量小的原则下不断剔除冗余场景把被删场景的概率合并到离它最近的保留场景上。3.2 同步回代削减的具体步骤和“快速”体现在哪实际工程中应用最广的快速概率距离削减法学名是同步回代削减Simultaneous Backward Reduction。它的核心是一个贪婪迭代过程每次删掉一个对总距离贡献最小的场景计算所有保留场景两两之间的距离 (D(i,j))距离可以用欧氏距离也可以根据调度目标用加权距离比如对晚高峰时段误差给更高权重。对每个场景 (i)找到它到最近的其他场景 (j) 的距离 (d_i \min_j D(i,j))乘上它的概率 (p_i)得到 (p_i d_i)。找出 (p_i d_i) 最小的那个场景把它删掉然后把这个场景的概率加到它最近的场景上。更新最近邻关系和概率重复直到保留场景数达到目标个数。“快速”体现在两个地方第一距离矩阵只需要在开始时完整计算一次后续删除场景后只更新与被删场景相关的行和列第二每次不需要重新全局搜索只需要维护每个场景的“当前最近邻居”候选集合。实测环境里从1000个场景削减到10个场景单次计算用时大概在几十毫秒到一两秒之间比重新做聚类快得多而且不需要预机器数。另一个替代方案是K-means聚类但不是所有项目都适合。聚类方法倾向于把场景按几何形态分组可能导致小概率但经济性冲击很大的极端场景被“平均”到普通场景里而基于概率距离的削减会概率加权平衡对稀有高损失场景保留得更完整。如果希望极端风险也能体现在模型中我推荐用同步回代削减而不是聚类。3.3 削减后场景质量的检验削减完之后不能直接拍板我一般会做两个检验。第一看削减前后集合的Kantorovich距离变化是否在可接受范围。如果削减到15个场景时距离已经明显抬升说明继续削减会破坏分布结构。第二把削减后各场景的关键统计量光伏总出力期望、负荷峰值期望、光伏低于某阈值的概率跟原始1000个场景的统计量对比偏差控制在1%-3%以内基本可用。一个常见的认知误区是削减后场景数量越少计算越快因此希望削减到5个甚至3个。但在含有切负荷惩罚、购售电价日内波动明显的模型里5个场景往往不足以表现“光伏高出力低谷价”和“光伏低出力高峰价”两种完全不同运行状态。我的经验是针对日前调度场景数不要少于10如果机组组合变量很多则15到20个是性价比最高的区间。4. 随机优化模型的数学骨架目标函数、约束与两阶段决策结构4.1 决策变量的两阶段划分哪些要事先承诺哪些可以事后调整建模前必须想清楚哪些决策属于“事前承诺”哪些决策属于“事后调整”。以我做的虚拟电厂微网为例第一阶段日前计划不随场景改变包括柴油机的开停机状态 (u_t)与上级电网的日前购售电计划 (P_{grid,t}^{day})储能系统是否参与调度的充放电状态二元变量如果需要。第二阶段日内响应依赖场景 (s)包括柴油机出力的调整量 (P_{diesel,t,s})储能实际充放电功率 (P_{ch,t,s}, P_{dis,t,s})可调负荷的削减量 (P_{shed,t,s})对联络线计划的小幅修正量。在申报场景中第一阶段的日前联络线计划是刚性的但在一些内部微网调节能力很强的场合联络线可能需要允许一定比例的范围调整。我处理的方式是引入一个松弛变量 (\Delta P_{grid,t,s})允许联络线在日前计划附近做有限调整同时对该调整设置高惩罚系数。这样模型不会死守一个可能在极端场景下根本不可行的计划也不会随意打破与上级电网的合约。4.2 目标函数期望成本而不是单一场景成本随机优化的目标函数是“所有场景运行成本的期望值最小化”一般可以写成[ \min \sum_{t1}^{T} \left( c_{diesel} P_{diesel,t}^{base} c_{grid,buy} P_{grid,t}^{day} - c_{grid,sell} P_{sell,t}^{day} \right) \sum_{s1}^{S} p_s \cdot \sum_{t1}^{T} C_{adjust,t,s} ]其中 (p_s) 是削减后各场景的概率注意这里不再是等概率后面求解代码里最容易出错的就是忘记用 (p_s) 加权。(C_{adjust,t,s}) 包含场景内调整成本柴油机额外燃料成本、储能充放电循环老化成本、切负荷惩罚成本、联络线调节惩罚成本等。对于所有场景都相同的基准项可以直接用预测期望值计算但是机组启停、储能决策如果涉及整数变量期望展开后会导致S倍大的MILP。实际工程中我会把第一阶段启停变量和第二阶段运行变量的耦合写清楚目标函数里所有成本项目都乘上场景概率否则结果会出现不同场景分配权重异常的情形。4.3 约束条件从功率平衡到储能SOC的时序耦合约束是随机优化模型最容易膨胀的地方。考虑到每个场景都要满足约束我这里列几个典型约束功率平衡约束每个场景、每个时刻 [ P_{pv,t,s} P_{diesel,t,s} P_{dis,t,s} P_{grid,t,s} P_{load,t,s} P_{ch,t,s} P_{shed,t,s} ] 其中 (P_{grid,t,s} P_{grid,t}^{day} \Delta P_{grid,t,s})。如果微网允许向电网倒送电则 (P_{grid,t,s}) 可正可负。柴油机出力约束出力在上下限之间并且爬坡约束也在每个场景内成立 [ P_{diesel,t-1,s} - P_{diesel,t,s} \le R_{down}, \quad P_{diesel,t,s} - P_{diesel,t-1,s} \le R_{up} ]储能SOC约束 [ SOC_{t,s} SOC_{t-1,s} \eta_{ch} P_{ch,t,s} - \frac{P_{dis,t,s}}{\eta_{dis}} ] 同时SOC要在0和容量之间。这里有个重要建模决策储能SOC是否允许在不同场景间“分叉”。如果SOC也是第二阶段的调整量那么不同场景下SOC会不同日内调度时储能可以根据实际源荷情况重新决策这更真实。但此时要在目标函数中加入SOC末态惩罚避免模型为了套利把电池无限充电或放空。可调负荷约束削减量有上限且削减总时长有上限避免用户侧被过度频繁地切除。联络线调整约束(\Delta P_{grid,t,s}) 的绝对值一般限制在日前计划的一定百分比内这一条约束决定了虚拟电厂对外表现的“可预测性”。这些约束按场景展开后一个24时段、15场景、含储能的模型可能产生几千个连续变量和上百个整数变量。只要用MILP求解器配合线性化求解时间一般可控制在几分钟内。5. 一个能直接跑的简化版本建模工具选择与核心代码骨架5.1 工具选择PythonPyomo/Gurobi是目前最省事的组合研究阶段我用过MATLABYalmip部署阶段则更倾向PythonPyomo。Pyomo的好处是建模语法跟代数表达式非常接近方便把两阶段随机模型写清楚而且求解器可以无缝切到Gurobi、Cplex或开源求解器HiGHS。如果项目只允许开源工具HiGHS对MILP的支持已经很让人意外处理中小规模日前调度模型够用。Gurobi的优势是求解MIP的MIP gap收敛明显更快特别是整数变量多的时候。场景生成和削减建议用独立的Python函数写好不要混进优化模型代码。一方面方便单独测试另一方面削减算法和优化模型解耦后后续如果加入CVaR约束或鲁棒约束只需要替换场景集输入不需要动模型主体。5.2 场景生成与削减的代码骨架场景生成部分的核心代码框架大概是import numpy as np np.random.seed(42) n_scenarios 1000 n_periods 24 # 光伏出力预测值 beta分布抽样 pv_samples np.random.beta(a2.5, b2.0, size(n_scenarios, n_periods)) pv_scenarios pv_forecast * pv_samples # 负荷误差正态分布每个时段单独标准差 load_error np.random.normal(loc0, scaleload_std, size(n_scenarios, n_periods)) load_error np.clip(load_error, -3 * load_std, 3 * load_std) load_scenarios load_forecast load_error削减函数如果用同步回代削减需要自己实现距离矩阵循环删除。实际写的时候我建议先用小规模数据验证删除逻辑比如先用100个场景削减到5个画出削减前后光伏出力分布直方图确认概率没有发生异常偏移后再上1000个场景。另外要强调削减后每个场景的概率不是1/目标场景数。这是新手最容易错的地方。例如原始1000个等概率场景削减到15个场景后每个场景概率是从多个原始场景合并来的最大概率可能超过0.2最小可能只有0.02。在随机优化模型里权重一定用scenario_prob[s]不能默认为1/S。5.3 Pyomo模型骨架示例把模型写成一个确定性等价MILP核心结构如下import pyomo.environ as pyo model pyo.ConcreteModel() model.T pyo.Set(initializerange(24)) model.S pyo.Set(initializerange(n_reduced_scenarios)) # 第一阶段变量 model.u_diesel pyo.Var(model.T, withinpyo.Binary) # 柴油机启停 model.grid_day pyo.Var(model.T, bounds(-max_export, max_import)) # 第二阶段变量 model.p_diesel pyo.Var(model.T, model.S, bounds(0, diesel_max)) model.p_ch pyo.Var(model.T, model.S, bounds(0, ch_max)) model.p_dis pyo.Var(model.T, model.S, bounds(0, dis_max)) model.soc pyo.Var(model.T, model.S, bounds(soc_min, soc_max)) model.shed pyo.Var(model.T, model.S, bounds(0, shed_max)) model.grid_adj pyo.Var(model.T, model.S, bounds(-adj_limit, adj_limit)) def objective_rule(m): base_cost sum(...) # 第一阶段成本 scen_cost sum(scenario_prob[s] * sum(...) for s in m.S for t in m.T) return base_cost scen_cost model.obj pyo.Objective(ruleobjective_rule, sensepyo.minimize) def power_balance_rule(m, t, s): return (pv[t, s] m.p_diesel[t, s] m.p_dis[t, s] m.grid_day[t] m.grid_adj[t, s] load[t, s] m.p_ch[t, s] m.shed[t, s]) model.power_balance pyo.Constraint(model.T, model.S, rulepower_balance_rule) def soc_update_rule(m, t, s): if t 0: return m.soc[t, s] soc_init ch_eff * m.p_ch[t, s] - m.p_dis[t, s] / dis_eff return m.soc[t, s] m.soc[t-1, s] ch_eff * m.p_ch[t, s] - m.p_dis[t, s] / dis_eff model.soc_update pyo.Constraint(model.T, model.S, rulesoc_update_rule)如果你的模型量级更大比如包含多台机组、多个储能并且需要更长的调度周期建议把场景同时切块并行建模型最后汇总结果。Pyomo中可以用ConcreteModel的组件高效索引不必写几十个重复变量使用Constraints和Objective时注意时间索引的连续性避免因为某个时段缺约束导致SOC“断点”。6. 算例对比结果与六个容易忽视的工程细节6.1 确定性模型与随机优化模型的算例对比我手头一个典型算例光伏装机10MW负荷峰值12MW储能容量5MWh柴油机容量4MW联络线功率上限6MW光伏预测误差标准差0.12pu负荷预测误差标准差0.08pu。先用单场景确定性模型算一次备用按光伏20%加负荷10%再用蒙特卡洛生成1000个场景削减到15个场景做随机优化。两种方案都用同一组历史场景回验。指标确定性模型随机优化模型15场景日前计划运行成本期望值基准降低4.2%实际运行切负荷量期望18.6 MWh7.3 MWh联络线调整次数11次5次求解时间4秒137秒这个结果很能说明问题确定性模型的日前计划看着便宜但因为没考虑光伏实际高出力和负荷低期望的“错位”实际运行时需要频繁调整切负荷量和联络线调整次数都明显偏高。随机优化模型虽然求解慢了两分钟但实际运行成本更低且对外表现更平稳。这就是从“点预测最优”到“期望最优”的价值。6.2 坑1削减后概率权重用错这个问题我在测试阶段就踩过。代码里场景概率最开始写成了1 / n_reduced结果每个场景平均概率只有0.067但实际削减后最大的场景概率接近0.19。目标函数被严重低估调度结果偏好削峰填谷力度不够。修正成scenario_prob[s]后模型才会知道高概率场景必须在功率平衡里被优先满足。6.3 坑2储能的SOC在场景间能否独立如果不小心把SOC设成所有场景共享同一个变量序列模型会把电池调度“平均化”导致某些场景过度充电、某些场景电量不够。正确做法是让SOC跟随场景独立变化但在目标函数里给末时段SOC设定目标值比如强制回到初始SOC的95%以上否则模型会为了贪便宜把所有电能提前放完。另外如果储能参与的是日前申报的固定计划那储能充放电本身变成第一阶段变量这时另一个问题是日内偏差无法被储能消化备用压力会更高。6.4 坑3场景削减数量与风险损失的关系建议做一次敏感分析分别削减到5、10、15、20、30个场景观察期望成本和最大场景下切负荷惩罚的变化。很多模型在10个场景时成本已经收敛但切负荷惩罚在极端场景下仍然偏高这种情况说明需要增加场景数或引入风险约束比如CVaR来压制极端场景。我最终选了15个场景就是综合了成本收敛和极端场景惩罚两个指标的结果。6.5 坑4随机数种子和分布参数的科学性蒙特卡洛模拟如果每次换种子结果都变化说明场景样本容量不足或者削减后场景代表性不够。至少对比3个不同种子如果调度成本波动超过5%就要增加原始场景数或调整削减策略。另外分布参数不能拍脑袋光伏贝塔分布的形状参数建议从历史出力数据用最大似然估计拟合负荷标准差要分时段统计早高峰和午间的误差特性完全不同。6.6 坑5极端场景被“平均”掉快速概率距离削减法虽然有概率加权但仍然有人为设置最大保留概率或最小保留概率的空间。有些真实场景很极端但发生概率只有0.001对期望成本影响小却可能让调度计划在某些几乎不发生的时刻失效。如果你所在地区对考核电量很敏感建议在削减后手动补充几个极端场景比如“光伏低出力和负荷高峰同时出现”把它们的概率设得保守一点。这种操作在工程上很常见本质是把概率约束换成鲁棒性兜底。6.7 坑6不要把随机优化当成万能药最后说一个认知层面的坑。随机优化解决的是“在已知概率分布下如何做期望最优决策”但如果分布本身估计错误模型一样会翻车。尤其光伏预测误差分布会随季节变化冬季和夏季的贝塔分布形状参数差异很大。实际部署时最好按月份重新拟合分布参数模型参数不要一套用一年。我能给的最终建议是做日前调度先用历史回验跑三个月观察随机优化模型的切负荷率和成本波动确认稳定后再切正式运行不要第一天上线就指望它比老调度员凭经验做出来的计划更稳。回头复盘这个项目最大的收获不是把模型跑通而是在处理源荷双重不确定性时找到了“场景生成—场景削减—随机优化”这一套可复用的方法论。蒙特卡洛模拟负责提供足够丰富的未来可能性快速概率距离削减法负责把可能性压缩成可计算的代表性方案随机优化模型则把风险成本量化到了每个时段的决策里。这套思路不仅适用于虚拟电厂微网日前调度也适用于园区级综合能源系统、储能电站日前充放电策略优化等场景。如果你正在被预测误差折磨不妨先从小规模算例开始把场景削减和概率权重这两件事做扎实再逐步扩展模型复杂度。