做微电网优化调度这个方向有一段时间了最近把一个比较完整的算例跑通了基于风光储能和需求响应的微电网日前经济调度用Matlab结合Yalmip建模求解器走的是商用MILP求解器。整个过程从搭参数表、写约束、调数值到最终把调度结果画出来前前后后花了不少精力。这篇文章就把这个算例完整拆开从物理问题到数学模型再到代码实现和调试经验一次说清楚。适合正在做微电网能量管理、综合能源系统优化或者准备写相关课程设计和论文的同学参考。微电网这东西说大不大但麻雀虽小五脏俱全。风机、光伏、储能、可调负荷、上级电网每个部分都有自己的运行约束和经济特性。把这些东西塞进一个优化模型里让它在“满足所有运行约束”的前提下“一天的总运行成本最低”就是日前经济调度的核心。真正动手写代码之前如果不把物理过程理解透后面调模型的时候会被各种莫名其妙的数值问题折磨到怀疑人生。1. 先把问题拆清楚微电网经济调度到底在算谁的钱1.1 多源协同调度的目标拆解这个问题的本质是一个多时段的资源分配问题。微电网里面有几种可用的“供电来源”上级电网购电、风电、光伏、储能放电、微型燃气轮机发电还有一面比较特殊的“柔性资源”——需求响应它可以通过削减负荷来降低供电压力。有来源就得有支出。购电要花钱电费按照分时电价走微燃机发电要烧天然气有燃料成本储能充放电会造成电池寿命损耗可以折算成一个退化成本需求响应需要给用户补偿如果风机和光伏的预测出力没有被完全消纳产生了弃风弃光也要在经济性上受到惩罚。把这些成本全部加起来作为目标函数然后让优化器去决定每个时段储能到底该充电还是放电、充放多少微燃机发多少电从电网买多少电需求响应削减多少负荷允许弃多少风、弃多少光。这个决策序列就是“日前调度计划”。1.2 为什么调度周期选“日前”而不是实时很多刚接触这个方向的人会问为什么一定要做“日前”调度实时调整不行吗答案是行但那是另一个层面的问题。日前调度解决的是“计划”问题实时调度解决的是“纠偏”问题。储能的充放电决策、微燃机的出力调整、需求响应的提前通知这些动作都有很强的时间惯性。你不可能等到电价涨了才发现需要储能放电储能前一天晚上就得充满。分时电价、负荷预测、风电光伏预测这些信息在日前尺度上是可以拿到且相对可靠的。所以工程上一般的做法是先做日前计划确定各个设备的运行基点日内再根据实际出力偏差做滚动修正。这篇博文讨论的就是第一层日前计划怎么做。2. 四个参与主体的建模每一步都不能含糊2.1 风光出力预测值加弃用惩罚怎么处理风电和光伏在微电网里的角色比较特殊它们的特点是发电成本几乎为零但有出力上限而且上限值取决于天气。在日前调度里我们把预测的风电和光伏出力当成“额定可用功率”优化器只能在0到这个可用功率之间决定实际消纳多少。那消纳不完的部分怎么办工程上的做法是引入一个“弃用变量”。比如风电可用功率是 600 kW优化器决定只消纳 550 kW剩下的 50 kW 就是弃风量目标函数里加上一个惩罚系数。这个系数不能随便设它代表“浪费了本来可以白用的电”的机会成本。我之前算例里把弃风弃光惩罚设成 1.5 元/kWh比峰时购电价1.2 元/kWh还高一点这样优化器会优先想办法消纳风电比如给储能充电而不是直接弃掉。注意弃风弃光变量必须是非负的而且不能超过该时段的可用功率。约束写的是 (0 \le P_{wc}(t) \le P_{w}^{avail}(t))等会儿代码里会体现。2.2 储能SOC递推与充放互斥建模储能是微电网里最“聪明”的元件因为它有跨时段转移能量的能力。低电价时段充电高电价时段放电赚的就是峰谷价差。建模的时候最关键的是一个递推公式[ SOC(t1) SOC(t) \eta_{ch} P_{ch}(t) - \frac{P_{dis}(t)}{\eta_{dis}} ]SOC是荷电状态可以理解成储能的“油量表”。(\eta_{ch}) 和 (\eta_{dis}) 是充放电效率。注意这里有个小细节充电时能量进去要打折扣放电时能量出来也要打折扣两个都不能省。很多初写模型的人只写了一个效率算出来的经济性会偏乐观实际工程上是不对的。储能的充放电功率是有限制的SOC也不能超过上下限。还有一条同一个时段不能同时充电和放电。这个约束在数学上可以用两个二进制变量配合实现也可以用互补约束但互补约束不是线性约束求解器处理起来麻烦。工程上最直接、最稳的方式就是引入两个二进制变量[ 0 \le P_{ch}(t) \le u_{ch}(t) P_{ch}^{max}, \quad 0 \le P_{dis}(t) \le u_{dis}(t) P_{dis}^{max}, \quad u_{ch}(t) u_{dis}(t) \le 1 ]这样一来模型就变成混合整数线性规划MILP了需要用支持0/1变量的求解器来解。还有个容易被忽略的约束调度周期结束时的SOC一般要求回到初始值。为什么因为微电网是天天连续运行的如果今天结束时SOC是95%明天结束时也是95%那相当于每天都在从高荷电状态起步长期看储能并没有“透支”能量这个调度方案才具有周期可重复性。2.3 需求响应激励型DR的线性化表达需求响应的建模要分两类价格型DR和激励型DR。价格型DR需要描述电价变化如何影响用户用电行为常用的工具是价格弹性矩阵建模比较复杂需要额外的历史数据支撑。我在这个算例里用的是激励型DR也就是直接假设有一部分负荷是“可削减的”优化器可以决定每个时段削减多少削减了就要付补偿。这个简化非常实用而且完全线性。模型里只需要一个变量 (\Delta P_{DR}(t))表示 t 时段削减的负荷量约束是[ 0 \le \Delta P_{DR}(t) \le \delta_{DR} \cdot L(t) ]其中 (\delta_{DR}) 是可削减负荷比例一般取10%~20%。目标函数里加一项补偿成本用固定的补偿单价乘以削减量。补偿单价的选取很有讲究。如果补偿单价高于某些时段的购电价比如平段电价0.75元/kWh你把补偿设为0.88元/kWh那么优化器永远不会在平段调用DR因为削减负荷省下的购电费还没有付出的补偿多不划算。但在峰段电价1.2元/kWh时削减1 kWh负荷只花0.88元却省了1.2元的购电费净省0.32元优化器就会积极调用DR。这个机制本质上是在让模型自己做“成本效益分析”不需要我们人工指定每个时段削减多少。2.4 与上级电网的交互购电成本与联络线约束微电网不是孤岛运行时和上级电网之间有功率交换。这个交互环节的建模其实不难就是给定一条联络线功率上限以及分时购电价然后把购电成本写进目标函数。不过有两个地方值得提醒。第一如果微电网允许向电网反送电也就是售电那就需要增加一个售电价格和售电变量如果只是想研究最简单的情况建议先只考虑购电不让微电网倒送功率模型会简单很多。第二联络线功率约束要覆盖到每个时段不能只在总体上有约束否则模型可能在某个时段买电买超了。这是物理设备的硬限制。现实工程里微电网并网点的变压器容量是固定的联络线功率一旦越限可能直接触发保护动作。别问我怎么知道的真实项目里这种问题很常见。3. 目标函数和约束条件是怎么拧成一个数学模型的3.1 目标函数从总成本到每个分项把前面的分析汇总日前经济调度的目标函数可以写成这样[ \min \sum_{t1}^{24} \left( \pi^{grid}(t) P_{grid}(t)c^{mt} P_{mt}(t)c^{dr} \Delta P_{DR}(t)c^{cur} (P_{wc}(t) P_{pc}(t))c^{bat} (P_{ch}(t) P_{dis}(t)) \right) ]逐项解释一下(\pi^{grid}(t) P_{grid}(t))向电网购电的费用电费随分时电价变化。(c^{mt} P_{mt}(t))微燃机发电的燃料成本按发电量线性计费。严格说还要考虑启停成本但日前调度如果假设微燃机全天运行或者启停次数有限制可以先忽略启停成本简化模型。(c^{dr} \Delta P_{DR}(t))需求响应补偿每削减1 kWh发多少补贴。(c^{cur} (P_{wc}(t) P_{pc}(t)))弃风弃光惩罚体现新能源消纳的优先级。(c^{bat} (P_{ch}(t) P_{dis}(t)))储能损耗成本。只要储能动作就按吞吐量计损耗这个模型对储能“无意义的来回充放”有天然的抑制作用。3.2 约束条件每一类都对应什么物理规律完整模型里必须写的约束至少有五类。第一类是功率平衡约束这是模型中地位最高的等式约束。它表示任意时段微电网内部所有电源出力之和等于所有负荷消耗之和。用数学式子表达就是[ P_{grid}(t) P_{mt}(t) P_{w}^{use}(t) P_{pv}^{use}(t) P_{dis}(t) L(t) - \Delta P_{DR}(t) P_{ch}(t) ]注意等式右边的含义原始负荷减去需求响应削减量再加上储能充电功率。储能充电本质上也是一个“负荷”它消耗电能要放在等式的右侧。第二类是储能约束包括SOC递推式、SOC上下限、充放电功率上限、充放互斥约束以及始末SOC相等的循环约束。第三类是微燃机约束包括出力上下限和爬坡约束。爬坡约束很关键它限制了微燃机相邻时段出力的变化速率。如果忽略爬坡约束优化器可能会给出一个每5分钟跳变200 kW的出力计划现实中任何机组都做不到。第四类是联络线功率约束写成交互功率上限。第五类是需求响应约束削减量比例限制。这五类约束全部是线性的加上目标函数里也都是线性项整个模型就是一个标准的MILP可以交给求解器去求解。3.3 量纲与标度求解器喜欢什么样的模型这里分享一个非常容易踩坑的点量纲和数值标度。我之前遇到过一种情况模型写出来逻辑上完全正确但求解器一跑就是数值问题或者解出来的结果里有非常微小的“抖动”。后来排查发现问题出在数值量级差异太大。比如功率用的是kW几千、电费单价是元/kWh零点几、惩罚系数是1.5而储能退化成本是0.2它们之间没有数量级的问题。但如果你把负荷定义成了W数值是几百万而电费还是元/kWh目标函数里几百亿的数加上零点几的数求解器的容差很容易出问题。建议做法是统一用kW作为功率单位kWh作为能量单位元作为货币单位。目标函数的量级控制在 (10^2) 到 (10^5) 之间求解器跑起来会非常舒服。这个细节看起来不起眼但能帮你省掉大量的调试时间。4. Matlab代码实战从参数表到调度结果4.1 代码总体结构与数据准备Matlab代码的整体结构分四块参数定义、决策变量声明、约束与目标函数构建、求解与结果输出。这里先给出一组完整的算例参数建议直接照着敲进去。参数数值说明负荷峰值1820 kW24小时负荷曲线风机额定功率600 kW预测出力曲线光伏额定功率500 kW预测出力曲线储能容量1000 kWhSOC上下限0.2~0.9储能最大充放电功率200 kW充放电效率0.95微燃机额定功率500 kW出力下限50 kW爬坡100 kW/h购电电价谷0.35平0.75峰1.2 元/kWh峰平谷时段分开设定需求响应比例15%补偿单价0.88元/kWh弃风弃光惩罚1.5元/kWh高于峰时电价储能退化成本0.2元/kWh按充放电吞吐量计数据准备部分核心代码如下%% 基础数据定义 T 24; % 调度时段数单位小时 dt 1; % 时段长度 % 24小时负荷预测曲线kW Load [510 480 460 440 430 450 520 650 880 1200 1450 1520 ... 1480 1400 1420 1510 1580 1700 1820 1760 1600 1350 980 720]; % 24小时风电出力预测曲线kW Wind [320 350 380 360 340 300 280 260 240 220 200 180 ... 160 170 190 210 230 250 260 250 240 220 200 180]; % 24小时光伏出力预测曲线kW PV [0 0 0 0 0 20 80 180 300 400 450 480 ... 460 420 350 250 150 50 0 0 0 0 0 0]; % 分时购电价元/kWh price zeros(1, T); price(1:7) 0.35; % 谷段 23:00-06:00 price(8:17) 0.75; % 平段 07:00-16:00 price(18:21) 1.2; % 峰段 17:00-20:00 price(22:24) 0.35; % 谷段 21:00-22:00这里有意识地让风电出力集中在中夜时段光伏集中在午间负荷高峰出现在傍晚三个特征叠加在一起就会形成一个非常典型的调度场景白天光伏多中午负荷相对低晚上负荷冲高光伏归零这时候就是储能放电和微燃机出力的时候。这样的数据设置能让结果图看起来非常有“戏”。4.2 基于Yalmip的建模代码实战变量定义是建模的第一步。决策变量有联络线购电功率、微燃机出力、储能充放电功率、两个二进制变量、SOC、需求响应削减量、弃风弃光量。%% 决策变量声明 P_grid sdpvar(1, T); % 购电功率 P_mt sdpvar(1, T); % 微燃机出力 P_ch sdpvar(1, T); % 储能充电功率 P_dis sdpvar(1, T); % 储能放电功率 u_ch binvar(1, T); % 充电状态 u_dis binvar(1, T); % 放电状态 SOC sdpvar(1, T1); % 荷电状态多一个末端点 P_dr sdpvar(1, T); % 需求响应削减量 P_wc sdpvar(1, T); % 弃风量 P_pc sdpvar(1, T); % 弃光量接下来构建约束条件。这里有个小技巧SOC变量定义成 T1 维让 SOC(t1) 有位置可放递推式写起来更顺畅。SOC(1)就是初始值SOC(25)是调度周期结束时的值。%% 约束条件 Constraints []; % 功率平衡约束 P_w_use Wind - P_wc; P_pv_use PV - P_pc; for t 1:T Constraints [Constraints, ... P_grid(t) P_mt(t) P_w_use(t) P_pv_use(t) P_dis(t) ... Load(t) - P_dr(t) P_ch(t)]; end % 储能SOC递推与约束 eta 0.95; SOC_max 0.9; SOC_min 0.2; SOC_init 0.45; P_bat_max 200; Constraints [Constraints, SOC(1) SOC_init]; for t 1:T Constraints [Constraints, ... SOC(t1) SOC(t) eta * P_ch(t) - P_dis(t) / eta]; Constraints [Constraints, ... SOC(t1) SOC_min, SOC(t1) SOC_max]; Constraints [Constraints, ... 0 P_ch(t) P_bat_max * u_ch(t)]; Constraints [Constraints, ... 0 P_dis(t) P_bat_max * u_dis(t)]; Constraints [Constraints, ... u_ch(t) u_dis(t) 1]; end % 调度周期末SOC回到初始值保证方案可循环执行 Constraints [Constraints, SOC(T1) SOC_init]; % 微燃机约束 P_mt_min 50; P_mt_max 500; P_mt_ramp 100; Constraints [Constraints, P_mt P_mt_min * ones(1,T), ... P_mt P_mt_max * ones(1,T)]; for t 2:T Constraints [Constraints, ... abs(P_mt(t) - P_mt(t-1)) P_mt_ramp]; end % 联络线功率约束 P_grid_max 1200; Constraints [Constraints, 0 P_grid P_grid_max]; % 需求响应约束 dr_ratio 0.15; Constraints [Constraints, 0 P_dr dr_ratio * Load]; % 弃风弃光约束 Constraints [Constraints, 0 P_wc Wind, 0 P_pc PV];这里有两个运算细节说明一下。第一sdpvar 是 Yalmip 工具箱定义的变量类型本质上是一个“符号变量”的容器配合 optimize 函数使用如果你没用过 Yalmip先花半个小时过一遍基础语法后面就顺了。第二微燃机的爬坡约束直接写了不等式因为 P_mt 本身可正可负abs 函数在 Yalmip 里会自动展开成线性不等式约束不需要手动拆成两个约束。目标函数的构建也很直接%% 目标函数 c_mt 0.45; % 微燃机发电成本 元/kWh c_dr 0.88; % 需求响应补偿单价 c_cur 1.5; % 弃风弃光惩罚 c_bat 0.2; % 储能退化成本 obj sum(price .* P_grid) c_mt * sum(P_mt) c_dr * sum(P_dr) ... c_cur * (sum(P_wc) sum(P_pc)) c_bat * (sum(P_ch) sum(P_dis));4.3 求解器配置与结果输出求解器配置很简单但这一步对求解速度影响很大。我一般用某个商业MILP求解器在学术授权下可以正常使用。如果你没有商业求解器Yalmip 也可以调用开源的 Cbc 求解器面对这个规模的算例24时段、几十个0/1变量完全扛得住只是求解时间会略长一点。%% 求解 ops sdpsettings(solver, gurobi, verbose, 2, ... showprogress, 0, savesolveroutput, 1); result optimize(Constraints, obj, ops); if result.problem 0 % 求解成功 P_grid_opt value(P_grid); P_mt_opt value(P_mt); P_ch_opt value(P_ch); P_dis_opt value(P_dis); SOC_opt value(SOC); P_dr_opt value(P_dr); P_wc_opt value(P_wc); P_pc_opt value(P_pc); else disp(求解失败错误信息); disp(result.info); end这里 result.problem 0 表示求解成功这是 Yalmip 的一个约定。如果返回的是别的数字可以去查 Yalmip 的文档对照错误类型后面调试部分会细说。4.4 结果图表解读从调度曲线到运行经济性求解完成之后画图这个环节千万不要偷懒。调度结果的合理性很多时候一眼看图就能看出来光看数值矩阵反而容易看花眼。我建议至少画三张图。第一张是功率平衡堆叠图横轴24小时纵轴是功率把购电、微燃机、风电消纳、光伏消纳、储能放电放在正值区负荷曲线叠加在图上能直观看到每个时段是哪几个电源在“顶”负荷。第二张是SOC曲线图重点看SOC的变化是否平滑有没有突然的锯齿是否保持在上下限范围内是否回到了初始值。第三张是需求响应削减量柱状图和原始负荷、响应后负荷放一起对比。从图里能读出什么举我这次跑出来的结果为例凌晨0点到5点负荷处于低谷电价0.35元/kWh风电出力却很强此时优化器选择给储能充电把多余的风电吸收掉同时少买电网电早上6点到8点光伏开始爬坡电价进入平段储能转成放电配合光伏把早高峰顶下来到了傍晚17点到20点光伏归零、负荷冲到峰值此时电价1.2元/kWh优化器让储能全功率放电同时调用了一大块需求响应把峰时购电需求压到联络线限值以内。这套调度逻辑和“谷充峰放、削峰填谷”的工程直觉完全一致说明模型建对了。提醒一点如果最终的SOC曲线在一天内几乎没怎么动优化器给出的方案可能是“储能闲置”。这种情况要么是峰谷价差太小要么是储能退化成本定得太高导致套利收益覆盖不了损耗。不要急着改模型先检查参数设置的合理性。5. 跑模型踩过的坑无解、锯齿和DR泛滥5.1 模型无解的排查路径初学者遇到“Infeasible problem”是最崩溃的。明明感觉所有约束都是对的为什么无解我的经验是按下面这个顺序排查第一先查等式约束是否存在矛盾。功率平衡等式如果没有写对比如储能充电负荷忘了放在等式右边很可能在某个时段出现“发电机加起来都不够用”的局面。第二查SOC初值和上下限是否矛盾。比如你把SOC初值设成0.1但下限是0.2那第一个约束就无解这个错误在模型初始化阶段就能找出来。第三查二进制互斥约束。u_ch加u_dis小于等于1这个约束本身不会导致无解但如果同时规定了充放电功率的上下限而上限设置得太小也可能让某些时段没有可行动作。第四用“逐步卸载法”做测试。先把目标函数去掉只求一个可行解如果可行解都没有再逐个删约束删到哪一条模型突然可解了问题就出在哪一条上。这个方法虽然笨但非常有效。5.2 SOC曲线锯齿与数值振荡另一种常见的异常是SOC曲线呈锯齿状比如相邻时段SOC在0.45和0.5之间来回跳。出现这种情况往往不是模型错了而是最优解本身不唯一目标函数对储能充放电的某些细节不敏感。比如峰谷价差比较小的时候充多少电、放多少电对总成本差异不大优化器就在多种等价方案里随意挑了一个结果看起来就“抖”。解决办法是给目标函数加一个小的正则项比如对储能吞吐量加一个极小的惩罚系数或者对储能出力变化率加一个小惩罚让优化器偏向“少动”的方案。我一般把这种正则系数设成1e-4元/kWh对总成本影响可以忽略但能有效抑制锯齿。5.3 需求响应“过度调用”的陷阱需求响应补偿单价设置不当会出现一个很有意思的现象DR削减量几乎每个时段都顶到上限优化器“逢峰必削”。如果补偿单价高于所有时段的购电价模型当然会疯狂调用DR因为它把DR当成“免费的负负荷”。但现实中用户不可能无条件接受削负荷。所以除了给DR设置比例上限之外还要检查经济参数上DR是否被过度倾斜。一个工程上常用的做法是把DR补偿单价定在峰段电价和平段电价之间比如我算例里的0.88元/kWh。这样设计后模型只会在电价超过0.88的时段调用DR自动体现了经济选择性。5.4 求解器报错与Yalmip版本问题Matlab这套生态里Yalmip本身是不熟的人容易出问题的点。有些老版本的Yalmip对二进制变量和sdpvar的运算规则支持不完整建议直接装最新版。求解器方面如果报“No appropriate solver for problem type”说明求解器没装好或者没有被Yalmip识别。在Matlab命令行里跑一下yalmiptest它会自动检测当前环境里有哪些可用的求解器。另外目标函数里如果把sum(price .* P_grid)写成了price * P_grid虽然Yalmip多数时候能正确处理但还是建议保持维度习惯一致都是1x24的行向量直接用元素乘和sum清楚又不容易出错。6. 做完这个算例还能往哪些方向继续深挖6.1 不确定性建模从确定性到场景法日前调度用的是预测值但风电和光伏预测不可能完全准确。工程上更严谨的做法是考虑不确定性比如用场景法生成一组典型的风光出力场景每个场景带一个发生概率目标函数变成所有场景的期望成本。这样模型规模会成倍增长但调度方案对预测误差的适应性会强很多。如果要进一步做鲁棒优化那就是另一套建模思路了需要把约束条件改造成“在最坏场景下也满足”的形式。鲁棒优化的好处是不需要假设概率分布但计算复杂度更高而且最优解可能偏保守。对初学者来说我建议先把确定性模型吃透再往场景法方向过渡。6.2 碳交易与低碳调度双碳目标下碳排放约束已经成了微电网调度绕不开的话题。可以在目标函数中增加碳交易成本项超过免费配额的部分要去碳市场购买低于配额的部分可以卖出。碳排放主要来自上级电网购电火电为主和微燃机发电风电光伏和储能不产生直接碳排放。碳价给到一个合理区间后优化器会自发增加新能源消纳、减少峰时购电结果往往比单纯的经济调度更符合政策导向。6.3 考虑EV充电负荷与多微网协同电动汽车大规模接入之后微电网里的负荷结构会发生质变。EV充电桩既是负荷也是潜在的移动储能资源有序充电和无序充电对微电网运行成本的影响差异非常大。把EV充电功率作为可控变量加进模型一共也就多几个约束和一个控制变量但研究价值一下子提升不少。多条微电网互联之后彼此之间可以共享功率、互为备用这时候就涉及到多主体协同优化的问题了目标函数从单个微电网成本最小变成系统级社会福利最大化求解算法也要升级成分解算法或者分布式优化方法。做完单个微电网的日前调度之后往这个方向走基本就是当前微电网领域比较活跃的研究前沿了。我个人跑完这一整套算例最大的体会是模型的价值不在于用了多高深的算法而在于每一处细节是否贴合物理实际。一个SOC递推公式的单位写错能让整个优化结果失去工程意义一个需求响应的补偿单价设得不合理能让DR在结果里“疯狂打工”。如果你现在也正在啃微电网经济调度这个方向建议先把这条文章里的基础算例完整复现一遍再慢慢往里加自己的东西。调度曲线跑通的那一刻你对整个微电网系统的理解会有一次质的提升。
