源荷双重不确定性下微网/虚拟电厂日前随机优化调度MATLAB实现

1. 为什么虚拟电厂/微网调度要逼自己正视双重不确定性虚拟电厂和微网这类多能互补系统的运行决策最磨人的往往不是设备选型或拓扑搭建而是明天到底怎么过。风电出力看气象脸色光伏曲线被云层牵制负荷端更是没有哪一天完全复刻前一天。我早前写调度代码时习惯直接抓一条预测曲线塞进优化模型跑出一份又漂亮又整齐的机组组合方案。但真正拿到运行日去对照风电出力偏出预测值二三十个百分点都不稀奇负荷侧的偏差也常年盘踞在5%到10%。预测时有多自信,执行时就有多狼狈。这种纸面完美、实际翻车的体验让我把研究重心从确定性优化挪到了随机优化。所谓基于源-荷双重不确定性的日前随机优化调度核心思路并不复杂在建模阶段把源侧风电、光伏和荷侧负荷的不确定性显式表达出来用一批可能存在的情景代替单条预测曲线让日前计划在未来可能出现的多种情形下综合表现最优。这篇文章就是围绕一套可直接复现的MATLAB代码把我踩过的坑、建模时的取舍、代码实现的关键节点全部梳理出来。这篇笔记适合正在做微网、虚拟电厂、园区综合能源调度方向的同学或工程师需要具备线性规划基础懂得基本的机组组合和储能充放电模型。代码层面用的是YALMIP加CPLEX/Gurobi的组合如果你手头没有商业求解器我也会在后面给出基于MATLAB内置求解器的替代思路。2. 为什么不算清账日前计划就只能当废纸看2.1 确定性模型的最后一公里总在关键时刻掉链子先把确定性模型的毛病说透。传统日前调度是这么玩的把风电、光伏、负荷的预测值当作已知参数建立一个确定性的优化问题求解出一组机组出力和储能充放电计划。在预测精度很高的极端理想场景下这样做没错。可实际问题在于源侧和荷侧的不确定性是叠加放大的不是简单相加。举个例子某天光伏预测出力是100MW实际因为云层遮挡只发出来70MW。与此同时负荷预测是80MW实际因为气温骤降升到88MW。两边一夹击系统的净负荷缺口就是100-7088-8038MW而不是你计划里那20MW的裕度。如果调度方案里没有足够的灵活调节能力38MW的缺口就得靠切负荷或紧急购电来填成本直接翻倍。我见过很多初学调度的朋友犯同样的错误在目标函数里加了弃风弃光惩罚却完全没考虑负荷侧的随机波动觉得风电光伏才是变量负荷预测挺准的。实际上负荷预测误差在温度剧烈变化的季节非常可观尤其在节假日或大型赛事活动期间误差可以轻松突破15%。只对源侧做鲁棒优化、忽略荷侧波动本质上还是在赌。2.2 源荷相关性这个隐藏坑比你想的更值得留意再往上走一步源侧和荷侧并不是互相独立的随机变量。光伏出力和空调负荷在夏季午后呈现很强的正相关——太阳辐射越强光伏出力越大同时气温升高带动制冷负荷爬升。这在数学上意味着如果建模时假设源、荷完全独立生成场景集合的方式就是有偏的最终求出来的调度方案也会偏离真实最优。比较讲究的做法是在场景生成阶段引入相关系数矩阵或者采用Copula函数把边缘分布黏合起来。对多数工程应用来说用相关系数矩阵修正拉丁超立方采样LHS就足够了。我在实际项目中常用经验Copula加上Pearson相关系数修正计算量可控效果也稳定。代码里我会固定一个随机种子方便复现和对比实验。2.3 从一根曲线到一棵场景树的认知转变随机优化的本质是把单点预测变成概率场景集合。想象你从明天风电一定是80MW切换到明天风电有30%可能是60MW40%可能是80MW30%可能是100MW——决策就从找一个固定最优解变成在所有可能出现的情况下找一个期望意义上的最优解。这正是场景树或者说多场景随机规划的基本形态。这种转变带来的直接收益是调度方案不再在某一个特定预测值下最优而是在一整片概率空间里稳健。代价是模型规模膨胀、求解时间变长、代码复杂度上升。这两者的权衡正是做随机优化调度最需要拿捏的地方。3. 源、荷不确定性怎么建模从概率分布到场景集的完整链路3.1 风电与光伏的概率模型到底怎么选风电出力的随机性来源于风速波动。工程上最常用的风速模型是两参数Weibull分布概率密度函数长这样f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k是形状参数c是尺度参数。实际应用中一般通过历史风速数据拟合得到k和c然后借助风电机组的功率曲线把风速映射成出力。功率曲线通常是一个分段函数切入风速以下出力为0额定风速以上出力封顶中间段近似线性或三次曲线关系。光伏出力模型比风电稍微复杂一点。常见的做法是先对光照辐照度建模——多数文献采用Beta分布因为它在[0,1]区间上有良好的灵活性能够逼近各种偏态分布。然后通过光伏面板的等效模型把辐照度换算成功率。温度对光伏转换效率的影响也不能无视工程上常用温度修正系数组件温度每上升1摄氏度效率大约下降0.4%到0.5%。3.2 负荷侧预测误差的分布假设负荷预测误差的经验分布虽然存在尖峰厚尾特征但在大多数日前调度研究中为了保持模型线性大家默认接受正态分布假设。设预测负荷为P_load_fc实际负荷P_load可以表示为P_load P_load_fc ε_load其中ε_load服从均值为0、标准差为σ_load的正态分布。σ_load的取值通常为预测负荷的5%到10%。需要提醒的是这种正态假设在极端天气事件下会失效如果要做极端场景分析建议用分位数或非参数估计来替代。3.3 场景生成和缩减蒙特卡洛与同步回代消除有了每个随机变量的分布形式接下来就要生成足量的离散场景。初始场景生成我一般用两种方法。第一种是朴素蒙特卡洛直接按照各随机变量的分布随机抽样每个样本点组成一个场景。优点是实现简单缺点是需要的样本量很大收敛速度只有O(1/√N)而且容易出现样本聚堆的情况。第二种是拉丁超立方采样它的思路是把每个随机变量的累积分布函数等分成N个区间在每个区间内各抽取一个样本然后通过随机洗牌把各个变量的样本组合成N个场景。LHS的优点是覆盖均匀同样的场景数下比朴素蒙特卡洛精度高不少。具体实现核心代码如下% 拉丁超立方采样生成源-荷初始场景 n_scen 500; % 初始场景数 n_var 3; % 变量数风电、光伏、负荷 rand(seed, 2024); % 固定随机种子方便复现 for i 1:n_var u (rand(n_scen, 1) (0:n_scen-1)) / n_scen; % 分层抽样 X(:, i) icdf(dist_obj{i}, u); % 逆变换采样 end % 对X的列进行随机重排破坏分层带来的相关性 % 再按相关系数矩阵做Cholesky修正场景生成之后紧接着是场景缩减。500个场景直接丢进优化模型决策变量和约束数量会爆炸求解时间动不动几个小时这在日前调度里是不可接受的。标准的做法是用同步回代消除也叫快速前代消除每次迭代删除一个与原场景概率距离最小的场景并把被删场景的概率累加到距离它最近的场景上。最终把场景数从几百个压到10到20个既保留概率分布的主要形态又把模型规模控制在可解范围内。% 同步回代消除的核心思路 % 1. 计算场景两两之间的距离矩阵 dist(i,j) % 2. 对每个场景i找到最近邻 j % 3. 寻找最小化 概率加权距离 的场景对 % 4. 删除其中一个场景概率累加到另一个 % 5. 重复直到达到目标场景数这一步实操时我通常会多看两眼缩减后的场景集合如果某些场景概率过高比如超过40%说明缩减过头丢失了尾部风险信息如果概率分布过于平均说明缩减程度不够还可以再压一压。场景数和概率分布的合理性直接决定后续优化结果的可靠性。4. 随机优化调度的数学模型目标函数、约束体系与决策变量怎么定4.1 目标函数里各项成本的博弈日前随机调度的目标函数我习惯写成最小化期望总运行成本。它由四个部分构成从主网购电/售电的功率交互成本、微型燃气轮机的燃料成本、储能充放电的折旧成本以及弃风弃光和失负荷的惩罚成本。网电成本和燃料成本是显性成本数值上可以用分时电价和燃料单价直接算出。储能折旧成本容易被初学者漏掉但它其实非常关键——如果储能每次充放电不计损耗优化模型会拼命让储能频繁动作得到一份理论上漂亮但根本跑不动的方案。我给储能折旧设一个单位功率成本大概在每兆瓦时几十块钱的量级。弃风弃光和失负荷惩罚成本的设置本质上是给优化模型一个不得已而为之的出口。惩罚系数不能太高太高会让模型为了规避极小概率的极端场景过度保守也不能太低太低会让模型在概率很小的场景里随意切负荷。我的经验值是失负荷惩罚取实时电价的五到十倍弃风弃光惩罚取两到三倍。这个比例可以在敏感性分析时再细化。目标函数用数学式子表达就是这样min Σ_t ( Σ_i C_i·P_i(t) C_grid(t)·P_grid(t) C_ESS·|P_ess(t)| λ_w·P_curt_w(t) λ_l·P_loss(t) )注意这里所有带不确定性的量都按场景取了期望所以严格写法是在每个场景下算成本、再对各场景概率加权求和。4.2 约束条件逐条拆解约束体系是整个优化模型的大头拆开来看有这么几类。功率平衡约束是最核心的一条在任意时间段微网内部的电源出力加上从主网购买的电量必须等于负荷加上储能充电功率。在场景s下它写成P_wind(s,t) P_pv(s,t) Σ P_i(s,t) P_grid(s,t) P_load(s,t) P_ess_ch(s,t) - P_ess_dis(s,t) P_curt_w(s,t)注意等式两边的每条量都带场景下标这意味着约束也要随场景扩展。机组运行约束包括出力上下限、爬坡速率约束、最小启停时间约束。燃气轮机这类可控机组还需要考虑启停状态变量这会把模型推入混合整数线性规划的范畴。储能约束是另一块大头。储能SOC动态方程、SOC上下限、充放电功率上下限、以及同一时刻不能同时充放。SOC动态方程写起来要小心时间步长的换算SOC(t1) SOC(t) η_ch·P_ch(t)·Δt - P_dis(t)·Δt/η_disη是充放电效率一般取0.9到0.95。初看没什么但如果你把时间步长设置成15分钟Δt就是0.25功率单位是MW能量单位是MWh单位换算错一位整条SOC曲线就会歪掉这是我早年踩过的实打实的坑。电网交互约束也不可忽视。虚拟电厂和配电网之间的联络线功率有上限买卖电的功率方向约束通常用一个大M法转化成线性约束- P_line_max ≤ P_grid(t) ≤ P_line_max4.3 为什么选用两阶段随机规划结构我把机组启停这类今天就得拍板的决策放到第一阶段把各场景下的出力调整、储能充放电功率放到第二阶段。原因很实际日前市场要求你在运行日前一天就把机组组合确定下来但具体出力还可以等到实际运行日根据风光和负荷的实时变化灵活调节。这种两阶段结构用数学语言写出来就是第一阶段的决策变量对所有场景保持一致第二阶段的决策变量带场景下标、允许随场景变化。在YALMIP里实现时第一阶段变量定义为普通sdpvar第二阶段变量定义为维度场景数×时段数的sdpvar矩阵约束用for循环逐场景写。它比单阶段鲁棒优化灵活又不至于像全场景随机规划那样臃肿。5. MATLAB代码实现从场景生成到YALMIP建模的完整流程5.1 代码整体框架与模块拆分再看回这套MATLAB代码整体结构分四块参数初始化、场景生成与缩减、优化建模与求解、结果分析与绘图。模块拆分要清晰尤其场景生成和优化求解之间通过一个标准化数据结构衔接。我通常用结构体数组来存场景数据每个场景包含风电、光伏、负荷三条时序曲线以及它的概率权重。% 主流程示意 params load_case(case_study); % 加载基础参数 [scen, prob] generate_scen(params); % 生成并缩减场景 result run_stochastic_dispatch(params, scen, prob); plot_dispatch_result(result);这种模块化设计的最大好处是后面想换场景生成算法、换目标函数结构都不用推倒重来。5.2 场景生成模块的核心代码场景生成模块我用蒙特卡洛加Cholesky相关修正来生成初始场景然后用同步回代消除做缩减。这里给一个场景生成的简化版本function [ScenAll, ProbAll] generate_scen(params) n_init params.n_init; % 初始场景数例如500 n_keep params.n_keep; % 缩减后场景数例如15 T params.T; % 调度时段数例如24 % 1. 拉丁超立方采样生成初始样本 wind_scen lhsnorm(params.wind_mu, params.wind_sigma, n_init); pv_scen lhsnorm(params.pv_mu, params.pv_sigma, n_init); load_scen lhsnorm(params.load_mu, params.load_sigma, n_init); % 2. 相关性修正源-荷相关系数矩阵 R [1, params.rho_wp, params.rho_wl; params.rho_wp, 1, params.rho_pl; params.rho_wl, params.rho_pl, 1]; L chol(R, lower); Z [wind_scen; pv_scen; load_scen]; Z_corr L * Z; ScenAll Z_corr; % 3. 概率均等赋给每个初始场景 ProbAll ones(n_init, 1) / n_init; % 4. 同步回代消除从n_init缩减到n_keep [ScenAll, ProbAll] scen_reduction(ScenAll, ProbAll, n_keep); end场景缩减函数scen_reduction是基于概率距离的迭代消除这部分代码稍微长一点核心逻辑我上面已经说过。需要留意的是缩减后要重新归一化概率确保所有场景概率之和为1否则优化求解时可能出现很隐蔽的数值问题。5.3 YALMIP建模与求解器配置优化建模部分我用YALMIP它让约束表达非常接近数学语言可读性好、不容易出错。以下是功率平衡约束和储能约束的示例% 定义决策变量 P_grid sdpvar(3, T, 24); % 网电交互功率 P_ess sdpvar(3, T, 24); % 储能充放电正为放电负为充电 SOC sdpvar(3, T, 24); % 储能SOC u_gas binvar(T, 1); % 燃气轮机启停状态第一阶段变量 % 功率平衡约束对每个场景、每个时段 Constraints []; for s 1:n_keep for t 1:T Constraints [Constraints, ... P_wind(s,t) P_pv(s,t) P_gas(s,t) P_grid(s,t) ... P_load(s,t) P_ess_ch(s,t) - P_ess_dis(s,t)]; end end求解器配置方面CPLEX是我用得最多的选择Gurobi也不错两者对MILP的求解性能都很强。YALMIP里指定求解器只需要一行ops sdpsettings(solver, cplex, verbose, 1, showprogress, 1); sol optimize(Constraints, Objective, ops);如果你没有CPLEX或Gurobi的授权替代方案是用MATLAB自带的intlinprog。YALMIP也支持intlinprog作为底层求解器。差别主要在求解速度上场景数少比如10个以内时intlinprog完全能顶住再往上就推荐商业求解器了。5.4 结果输出与调度方案解析求解完成后的输出我用四张图来呈现第一张是各场景下的净负荷曲线和日前调度出力的堆叠图重点看功率平衡是否在每个时段都闭合。第二张是储能SOC曲线检查SOC是否一直保持在合理范围内、有没有违反约束的跳变。第三张是各场景概率加权后的期望成本分解直观看出钱花在网电、燃气还是惩罚项上。第四张是随机优化和确定性优化的调度方案对比。绘图时有一个小技巧先判断sol.problem是否为0YALMIP返回值非0说明求解失败或数值异常。遇到这种情况优先检查场景概率是否归一化、变量维度是否匹配、有没有意外的Inf或NaN。6. 调试运行中避不开的坑与改进思路6.1 初始场景数太大造成的维度灾难我最开始做场景生成时初始场景数拍脑袋选了2000个缩减到50个结果YALMIP建模时约束矩阵直接撑爆内存MATLAB卡死。后来把初始场景降到500、缩减目标调到15问题迎刃而解。这里面有个基本权衡场景太多模型规模指数膨胀场景太少不确定性刻画失真。我的经验值是初始300到500、最终8到15具体取决于你的机组数量和时间粒度。如果时间分辨率是15分钟96个时段场景数再翻倍模型会非常吃力这时候优先缩减场景而不是压缩时段。6.2 收敛性与求解时间的平衡策略MILP求解时间对场景数和整数变量数非常敏感。燃气轮机启停变量是整数变量每个机组每个时段一个二进制变量场景再多整数变量的规模就是T乘机组数。如果你发现求解时间超过半小时可以尝试这几种手段一是给求解器设置MIP相对间隙MIP gap容忍度CPLEX里设置成0.01或0.02允许1%到2%的次优性求解时间通常能缩短一半以上。二是启用YALMIP的冗余约束消除或边界收紧选项有些时候模型本身约束写冗余了求解器预处理会花掉大量时间。三是调整场景缩减策略把概率权重极端化的场景剔除。6.3 从代码到实际工程场景的移植扩展这套代码的逻辑核心是场景树加两阶段随机规划换应用场景只需要替换三块参数文件、场景生成模块的概率模型、目标函数里的成本系数。比如从微网扩展到虚拟电厂聚合多个分布式资源只要把各个分布式电源的出力变量和成本项加进去约束体系跟着扩展就行。再比如加入电动汽车充电负荷弹性只需要把负荷变量从固定参数改成可调变量加上充电桩的功率约束。6.4 我个人的一点体会最后说点实在的。做随机优化调度最关键的往往不是数学公式多漂亮、求解器多高端而是你对不确定性本身有没有敬畏心。一套调度方案如果只在预测曲线下表现完美在真实场景里就大概率翻车。反过来一味增加场景数量追求更优而不考虑计算代价代码也落不了地。我从这套代码里学到的最重要一课是先跑通10个场景的简化版把逻辑完全理顺再逐步扩展场景规模。整个过程里固定随机种子、保留每次实验的场景集、记录求解时间这三件事让我少走了一半弯路。如果你想拿这套代码做研究对比强烈建议把确定性模型和随机模型放在同一个测试集上跑分别统计期望成本、最差场景成本和弃风弃光率那一瞬间你会直观感受到不确定性四个字的分量。