微电网日前经济调度：风光储与需求响应建模及Python求解

前段时间我刚把一套微电网日前经济调度程序跑通。项目不大但涉及的模块不少光伏、风电、储能再加上需求响应最后用 Python 把整套日前优化模型搭出来。项目标题挂的就是“基于风光储能和需求响应的微电网日前经济调度”写的也是这些内容。今天不打算念 PPT就按我实际开发的过程把模型怎么建、代码怎么写、坑在哪一条条说清楚。如果你准备做微电网课程设计、毕业论文或者刚开始接手园区微电网能量管理系统这篇可以直接当参考。我会尽量说人话数学公式点到为止重点放在可复现的代码逻辑和实操注意点上。1. 先把调度问题的边界想清楚1.1 这个优化到底在优化什么日前经济调度的“日前”两个字意思是提前一天给未来 24 小时做运行计划。常见的时间粒度是一小时一段也就是 24 个时段也有用 15 分钟一个点、总共 96 个时段的做法。时间粒度越小控制精度越高但模型变量和求解时间都会涨。优化对象是微电网未来一天的总运行成本。成本来源主要有三块向大电网买电的费用、卖电收入如果允许余电上网、调用需求响应资源付出的补偿费用另外还可以叠加储能损耗或折旧成本。在满足负荷需求、设备物理约束的前提下求解器会帮你找到一组“明天怎么运行最省钱”的指令包括储能什么时候充、什么时候放、需求响应在哪个时段减多少负荷、什么时候从电网买电或卖电。这就像你前一天晚上根据天气预报和峰谷电价决定第二天中午太阳能多的时候给家里储能充满、晚上用电高峰时再放出来用。微电网只是把这个过程变成了数学规划问题交给求解器算。1.2 风光储和需求响应各自扮演什么角色我接手这个项目时首先把系统里的电源和灵活性资源分了个类因为不同类型资源的建模方式完全不一样。资源是否可控调度变量特点光伏仅可弃光实际出力/弃光量只能向下调不可跨时段搬移风电仅可弃风实际出力/弃风量随机性强预测误差大储能完全可控充、放功率SOC能搬移电量是灵活性核心需求响应部分可控削减/转移负荷量用户侧资源需支付补偿风光出力本质上是“预测值 上限值”。在日前优化里光伏和风机出力不需要作为连续可控电源去优化只需约束实际出力不大于预测值允许弃光弃风即可因为如果电价很高或本地缺电一定满发电价很低甚至负电价时可能主动弃光。储能是唯一能在时间维度上搬移能量的设备能量守恒约束决定了它不能凭空造电。而需求响应相当于把一部分负荷变成了柔性量系统可以在高电价或功率紧张时削减它代价是给用户补偿。这个分类想清楚后后面的数学建模就很自然。2. 数学模型怎么搭目标函数与约束条件2.1 目标函数成本最小但不是越小越“好”这个模型是典型的混合整数线性规划MILP。目标函数我按下面的方式写假设共有 T 个时段记时段 t 0,1,...,T-1Δt 为每个时段的长度小时。目标为最小化总运行成本min Σ [ λbuy(t)·Pbuy(t) - λsell(t)·Psell(t) λdr·Pdr(t) λbat·(Pch(t)Pdis(t)) ]·Δt其中λbuy(t) 是 t 时段微电网从大电网购电的电价Pbuy(t) 是 t 时段购电功率λsell(t) 是 t 时段余电上网电价Psell(t) 是 t 时段售电功率λdr 是需求响应单位补偿成本Pdr(t) 是 t 时段削减的负荷功率λbat 是储能充放电的单位损耗惩罚系数用来抑制储能不必要的频繁动作。这里有个容易被忽略的细节λbat 不一定要很大我一般取电价均值的 1%~5%目的不是把储能“用废”而是让模型不会为了极小成本差让储能来回充放。另外购电和售电两个变量同时出现时必须加上不同时发生的约束否则模型可能“左手买右手卖”套利结果一点都不物理。需求响应成本项很多人会漏掉。需求响应不是免费的模型之所以削减负荷是因为削减这部分负荷省下的购电成本比补偿费用高。如果只加 Pdr(t) 变量而不在目标函数里放惩罚成本那么优化结果一定会把能减的负荷全减光因为不花钱还省成本这点一定注意。2.2 功率平衡与非同时充放电约束微电网不管怎么调度每一时刻都必须满足功率平衡。我用的功率平衡关系是P_pv(t) P_wt(t) P_buy(t) P_dis(t) P_load(t) - P_dr(t) P_ch(t) P_sell(t)其中 P_pv 和 P_wt 是光伏/风电实际出力P_dis、P_ch 是储能放电/充电功率P_load 是原始负荷预测值。左侧是所有进入系统的功率右侧是流出和消耗的功率需求响应削减后负荷变成 P_load-P_dr。储能建模的核心是 SOC荷电状态递推SOC(t) SOC(t-1) (η_ch·P_ch(t) - P_dis(t)/η_dis) · Δt / E_max其中 η_ch 和 η_dis 分别是充、放电效率E_max 是储能额定容量。SOC 要限制在安全范围内比如 10%~90%充放电功率要设上限还有一个很多人会漏掉的约束同一时段不能既充电又放电。虽然目标函数一般不会故意同时充放但没有二进制约束约束住数值求解时有可能出现“一边充、一边放”的伪平衡解。解决办法是引入两个二进制变量 u_ch(t)、u_dis(t)满足P_ch(t) ≤ P_ch_max · u_ch(t) P_dis(t) ≤ P_dis_max · u_dis(t) u_ch(t) u_dis(t) ≤ 1需求响应建模要分类型。这个项目里我用的是可削减负荷模型即每个时段最多只能削减原始负荷的一定比例0 ≤ P_dr(t) ≤ δ·P_load(t)除了单时段比例限制我还会加一个总需求量约束例如全天削减量不超过总负荷的一定比例Σ P_dr(t)·Δt ≤ ε·Σ P_load(t)·Δt这样才不会出现模型为了省钱把某几个时段负荷削到零的极端方案。购售电之间也存在类似的互斥约束。引入二进制变量 u_buy、u_sell分别控制购电和售电状态保证两者不同时发生。电网交换功率也要设上限比如联络线容量。2.3 调度周期、时间粒度与参数表做日前调度前先列参数表再写代码。表里每项都要有明确单位和来源。我一般用下面这张表参数符号单位说明分时购电价λbuy(t)元/kWh日前公布的分时电价光伏预测出力P_pv_fc(t)kW根据气象预测生成风电预测出力P_wt_fc(t)kW根据风速预测生成负荷预测P_load(t)kW历史负荷预测算法储能容量E_maxkWh电芯可用容量储能最大充电功率P_ch_maxkW受 PCS 限制储能最大放电功率P_dis_maxkW受 PCS 限制储能充电效率η_ch/一般 0.90~0.98储能放电效率η_dis/一般 0.90~0.98SOC 上下限SOC_min, SOC_max/保护电池需求响应补偿λ_dr元/kWh合同约定需求响应比例上限δ/单时段上限全日需求响应比例上限ε/总量上限时间粒度为 1 小时时 T 就是 24如果负荷波动大建议按 15 分钟取 96 个点。时间粒度变小后储能递推公式里的 Δt 必须相应调整这是新手最容易栽的地方。3. Python 代码实现的完整套路3.1 数据准备把 CSV 读进来再转成字典我习惯先建一个 load_data.py把原始 CSV/Excel 读进来处理缺失值和单位然后统一转成字典方便传给 Pyomo。import pandas as pd df pd.read_csv(input_data.csv, parse_dates[time]) df df.set_index(time) df df.resample(1H).mean().interpolate() load df[load].to_dict() pv_fc df[pv_forecast].to_dict() wt_fc df[wind_forecast].to_dict() buy_price df[buy_price].to_dict() sell_price df[sell_price].to_dict()这段代码里做了“重采样 插值”因为原始数据可能是 15 分钟一个点而日前调度想用 1 小时一个点就必须先统一时间粒度。如果你反过来把 1 小时的原始数据凑成 15 分钟插值时会有未来数据参与造成“偷看未来”的问题应该用向后填充而不是插值。这点我在第四节还会提。注意字典的 key 是 pandas 的 Timestamp传到 Pyomo 里没问题但为了方便索引我会把 key 统一改成 0 到 23 的整数。3.2 用 Pyomo 搭模型变量、约束、目标函数下面是我模型骨架里最核心的一部分。变量全部用 Pyomo 定义约束用 rule 函数写。import pyomo.environ as pyo T range(24) dt 1.0 E_max 500 P_ch_max 100 P_dis_max 100 eta_ch 0.95 eta_dis 0.95 soc_init 0.5 model pyo.ConcreteModel() model.T pyo.Set(initializeT) model.p_buy pyo.Var(model.T, withinpyo.NonNegativeReals) model.p_sell pyo.Var(model.T, withinpyo.NonNegativeReals) model.p_ch pyo.Var(model.T, withinpyo.NonNegativeReals, bounds(0, P_ch_max)) model.p_dis pyo.Var(model.T, withinpyo.NonNegativeReals, bounds(0, P_dis_max)) model.dr pyo.Var(model.T, withinpyo.NonNegativeReals) model.soc pyo.Var(model.T, bounds(0.1, 0.9)) model.u_ch pyo.Var(model.T, withinpyo.Binary) model.u_dis pyo.Var(model.T, withinpyo.Binary) model.u_buy pyo.Var(model.T, withinpyo.Binary) model.u_sell pyo.Var(model.T, withinpyo.Binary) def balance_rule(m, t): return (pv_fc[t] wt_fc[t] m.p_buy[t] m.p_dis[t] load[t] - m.dr[t] m.p_ch[t] m.p_sell[t]) model.balance_con pyo.Constraint(model.T, rulebalance_rule) def soc_rule(m, t): if t 0: prev_soc soc_init else: prev_soc m.soc[t-1] return m.soc[t] prev_soc ( eta_ch * m.p_ch[t] - m.p_dis[t] / eta_dis) * dt / E_max model.soc_con pyo.Constraint(model.T, rulesoc_rule) def dr_rule(m, t): return m.dr[t] dr_ratio_max * load[t] model.dr_con pyo.Constraint(model.T, ruledr_rule)这段代码里最值得说的是soc_rule对 t0 的特殊处理。Pyomo 的集合索引从 0 开始如果直接用m.soc[t-1]t0 时会取出m.soc[-1]求解器直接报错。所以要么从 t1 开始建约束要么像我这样在函数里加一个 if。下面把目标函数和互斥约束补全def objective_rule(m): return sum( (buy_price[t] * m.p_buy[t] - sell_price[t] * m.p_sell[t] dr_price * m.dr[t] bat_cost * (m.p_ch[t] m.p_dis[t])) * dt for t in m.T) model.cost pyo.Objective(ruleobjective_rule, sensepyo.minimize) def ch_dis_exclusive_rule(m, t): return m.u_ch[t] m.u_dis[t] 1 model.ch_dis_exclusive pyo.Constraint(model.T, rulech_dis_exclusive_rule) def pch_link_rule(m, t): return m.p_ch[t] P_ch_max * m.u_ch[t] model.pch_link pyo.Constraint(model.T, rulepch_link_rule) def pdis_link_rule(m, t): return m.p_dis[t] P_dis_max * m.u_dis[t] model.pdis_link pyo.Constraint(model.T, rulepdis_link_rule)这里还要注意p_buy和p_sell也需要互斥否则可能出现同时购售。加法和上面一样用u_buy[t] u_sell[t] 1以及两个上限约束即可。容量上限我一般在 Var 定义时直接写 bounds例如bounds(0, grid_max)。3.3 求解与结果导出模型建好后调用求解器求解。我一般用开源分支定界求解器因为模型规模不大24 时段混合整数线性规划用这类求解器足够。solver pyo.SolverFactory(cbc) results solver.solve(model, teeFalse) opt_cost pyo.value(model.cost) res pd.DataFrame({ time: list(model.T), buy: [pyo.value(model.p_buy[t]) for t in model.T], sell: [pyo.value(model.p_sell[t]) for t in model.T], ch: [pyo.value(model.p_ch[t]) for t in model.T], dis: [pyo.value(model.p_dis[t]) for t in model.T], soc: [pyo.value(model.soc[t]) for t in model.T], dr: [pyo.value(model.dr[t]) for t in model.T], }) res.to_csv(schedule_result.csv, indexFalse)求解器如果没有安装好Pyomo 会直接抛异常所以跑之前先用一个最简单的线性规划试试求解器是否可用。输出我用 DataFrame 保存后面画图、生成报表都会方便。4. 结果怎么看调度曲线、SOC 和需求响应量4.1 典型日的调度过程我用一组模拟数据跑了一遍整体趋势非常典型。白天 9 点到 15 点光伏出力大负荷相对较低所以储能在这段时间以恒功率充电把多余的光伏电量存起来下午 14 点前后还有一部分余电上网因为午间负荷低谷时段从电网买入电不划算但卖电收益有限。到了 18 点到 21 点晚高峰光伏已经掉下去电价走高储能开始放电同时需求响应也在峰荷时段削减了一部分负荷。下面这张表截取的是典型时段的调度结果正数表示放电负数表示充电时段光伏/kW负荷/kW储能充放/kW购电/kW需求响应/kW10:00180120-600012:00200160-1000014:0016019000015:00120200-400019:00303201001403020:0010360802204021:00034002903022:000280602200可以看到储能基本遵守“白天充、晚上放”的规律晚高峰需求响应削减了 30~40 kW这个削减量刚好让购电功率压在联络线限制以内。我从这个结果里还发现一个有用的判断最优 SOC 曲线应该是一整天内在上下限之间平滑变化如果看到 SOC 剧烈震荡基本是约束或时间粒度出了问题。4.2 成本构成与需求响应的价值验证为了证明需求响应确实带来经济收益我又跑了一组对比去掉需求响应只保留储能调度再跑一个含需求响应的完整模型。结果如下场景购电成本/元售电收入/元DR补偿/元总运行成本/元无DR528631004976含DR49283154305043上面的数字是示意但结论方向是对的削减晚高峰负荷以后从电网高价购电的量减少购电成本降低但因为要给用户补偿总成本可能不一定降低很多甚至略高。这时候就要看 DR 的合同报价和削减潜力是否划算。如果你发现加入需求响应后总成本反而上升不一定代表模型有问题可能是补偿单价设得偏高或者削减时段选得不合理。我在项目中一般会跑一个“DR 灵敏度分析”看不同补偿单价下的最优削减量帮助业务端制定报价策略。这一点给做工程的人提了个醒经济调度算出来的结果不只是“一套计划”还可以用来评估需求响应资源的应调用水平甚至反推合同电价。5. 调试过程中踩过的坑与排查技巧5.1 最容易被忽略的条件调试这套模型时我踩了几个比较典型的坑。第一个是起始 SOC 和终止 SOC 都设成固定值。如果只设起始 SOC不设终止 SOC模型很可能把储能里的电在一天最后几个时段全放光因为这样成本最低。但放到真实系统里第二天还要继续运行所以必须加一个终止 SOC 约束我一般让它不小于起始值。如果要求更严格可以在最后时段设置 SOC 等于给定值不过这样会损失一点灵活性。第二个是购售电互斥。没有这个约束时模型可能在某个时段同时买电和卖电只要卖电价高于对应时刻的边际价格两边差价看起来有收益。这是纯数学套利物理上不允许。我排查过很多这类问题基本都是二进制变量漏了。第三个是储能同时充放电。模型没有加互斥约束时由于充放电在功率平衡方程里符号相反会出现“充 50 kW 放 50 kW”但 SOC 不变的操作相当于白浪费能量。加了二进制约束后这种现象立刻消失。第四个是需求响应的单时段比例上限。如果只加目标函数里补偿项不给单时段削减比例设上限模型很可能在某个光伏午间低谷时段把负荷削减到接近 0因为这时候省下的购电成本可能大于补偿。但真实用户不可能完全接受所以上限约束必须作为硬约束保留。5.2 排查不可行问题的几条路径我调试时最常听到的报错就是“infeasible”不可行。碰到这个报错我一般按下面的顺序排查先去掉终止 SOC 约束看系统是否可行。如果可行说明问题出在“电量不够”或“SOC 范围太小”上。把储能充放电功率上限调大 10 倍或者把 SOC 范围放开做一个“宽松模型”测试确认问题是不是充放电能力不足导致。检查功率平衡约束有没有写错符号特别是需求响应 P_dr 到底是加在负荷侧还是电源侧。我见过不少新手把符号写错导致负荷永远是负的。检查每个变量 bounds 和耦合约束是否冲突。比如 p_dis 上限是 100 kW但同时又要求 p_ch p_dis ≤ 150 kW这种约束如果共存可能在某些时段无解。最后一个办法是把所有整数变量全改成连续变量放松成 LP 看是否可行。如果 LP 可行再逐个恢复整数变量往往就能定位到哪条整数约束太紧。定位不可行问题最好的工具不是猜而是看求解器返回的不可行约束报告。Pyomo 配合某些求解器可以输出约束违反量能直接找出是哪条约束卡住了。5.3 求解器选型与数值稳定性建议我这个模型有 24 个时段整数变量不多开源求解器完全可以。但如果时间粒度变成 96 点或者加入可平移负荷整数变量数量会上涨很多这时候开源求解器可能会变慢。实际工程中我一般这样选场景推荐思路课程设计/演示开源求解器简单易装小规模园区周期日调度开源求解器足够数百节点、多微网、日内滚动商用求解器性能更稳定嵌入实时控制系统商用求解器 热启动数值方面可以把功率单位统一成 kW容量单位 kWh电价单位元/kWh。不同数量级差距太大的系数会影响求解器数值稳定性我通常对信息化简到几百以内的量级。如果出现奇怪的最优解可以先检查是不是 SOC 递推公式里 Δt 的单位写错了。Δt 单位是小时充电功率 kW 乘以小时就是 kWh再除以 E_max(kWh)SOC 才是无量纲的比例值。6. 还可以怎么扩展从日前模型到工程应用6.1 可平移负荷怎么建模前面只做了可削减负荷。实际项目中生产线、洗衣机这类负荷更常用可平移负荷模型。可平移负荷的关键是用一个整数变量表示转移后的开始时段然后锁定额外约束。例如某台设备需要在 [a, b] 窗口内启动一旦启动就连续运行 L 小时功率恒定 P_shift。建模时引入二进制变量 y_s表示该设备在 s 时段开始运行满足 Σ y_s 1同时该设备在时段 t 的功率等于 Σ_{s: t-L1 ≤ s ≤ t} y_s·P_shift。这种模型会让问题变成真正的 MILP约束数量增加不少。我建议在可削减负荷模型跑通后再加可平移负荷否则调试难度会一次拉满。加入之后需求响应约束从“每个时段最多削多少”变成“哪些时段可以平移”更接近工业现场。6.2 加一层滚动修正日前调度用的是前一天的预测数据实际情况中光伏预测误差经常超过 20%。为了应对这一点我通常会在日前计划外面再加一层日内滚动修正类似模型预测控制MPC每 1 小时或每 15 分钟更新一次预测只优化未来 4 到 6 小时但只执行第一个控制动作。这样既保留日前调度的经济性又能用最新预测修正误差。如果在论文或项目里写我会把日前调度作为上层经济优化日内修正作为下层实时控制两者形成双层框架。不过要注意日内修正模型如果优化频率很高求解时间必须压得很低所以通常会把需求响应变量保持白天固定只让储能动态调节降低在线求解压力。这个扩展方向还可以继续做下去比如引入风光出力预测误差的场景集、做两阶段随机优化或者把储能寿命模型加进去做长期规划。但所有扩展都建议先把基础 MILP 模型吃透再一层层往上加。最后分享一个我实操中的小习惯模型结果出来之后我会先画一张“功率平衡堆叠图”把购电、光伏、风电、储能放电、需求响应削减量堆起来看看叠出来的总值和实际负荷是否对齐。如果图对不上说明功率平衡或符号处理有问题先别急着分析经济性。这个检查只要几行代码但能省掉大量排查时间。用 Python 做微电网日前经济调度本身并不是一个高不可攀的任务。难的地方在于把物理约束、经济逻辑和代码堆叠到一起后还能保持模型的解释性和可维护性。按照本文这套“先建模型、再写代码、最后逼自己画图验证”的流程走绝大多数问题都能提前暴露。希望这篇内容能帮你少走弯路。