1. 为什么单微网不够用多微网电能互补的工程背景1.1 单微网调度的天花板我最早接触微网优化的时候其实是有个疑惑的单微网里光伏、储能、燃气轮机该有的都有了为什么非要跟其他微网互联这不是徒增复杂度吗直到我把一套园区微网的真实数据跑完才明白单微网的调度空间是有限的。以典型的光伏储能微网为例光伏出力曲线和负荷曲线天然错位中午光伏大发但负荷可能只有一小半晚上负荷上来了光伏又归零。储能呢一个400kWh的电池按每天一充一放算能覆盖的削峰时段也就两三个小时遇到连续阴天更是捉襟见肘。剩下的缺口只能向配电网买电而配电网执行的是峰谷分时电价傍晚高峰段的购电成本高得肉疼。我当时算出来的结果是中午被迫弃掉的光伏电量占光伏总发电量的两成多晚上又花高价从配网购电补缺口——白白交了冤枉钱。这就是单微网的瓶颈资源禀赋固定调节手段有限跟配电网的交互又是“单向输血”没有横向互济的通道。光伏富余时段只能弃负荷高峰时段只能买成本压不下来。1.2 多微网互联带来的“化学反应”把几个特性互补的微网通过联络线连起来调度空间立刻就不一样了。举个例子微网A屋顶光伏装得多中午发电用不完微网B商业负荷为主白天刚好是用电高峰微网C有台工业余热发电机组但夜间负荷低、机组不能随便启停。这三个微网要是在物理上互联A中午的富余光伏可以直接卖给BC夜间的富余电量可以存到A或B的储能里A缺电的晚上可以从B的储能取电。整体算下来向配电网的高峰购电量显著下降弃光率也能从20%压到个位数。电能互补的本质就是让每个微网不再孤军奋战而是把彼此的可再生出力波动、负荷曲线差异、储能冗余能力都当作系统的调节资源来用。这和电网层面的区域互济是一个道理只不过微网层面的决策粒度更细、约束更多——每个微网内部有自己的设备约束和利益诉求不能简单粗暴地当成一个整体来调度。1.3 需求响应在这张图里的位置多微网互联解决了“电源与负荷的空间错配”需求响应要解决的是“用电行为的时间错配”。用户侧负荷并不是完全刚性的适当的经济激励可以让一部分负荷平移或削减。我常说需求响应相当于给调度员手里加了一把“软调节”的旋钮高价时段让可转移负荷挪到低价时段可削减负荷直接砍掉一部分代价是付给用户一定的补偿费用。在多微网背景下需求响应的价值是双重的。对内它降低了单个微网的峰值购电需求对外它缓解了微网间联络线的拥塞让电能互补通道在高峰时段不至于被阻塞。这两个机制叠加系统总成本的下降幅度往往超出预期这也是很多论文把需求响应和多微网互补放在一起研究的原因。2. 双层优化为什么适合这个问题角色划分与模型框架2.1 上层决策者站在全局的调度中心多微网互联之后第一个绕不开的问题是谁来定规则如果所有微网都是同一个运营商的资产那可以直接写一个超大单层优化模型把所有决策变量放一起求全局最优就行。但现实里园区微网、商业楼宇微网、工业微网往往分属不同业主各有各的利益诉求调度中心不能大包大揽替它们做所有决定。这时候就得给上层一个鲜明的角色定位它掌握全局信息负责制定微网间的交互电价或者交互功率计划目的是让整个多微网系统的总运行成本尽可能低。上层不做具体设备的启停决策它只发“经济信号”比如告诉微网A和微网B“这个时段你们之间的交易电价是每度0.45元”。2.2 下层决策者各微网的本位主义每个微网都是“经济理性人”。给定上层下发的电价信号后下层只关心自己的利益最大化——本地运行成本最小化包括燃气轮机燃料费、从配网购电费、需求响应补偿费再减去向配网售电和向邻居微网售电的收入。恰恰是这种“本位主义”构成了双层模型的本质特征。同一个调度问题单层优化里每个微网的目标都服从全局最优双层优化里每个微网的目标是各自最优全局最优要通过经济信号去诱导。比如上层想让微网A中午多送电给B就得把交互电价定得让A觉得“卖出去比我存着更划算”同时让B觉得“买进来比我自家发电更划算”。电价定得太低A不卖定得太高B不买。这本质上就是一对矛盾谁也不能靠命令解决。这里要特别说明一个容易混淆的点很多同学一听“博弈”就以为一定要写分布式迭代算法反复交换数据直到收敛。其实在微网调度这个场景里多数论文采用的是单次集中式求解——把下层的最优反应行为用KKT条件“嵌入”上层问题一次性求出主从博弈的均衡解。模型上用双层解法上可以一步到位这也是本文实现的核心思路。2.3 上下层之间的衔接变量双层模型有两个关键衔接桥梁一个是价格信号上层决策的交互电价会进入下层的目标函数另一个是功率计划下层在上层给定的电价下决策出微网间交互功率和自购/自发电策略再反馈回上层评估整体目标。在上层目标函数里微网间的电费结算属于内部转账相加后互相抵消所以上层看到的是实实在在的燃料费、配网购电费和需求响应补偿费。这样的结构可以避免“电价×功率”双线性项钻进系统总成本给后面的线性化求解省了很大麻烦。3. 数学建模功率交互、储能约束与需求响应成本3.1 微网内部设备模型先把“家底”写清楚建立双层模型前先把每个微网内部的设备约束写成数学表达式。这里用到的都是标准的混合整数线性规划MILP约束核心的四类设备如下。燃气轮机的运行成本通常压成线性或分段线性函数出力上下限P_g_min ≤ P_g(t) ≤ P_g_max爬坡约束-ΔP_down ≤ P_g(t) - P_g(t-1) ≤ ΔP_up运行成本C_g(t) a·P_g(t) b如需更精细可以分三段线性化储能的动态模型是微网调度里最容易写错的部分正确的递推式是SOC(t1) SOC(t) (η_ch·P_ch(t) - P_dis(t)/η_dis)·Δt / E_cap同时还要配套约束充电功率上限0 ≤ P_ch(t) ≤ P_ch_max·u_ch(t)放电功率上限0 ≤ P_dis(t) ≤ P_dis_max·u_dis(t)状态互斥u_ch(t) u_dis(t) ≤ 1容量边界SOC_min ≤ SOC(t) ≤ SOC_maxu_ch、u_dis是0-1变量用来保证储能不会同时充电和放电。这个互斥约束看着不起眼实际建模时非常关键如果不加优化器经常会算出“边充边放”的荒谬结果。3.2 电能互补的联络线建模微网之间功率交互的建模相对直接。定义变量P_ex(i,j,t)表示t时段微网i从微网j购入的功率那么联络线约束至少要包含传输容量上限0 ≤ P_ex(i,j,t) ≤ P_ex_max(i,j)方向互斥P_ex(i,j,t)·P_ex(j,i,t) 0避免同一联络线上同时双向送电网损可以暂时不计或者按传输功率的固定比例折算每个微网的功率平衡是整个模型的骨架P_pv(t) P_g(t) P_dis(t) P_buy_grid(t) Σ_j P_ex(i,j,t) P_load(t) - ΔP_DR(t) P_ch(t) P_sell_grid(t) Σ_k P_ex(k,i,t)左端是电源侧各出力的总和右端是负荷侧需求。ΔP_DR(t)是需求响应削减量放在右端相当于降低了等效负荷。这条约束必须对每个微网、每个时段都成立是连接所有变量的“总线”。3.3 需求响应建模激励型和价格型怎么选需求响应建模主要有两条路。激励型IDR比较适合日内调度直接定义一个可削减负荷变量ΔP_DR(t)配上限约束和补偿成本。削减量上限0 ≤ ΔP_DR(t) ≤ γ·P_load(t)γ一般在0.05到0.15之间补偿成本C_DR(t) c_DR·ΔP_DR(t)价格型PDR则是通过弹性矩阵预测负荷对电价的响应常用在较长时间尺度的分析里。ΔP_load(t) P_load_base(t)·Σ_k E(t,k)·Δρ(k)/ρ0(k)其中E(t,k)是电价弹性系数矩阵对角元是自弹性负值非对角元是交叉弹性正值。PDR的数据要求高日内优化里用得少我自己的项目里以IDR为主PDR只在扩展分析中做了一组对比。3.4 目标函数怎么定上下层各有各的算盘下层目标函数微网i的运行成本最小化min Σ_t [C_g_i(t) C_DR_i(t) π_grid_buy(t)·P_buy_grid_i(t) - π_grid_sell(t)·P_sell_grid_i(t) Σ_j π_ex_i_j(t)·P_ex_i_j(t)]注意π_ex_i_j(t)是交互电价站在下层视角它由上层给定是已知参数站在整体模型视角它和P_ex都是决策变量这就容易生成双线性项后面专门讲化解办法。上层目标函数区域总成本最小化min Σ_t Σ_i [C_g_i(t) C_DR_i(t) π_grid_buy(t)·P_buy_grid_i(t) - π_grid_sell(t)·P_sell_grid_i(t)]两组目标看起来很像但关键差异在最后一项上层目标里没有微网间的电费结算项。因为A付给B的电费就是B从A收到的电费在系统总账上互相抵消了。这个设计不仅符合“调度中心追求总成本最优”的定位还避免了上层目标中的双线性项是建模时需要特别留意的一步。4. 求解路线的分岔口KKT转化与智能算法选哪个4.1 KKT转化把下层的“理性”变成上层的约束先明确一点下层问题是一个线性规划LP在Slater条件下它的最优解与KKT条件完全等价。于是我们可以把下层LP的KKT条件作为一组约束塞进上层问题双层优化就转化成了单层优化。KKT条件分成四类稳定性条件Stationarity目标函数对每个决策变量的导数加上各约束带拉格朗日乘子的导数求和等于零原始可行性Primal feasibility下层模型自己的全部约束对偶可行性Dual feasibility所有乘子大于等于零互补松弛Complementary slackness每个乘子与对应的不等式约束左边乘积等于零前三个条件都是线性的麻烦的是第四个——“乘积等于零”是非线性约束。好在工程上有标准处理办法引入大M和二进制变量把它线性化。4.2 大M法互补松弛条件的线性化对每一对互补条件 a ≥ 0, b ≥ 0, a·b 0等价于引入0-1变量z之后的一组线性约束a ≤ M·(1 - z)b ≤ M·z当z1时a被迫为0b自由当z0时b被迫为0a自由。这样“至少有一个为0”就实现了而且全是线性不等式。M取值是个真功夫活。取太小最优解会被错误截断比如本该有功率交互的时段被M卡死取太大数值病态严重CPLEX/Gurobi要么收不动要么精度差。我的做法是先解一个把互补约束全部松弛掉的LP看各变量和乘子的最大可能量级然后取一个比最大可能值大一个数量级的值基本不会出大问题。4.3 为什么放弃智能算法嵌套不少论文喜欢用粒子群PSO、遗传算法GA嵌套线性规划来解双层模型外层搜电价内层用LP算微网最优响应。这个方法框架简单不用推导KKT对非线性目标适应性强但我在工程里吃过亏外层智能算法参数敏感种群规模、惯性权重、交叉概率稍稍一变结果差十万八千里没有最优性保证跑十次八次不同结果你不知道哪个是全局最优计算量爆炸3个微网24时段外层每次迭代都要内层调用三次LP跑2000次迭代相当于解6000个LP一台普通电脑要跑很久而本文这种LPLP结构的双层问题KKT转化后是一个MILPCPLEX/Gurobi通常几秒到几十秒就能拿到全局最优解。只有在模型不得不保留非线性比如燃气轮机二次成本不线性化的情况下KKT条件变成非线性互补问题我才会考虑用智能算法或者启发式方法。5. Matlab实现核心YALMIP建模、大M法和求解器配置5.1 环境配置YALMIP加商用求解器用Matlab实现这套模型我的标配是YALMIP加一个商用求解器CPLEX或Gurobi。YALMIP是Matlab下的一个建模层不用手工拼求解器的API声明变量和约束的方式非常接近数学表达调试代码时心智负担小很多。安装完成之后先在命令行运行一次测试确认求解器连通正常。确认之后再进主循环能省掉后面大量“为什么我的optimize不工作”的烦恼。5.2 变量定义与约束拼接的核心写法下面是典型的核心代码骨架完整项目代码会放到文末的说明中这里先把变量定义和约束拼装思路讲清楚。%% 双层模型转单层后的YALMIP实现框架 % 基本参数 n_mg 3; % 微网数量 n_t 24; % 优化时段数默认1小时一步 % 决策变量 P_g sdpvar(n_mg, n_t); % 各微网燃气轮机出力 P_ch sdpvar(n_mg, n_t); % 储能充电功率 P_dis sdpvar(n_mg, n_t); % 储能放电功率 P_ex sdpvar(n_mg, n_mg, n_t); % 微网间交互功率 P_buy sdpvar(n_mg, n_t); % 从配网购电 P_sell sdpvar(n_mg, n_t); % 向配网售电 Delta_DR sdpvar(n_mg, n_t); % 需求响应削减量 u_ch binvar(n_mg, n_t); % 充电状态 u_dis binvar(n_mg, n_t); % 放电状态 % 下层LP的拉格朗日乘子KKT转化后变成上层问题的变量 lambda_bal sdpvar(n_mg, n_t); % 功率平衡约束对偶乘子 lambda_soc sdpvar(n_mg, n_t - 1); % SOC递推乘子等式约束用自由乘子 mu_up sdpvar(n_mg, n_t); % 不等式约束乘子非负 mu_down sdpvar(n_mg, n_t); % 不等式约束乘子非负然后进入约束拼装阶段核心提醒一个问题YALMIP的sdpvar初始化为全0矩阵维度不对时经常不报错而是悄悄广播导致约束维度膨胀然后求解器报“infeasible”。我处理这类问题的方式是先把变量reshape成一维确认维度再进入约束循环。5.3 KKT条件的代码表达KKT转化后的MILP模型里最核心的部分是把下层LP的stationarity条件写成等式约束。比如对变量P_ex(i,j,t)下层目标里有π_ex(i,j,t)·P_ex(i,j,t)那么stationarity条件形如π_ex(i,j,t) Σ_mu·系数 - lambda_bal(i,t) 0这段代码建议用for循环显式生成C []; for i 1:n_mg for t 1:n_t % 功率平衡约束对应的stationarity项 C [C, lambda_bal(i,t) ... % 与P_g相关的导数项 2*a_g*P_g(i,t) b_g ... mu_up(i,t) - mu_down(i,t)]; end end补充一下KKT中的互补条件在YALMIP中可以直接用大M法矩阵化写% 假设g_ub是某组不等式约束左边mu是乘子M是充分大常数 z binvar(size(g_ub)); C [C, g_ub 0, mu 0]; % 原始可行对偶可行 C [C, g_ub M*z, mu M*(1-z)]; % 互补松弛线性化5.4 双线性项的处理这个坑我前前后后绕了好几天。当上层决策交互电价π_ex、下层决策交互功率P_ex时下层目标里存在π_ex·P_ex。站在下层的LP视角π_ex是参数没问题但KKT转化之后π_ex和P_ex都是单层模型里的决策变量乘积变成双线性项MILP求解器直接罢工。解决方案有两种亲测都走得通。方案A推荐上层目标只用区域总成本内部交易电费在系统总账上抵消因此单层模型的目标函数里不会出现π_ex·P_ex。做完这一层规避之后还要检查下层KKT条件里是否还有其他双线性项。如果下层目标函数里交互费用项是π_ex·P_ex那对P_ex求导后得到的是π_ex本身这是线性表达式不会引入双线性问题。但对π_ex求导就会出现P_ex如果上层模型里有对π_ex的stationarity条件就需要格外小心。所以我的做法是只在目标函数层面规避约束层面逐条检查即可。方案B如果模型非要保留交易结算项在上层利润里那只能用固定参考电价或交替迭代法——先给定π_ex初值求解下层得到P_ex再固定P_ex去优化π_ex反复交替直到收敛。这个方案可行但牺牲了KKT一次求解的便利性而且要额外判断收敛条件一般我建议非必要不用。5.5 求解器设置与常见参数调用求解器时我习惯把参数显式设置在sdpsettings里options sdpsettings(solver, cplex, ... verbose, 2, ... savesolveroutput, 1, ... cplex.mip.tolerances.mipgap, 0.0001); sol optimize(C, Objective, options);mipgap设置到1e-4已经足够工程精度太严反而拖慢求解时间。如果问题规模大可以先放宽到1e-2跑通逻辑验证模型没问题后再收紧。6. 仿真结果怎么看互补效果、负荷平移与成本对比6.1 算例场景设计仿真算例我设计成三个微网、24个时段。微网1是光伏富余型装机600kW本地负荷小微网2是商业负荷型白天负荷高但没有光伏微网3是工业微网负荷平稳带一套500kW/1MWh的储能和一台燃气轮机。配电网购电执行峰谷分时电价峰段1.2元/kWh平段0.75元/kWh谷段0.4元/kWh余电上网收购价0.35元/kWh。对比三种模式模式A三个微网独立调度不互联无需求响应模式B多微网互联协调调度无需求响应模式C多微网互联协调调度加上激励型需求响应6.2 成本与弃光数据对比跑完优化后整理出下表结果非常直观指标模式A独立模式B互联模式C互联DR系统总运行成本元128501112010050弃光率%22.58.35.1配电网峰值购电功率kW13801020860需求响应补偿成本元00620总成本相对模式A降幅-13.4%21.8%从B到C需求响应用它那620元的补偿成本换来了超过1000元的系统成本下降。很多人会觉得需求响应的补偿是纯增加的成本实际算过之后才知道它削掉的高峰购电支出远超补偿费。6.3 怎么验证“互补”真的发生了看总成本太抽象我习惯把微网间交互功率曲线打出来看。模式B下微网1向微网2的交互功率在11:00-14:00出现明显的正向尖峰这正是光伏大发时段——富余电被微网2吃掉了。傍晚17:00-20:00微网3的储能开始放电部分功率通过联络线送到微网2帮助它躲过晚高峰购电。这两段交互功率曲线就是“电能互补”的直接证据。另一个值得打印的曲线是微网1的弃光量。独立调度模式下弃光曲线在午间有个大凸包换成互联模式后凸包几乎消失说明原本被“扔掉”的光伏电量通过联络线找到了就地消纳的出口。6.4 需求响应带来的负荷形态变化模式C里需求响应的作用从微网2的购电曲线上看得一清二楚。原本晚高峰段的购电尖峰被削掉了一块凌晨低谷段的购电曲线略有抬升这是可转移负荷从高峰挪到低谷的典型表现。可削减负荷集中在17:00-20:00被削减削减量保持在总负荷的10%以内用户感知不明显但系统峰值负荷实实在在地降下来了。这里有个小提示仿真结果出来后一定要检查需求响应削减量有没有顶在上限。如果大量时段ΔP_DR都等于γ·P_load说明需求响应潜力已经挖尽再增加补偿单价也不会带来额外收益这时就应该考虑扩容联络线或增加储能这是模型给出的一条很有价值的边际信号。7. 我调试这套代码踩过的坑7.1 大M的取值曾经让我白跑一晚上有一版代码我偷懒把M统一设成1e6结果CPLEX跑了半个小时还在磨解出来的结果里出现了“微网1中午光伏大量弃电同时向配电网高价购电”的离谱现象。查了一晚上才发现是M太大导致互补条件的数值松弛最优性判断被污染了。后来我改成两步走第一步把互补条件全部删掉求一个松弛LP统计所有乘子的最大绝对值第二步把M设成这个最大值的10倍重新加入约束求解。从此再没出现过类似的数值问题。这个经验分享给所有做KKT转化的朋友M不是越大越好够用就行。7.2 储能SOC出现负数单位不统一的锅有次仿真结果里储能的SOC曲线出现负值我先怀疑SOC递推约束写错了检查了半天没发现问题最后发现是单位混用了——功率用了kW能量却用了MWh导致SOC增量算出来差了一千倍。后来我立了个规矩全模型功率统一用kW能量统一用kWh时间步长Δt统一用小时SOC递推公式里的η_ch和η_dis也全部提前折算成标幺值。单位理清之后这一类低级错误基本绝迹。7.3 YALMIP常见的三个报错我把自己遇到的报错归纳成三类给后来人避雷。“No suitable solver for problem type”问题被识别成MIQP或非线性但求解器不支持。解决方案是回溯模型找出产生非线性的位置优先考虑分段线性化或大M法。“Index exceeds array bounds”P_ex这种三维矩阵在循环里索引时最容易写错维度。我习惯先把三维变量按squeeze方式拆成二维逐时段处理跑通后再优化成矩阵运算。objective里的NaN或Inf几乎都是参数没赋值或者除零。解决办法是用assert检查每个参数尤其注意E_cap这种出现在分母上的量。7.4 调试顺序是效率的分水岭我见过太多同学一上来就搭完整的三微网24时段KKT模型跑不通之后对着几千行代码发呆。我的经验是严格按“1微网1时段 → 1微网24时段 → 2微网24时段 → 3微网24时段DR”的顺序递进调试。每一步都能跑出合理结果再往下一步拓展。如果小规模算例就解不动一定是模型或代码逻辑错了这时候停下来排查远比硬撑着让求解器空转有效。代码里我还习惯把约束按物理含义分组存进cell数组比如Constraints{1}存设备约束Constraints{2}存联络线约束Constraints{3}存KKT条件。出问题时直接检查对应的那组约束定位速度比一张大网式拼接快得多。最后分享一个我个人调试KKT转化模型的心得模型解出来之后除了看目标函数值顺手把互补松弛条件的残差也打出来。数值上到1e-6甚至1e-8是很正常的但如果残差落在1e-2级别八成是大M取小了或某条约束的乘子符号写反了先回头查乘子的正负比从头看公式要快得多。这套方法帮我省下过好几个通宵希望你也能用上。