Zonotope方法在虚拟电厂广域聚合调控中的鲁棒调度应用

简介面向智能电网虚拟电厂VPP场景这份资料复现了基于奇诺多面体Zonotope的分布式资源广域聚合调控方法适合具有电力系统背景和Python编程基础的研究人员与开发人员。内容涵盖空调负荷、储能设备、柴油发电机三类关键资源的可行域刻画并给出从半空间约束到Zonotope描述、Minkowski求和聚合以及线性规划调度等核心Python代码与逐段解释可支撑大型园区或多用户场景下的精细化调度与鲁棒性分析。压缩包仅含1个docx文档体积22KB便于快速查阅目前已有237人学习下载。通过具体时间段功率输出预测及展示图读者可直观理解资源聚合与调控流程并可将参数调整应用于实际项目值得说明的是预览中提到复现提供的是核心代码框架完整实现需结合实际数据进一步丰富适合作为算法复现与二次开发的基础。1. 用Zonotope做广域聚合调控先算清虚拟电厂“不确定性”这笔账做虚拟电厂的人都知道分布式资源一变多广域聚合调控就开始不好伺候光伏骤降、负荷毛刺、储能SOC标称和实际对不上每一版预测误差都在给调度中心添乱。Zonotope方法就是拿来把这笔账算清楚的——给每个分布式资源的预测误差画一个紧致凸集合几百个资源用Minkowski加法聚成一个整体再把鲁棒调度写成线性规划整套流程在Python里从建模到求解都能跑起来。这篇会把Zonotope的原理、资源级建模、聚合后的鲁棒调度代码和避坑细节一次讲完适合正在做虚拟电厂调度、聚合商报价或者鲁棒优化入门的人。2. 广域聚合为什么选Zonotope从区间盒子、椭球到凸多面体的取舍2.1 三种不确定性集合在聚合调控里的“手感和账本”先把候选工具摆出来区间盒子、椭球、Zonotope。它们在描述“预测误差到底可能落在哪”这件事上手感差别很大。区间盒子是默认起点。每个分布式资源给一个预测区间联合起来就是一个轴对齐的立方体。优点是直观缺点是它默认所有资源的误差互相独立、并且最坏情况会在同一个方向上同时发生。实际场景里区域光伏和本地负荷在时序上往往有互补关系一个涨的时候另一个多半在跌盒子完全看不见这种结构。更麻烦的是Minkowski和两个盒子就是把边长直接相加资源数量一上来可行域被撑得特别大调度结果会显得非常“笨”——备用容量永远留过头聚合商报出去的曲线没有竞争力。椭球能刻画相关性协方差矩阵天生就是干这个的。但椭球在Minkowski和之后不再是一个椭球工程上只能再做一层外层近似这个近似本身又会把形状撑大一圈线性映射虽然保持椭球形式但协方差被左右各乘一次之后膨胀得很厉害求解还要上二阶锥规划。对一线的调度系统来说能用线性规划解决的事尽量别上SOCP这是很现实的理由。Zonotope是凸多边形的推广定义式是Z {c G·x : ‖x‖∞ ≤ 1}其中c是中心向量G是生成器矩阵。它的核心优势是两条对Minkowski和封闭对线性映射封闭。前者意味着几百个资源的聚合就是生成器矩阵水平拼接后者意味着功率平衡、线路潮流这些线性变换不会破坏集合的数学形式。而且生成器的列数可以截断——想要更紧就多留几个方向想跑得快就少留几个谱系上有明确的操作旋钮。我一般做选型时会用下面这张表快速过一遍集合类型Minkowski和线性映射相关性刻画典型求解保守性控制区间盒子边长逐维相加逐维绝对值和不能LP只能整体缩放椭球需要外层近似协方差左右膨胀能SOCP靠缩放因子Zonotope生成器水平拼接左乘矩阵即可能生成器方向体现LP截断生成器列数2.2 Zonotope的支撑函数把“任意不确定性”变成一条线性不等式Zonotope在鲁棒调度里最值钱的不是它的几何形状而是它的支撑函数可以被解析地算出来。对任意方向向量a支撑函数定义为hZ(a) max aᵀzz∈Z对Zonotope来说结果很简单hZ(a) aᵀc ‖Gᵀa‖₁也就是说Zonotope在任意方向上的最大投影等于中心投影加上生成器矩阵在该方向上投影的1-范数。工程含义非常直接如果你有一条安全约束“aᵀ(p u) ≤ b”其中u落在零中心的Zonotope里那么这条约束对所有可能的u成立当且仅当aᵀp ‖Gᵀa‖₁ ≤ b。无穷多个不确定场景被压缩成了一个1-范数项这是Zonotope方法能在优化求解器里落地的根本原因。写代码验证这个公式时我一直用的方法是拿随机采样对拍防止自己实现错了方向还浑然不觉import numpy as np rng np.random.default_rng(42) # 一个 2 维、3 个生成器的 Zonotope故意取不对称形状 center np.array([0.0, 0.0]) G np.array([ [1.2, 0.5, 0.3], [0.4, 1.0, -0.6], ]) def support(direction): a np.asarray(direction, dtypefloat) # 支撑函数解析式a^T c sum_j |a^T g_j| return float(a center np.abs(a G).sum()) # 随机采样 50 万个 x用定义式生成 z 点再取方向投影的最大值 xs rng.uniform(-1, 1, size(500_000, G.shape[1])) z_points center xs G.T a np.array([1.0, 0.0]) # 目标方向 max_by_sampling float((z_points a).max()) print(解析支撑函数值:, support(a)) print(采样最大值:, max_by_sampling)这段代码里的support函数实现了解析式采样部分只是用来对拍——在真实的优化求解环境里你不可能采样五十万个点再去找最值支撑函数的价值就是把这个最值用一条闭式表达式算出来。参数方面采样数量50万是我的习惯足够多又不至于让演示脚本跑太久固定随机种子42保证每次对拍结果可复现。如果发现解析值和采样最大值对不上优先检查生成器矩阵的维度方向我习惯把G存成(n_dim, n_gen)也就是每列是一个生成器向量支撑函数里要转置后和方向向量a做矩阵乘。2.3 拿历史数据生成生成器矩阵PCA主方向加经验分位数残差盒Zonotope的形状不是拍脑袋定的而是从“预测值-实际值”的历史误差序列里挤出来的。我常用的构造流程分三步零中心化、取主方向、残差用分位数兜底。第一步对每条分布式资源的历史误差序列按维度减去中位数而不是均值。预测误差经常有偏比如光伏预测在阴天系统性低估用中位数能避免一个偏离很大的样本把中心拉走。第二步对零中心化后的误差矩阵做协方差分解取前k个特征向量和对应的特征值平方根作为主生成器。第三步把原始误差投影到主生成器张成的子空间上投影后的残差逐维取经验分位数比如90%分位数的半区间作为若干条对角生成器兜住那些主方向没解释掉的尾巴。这个流程对应的Python函数下面这段可以直接抄依赖就是numpydef build_zonotope_from_errors(errors, n_keep3, quantile0.9, alpha1.1): # errors: (n_samples, n_dim)每一行是一条历史预测误差样本 err errors - np.median(errors, axis0, keepdimsTrue) cov np.cov(err, rowvarFalse) w, V np.linalg.eigh(cov) idx np.argsort(w)[::-1][:n_keep] # 取特征值最大的 n_keep 个方向 V_k V[:, idx] proj (err V_k) V_k.T # 向主方向子空间投影 resid err - proj G_main V_k * np.sqrt(w[idx])[None, :] * alpha G_res np.diag(np.quantile(np.abs(resid), quantile, axis0)) center np.zeros(err.shape[1]) return center, np.hstack([G_main, G_res])这里三个参数需要根据资源特性调。n_keep是主方向数量一般取3到5光伏这种强时序相关的取3就够负荷形态复杂的可以加到5取多了生成器列数膨胀取少了残差盒子的半区间会很大整体形状反而松。quantile是残差盒子的分位数我一般取0.9到0.95它直接决定鲁棒约束的覆盖概率取太大保守性明显上升。alpha是校准系数作用是补偿有限样本带来的估计偏差取值在1.0到1.2之间工程上通常保留默认1.1等回测结果出来再微调。这个构造方式不假设误差服从高斯分布它只要求历史数据覆盖了实际会发生的情况。天气模式一变误差形态可能整体漂移所以这套生成器矩阵不是一劳永逸的——后面第6章会讲怎么定期更新。3. 把分布式资源“包”成Zonotope光伏、负荷、储能的建模与Python实现3.1 一个够用的Zonotope类中心、生成器矩阵与三个核心方法建模之前先把基础设施写好。一个最小的Zonotope类只需要四个东西中心c、生成器矩阵G、支撑函数、Minkowski加法。线性映射方法也顺手加上后面做网络约束时会用到。class Zonotope: def __init__(self, center, generators): # generators 形状约定为 (n_dim, n_gen)每列是一个生成器向量 self.c np.asarray(center, dtypefloat).reshape(-1) self.G np.asarray(generators, dtypefloat) self.n_dim self.G.shape[0] self.n_gen self.G.shape[1] def support(self, direction): # 方向向量 shape 必须与 self.c 一致 a np.asarray(direction, dtypefloat).reshape(-1) return float(a self.c np.abs(a self.G).sum()) def __add__(self, other): # Minkowski 和中心相加生成器水平拼接 return Zonotope(self.c other.c, np.hstack([self.G, other.G])) def linear_map(self, A): # 线性映射A 作用于集合中的每个点 return Zonotope(A self.c, A self.G)这个类里的__add__是整个广域聚合的核心两个Zonotope相加中心直接相加生成器矩阵水平拼接。注意不是求凸包凸包会丢掉内部结构、把集合撑大Minkowski和保留了每个资源独立的误差维度代价仅仅是生成器列数变多。linear_map对应功率平衡方程里的系数矩阵约束变换时用得上。有一点要提醒两个Zonotope相加时n_dim必须一致否则np.hstack会拼出形状错乱的矩阵报错信息看着会很懵。3.2 资源级建模光伏预测残差、负荷基线偏差、储能执行误差不同类型的分布式资源在聚合里的角色不一样。光伏和负荷是不可控的扰动量它们的预测误差要进Zonotope储能是可控量功率大小是调度变量不需要用Zonotope描述可调范围但可以考虑通信延迟或功率跟踪偏差带来的执行误差——这个误差通常很小一条生成器就够了。下面用模拟数据走一遍完整流程。光伏预测误差用AR(1)过程生成因为它有明显的时序相关性负荷误差同样带自相关但波动幅度比光伏小。储能执行误差直接给定一个比例常数。rng np.random.default_rng(7) T 24 # 日内调度时段数 n_days 800 # 历史样本天数 def make_ar_errors(sigma, rho, seed): r np.random.default_rng(seed) e np.zeros((n_days, T)) for t in range(1, T): e[:, t] rho * e[:, t-1] r.normal(0, sigma, n_days) return e err_pv make_ar_errors(sigma0.15, rho0.75, seed1) err_load make_ar_errors(sigma0.10, rho0.90, seed2) z_pv Zonotope(*build_zonotope_from_errors(err_pv, n_keep3, quantile0.92)) z_load Zonotope(*build_zonotope_from_errors(err_load, n_keep2, quantile0.92)) z_bess Zonotope(np.zeros(T), 0.02 * np.eye(T)[:, :1]) print(光伏 Zonotope 生成器数量:, z_pv.n_gen) print(负荷 Zonotope 生成器数量:, z_load.n_gen) print(储能 Zonotope 生成器数量:, z_bess.n_gen)光伏的G是3条主生成器加24条残差对角线一共27列负荷是2条主生成器加24条残差对角线共26列储能只有1列。这里的量纲都归一到标幺值光伏sigma0.15意味着预测误差标准差约为装机容量的15%负荷sigma0.10类似。rho控制时序相关性rho越大相邻时段误差越容易出现同向的“长跑”后面验证跨时段行为时要留意。储能执行误差0.02看着小但它每天充放上百个循环累积下来对SOC边界的影响不小不能直接省略。3.3 广域聚合上百个资源用Minkowski加法合并成一个Zonotope资源级Zonotope构造好之后聚合就是一次循环叠加。这一步之所以适合“广域”是因为每个资源只需要向聚合方上传一个中心向量和生成器矩阵不要求所有资源在同一时刻高频通信聚合方收到的是一组低维参数而不是成千上万个场景样本。def aggregate(resources): z_agg resources[0] for z in resources[1:]: z_agg z_agg z return z_agg # 举例3 个光伏厂站、2 个负荷聚合商、1 个储能电站 pv_fleet [z_pv] * 3 load_fleet [z_load] * 2 resources_all pv_fleet load_fleet [z_bess] z_agg aggregate(resources_all) print(聚合后生成器总列数:, z_agg.n_gen)聚合后的生成器列数等于各资源列数之和计算量和资源数量线性关系。实际工程里如果资源数量上百可以用concurrent.futures或者协程并行构造资源级Zonotope聚合这段只有矩阵拼接瓶颈不在聚合本身而在每个资源的误差数据清洗上。有一点要明确Minkowski和的这种操作假设各资源的误差变量相互独立地落在各自范围内如果两个资源实际上高度相关——比如同一片云团下的两座光伏电站——聚合结果会偏保守因为集合允许一个出现正误差而另一个出现负误差的组合而现实里它俩总是同涨同跌。所以我在做聚合前会先算一下资源间的误差相关系数超过0.6的会先合并成一个资源再进聚合列表。聚合Zonotope的中心是所有资源中心的向量和。上面示例里每个资源的中心都是零向量所以聚合中心也是零向量表示“预测误差的期望是零”。实际数据里如果某个资源有系统性偏差中心就会偏移调度时必须把这个偏移叠加到预测曲线上不然名义功率平衡方程会整体失真。4. 聚合后的鲁棒调度把Zonotope推进cvxpy求解虚拟电厂经济调控4.1 调度问题怎么描述目标函数、设备模型和三类约束聚合Zonotope就位后调度问题可以写成一个小型的鲁棒经济调度。目标函数是主网购电成本、弃光惩罚、负荷削减补偿、储能损耗四项之和设备模型包括储能SOC递推和功率上下限约束包括功率平衡、联络线容量、SOC边界。不确定性通过光伏和负荷的预测误差进入功率平衡方程鲁棒性的体现方式是联络线功率在最坏误差组合下也不能越限。这里有一个关键约定我们调度的是名义变量实际功率等于名义功率减去误差扰动。光伏出力实际值 预测值 uu落在聚合Zonotope里平衡方程里主网购电随u被动变化所以联络线约束必须对集合内所有u成立。目标函数用名义变量的成本因为误差的期望是零线性目标下的期望成本就等于名义成本——这个等价关系让问题保持线性。储能SOC模型按15分钟或者1小时的采样周期离散化我这里用1小时时段共24个时段。SOC递推关系是soc[t1] soc[t] 充电功率×充电效率 - 放电功率/放电效率起止SOC都固定在50%防止调度结果把储能“掏空”。4.2 用支撑函数把鲁棒约束写成确定性线性约束给定时段t聚合Zonotope在该时段维度上的投影范围用支撑函数算。联络线功率P_grid_t的实际值等于名义值减去误差扰动在t时刻的分量u_t。要求P_grid_t落在[P_min, P_max]内对任意u_t成立等价于给名义联络线功率两边各让出支撑函数值。具体推导是实际功率最大出现在u_t取最小值时而u_t的最小值是 -hZ(-e_t)所以名义功率上限要减去hZ(-e_t)实际功率最小出现在u_t取最大值时即hZ(e_t)所以名义功率下限要加上hZ(e_t)。用代码表述就是s_plus np.array([z_agg.support(np.eye(T)[t]) for t in range(T)]) s_minus np.array([z_agg.support(-np.eye(T)[t]) for t in range(T)])如果聚合Zonotope中心确实是零向量s_plus和s_minus数值上相等分开算的好处是万一某个资源中心有偏移约束仍然不会漏。这两条向量就代表着“这个时段调度中心必须为不确定性预留的联络线容量”值越大说明该时段的预测误差集合越宽——早晨光伏爬坡和傍晚负荷高峰时段它们通常最大。4.3 完整可运行调度代码与参数设置下面这段cvxpy模型可以直接跑通依赖前面定义的Zonotope类和aggregate函数还需要一套预测曲线数据。为避免读者被无关细节干扰光伏预测和负荷需求的曲线我用数组直接给实际部署时换成SCADA或气象预报接口的输出即可。import cvxpy as cp import numpy as np T 24 eta_ch, eta_dis 0.95, 0.95 # 储能充放电效率 cap 3.0 # 储能容量MWh P_b_max 1.0 # 储能功率上限MW P_grid_min, P_grid_max -2.0, 5.0 # 联络线功率界限MW price 0.4 0.1 * np.sin(np.linspace(0, np.pi, T)) # 模拟分时电价 pv_hat 0.5 0.3 * np.sin(np.linspace(0, np.pi, T)) # 光伏预测曲线 demand 1.0 0.4 * np.sin(np.linspace(0, np.pi, T) 0.5) # 负荷需求曲线 p_b_ch cp.Variable(T, nonnegTrue) p_b_dis cp.Variable(T, nonnegTrue) p_shed cp.Variable(T, nonnegTrue) p_grid cp.Variable(T) soc cp.Variable(T 1) constraints [] constraints [soc[0] 0.5 * cap, soc[T] 0.5 * cap] constraints [soc[t 1] soc[t] p_b_ch[t] * eta_ch - p_b_dis[t] / eta_dis for t in range(T)] constraints [0 p_b_ch P_b_max, 0 p_b_dis P_b_max, 0 soc cap] constraints [p_grid[t] pv_hat[t] p_b_dis[t] - p_b_ch[t] demand[t] - p_shed[t] for t in range(T)] constraints [0 p_shed 0.2 * demand[t] for t in range(T)] constraints [P_grid_min s_plus[t] p_grid[t] for t in range(T)] constraints [p_grid[t] P_grid_max - s_minus[t] for t in range(T)] objective cp.Minimize(cp.sum(price p_grid 5.0 * p_shed 0.01 * (p_b_ch p_b_dis))) prob cp.Problem(objective, constraints) prob.solve(solvercp.CLARABEL) print(总成本:, prob.value) print(联络线功率前6时段:, p_grid.value[:6]) print(储能SOC末值:, soc.value[-1])solver我常用CLARABEL因为它对中等问题求解快没有的话换成OSQP或者ECOS也行这些求解器处理24时段的线性规划都是秒级。几个参数值得说p_shed的惩罚系数5.0是电价量级的十几倍保证只有极端情况下才削减负荷储能损耗系数0.01是虚拟的小惩罚避免充电和放电在同一个时段内无意义地同时发生。如果你把惩罚系数调得太小会出现“弃掉光伏却同时从主网购电”的反常结果那不是模型错了而是目标函数里光储互斥的指引力不够。运行后检查s_plus较大的时段对应的p_grid如果某个时段上下界挤压后可行域为空模型会直接报错或者给一个infeasible状态。这种情况往往意味着聚合Zonotope在该时段投影太宽需要用第6章讲的生成器约简把保守性降下来或者放宽联络线容量。5. 避坑Zonotope从建模到求解最容易翻车的5个细节5.1 生成器矩阵“爆列”求解器从秒级拖到分钟级现象资源数量一多聚合后的G列数可能是几百上千cvxpy建模时约束里的1-范数项逐个展开求解时间肉眼可见地变慢。原因Minkowski和的代价就是生成器列数线性累加每个生成器在支撑函数里会变成vstack里的一个变量组列数膨胀直接拉大问题尺寸。解决对聚合后的Zonotope做约简保留列范数最大的前m条生成器把被丢弃方向的宽度吸收进一条对角残差盒子里。检查方法很简单打印一下z_agg.G.shape如果列数超过30求解时间通常已经不可接受了。5.2 中心偏移被忽略调度基准整体偏低现象蒙特卡洛回测里实际购电成本总是比调度目标值高出一截而且偏差方向一致。原因某个资源的预测误差有系统性偏差但构造Zonotope时用了均值中心化或者聚合后没有检查z_agg.c是否为零向量。名义功率平衡方程里若没有把这个偏移补回去调度中心会以为自己买电买少了。解决构造生成器矩阵前先对误差序列减中位数聚合后打印z_agg.c的范数非零就说明需要把中心向量叠加到pv_hat或者demand上再进调度模型。5.3 把Zonotope和Box混用鲁棒约束看似宽松实则漏风现象调度结果通过了Zonotope鲁棒约束但蒙特卡洛回测发现联络线功率偶尔越限。原因如果验证用的不确定性集合是独立盒子逐维区间而实际误差集合的形状是倾斜的Zonotope盒子的四个角覆盖了Zonotope没有覆盖的方向就会出现约束评估不一致。反过来也成立用Zonotope做约束但用盒子做评估会高估保守性。解决约束和验证必须用同一个集合对象。我的习惯是构造完Zonotope后单独跑一段采样脚本统计历史误差样本落在Zonotope内的比例保证覆盖率在95%以上然后再把调度约束挂到这个Zonotope上。5.4 逐时段鲁棒约束忽略跨时段相关性备用容量白白多留现象调度结果比用块对角联合Zonotope算出来的成本高出一截尤其在连续阴天或者持续高温时段。原因每个时段单独建一个Zonotope等于允许相邻时段误差任意反向组合。但光伏误差有强自相关傍晚时段和前一小时的误差高度同向逐时段约束会把这种“本来不会发生”的组合也兜进去。解决对自相关强的时段窗口把连续4到6小时的误差样本拼成一个高维误差矩阵构造联合Zonotope再做支撑函数。这样生成器列数会多一些但约束更紧成本能降下来不少。5.5 验证图横坐标太密集曲线边界看不出问题现象画出24时段的名义功率和鲁棒边界后x轴刻度全部挤在一起边界有没有越限肉眼根本判断不了。原因matplotlib默认在x轴上给每个数据点标一个刻度24个数字挤在七八厘米宽的图里就成了黑疙瘩。解决用plt.xticks(range(0, 24, 4), rotation45)控制每隔4小时标一个刻度并旋转45度。这个习惯我是在第一次画“越限点在图上搅成一团”之后养成的看似小问题实际会浪费很多排查时间。6. 验证与进阶蒙特卡洛回测、生成器约简和两个实用技巧6.1 蒙特卡洛回测采样x轴点看边际约束有没有“漏”鲁棒约束是否真的有效需要回测确认。做法是保持调度决策变量不变从[-1,1]^k采样一组x用z c G·x生成一个误差场景把它注入平衡方程算出实际联络线功率统计越限比例。rng np.random.default_rng(2024) n_sim 2000 xs rng.uniform(-1, 1, size(n_sim, z_agg.n_gen)) violations 0 for x in xs: u z_agg.G x p_grid_actual p_grid.value pv_hat - (pv_hat u) # 实际购电随扰动变化 if (p_grid_actual P_grid_min).any() or (p_grid_actual P_grid_max).any(): violations 1 print(越限比例:, violations / n_sim)这段代码的物理含义是把聚合Zonotope定义的每个极端组合都试一遍。越限比例超过2%先回去检查分位数参数和中心偏移而不是急着调调度权重。回测时把样本数固定在2000避免画图横坐标过于密集也是这个套路顺便解决的。6.2 生成器约简拓扑成型后要控制列数总量约简的常见做法是保留列范数最大的m条生成器把剩下的列对每个维度的贡献取绝对值求和作为一条新的对角盒子生成器。这个技巧在Zonotope字面量上做减维效果相当于把次要方向的不确定性膨胀成一个略大的盒子兜底能显著压缩求解规模代价是形状略微外扩。m取10到15之间通常够用超过20求解速度收益就开始变弱。6.3 两个实用技巧与最终收尾第一个技巧是把Zonotope调度结果和随机场景法对拍。随机采样1000个场景直接进调度模型求期望成本如果两者成本差距在5%以内说明Zonotope的鲁棒约束没有过度保守如果差距过大检查主方向数量是不是取少了。第二个技巧是定期更新生成器矩阵我一般每个季度重新跑一次build_zonotope_from_errors因为光伏组件衰减、负荷结构变化都会让历史误差的形态漂移一套参数用到底会逐渐失真。这些方法组合起来就是一套能从数据清洗一路做到调度求解再到验证闭环的Zonotope实施路径。我最早实现的时候光是在零中心化和分位数量纲上就折腾了两天后来把G的每一列打印出来看范数才找到问题——此后每次换数据源第一件事就是检查c和G的统计量。希望帮到你。本文还有配套的精品资源点击获取