Python实现风光储联合优化调度:MILP建模与废弃矿井抽蓄协同
发布时间:2026/10/2 18:35:12来源:尧图网络
干过电力系统调度优化的人应该都有体会风电、光伏和储能放在一起做联合优化调度难度不是简单叠加。风光的随机性、储能的多时间尺度特性、不同储能形式的性能差异任何一个环节没处理好优化结果就会明显失真。我最近用Python把这套“风电光伏电池储能废弃矿井小型抽水蓄能”的互补调度模型完整跑了一遍从数据生成、数学建模、求解到可视化都做了落地今天把过程中的思路、关键代码和踩过的坑一并分享出来。这篇文章适合两类人一是正在做新能源电力系统优化、储能调度方向研究的学生论文里需要一套能跑通、能复现的基准模型二是实际搞微电网或园区综合能源运行的工程师想用Python把风光储联合调度从概念验证推向工程原型。我会把目标函数、约束条件、代码实现和调试经验都讲透特别是废弃矿井小型抽水蓄能这个比较少见的储能形式它在建模上和电池有本质区别我会单独拿出来讲。1. 先理清思路三种电源的“脾气”和互补逻辑1.1 风电和光伏的出力特征决定了调度为什么难风电出力的核心特点是随机性和反调峰性。风速本身服从威布尔分布折算到功率之后波动幅度会被放大——风速增加一倍理论上功率可能增加八倍。尤其在我国北方地区夜间风速往往高于白天风电出力高峰恰好落在负荷低谷时段这就是典型的反调峰现象。光伏出力则完全不同它严格依赖太阳辐射白天出现高峰午间甚至可能超过负荷需求但夜间出力为零而且受云层遮挡影响分钟级波动非常剧烈。这两种电源放到一起单纯从“量”上看确实有互补性——白天靠光伏夜间靠风电但“互补”不等于“稳定”。调度的第一个难点就是风光出力曲线和负荷曲线在时间尺度上错位在波动尺度上叠加。负荷高峰通常在上午和晚间而光伏高峰在正午风电高峰在深夜三个曲线凑在一起如果不做主动的储能调节要么大量弃风弃光要么系统需要频繁从电网买电来平衡短时功率缺口。1.2 储能不是“一个东西”电池和抽水蓄能扮演的角色完全不同很多人把储能当成一个统一的概念建模实际上电池储能和小型抽水蓄能的特性差异非常大这个差异直接决定了调度模型的约束写法。电池储能是典型的功率型储能响应速度快秒级甚至毫秒级就能完成功率调节往返效率高目前锂电系统交流侧往返效率能做到90%以上但它的能量容量和功率是耦合的——想多存电就得加电芯成本线性上涨。所以电池适合处理小时级以内的高频波动相当于系统里的“短期缓冲器”。抽水蓄能则是典型的能量型储能它靠上下水库之间的水位差存储势能能量容量由库容决定做大容量的边际成本很低但响应速度慢从机组启动到满功率通常需要几分钟到十几分钟往返效率也偏低综合效率一般70%到80%。它适合处理跨小时的日内能量搬移比如把午间过剩的光伏电能转移到晚高峰使用。废弃矿井小型抽蓄就是这个思路的延伸——利用废弃矿井的竖井和巷道高差建设数十米到数百米水头的分布式抽蓄电站单机功率一般在0.5MW到10MW量级。1.3 为什么我选了Python而不是其他工具调度优化在工业界常用GAMS、MATLAB加YALMIP这些工具但我在这个项目里坚持用Python核心原因有三个。第一Python的数据生态和优化生态都极其完整pandas处理时间序列、numpy做矩阵运算、scipy直接提供线性规划和混合整数线性规划求解器matplotlib出图整个链路不需要换工具。第二可复现性强论文审稿或团队协作时工程的上一段代码在任何机器上能跑出一致结果这是很多图形化工具做不到的。第三后续如果要接入机器学习功率预测模型比如LSTM、TransformerPython天然无缝衔接预测输出可以直接作为调度模型的输入不用做跨语言的接口转换。当然Python自带的求解器在大规模问题上性能确实不如商业求解器Gurobi、CPLEX但在教学演示、方案验证、中小规模微电网调度场景下已经足够。这篇文章的所有代码只要安装numpy、scipy、matplotlib三个库就能跑通。2. 调度模型怎么搭从目标函数到约束条件的完整推导2.1 数据准备预测曲线从哪里来怎么生成测试场景调度模型的前提是已知未来一段时间本文以24小时为例的风电、光伏可用出力和负荷曲线。如果有数值天气预报数据可以用功率预测模型生成这些序列没有实测数据时用统计分布模拟典型场景是学术界和工程上通用的做法。风电可用出力我用威布尔分布随机生成风速序列再通过风速-功率特性曲线转成功率。公式是切入风速以下出力为0额定风速以上封顶中间段按三次方关系计算。光伏出力用晴空模型乘以随机云遮系数模拟负荷曲线则在基础负荷上叠加早晚两个高峰。这样生成的场景虽然不是某一天的真实数据但保留了风光出力的典型统计特征足以验证调度模型的正确性。时间分辨率方面本文取1小时、24个时段这对概念验证是合适的。实际工程中建议加密到15分钟甚至5分钟因为风光分钟级波动对电池充放电策略影响很大但分辨率越高变量数量越大求解时间越长需要根据场景规模权衡。2.2 目标函数运行成本最小化弃电惩罚必须加调度模型的目标函数我设计为购电成本加上弃风弃光惩罚的最小化。公式可以写成min Σ [ price(t) × P_import(t) ] λ × Σ [ (P_wt_fc(t) − P_wt(t)) (P_pv_fc(t) − P_pv(t)) ]其中P_import(t)是t时段从电网购入的功率price(t)是分时电价P_wt(t)和P_pv(t)是实际消纳的风电和光伏功率P_wt_fc(t)和P_pv_fc(t)是预测可用出力λ是弃电惩罚系数。这个惩罚项非常关键。如果不加它优化器为了减少购电成本会在某些时段故意不消纳风电用网购电代替——这在经济上可能成立但完全违背了新能源调度的初衷。λ的取值一般要高于最低电价我在这组参数里取200元/MWh目的是让优化器只有在“实在消纳不了”的情况下才弃电。这里有一个实现细节值得注意由于弃电量等于预测值减去实际消纳值目标函数展开后会出现常数项λ×Σ(P_wt_fcP_pv_fc)。常数项不影响优化结果我在代码里直接把它剥离了所以在目标函数系数中P_wt(t)和P_pv(t)对应的系数是负的λ。第一次写的时候如果不理解这一点很容易把自己绕晕。2.3 约束条件功率平衡、储能递推、互斥逻辑约束条件是调度模型的骨架我把它分成四组。第一组是系统功率平衡约束任意时刻所有电源出力风电光伏电池放电抽蓄发电购电必须等于负荷加上所有充电功率电池充电抽蓄泵水。第二组是储能能量状态递推。电池的模型是E_bat(t1) E_bat(t) η_c × P_bat_c(t) × Δt − P_bat_d(t) / η_d × Δt其中η_c和η_d分别是充放电效率。抽蓄的模型形式上完全一样只是把充放电换成泵水和发电工况效率换成抽蓄的泵水效率和发电效率。第三组是储能功率约束这里有个常见的坑电池不能同时充电和放电。这个“互斥逻辑”是0-1整数约束必须引入二进制变量才能严格表达。如果只用线性规划而忽略互斥优化器很可能给出“一边充电一边放电”的荒谬结果——因为充电和放电在功率平衡方程里符号相反同时进行意味着无意义的能量循环成本上看似没增加实际是物理上不可能的。第四组是上下限约束电池SOC维持在10%到90%之间抽蓄等效能量状态维持在10%到90%之间购电功率不超过联络线容量上限。3. Python实现数据生成、MILP建模与结果可视化全流程3.1 工程结构模块化拆分别把代码都堆在一个文件里项目我拆成四个部分数据生成模块、参数配置模块、优化求解模块、可视化模块。这样做的好处是改参数时不用翻一长串代码跑实验时也可以单独替换数据生成逻辑。下面这段代码是数据生成和参数配置的核心部分直接附上。import numpy as np from scipy.optimize import milp, LinearConstraint, Bounds import matplotlib.pyplot as plt # ---------- 1. 生成测试场景数据 ---------- rng np.random.default_rng(42) T 24 # 调度周期24小时 dt 1 # 时段时长1小时 # 风电功率: 威布尔风速 - 功率曲线 v rng.weibull(2.2, T) * 8 v_in, v_r, v_out 3.0, 12.0, 25.0 P_wt_fc np.zeros(T) mask1 (v v_in) (v v_r) mask2 (v v_r) (v v_out) P_wt_fc[mask1] 150 * (v[mask1]**3 - v_in**3) / (v_r**3 - v_in**3) P_wt_fc[mask2] 150 # 光伏功率: 晴空模型 x 云遮系数 h np.arange(T) clear np.clip(np.sin(np.pi * (h - 6) / 12), 0, 1) cloud np.clip(rng.normal(0.8, 0.15, T), 0, 1.2) P_pv_fc 100 * clear * cloud # 负荷: 基础负荷 早晚高峰 随机扰动 load_base 200 load_peak1 30 * np.exp(-((h - 8) / 3)**2) load_peak2 40 * np.exp(-((h - 19) / 4)**2) P_load load_base load_peak1 load_peak2 rng.normal(0, 5, T) # 分时电价(元/MWh): 峰平谷三段 price np.where((h 10) (h 15), 900, np.where((h 18) (h 22), 800, 400))3.2 核心模型scipy的MILP怎么搭矩阵这里我直接给出混合整数线性规划的完整建模代码。变量排列是关键我按每个时段11个变量组织风电消纳、光伏消纳、电池充电、电池放电、抽蓄泵水、抽蓄发电、购电功率、电池SOC、抽蓄能量状态、电池互斥标志、抽蓄互斥标志。这样在构建约束矩阵时通过索引偏移就能精确定位变量。# ---------- 2. 系统参数 ---------- lambda_curtail 200.0 # 弃电惩罚系数 P_import_max 50.0 # 联络线上限 P_bat_rated 20.0 # 电池额定功率MW E_bat_max 100.0 # 电池容量MWh E_bat_init 50.0 # 初始SOC对应能量 eta_bat_c, eta_bat_d 0.95, 0.95 P_ph_rated 30.0 # 抽蓄额定功率MW E_ph_max 300.0 # 抽蓄等效能量容量MWh E_ph_init 150.0 eta_ph_pump, eta_ph_gen 0.85, 0.90 # ---------- 3. 构建优化模型 ---------- NVAR T * 11 c np.zeros(NVAR) for t in range(T): base t * 11 # 剥离常数项后, 多消纳风电/光伏相当于减少罚金, 系数为负 c[base] -lambda_curtail # 风电 c[base 1] -lambda_curtail # 光伏 c[base 6] price[t] # 购电 rows [] lb [] ub [] for t in range(T): base t * 11 idx_wt, idx_pv base, base 1 idx_bc, idx_bd base 2, base 3 idx_pp, idx_pg base 4, base 5 idx_imp base 6 idx_eb, idx_ep base 7, base 8 idx_yb, idx_yp base 9, base 10 # 功率平衡: 源端 荷端 row np.zeros(NVAR) row[[idx_wt, idx_pv, idx_bd, idx_pg, idx_imp]] 1 row[[idx_bc, idx_pp]] -1 rows.append(row) lb.append(P_load[t]) ub.append(P_load[t]) # 电池互斥: 充电时不能放电 row np.zeros(NVAR) row[idx_bc] 1 row[idx_yb] -P_bat_rated rows.append(row) lb.append(-np.inf) ub.append(0) row np.zeros(NVAR) row[idx_bd] 1 row[idx_yb] P_bat_rated rows.append(row) lb.append(-np.inf) ub.append(P_bat_rated) # 抽蓄互斥: 泵水工况与发电工况互斥 row np.zeros(NVAR) row[idx_pp] 1 row[idx_yp] -P_ph_rated rows.append(row) lb.append(-np.inf) ub.append(0) row np.zeros(NVAR) row[idx_pg] 1 row[idx_yp] P_ph_rated rows.append(row) lb.append(-np.inf) ub.append(P_ph_rated) # 电池SOC递推: E(t) - E_start - eta_c*P_c P_d/eta_d 0 row np.zeros(NVAR) row[idx_eb] 1 row[idx_bc] -eta_bat_c * dt row[idx_bd] dt / eta_bat_d if t 0: rhs E_bat_init else: row[(t - 1) * 11 7] -1 rhs 0 rows.append(row) lb.append(rhs) ub.append(rhs) # 抽蓄能量递推 row np.zeros(NVAR) row[idx_ep] 1 row[idx_pp] -eta_ph_pump * dt row[idx_pg] dt / eta_ph_gen if t 0: rhs E_ph_init else: row[(t - 1) * 11 8] -1 rhs 0 rows.append(row) lb.append(rhs) ub.append(rhs) # 变量边界 bounds [] for t in range(T): base t * 11 bounds [(0, P_wt_fc[t]), (0, P_pv_fc[t]), (0, P_bat_rated), (0, P_bat_rated), (0, P_ph_rated), (0, P_ph_rated), (0, P_import_max), (0.1 * E_bat_max, 0.9 * E_bat_max), (0.1 * E_ph_max, 0.9 * E_ph_max), (0, 1), (0, 1)] # 固定初始能量状态 bounds[0 * 11 7] (E_bat_init, E_bat_init) bounds[0 * 11 8] (E_ph_init, E_ph_init) # 整数变量位置: 每个时段的第10、11个变量 integrality np.zeros(NVAR) integrality[9::11] 1 integrality[10::11] 1 # 求解 constraints LinearConstraint(np.array(rows), np.array(lb), np.array(ub)) res milp(cc, integralityintegrality, boundsBounds(np.array([b[0] for b in bounds]), np.array([b[1] for b in bounds])), constraintsconstraints, options{time_limit: 10}) assert res.success, f求解失败: {res.message} x_opt res.x.reshape(T, 11) # 结果指标 cost_total sum(price[t] * x_opt[t, 6] for t in range(T)) curtail_wt sum(P_wt_fc[t] - x_opt[t, 0] for t in range(T)) curtail_pv sum(P_pv_fc[t] - x_opt[t, 1] for t in range(T)) print(f总购电成本: {cost_total:.1f} 元) print(f弃风量: {curtail_wt:.2f} MWh, 弃光量: {curtail_pv:.2f} MWh)这一段代码有几个地方是新手最容易写错的。第一个是SOC递推公式里系数的符号充电效率乘在充电功率上放电效率除在放电功率上一旦写反结果里储能会像“永动机”一样越充越多。第二个是互斥约束的方向p_bd P_rated × y_bat ≤ P_rated 的含义是当y_bat1充电状态时放电功率上限为0当y_bat0非充电时放电功率上限为额定值。写反了会导致约束完全失效。3.3 结果可视化光看数字不够调度曲线才能暴露问题优化结果最终要落到曲线上看不然很难发现异常。我习惯用三个子图展示上面是源荷功率曲线中间是储能功率与能量状态下面是电网交互与弃电情况。fig, axes plt.subplots(3, 1, figsize(11, 10), sharexTrue) ax axes[0] ax.plot(h, P_wt_fc, label风电可用出力) ax.plot(h, P_pv_fc, label光伏可用出力) ax.plot(h, x_opt[:, 0] x_opt[:, 1], --, label实际消纳) ax.plot(h, P_load, colorgray, label负荷) ax.set_ylabel(功率(MW)) ax.legend() ax axes[1] ax.step(h, x_opt[:, 2] - x_opt[:, 3], wheremid, label电池充放(充/-放), colortab:blue) ax.step(h, x_opt[:, 4] - x_opt[:, 5], wheremid, label抽蓄泵水/发电, colortab:green) ax.plot(h, x_opt[:, 7], colortab:blue, linestyle--, label电池SOC) ax.plot(h, x_opt[:, 8], colortab:green, linestyle--, label抽蓄能量状态) ax.set_ylabel(功率(MW)/能量(MWh)) ax.legend() ax axes[2] ax.bar(h, x_opt[:, 6], label购电功率, colortab:orange) ax.plot(h, P_wt_fc - x_opt[:, 0], colortab:red, label弃风) ax.plot(h, P_pv_fc - x_opt[:, 1], colortab:purple, label弃光) ax.set_xlabel(时刻(h)) ax.set_ylabel(功率(MW)) ax.legend() fig.tight_layout() plt.show()正常的调度结果长什么样简单说低谷时段电价低电池和抽蓄都倾向于充电购电做能量搬移高峰时段优先用储能放电而不是高价购买电网电力弃风弃光只在风光大发且储能满仓、电网无法消纳的极端时段出现。如果曲线完全违背这些直觉那大概率是某个约束写错了。4. 废弃矿井抽水蓄能特性、建模与协同要点4.1 为什么废弃矿井能重新变成储能电站废弃矿井做抽水蓄能本质上利用的是“高度差”和“地下空间”。煤矿关停后竖井的深度一般在200米到800米之间巷道网络形成巨大的可储水空间。工作原理很简单上水库放在井口附近的地表蓄水池下水库利用井底巷道空间需要储能时把上水库的水放到井下需要发电时把井下蓄水抽回井上水流经水轮机带动发电机。相比传统抽水蓄能电站动辄上千万立方米库容废弃矿井小型抽蓄的优势在于分布式和灵活性。单机功率小选址限制少而且把一个废弃矿山改造成储能设施本身就实现了资源再利用减少地面沉降和水位变化带来的环境隐患。更重要的是这种储能形式能量容量由地下空间决定扩容成本低——矿井巷道挖好了就是现成的库容不需要像电池一样每增加一度电容量都要采购电芯。4.2 建模差异效率、规模和运行约束都和电池不一样从调度模型的视角电池和抽蓄虽然都可以抽象成“能量状态功率输入输出”但参数和约束的差异很大。第一是效率结构电池充放电效率可以分别设置而抽蓄的泵水效率和发电效率也分开设置但综合往返效率明显更低——我做参数设置时电池往返效率约90%抽蓄综合往返效率约76%这意味着抽蓄不适合频繁的短时充放循环否则效率损耗会把经济性吃掉。第二是能量容量规模本文中抽蓄等效能量容量设为300MWh是电池的3倍这是刻意体现两者定位的差异——电池负责短时功率调节抽蓄负责长时能量搬移。第三是运行约束的颗粒度真实的抽蓄机组有最小技术出力、启停次数限制、工况转换时间等约束在这些实验代码中没有全加但做工程落地时必须考虑否则调度指令下发到现场根本执行不了。4.3 与电池储能的协同一个管高频一个管搬运两种储能协同调度的价值在结果曲线上体现得非常直观。看优化结果时你会发现电池的SOC曲线波动频繁充电和放电的切换次数多明显在跟踪分钟级到小时级的功率波动充当“高频缓冲器”的角色而抽蓄的能量状态呈现平滑的长周期摆动——通常是夜间泵水蓄能、午间也吸收多余光伏、早晚高峰发电一天只切换一两个完整循环。这种协同模式背后有明确的物理和经济逻辑电池响应快但容量有限让它承担基荷级的长时间能量搬移会过早耗尽容量且循环次数激增影响寿命抽蓄效率低但容量大、每度电的存储成本便宜让它频繁启停去应对短时波动不仅响应跟不上效率损失也不划算。MILP模型天然就能捕捉到这个分工因为优化器会权衡两种储能在效率、容量和功率上的差异自动把任务分配给“更擅长”的一方。这就是为什么模型里必须把两种储能分开建模而不是合并成一个“总储能”变量。5. 常见报错与结果异常排查我的调试实录5.1 求解失败和矩阵装配错误我调试过程中遇到最多的报错是维度和边界不匹配。scipy的milp函数要求约束矩阵的维度、lb/ub数组的长度、bounds的个数必须严格一致任何一处少了一个变量都会抛出ValueError。排查方法很简单在求解前打印NVAR、矩阵行数、lb长度这三个数值核对是否满足“行数等于lb长度”这个基本关系。更难定位的是逻辑错误而不是语法错误。比如SOC递推约束中如果t0时忘记把上一时段的e_bat[t-1]加入矩阵约束就变成每一时段只与初始值相关储能能量会在一天结束后出现明显漂移。这类错误不会导致求解失败但结果会非常离谱需要靠绘图才能发现。5.2 结果“看起来不对”的几个典型症状症状一电池和抽蓄同时大功率充放电。这几乎可以肯定是互斥约束写错了比如互斥方向弄反或者用了近似松弛却没加功率合限。症状二SOC长期贴着上限或下限。这说明惩罚权重或电价结构有问题——如果目标函数里缺少对“储能被用完”的约束优化器会在最后几个时段把储能全部放空以便少购电单日调度勉强可以但滚动调度次日就无电可用。我之前在单日测试里没加末态约束结果最后两小时抽蓄强制满功率发电电量全部清零。加一条E(T) ≥ E_init的约束这个问题立刻解决。症状三弃风弃光量为全零但购电成本异常高。这通常是λ取值过大优化器宁可高价购电也不弃新能源逻辑上没毛病但违背了经济调度的常识。λ需要标定我建议和最低电价值对比λ略高于低谷电价系统才会优先消纳新能源。5.3 参数敏感性效率和容量的影响比想象中大最后说一个经验在做参数分析时储能效率对调度结果的影响远大于容量。我把电池效率从0.95降到0.85整个系统的购电成本立刻上升十几个百分点因为每一次储能循环都在“亏电”。而把电池容量翻倍对成本的影响反而不明显因为成本瓶颈在系统高峰时段的功率支撑能力而不是总能量。这引出一个更深的体会调度优化的重点不只是“让模型跑通”而是理解参数变化时系统行为的迁移规律。比如废弃矿井抽蓄的综合效率如果实测下来只有65%而不是模型的76%它的调度优先级会大幅下降很多原本分配给它的能量搬移任务会被电池接管。建模时留出参数接口运行后做敏感性分析这个习惯能帮你避开很多“模型做完才发现参数离谱”的坑。这套代码我目前还在持续迭代最近在尝试把预测模块从统计分布换成长短期记忆网络并引入滚动时域控制框架让每一天的调度决策都能根据最新的预测信息修正。从个人体验来看先把线性规划和混合整数规划版本的调度闭环跑通把约束条件、数据流和参数标定这些基础打牢再逐步加复杂度是最高效的路径。希望这篇分享能帮你少走几步弯路。
网站建设高端定制企业官网