简介MATLAB复现的虚拟电厂与配电公司协调调度模型中的场景削减模块面向电力系统优化调度、新能源不确定性建模方向的研究者与工程师尤其适合处理含高比例新能源的复杂不确定性场景。代码聚焦风电、光伏及电价场景模拟由一组确定性方案生成50个光伏场景再采用基于概率距离的快速削减算法将场景压缩至5个运行后直接输出削减前后的场景及对应概率既降低大规模场景带来的内存与计算负担又保留主要概率特征可直接嵌入含风光水的协调调度、备用优化及配电公司决策等模型。资源共5个文件以4个m脚本为主体分别对应风电、光伏、电价等场景生成与削减流程另有1个txt参考文献说明压缩包仅5KB轻量易读、结构清晰便于按需修改参数后复现或二次开发。该资源已有761人学习适合具备一定MATLAB基础、希望快速落地场景削减算法或开展虚拟电厂协调调度研究的研究生与工程技术人员。1. 含风光水的虚拟电厂与配电公司协调调度为什么偏偏要谈场景削减“场景削减”四个字看着像数据预处理的小步骤实际是含风光水的虚拟电厂VPP与配电公司DSO协调调度能否落地的胜负手。风电、光伏的出力曲线本质上是从随机分布里采出来的一批“可能明天长这样”的样本水电又受来水和库容双重约束三者叠加之后优化模型的约束矩阵会被几十上百条时序场景撑到难以直接求解。越想把随机性刻画得细场景数越多混合整数规划的求解时间就越失控——这就是为什么任何讲VPP与DSO协调调度的复现项目都要单拎出“场景削减部分”作为独立模块来谈。本文按照“先立模型、再造场景、再谈削减、最后校验”的顺序拆解。读者如果手上有一份还没跑通的MATLAB复现工程或者正在决定削减算法用K-means还是快速前向选择应该能从里面对上号。至于为什么不让场景削减走“直接抽样取均值”这种简单路第三、四章会给出可量化的答案。先说明一个基本立场下面所有代码按MATLAB R2022b及以上版本组织优化求解部分依赖YALMIP配cplex或gurobi场景生成与削减仅用MATLAB基础工具箱不需要额外安装包。2. 协调调度模型的目标函数与约束体系先把边界定清楚2.1 VPP与配电公司的交易框架谁向谁买电按什么价格买虚拟电厂聚合了风、光、水三种分布式资源外加储能和可调负荷对外表现为一个可以统一调度的发电主体。配电公司在这个框架里扮演两个角色一是网络运营商保证配电网节点电压和支路潮流不越限二是购电方从VPP收购电力以满足负荷需求。协调调度模型要做的事就是在给定的调度周期内通常取24小时分辨率1小时或15分钟算出VPP内部各机组的最优出力计划和VPP与DSO之间的交换功率使系统总运行成本最小同时保证网络安全约束。目标函数可以分解成三块VPP内部发电成本含水电机组的启停成本和弃风弃光惩罚、向配网购电的成本、调度周期内储能充放电的折旧成本。写成标准形式% 目标函数min Cost Cost sum(sum( ... % 对所有时段和场景求和 c_w * P_w_scen ... % 风电出力成本(弃风惩罚计入) c_pv * P_pv_scen ... % 光伏出力成本 c_h * (P_h_on P_h_off) ... % 水电运行成本(含启停) c_grid * P_buy_scen ... % 从配网购电成本 c_es * (P_ch P_dis) ... % 储能充放电折旧 ));注意这里P_w_scen、P_pv_scen是场景依赖的决策变量意味着模型是一个两阶段随机优化第一阶段决定VPP与DSO的交换功率基线非场景依赖第二阶段根据各场景的实际情况调整内部机组出力。c_w和c_pv设为很小的正值而不是零是为了让优化器在“弃掉可再生能源”和“以极小成本消纳”之间选择后者否则会出现退化情况。2.2 风光水的出力模型差异三种随机性不能套同一个抽样函数风电、光伏、水电的随机特性完全不同这是场景生成之前必须先看清的前提。风电出力跟风速的三次方近似成正比切入风速3m/s、额定风速12m/s、切出风速25m/s这一段曲线是分段函数概率分布上呈偏态。光伏出力主要受光照强度影响近似Beta分布且存在明显的白天时段约束——深夜光伏出力恒为0如果按全时段统一抽样会浪费大量场景样本在“不可能区间”。水电出力则受来水量的季节性和日调节约束存在最小出力和最大出力的硬边界且水库库容带来跨时段耦合。在MATLAB里分别构造这三种随机源时推荐的做法是这样% 风速场景Weibull分布抽样k2.1, c8.5为典型参数 wind_speed wblrnd(8.5, 2.1, n_scen, n_horizon); % 光照场景Beta分布抽样alpha2, beta5仅对日出时段 solar_irr betarnd(2, 5, n_scen, n_horizon); solar_irr(:, ~daylight_mask) 0; % 径流场景正态分布加约束裁剪 inflow normrnd(mu_inflow, sigma_inflow, n_scen, n_horizon); inflow max(inflow, inflow_min); inflow min(inflow, inflow_max);然后风电出力按功率曲线从speed映射到P_w光伏出力乘以装机容量和效率系数。这三组场景在进入模型前要各自归一化到0-1区间因为后续场景削减计算概率距离时不同量纲之间的欧氏距离会被风电的千瓦级数值淹没。2.3 配电网安全约束和VPP内部约束的衔接方式协调调度和VPP独立调度最大的区别就在这里VPP内部优化只需要管功率平衡和设备约束但协调调度还要把配电网潮流考虑进来。在复现时如果直接用完整DistFlow方程模型会变成非凸的求解速度断崖式下降。常见做法是线性化——忽略网损或把网损当作固定值处理用灵敏度系数建立节点电压与注入功率的线性关系。我一般这样组织约束矩阵% 配电网电压约束V_min V0 K * P_inj V_max % K为电压灵敏度矩阵由配电网潮流在一次基准点线性化得到 constraints [constraints, ... V_min * ones(n_bus, n_horizon) ... V_base K * P_inj V_max * ones(n_bus, n_horizon)]; % VPP内部功率平衡每个场景每个时段 constraints [constraints, ... P_w_scen P_pv_scen P_h_scen P_dis - P_ch P_buy ... P_load_scen];P_load_scen也是场景依赖的——负荷预测误差同样存在随机性。有的复现版本只对风光出力做场景负荷取期望值这会让调度结果在负荷波动大的日子严重偏乐观。更稳妥的做法是把负荷也纳入场景生成只不过负荷的相关系数矩阵跟风光不同要用历史负荷数据估算协方差矩阵而不是套用风速的Weibull分布。3. 场景生成与初始场景集构造蒙特卡洛抽样和关联关系处理3.1 独立抽样与相关性处理的差距如果对风、光、水三种出力分别独立抽样看起来省事但生产出来的场景集是“假”的——同一地区风电和光伏出力往往存在负相关多云天气光伏低但风速可能高径流与降雨相关而降雨又跟云量相关这导致独立抽样生成的风光水联合场景会出现在现实中几乎不会出现的组合比如“大风强光照丰水”同时发生。拿这种场景集去做调度优化结果会偏保守因为优化器为了应对几乎不存在的极端组合会预留过多备用容量。处理相关性有两种常用路线。第一种是Cholesky分解法适用于多维正态分布的情形。先对风光水出力历史数据做正态变换可以用逆CDF变换把非正态边缘分布映射到标准正态空间在正态空间里估计相关系数矩阵R然后对R做Cholesky分解得到下三角矩阵L最后把独立标准正态样本乘以L即可得到带相关性的样本。第二种是Copula法可以保留各自边缘分布的非正态特性更精细但计算量也更大。对于复现场景削减部分的工程目的Cholesky法足够。3.2 基于拉丁超立方抽样的初始场景生成普通蒙特卡洛抽样的缺点是样本点可能会扎堆——比如抽1000个风速样本极端高风速区间的样本可能只有个位数导致尾部分布刻画不充分。拉丁超立方抽样LHS把每个维度的累积概率区间等分成N个小区间在每个小区间内各抽一个样本保证覆盖均匀。MATLAB内置的lhsdesign可以直接用但它的默认设置是按[0,1]均匀分布设计的要用于风速和光照抽样需要配合逆CDF变换。% 生成500个初始场景24小时 n_scen 500; n_horizon 24; % LHS在[0,1]均匀采样 u lhsdesign(n_scen, n_horizon); % 转成风速样本u - Weibull逆CDF wind_scen wblinv(u, 8.5, 2.1); % 同理光伏光照用Beta逆CDF solar_scen betainv(u, 2, 5); solar_scen(:, ~daylight_mask) 0; % 夜间置零 % 相关性处理把独立LHS样本映射到相关空间 R [1, rho_ws, rho_wh; rho_ws, 1, rho_sh; rho_wh, rho_sh, 1]; % 风光水相关系数矩阵 L chol(R, lower); corr_samples L * [wind_scen(:); solar_scen(:); hydro_scen(:)]; % 重新还原为场景矩阵 wind_corr reshape(corr_samples(1,:), n_scen, n_horizon);wblinv和betainv是MATLAB统计工具箱里的逆CDF函数它们接受[0,1]区间输入并输出服从指定分布的随机数。chol(R,lower)得到的是下三角矩阵左乘到独立样本矩阵上就能把独立样本的相关系数矩阵从单位阵变换成目标矩阵R。有一点容易踩坑LHS抽样后经过Cholesky变换样本的均匀分层特性会被破坏一部分这是统计上不可避免的。所以在样本量选择上要留余量——如果削减后目标场景数是100初始场景建议生成500以上的LHS样本否则变换后尾部可能还是不够平滑。3.3 从连续分布到离散时间序列时序耦合怎么保风光的时序特性除了日内分布外还有小时级的相关性——15点的风速跟14点的风速高度相关如果每个时段独立抽样生成的风速曲线会剧烈跳变完全不像真实风速的平滑过程。这个问题的标准解法是使用一阶自回归模型AR(1)来生成时序样本。% AR(1)模型x(t) phi * x(t-1) eps(t) phi_w 0.85; % 风电时序自相关系数(典型值0.7-0.9) sigma_w sqrt(1 - phi_w^2); u_ar zeros(n_scen, n_horizon); u_ar(:,1) randn(n_scen, 1); % 初始时段标准正态 for t 2:n_horizon u_ar(:,t) phi_w * u_ar(:,t-1) sigma_w * randn(n_scen, 1); end % 再通过正态CDF映射到均匀分布最后经逆CDF转成风速 wind_scen wblinv(normcdf(u_ar), 8.5, 2.1);这里normcdf把AR(1)生成的时序正态样本映射回[0,1]区间再交给wblinv变成风速逻辑链是从“相关正态空间”转回“物理量空间”。phi_w越接近1风速曲线的平滑性越强但过度平滑会低估风速骤变带来的调峰压力。实际工程中可以用历史风速数据拟合自相关系数而不是拍脑袋取0.85。4. 场景削减的MATLAB实现快速前向选择与聚类法逐行拆解4.1 为什么不能直接对场景取平均——期望场景的致命缺陷很多初学者试图用“把500个场景求平均得到一条期望曲线”来减少计算量结果调度结果的可靠性惨不忍睹。原因很直白风电、光伏的出力曲线是高度非线性的风机功率曲线有死区和饱和区期望风速对应的功率并不等于功率的期望。举个数值例子两种风速场景分别是5m/s和13m/s平均风速9m/s对应的出力是500kW但5m/s对应的出力可能只有100kW13m/s已经达到额定1000kW两种场景的真实平均出力是550kW跟“用平均风速算”差了50kW。场景削减要保留的是“分布的形态”而不是把分布压缩成一个点。4.2 快速前向选择Fast Forward Selection的算法原理与MATLAB代码快速前向选择是目前随机规划中应用最广泛的场景削减算法。它的思路是从N个初始场景中逐步挑出要被删除的场景每删除一个场景就把它的概率加到离它最近的保留场景上直到剩余场景数达到预设值。每轮选择“删除哪个场景”依据的是该场景与最近场景之间的概率距离——距离越小说明它越能被其他场景代表优先删掉。function [scen_reduced, prob_reduced] fast_forward(scen, prob, n_keep) % 快速前向选择场景削减 % 输入: scen - N x T 场景矩阵(N个场景, T个时段) % prob - N x 1 各场景初始概率 % n_keep- 目标保留场景数 % 输出: scen_reduced - n_keep x T 削减后场景 % prob_reduced - n_keep x 1 削减后概率 N size(scen, 1); J (1:N); % 当前保留场景索引 p prob(:); % 当前概率 while length(J) n_keep % 计算保留场景两两之间的欧氏距离矩阵 D pdist2(scen(J,:), scen(J,:)); D D diag(inf(length(J),1)); % 自身距离设为inf % 每个场景到最近其他场景的距离 d_min min(D, [], 2); % 找距离最小的场景对中概率质量较小的那个删除 [~, idx_min] min(d_min); del_idx idx_min; % 被删除场景在J中的位置 % 找del_idx最近的其他场景把概率转移过去 [~, nearest] min(D(del_idx, :)); p(J(nearest)) p(J(nearest)) p(J(del_idx)); % 从保留集中移除 J(del_idx) []; end scen_reduced scen(J,:); prob_reduced p(J); endpdist2是MATLAB统计工具箱的距离矩阵计算函数默认欧氏距离。D diag(inf(...))的技巧是让对角线上的自距离不参与最小值比较这样每个场景的“最近场景”一定是另一个场景。当删除场景的概率转移到其最近邻后被删场景从索引数组J中移除循环继续直到保留数量达标。这段代码的复杂度是O(n³)——每轮删除都要重算距离矩阵。当初始场景数超过2000时运行时间会显著拉长。优化方案有两个一是只在被删除场景附近做局部距离更新而非全量重算代码复杂度高很多但速度快一个量级二是先用聚类把初始场景缩到几百个再跑快速前向选择。4.3 基于K-means聚类的场景削减路线与评价指标K-means聚类做场景削减的思路完全不同直接把初始场景当作N个数据点让K-means把它们分成n_keep个簇然后用每个簇的质心centroid作为削减后的典型场景把簇内所有场景的概率之和作为该质心的概率。实现简单计算速度快在初始场景数大5000时优势明显。% 用K-means做场景削减 n_keep 20; [cluster_idx, centroid] kmeans(scen_all, n_keep, ... Distance, sqeuclidean, ... Replicates, 5, ... % 多次重复找最优 MaxIter, 500); % 计算各簇概率 prob_red zeros(n_keep, 1); for k 1:n_keep prob_red(k) sum(prob_all(cluster_idx k)); endReplicates参数设为5或10很关键——K-means对初始中心敏感单次运行可能陷入局部最优多次重复取最小损失函数的那次才是稳的。Distance选sqeuclidean的原因是与K-means的原始目标函数簇内平方误差和一致选用cosine或cityblock会得到形状不同的簇但用途上对时序数据没有特别优势除非场景之间存在明显的幅值差异需要消除。聚类法和快速前向选择的核心差别在于K-means生成的质心是“虚构场景”——可能任何真实场景都不等于质心本身快速前向选择保留的是“真实历史场景”——每个保留场景都是初始集中真实存在的一条曲线。对于配电网安全校验这种需要“最坏情况分析”的场景虚构质心会平滑掉极端情况保留真实场景更稳妥。4.4 场景削减的评价指标削减前后的概率分布距离削减完不能只盯着剩余场景数看必须量化“削掉了多少信息”。常用的指标是削减前后的场景集与原始场景集之间的Kantorovich距离简称KD它衡量的是两个概率分布之间的最小传输代价也就是把原始分布“搬运”成削减后分布需要付出的最小成本。% 计算削减前后的Kantorovich距离 % 先计算每个原始场景到最近保留场景的距离 D_orig pdist2(scen_orig, scen_reduced); d_min_orig min(D_orig, [], 2); % KD 每个原始场景的最近距离乘以其概率之和 KD sum(prob_orig .* d_min_orig); % 计算各时段平均绝对误差作为辅助指标 MAE mean(abs(scen_orig * prob_orig - scen_reduced * prob_reduced));KD值越小代表削减质量越高但它没有绝对的好坏标准需要跟削减数量挂钩看变化趋势随着保留场景数从5增加到50KD应该快速下降然后趋于平缓形成肘部曲线。如果肘部出现在保留场景远大于预期值的位置说明初始场景间的差异过大需要回去检查场景生成部分是否把相关性空间映射搞错了。4.5 场景削减的参数设置表与常见误区参数快速前向选择K-means聚类建议初始场景数500-20002000-10000小于500时削减优势不明显目标场景数10-5010-100目标数取调度模型的整数变量规模的1/10左右距离度量欧氏距离欧氏距离量纲不一致时先归一化停止条件保留数达标簇数达标迭代收敛快速前向无收敛性问题最大瓶颈距离矩阵O(n²)内存初始中心敏感性大规模用K-means多次重复注意场景削减的目标场景数并非越少越好。目标数过小时削减后的调度方案会低估系统运行成本——因为不确定性被过度压缩。目标是“用尽可能少的场景保持概率特征”而不是“用最少的场景完成任务”。一般以削减后的期望运行成本与削减前成本的偏差不超过5%为可接受阈值。5. 削减后的协调调度求解与结果校验一套可复现的闭合流程5.1 YALMIP建模时的场景循环处理技巧场景削减完成后原来的机会约束或期望值目标函数就需要改写成对削减后场景逐一遍历的确定性等价形式。在YALMIP里一个常见错误是试图用三维决策变量一次塞进所有场景约束导致YALMIP解析时把每个场景当成独立优化问题处理变量数爆炸。标准做法是用循环给每个场景加约束让求解器看到统一的稀疏结构% 削减后场景数 n_keep P_w sdpvar(n_keep, n_horizon, full); % 每个场景的风电出力 P_h sdpvar(n_keep, n_horizon, full); P_buy sdpvar(n_keep, n_horizon, full); P_ex sdpvar(1, n_horizon, full); % 一阶段交换功率(非场景依赖) constraints []; for s 1:n_keep constraints [constraints, ... P_w(s,:) P_pv(s,:) P_h(s,:) P_ex P_load(s,:), ... P_w(s,:) P_w_max * scen_w(s,:), ... % 场景数据作为系数 P_h(s,:) P_h_max, ... P_buy(s,:) 0]; end % 目标函数按场景概率加权 objective 0; for s 1:n_keep objective objective prob_red(s) * ... (c_w * sum(P_w(s,:)) c_buy * sum(P_buy(s,:)) ...); end optimize(constraints, objective, sdpsettings(solver, gurobi));sdpsettings里的solver字段要显式指定。MATLAB环境下默认求解器可能是sedumi或sdpt3它们对混合整数规划无能为力模型里只要带一个整数变量就会报错。指定gurobi或cplex后求解速度通常有10倍以上的差距。5.2 削减效果的验证维度成本偏差、风险价值与计算时间调度结果出来之后必须回答“削减后的方案跟削减前比差多少”这个问题。三个维度缺一不可第一是期望运行成本的偏差把削减后调度方案代回全部初始场景计算真实期望成本对比削减前模型的目标函数值第二是5%分位数的风险情况——最差场景下的运行成本是否被低估第三是求解时间这是场景削减的直接收益。% 回代验证把削减后求得的P_ex用于所有初始场景 cost_actual 0; for s 1:n_scen_orig cost_actual cost_actual prob_orig(s) * ... compute_cost(P_ex, P_w_opt(s,:), P_load_orig(s,:), price); end % 计算削减前后成本偏差 cost_reduced_model value(objective); % 削减后模型的目标函数值 deviation abs(cost_actual - cost_reduced_model) / cost_actual * 100; fprintf(成本偏差: %.2f%%\n, deviation);如果deviation超过8%说明削减后的场景没有忠实地代表原始场景分布。处理方式不是继续加目标场景数而是先检查KD肘部曲线——如果KD在目标场景数附近已经明显收敛偏差却依然很大就要考虑是场景生成阶段的时序相关性没有建模对导致削减算法在“结构错误”的场景集上做无用功。5.3 实战技巧用“迭代剪枝”快速定位最优保留场景数找最优保留场景数时不要一个个试。我一般用二分迭代先试50个场景计算期望成本再试10个场景计算期望成本如果两者偏差小于3%砍半到5个再试如果10个场景偏差超过5%就往25、35方向加密。这种方式通常只需要4-6次完整求解就能锁定最优范围比从1到50逐个遍历节省八成以上的时间。对于“含风光水的虚拟电厂与配电公司协调调度模型”这个框架来说场景削减的真正目的不是让模型变快而是让模型在面对真实分布时算得准。快速前向选择保真度高K-means计算快两者结合——先用K-means粗减到200个再用快速前向精减到20个——是我在复现类似项目时最常见的方案。初始场景500个以上、削减后20个左右、KD肘部收敛、回代成本偏差控制5%以内这条路径跑通整个协调调度模型的可信度就有了保障。本文还有配套的精品资源点击获取 含风光水的虚拟电厂与配电网协调调度中的场景削减方法解析
简介MATLAB复现的虚拟电厂与配电公司协调调度模型中的场景削减模块面向电力系统优化调度、新能源不确定性建模方向的研究者与工程师尤其适合处理含高比例新能源的复杂不确定性场景。代码聚焦风电、光伏及电价场景模拟由一组确定性方案生成50个光伏场景再采用基于概率距离的快速削减算法将场景压缩至5个运行后直接输出削减前后的场景及对应概率既降低大规模场景带来的内存与计算负担又保留主要概率特征可直接嵌入含风光水的协调调度、备用优化及配电公司决策等模型。资源共5个文件以4个m脚本为主体分别对应风电、光伏、电价等场景生成与削减流程另有1个txt参考文献说明压缩包仅5KB轻量易读、结构清晰便于按需修改参数后复现或二次开发。该资源已有761人学习适合具备一定MATLAB基础、希望快速落地场景削减算法或开展虚拟电厂协调调度研究的研究生与工程技术人员。1. 含风光水的虚拟电厂与配电公司协调调度为什么偏偏要谈场景削减“场景削减”四个字看着像数据预处理的小步骤实际是含风光水的虚拟电厂VPP与配电公司DSO协调调度能否落地的胜负手。风电、光伏的出力曲线本质上是从随机分布里采出来的一批“可能明天长这样”的样本水电又受来水和库容双重约束三者叠加之后优化模型的约束矩阵会被几十上百条时序场景撑到难以直接求解。越想把随机性刻画得细场景数越多混合整数规划的求解时间就越失控——这就是为什么任何讲VPP与DSO协调调度的复现项目都要单拎出“场景削减部分”作为独立模块来谈。本文按照“先立模型、再造场景、再谈削减、最后校验”的顺序拆解。读者如果手上有一份还没跑通的MATLAB复现工程或者正在决定削减算法用K-means还是快速前向选择应该能从里面对上号。至于为什么不让场景削减走“直接抽样取均值”这种简单路第三、四章会给出可量化的答案。先说明一个基本立场下面所有代码按MATLAB R2022b及以上版本组织优化求解部分依赖YALMIP配cplex或gurobi场景生成与削减仅用MATLAB基础工具箱不需要额外安装包。2. 协调调度模型的目标函数与约束体系先把边界定清楚2.1 VPP与配电公司的交易框架谁向谁买电按什么价格买虚拟电厂聚合了风、光、水三种分布式资源外加储能和可调负荷对外表现为一个可以统一调度的发电主体。配电公司在这个框架里扮演两个角色一是网络运营商保证配电网节点电压和支路潮流不越限二是购电方从VPP收购电力以满足负荷需求。协调调度模型要做的事就是在给定的调度周期内通常取24小时分辨率1小时或15分钟算出VPP内部各机组的最优出力计划和VPP与DSO之间的交换功率使系统总运行成本最小同时保证网络安全约束。目标函数可以分解成三块VPP内部发电成本含水电机组的启停成本和弃风弃光惩罚、向配网购电的成本、调度周期内储能充放电的折旧成本。写成标准形式% 目标函数min Cost Cost sum(sum( ... % 对所有时段和场景求和 c_w * P_w_scen ... % 风电出力成本(弃风惩罚计入) c_pv * P_pv_scen ... % 光伏出力成本 c_h * (P_h_on P_h_off) ... % 水电运行成本(含启停) c_grid * P_buy_scen ... % 从配网购电成本 c_es * (P_ch P_dis) ... % 储能充放电折旧 ));注意这里P_w_scen、P_pv_scen是场景依赖的决策变量意味着模型是一个两阶段随机优化第一阶段决定VPP与DSO的交换功率基线非场景依赖第二阶段根据各场景的实际情况调整内部机组出力。c_w和c_pv设为很小的正值而不是零是为了让优化器在“弃掉可再生能源”和“以极小成本消纳”之间选择后者否则会出现退化情况。2.2 风光水的出力模型差异三种随机性不能套同一个抽样函数风电、光伏、水电的随机特性完全不同这是场景生成之前必须先看清的前提。风电出力跟风速的三次方近似成正比切入风速3m/s、额定风速12m/s、切出风速25m/s这一段曲线是分段函数概率分布上呈偏态。光伏出力主要受光照强度影响近似Beta分布且存在明显的白天时段约束——深夜光伏出力恒为0如果按全时段统一抽样会浪费大量场景样本在“不可能区间”。水电出力则受来水量的季节性和日调节约束存在最小出力和最大出力的硬边界且水库库容带来跨时段耦合。在MATLAB里分别构造这三种随机源时推荐的做法是这样% 风速场景Weibull分布抽样k2.1, c8.5为典型参数 wind_speed wblrnd(8.5, 2.1, n_scen, n_horizon); % 光照场景Beta分布抽样alpha2, beta5仅对日出时段 solar_irr betarnd(2, 5, n_scen, n_horizon); solar_irr(:, ~daylight_mask) 0; % 径流场景正态分布加约束裁剪 inflow normrnd(mu_inflow, sigma_inflow, n_scen, n_horizon); inflow max(inflow, inflow_min); inflow min(inflow, inflow_max);然后风电出力按功率曲线从speed映射到P_w光伏出力乘以装机容量和效率系数。这三组场景在进入模型前要各自归一化到0-1区间因为后续场景削减计算概率距离时不同量纲之间的欧氏距离会被风电的千瓦级数值淹没。2.3 配电网安全约束和VPP内部约束的衔接方式协调调度和VPP独立调度最大的区别就在这里VPP内部优化只需要管功率平衡和设备约束但协调调度还要把配电网潮流考虑进来。在复现时如果直接用完整DistFlow方程模型会变成非凸的求解速度断崖式下降。常见做法是线性化——忽略网损或把网损当作固定值处理用灵敏度系数建立节点电压与注入功率的线性关系。我一般这样组织约束矩阵% 配电网电压约束V_min V0 K * P_inj V_max % K为电压灵敏度矩阵由配电网潮流在一次基准点线性化得到 constraints [constraints, ... V_min * ones(n_bus, n_horizon) ... V_base K * P_inj V_max * ones(n_bus, n_horizon)]; % VPP内部功率平衡每个场景每个时段 constraints [constraints, ... P_w_scen P_pv_scen P_h_scen P_dis - P_ch P_buy ... P_load_scen];P_load_scen也是场景依赖的——负荷预测误差同样存在随机性。有的复现版本只对风光出力做场景负荷取期望值这会让调度结果在负荷波动大的日子严重偏乐观。更稳妥的做法是把负荷也纳入场景生成只不过负荷的相关系数矩阵跟风光不同要用历史负荷数据估算协方差矩阵而不是套用风速的Weibull分布。3. 场景生成与初始场景集构造蒙特卡洛抽样和关联关系处理3.1 独立抽样与相关性处理的差距如果对风、光、水三种出力分别独立抽样看起来省事但生产出来的场景集是“假”的——同一地区风电和光伏出力往往存在负相关多云天气光伏低但风速可能高径流与降雨相关而降雨又跟云量相关这导致独立抽样生成的风光水联合场景会出现在现实中几乎不会出现的组合比如“大风强光照丰水”同时发生。拿这种场景集去做调度优化结果会偏保守因为优化器为了应对几乎不存在的极端组合会预留过多备用容量。处理相关性有两种常用路线。第一种是Cholesky分解法适用于多维正态分布的情形。先对风光水出力历史数据做正态变换可以用逆CDF变换把非正态边缘分布映射到标准正态空间在正态空间里估计相关系数矩阵R然后对R做Cholesky分解得到下三角矩阵L最后把独立标准正态样本乘以L即可得到带相关性的样本。第二种是Copula法可以保留各自边缘分布的非正态特性更精细但计算量也更大。对于复现场景削减部分的工程目的Cholesky法足够。3.2 基于拉丁超立方抽样的初始场景生成普通蒙特卡洛抽样的缺点是样本点可能会扎堆——比如抽1000个风速样本极端高风速区间的样本可能只有个位数导致尾部分布刻画不充分。拉丁超立方抽样LHS把每个维度的累积概率区间等分成N个小区间在每个小区间内各抽一个样本保证覆盖均匀。MATLAB内置的lhsdesign可以直接用但它的默认设置是按[0,1]均匀分布设计的要用于风速和光照抽样需要配合逆CDF变换。% 生成500个初始场景24小时 n_scen 500; n_horizon 24; % LHS在[0,1]均匀采样 u lhsdesign(n_scen, n_horizon); % 转成风速样本u - Weibull逆CDF wind_scen wblinv(u, 8.5, 2.1); % 同理光伏光照用Beta逆CDF solar_scen betainv(u, 2, 5); solar_scen(:, ~daylight_mask) 0; % 夜间置零 % 相关性处理把独立LHS样本映射到相关空间 R [1, rho_ws, rho_wh; rho_ws, 1, rho_sh; rho_wh, rho_sh, 1]; % 风光水相关系数矩阵 L chol(R, lower); corr_samples L * [wind_scen(:); solar_scen(:); hydro_scen(:)]; % 重新还原为场景矩阵 wind_corr reshape(corr_samples(1,:), n_scen, n_horizon);wblinv和betainv是MATLAB统计工具箱里的逆CDF函数它们接受[0,1]区间输入并输出服从指定分布的随机数。chol(R,lower)得到的是下三角矩阵左乘到独立样本矩阵上就能把独立样本的相关系数矩阵从单位阵变换成目标矩阵R。有一点容易踩坑LHS抽样后经过Cholesky变换样本的均匀分层特性会被破坏一部分这是统计上不可避免的。所以在样本量选择上要留余量——如果削减后目标场景数是100初始场景建议生成500以上的LHS样本否则变换后尾部可能还是不够平滑。3.3 从连续分布到离散时间序列时序耦合怎么保风光的时序特性除了日内分布外还有小时级的相关性——15点的风速跟14点的风速高度相关如果每个时段独立抽样生成的风速曲线会剧烈跳变完全不像真实风速的平滑过程。这个问题的标准解法是使用一阶自回归模型AR(1)来生成时序样本。% AR(1)模型x(t) phi * x(t-1) eps(t) phi_w 0.85; % 风电时序自相关系数(典型值0.7-0.9) sigma_w sqrt(1 - phi_w^2); u_ar zeros(n_scen, n_horizon); u_ar(:,1) randn(n_scen, 1); % 初始时段标准正态 for t 2:n_horizon u_ar(:,t) phi_w * u_ar(:,t-1) sigma_w * randn(n_scen, 1); end % 再通过正态CDF映射到均匀分布最后经逆CDF转成风速 wind_scen wblinv(normcdf(u_ar), 8.5, 2.1);这里normcdf把AR(1)生成的时序正态样本映射回[0,1]区间再交给wblinv变成风速逻辑链是从“相关正态空间”转回“物理量空间”。phi_w越接近1风速曲线的平滑性越强但过度平滑会低估风速骤变带来的调峰压力。实际工程中可以用历史风速数据拟合自相关系数而不是拍脑袋取0.85。4. 场景削减的MATLAB实现快速前向选择与聚类法逐行拆解4.1 为什么不能直接对场景取平均——期望场景的致命缺陷很多初学者试图用“把500个场景求平均得到一条期望曲线”来减少计算量结果调度结果的可靠性惨不忍睹。原因很直白风电、光伏的出力曲线是高度非线性的风机功率曲线有死区和饱和区期望风速对应的功率并不等于功率的期望。举个数值例子两种风速场景分别是5m/s和13m/s平均风速9m/s对应的出力是500kW但5m/s对应的出力可能只有100kW13m/s已经达到额定1000kW两种场景的真实平均出力是550kW跟“用平均风速算”差了50kW。场景削减要保留的是“分布的形态”而不是把分布压缩成一个点。4.2 快速前向选择Fast Forward Selection的算法原理与MATLAB代码快速前向选择是目前随机规划中应用最广泛的场景削减算法。它的思路是从N个初始场景中逐步挑出要被删除的场景每删除一个场景就把它的概率加到离它最近的保留场景上直到剩余场景数达到预设值。每轮选择“删除哪个场景”依据的是该场景与最近场景之间的概率距离——距离越小说明它越能被其他场景代表优先删掉。function [scen_reduced, prob_reduced] fast_forward(scen, prob, n_keep) % 快速前向选择场景削减 % 输入: scen - N x T 场景矩阵(N个场景, T个时段) % prob - N x 1 各场景初始概率 % n_keep- 目标保留场景数 % 输出: scen_reduced - n_keep x T 削减后场景 % prob_reduced - n_keep x 1 削减后概率 N size(scen, 1); J (1:N); % 当前保留场景索引 p prob(:); % 当前概率 while length(J) n_keep % 计算保留场景两两之间的欧氏距离矩阵 D pdist2(scen(J,:), scen(J,:)); D D diag(inf(length(J),1)); % 自身距离设为inf % 每个场景到最近其他场景的距离 d_min min(D, [], 2); % 找距离最小的场景对中概率质量较小的那个删除 [~, idx_min] min(d_min); del_idx idx_min; % 被删除场景在J中的位置 % 找del_idx最近的其他场景把概率转移过去 [~, nearest] min(D(del_idx, :)); p(J(nearest)) p(J(nearest)) p(J(del_idx)); % 从保留集中移除 J(del_idx) []; end scen_reduced scen(J,:); prob_reduced p(J); endpdist2是MATLAB统计工具箱的距离矩阵计算函数默认欧氏距离。D diag(inf(...))的技巧是让对角线上的自距离不参与最小值比较这样每个场景的“最近场景”一定是另一个场景。当删除场景的概率转移到其最近邻后被删场景从索引数组J中移除循环继续直到保留数量达标。这段代码的复杂度是O(n³)——每轮删除都要重算距离矩阵。当初始场景数超过2000时运行时间会显著拉长。优化方案有两个一是只在被删除场景附近做局部距离更新而非全量重算代码复杂度高很多但速度快一个量级二是先用聚类把初始场景缩到几百个再跑快速前向选择。4.3 基于K-means聚类的场景削减路线与评价指标K-means聚类做场景削减的思路完全不同直接把初始场景当作N个数据点让K-means把它们分成n_keep个簇然后用每个簇的质心centroid作为削减后的典型场景把簇内所有场景的概率之和作为该质心的概率。实现简单计算速度快在初始场景数大5000时优势明显。% 用K-means做场景削减 n_keep 20; [cluster_idx, centroid] kmeans(scen_all, n_keep, ... Distance, sqeuclidean, ... Replicates, 5, ... % 多次重复找最优 MaxIter, 500); % 计算各簇概率 prob_red zeros(n_keep, 1); for k 1:n_keep prob_red(k) sum(prob_all(cluster_idx k)); endReplicates参数设为5或10很关键——K-means对初始中心敏感单次运行可能陷入局部最优多次重复取最小损失函数的那次才是稳的。Distance选sqeuclidean的原因是与K-means的原始目标函数簇内平方误差和一致选用cosine或cityblock会得到形状不同的簇但用途上对时序数据没有特别优势除非场景之间存在明显的幅值差异需要消除。聚类法和快速前向选择的核心差别在于K-means生成的质心是“虚构场景”——可能任何真实场景都不等于质心本身快速前向选择保留的是“真实历史场景”——每个保留场景都是初始集中真实存在的一条曲线。对于配电网安全校验这种需要“最坏情况分析”的场景虚构质心会平滑掉极端情况保留真实场景更稳妥。4.4 场景削减的评价指标削减前后的概率分布距离削减完不能只盯着剩余场景数看必须量化“削掉了多少信息”。常用的指标是削减前后的场景集与原始场景集之间的Kantorovich距离简称KD它衡量的是两个概率分布之间的最小传输代价也就是把原始分布“搬运”成削减后分布需要付出的最小成本。% 计算削减前后的Kantorovich距离 % 先计算每个原始场景到最近保留场景的距离 D_orig pdist2(scen_orig, scen_reduced); d_min_orig min(D_orig, [], 2); % KD 每个原始场景的最近距离乘以其概率之和 KD sum(prob_orig .* d_min_orig); % 计算各时段平均绝对误差作为辅助指标 MAE mean(abs(scen_orig * prob_orig - scen_reduced * prob_reduced));KD值越小代表削减质量越高但它没有绝对的好坏标准需要跟削减数量挂钩看变化趋势随着保留场景数从5增加到50KD应该快速下降然后趋于平缓形成肘部曲线。如果肘部出现在保留场景远大于预期值的位置说明初始场景间的差异过大需要回去检查场景生成部分是否把相关性空间映射搞错了。4.5 场景削减的参数设置表与常见误区参数快速前向选择K-means聚类建议初始场景数500-20002000-10000小于500时削减优势不明显目标场景数10-5010-100目标数取调度模型的整数变量规模的1/10左右距离度量欧氏距离欧氏距离量纲不一致时先归一化停止条件保留数达标簇数达标迭代收敛快速前向无收敛性问题最大瓶颈距离矩阵O(n²)内存初始中心敏感性大规模用K-means多次重复注意场景削减的目标场景数并非越少越好。目标数过小时削减后的调度方案会低估系统运行成本——因为不确定性被过度压缩。目标是“用尽可能少的场景保持概率特征”而不是“用最少的场景完成任务”。一般以削减后的期望运行成本与削减前成本的偏差不超过5%为可接受阈值。5. 削减后的协调调度求解与结果校验一套可复现的闭合流程5.1 YALMIP建模时的场景循环处理技巧场景削减完成后原来的机会约束或期望值目标函数就需要改写成对削减后场景逐一遍历的确定性等价形式。在YALMIP里一个常见错误是试图用三维决策变量一次塞进所有场景约束导致YALMIP解析时把每个场景当成独立优化问题处理变量数爆炸。标准做法是用循环给每个场景加约束让求解器看到统一的稀疏结构% 削减后场景数 n_keep P_w sdpvar(n_keep, n_horizon, full); % 每个场景的风电出力 P_h sdpvar(n_keep, n_horizon, full); P_buy sdpvar(n_keep, n_horizon, full); P_ex sdpvar(1, n_horizon, full); % 一阶段交换功率(非场景依赖) constraints []; for s 1:n_keep constraints [constraints, ... P_w(s,:) P_pv(s,:) P_h(s,:) P_ex P_load(s,:), ... P_w(s,:) P_w_max * scen_w(s,:), ... % 场景数据作为系数 P_h(s,:) P_h_max, ... P_buy(s,:) 0]; end % 目标函数按场景概率加权 objective 0; for s 1:n_keep objective objective prob_red(s) * ... (c_w * sum(P_w(s,:)) c_buy * sum(P_buy(s,:)) ...); end optimize(constraints, objective, sdpsettings(solver, gurobi));sdpsettings里的solver字段要显式指定。MATLAB环境下默认求解器可能是sedumi或sdpt3它们对混合整数规划无能为力模型里只要带一个整数变量就会报错。指定gurobi或cplex后求解速度通常有10倍以上的差距。5.2 削减效果的验证维度成本偏差、风险价值与计算时间调度结果出来之后必须回答“削减后的方案跟削减前比差多少”这个问题。三个维度缺一不可第一是期望运行成本的偏差把削减后调度方案代回全部初始场景计算真实期望成本对比削减前模型的目标函数值第二是5%分位数的风险情况——最差场景下的运行成本是否被低估第三是求解时间这是场景削减的直接收益。% 回代验证把削减后求得的P_ex用于所有初始场景 cost_actual 0; for s 1:n_scen_orig cost_actual cost_actual prob_orig(s) * ... compute_cost(P_ex, P_w_opt(s,:), P_load_orig(s,:), price); end % 计算削减前后成本偏差 cost_reduced_model value(objective); % 削减后模型的目标函数值 deviation abs(cost_actual - cost_reduced_model) / cost_actual * 100; fprintf(成本偏差: %.2f%%\n, deviation);如果deviation超过8%说明削减后的场景没有忠实地代表原始场景分布。处理方式不是继续加目标场景数而是先检查KD肘部曲线——如果KD在目标场景数附近已经明显收敛偏差却依然很大就要考虑是场景生成阶段的时序相关性没有建模对导致削减算法在“结构错误”的场景集上做无用功。5.3 实战技巧用“迭代剪枝”快速定位最优保留场景数找最优保留场景数时不要一个个试。我一般用二分迭代先试50个场景计算期望成本再试10个场景计算期望成本如果两者偏差小于3%砍半到5个再试如果10个场景偏差超过5%就往25、35方向加密。这种方式通常只需要4-6次完整求解就能锁定最优范围比从1到50逐个遍历节省八成以上的时间。对于“含风光水的虚拟电厂与配电公司协调调度模型”这个框架来说场景削减的真正目的不是让模型变快而是让模型在面对真实分布时算得准。快速前向选择保真度高K-means计算快两者结合——先用K-means粗减到200个再用快速前向精减到20个——是我在复现类似项目时最常见的方案。初始场景500个以上、削减后20个左右、KD肘部收敛、回代成本偏差控制5%以内这条路径跑通整个协调调度模型的可信度就有了保障。本文还有配套的精品资源点击获取