这阵子正好在做一个三微网能量互联优化调度的课题把面向低碳经济运行目标的多微网建模方法和Matlab代码实现完整跑了一遍。说实话这类题目在学术论文里很常见但真正自己动手从零搭模型的时候坑远比想象中多。从目标函数怎么写、约束条件怎么列到Yalmip建模、求解器配置每一步都值得好好梳理一遍。这篇文章就把我这个项目的完整思路、数学模型、代码框架和踩过的坑全部整理出来给正在做多微网优化调度方向的同学或者准备做园区微网规划的朋友一个可以直接参考的模板。1. 项目到底在做什么一套模型的三个关键问题拿到“三微网能量互联优化调度”这个题目很多人第一反应是这不就是把三个微电网连起来嘛。但真要落地你需要先把三个关键问题想清楚为什么要互联、怎么体现低碳、系统用什么拓扑来建模。1.1 为什么单个微网不够用要搞多微网互联单个微网的痛点其实很直白新能源出力波动大负荷曲线和发电曲线经常对不上。比如白天光伏大发但负荷不高电用不完只能弃掉傍晚负荷上来了光伏又归零只能高价从主网购电。储能可以缓解一部分但容量有限成本也高。多微网互联之后逻辑就变了。微网A中午光伏有富余微网B刚好是工业负荷高峰那A直接通过联络线把电送给B双方都划算。这就像几个小区各自有屋顶光伏与其每家都装一块巨贵的储能电池不如小区之间拉一条线互通有无。从整个系统的角度看互联能显著降低弃风弃光率、减少主网购电压力同时提升供电可靠性。在这个项目里我用三个微网组成一个小的互联系统每个微网都配置了光伏、风电、微型燃气轮机和储能。三个微网的负荷特性和新能源装机容量故意设得不一样这样互联的调度优势才会在仿真结果里显现出来。如果三个微网长得一模一样那互联的意义就大打折扣。1.2 低碳和经济目标怎么在一个模型里统一多目标优化最怕的就是目标之间打架。运行成本最低的方案可能是凌晨低谷电价时疯狂从主网购电而这个方案折算下来碳排放往往很高。反过来碳排放最少的方案可能是尽量用光伏和风电但需要备用容量经济性又差了。常见的解决办法有三种。第一种是加权求和把碳排放量乘上一个惩罚系数加到成本里变成一个单目标问题。第二种是碳交易机制给系统一个免费碳排放配额实际排放超过配额就要花钱买少于配额就可以卖钱碳排放就有了价格信号。第三种是用多目标算法求帕累托前沿比如NSGA-II但结果是一组解集还得再做决策。我这个项目用的是第二种思路碳交易机制。原因很简单它最贴近实际政策环境而且能自然地把碳排放转换成经济成本让模型保持单目标形式可以直接用成熟的商业求解器求解。具体做法在第二章详细展开。1.3 三微网系统的组成与拓扑怎么定系统拓扑我选了最简单的辐射式结构也就是微网1和微网2之间有联络线微网2和微网3之间有联络线微网1和微网3不直接相连。别小看这个选择拓扑越复杂联络线变量就越多模型求解规模会成倍增加。对于一套用于方法验证的模型来说辐射式结构已经完全够用。每个微网内部的设备配置有一定差异这样互联调度的价值才能体现出来。微网1的负荷峰谷差大微网2的光伏装机特别多微网3的风电比较充足而且基础负荷平稳。这种“各有所长、各有所缺”的设定恰恰是实际工程中最常见的情况。2. 数学模型是灵魂目标函数和约束条件的完整拆解优化调度模型的核心就两块目标函数和约束条件。很多人卡在这一步不是因为看不懂公式而是不清楚每个变量、每条约束在工程上到底代表什么意思。我尽量用大白话把这部分讲透。2.1 目标函数钱和碳一起算我用的目标函数是所有微网总运行成本最小化同时通过碳交易机制体现低碳要求。总的表达式如下$$\min ; F \sum_{t1}^{T} \left[ C_{\text{buy}}(t) - C_{\text{sell}}(t) \sum_{i1}^{M} \left( C_{\text{fuel},i}(t) C_{\text{om},i}(t) \right) C_{\text{carbon}}(t) \right]$$逐项拆解一下$C_{\text{buy}}(t)$ 和 $C_{\text{sell}}(t)$ 是t时段系统从主网购电的费用和向主网售电的收入。购电价格采用分时电价峰时贵、谷时便宜这是驱动微网在低谷充电、高峰放电的关键信号。$C_{\text{fuel},i}(t)$ 是第i个微网微型燃气轮机的燃料成本。燃机出力-燃料消耗是二次函数关系但为了用线性求解器我做了分段线性化处理。$C_{\text{om},i}(t)$ 是设备运维成本包含燃机、储能、光伏逆变器的维护费用。运维成本通常按出力比例折算数值不大但必须算进去。$C_{\text{carbon}}(t)$ 是碳交易成本。它的表达式是碳价 ×实际碳排放量 − 免费配额。实际排放量包含两部分一部分是从主网购电对应的排放按电网平均排放因子折算另一部分是微燃机发电的直接排放。免费配额则根据系统总负荷水平乘以一个配额系数得到。整个目标函数的经济含义就是所有微网作为一个整体在各种技术约束和碳排放约束下找到一组最优的出力计划让总费用最低。碳排放并没有脱离经济目标单独优化而是通过价格信号内化到成本里。2.2 约束条件每条约束为什么非有不可约束条件决定了解空间的边界。少一条约束求解器就可能给出“物理上不可能”的方案。我这里列出了完整的约束体系每条都说明工程含义。功率平衡约束是硬性约束每个微网在每个时段都必须满足$$P_{\text{pv},i}(t) P_{\text{wt},i}(t) P_{\text{mt},i}(t) P_{\text{ess},i}^{\text{dis}}(t) P_{\text{buy},i}(t) \sum_{j \in \Omega_i} P_{\text{line},ij}(t) P_{\text{load},i}(t) P_{\text{ess},i}^{\text{ch}}(t) P_{\text{sell},i}(t)$$这个等式想表达的意思很朴素某个微网某一时刻的用电和充电需求必须由光伏、风电、燃机、储能放电、主网购电以及从相邻微网流入的功率加在一起满足。注意联络线功率 $P_{\text{line},ij}$ 是有正负方向的变量从j流向i为正。如果没有这条约束整个模型就没有意义了。微燃机出力约束包括上下限约束和爬坡约束$$P_{\text{mt},i}^{\min} \leq P_{\text{mt},i}(t) \leq P_{\text{mt},i}^{\max}$$$$-R_{i}^{\text{down}} \leq P_{\text{mt},i}(t) - P_{\text{mt},i}(t-1) \leq R_{i}^{\text{up}}$$很多初学者只写上下限漏掉爬坡约束结果调度方案里燃机出力从一个很小的值直接跳到满发这在工程上根本做不到。燃气轮机升负荷速率是有限制的必须写进模型。储能系统约束是模型里比较繁琐的部分包含充放电功率限制、SOC递推关系、以及同一时段不能同时充放电的互补约束$$0 \leq P_{\text{ess},i}^{\text{ch}}(t) \leq P_{\text{ess},i}^{\text{ch,max}} \cdot u_{\text{ch},i}(t)$$$$0 \leq P_{\text{ess},i}^{\text{dis}}(t) \leq P_{\text{ess},i}^{\text{dis,max}} \cdot u_{\text{dis},i}(t)$$$$u_{\text{ch},i}(t) u_{\text{dis},i}(t) \leq 1$$$$SOC_i(t1) SOC_i(t) \eta_{\text{ch}} P_{\text{ess},i}^{\text{ch}}(t) \Delta t - \frac{P_{\text{ess},i}^{\text{dis}}(t) \Delta t}{\eta_{\text{dis}}}$$这里的 $u_{\text{ch}}$ 和 $u_{\text{dis}}$ 是0-1二进制变量用来保证同一个储能系统不会在同一时段又充电又放电。这个约束很关键否则求解器很可能输出又充又放的“空转”方案白白增加损耗。联络线功率约束限制了两个微网之间交换功率的大小$$-P_{\text{line},ij}^{\max} \leq P_{\text{line},ij}(t) \leq P_{\text{line},ij}^{\max}$$同时还要注意能量守恒微网i流向微网j的功率等于微网j从微网i接收的功率即 $P_{\text{line},ij}(t) -P_{\text{line},ji}(t)$。在实际代码里我直接用一对变量正负表示方向避免了变量冗余。碳排放配额约束规定全系统在一个调度周期内的总碳排放量不得超过某个上限或者通过碳交易机制结算。如果采用纯惩罚形式也可以直接写进目标函数而不用单独约束。我的做法是保留碳交易项配额不达标的部分在目标函数里结算这样模型更灵活。2.3 建模细节线性化、二进制变量和量纲陷阱这一小节是我最想强调的部分。模型能不能被高效求解很大程度上取决于你的建模技巧。第一个细节是关于燃机燃料成本的非线性。燃机燃料成本与出力的关系通常在工程上用二次函数拟合但二次项会让模型变成二次约束规划求解难度大幅上升。我采用的是分段线性逼近把出力区间分成若干段每段用一条直线近似这样整个模型就变成混合整数线性规划MILPCplex和Gurobi求解MILP的效率远高于求解非线性规划。第二个细节是储能SOC的递推公式里充放电效率和时间的乘积必须统一量纲。如果功率用kW、时间用h那么能量就是kWh如果功率用MW、时间用min能量单位是MWh前后必须一致。我见过不少同学的模型结果非常离谱一查才发现是Δt写错了溶液量单位不对储能SOC一会儿超过100%一会儿变成负数。第三个细节是决策变量的边界要尽量收紧。Yalmip里定义sdpvar之后最好用约束把变量限定在一个合理的范围内。比如微燃机出力上界就是额定容量联络线功率上界就是线路允许的最大传输容量。边界越紧求解器的搜索空间越小求解速度越快。3. Matlab代码实现从数据到结果的完整流程模型写清楚之后代码实现就是按图索骥了。这个项目我全程用Matlab完成工具箱用的是Yalmip做建模求解器用Cplex。下面把整个实现流程拆开讲。3.1 工具链Yalmip搭配Cplex/Gurobi为什么是首选先说说求解器怎么选。Matlab内置的linprog和intlinprog都能解线性规划和整数规划但对于这种多微网调度模型内置求解器的性能和中大规模的通用求解器差距明显。Cplex和Gurobi是学术界和工业界公认的标杆对MILP问题有非常成熟的算法实现。学生可以通过学术许可免费申请这两个求解器这也是绝大部分论文和课题的标配组合。Yalmip是一个建模层工具它不属于求解器而是一个“翻译官”让你用接近数学表达的方式描述优化问题然后再自动转换成求解器需要的格式。用Yalmip写多微网调度模型比自己手工构造矩阵再用linprog调用要直观太多。而且Yalmip对不同求解器的兼容性好同一套代码换求解器只需改一行配置。使用Matlab版本建议在R2021b以上Yalmip在近几个版本兼容性都很好。安装Yalmip就是把解压后的文件夹加入Matlab路径然后用yalmiptest命令验证安装是否成功。Cplex安装后要注意命令行在Matlab里是否生效通常optimize命令执行时Yalmip会自动检测到可用的Cplex接口。3.2 代码框架与核心片段解读整个代码按照模块化思路组织大致分为数据定义、决策变量、约束生成、目标函数、求解和后处理六个模块。核心框架如下%% 参数定义 T 24; % 调度时段数 M 3; % 微网数量 dt 1; % 单时段时长单位小时 % 负荷数据: load_pv(i,t) 表示第i个微网t时段的光伏预测出力 % 分时电价: price_buy(t), price_sell(t) % 微燃机参数: mt_max(i), mt_min(i), ramp_up(i), ramp_down(i) % 储能参数: ess_cap(i), ess_ch_max(i), ess_dis_max(i), soc_init(i) %% 决策变量 P_mt sdpvar(M, T, full); P_ess_ch sdpvar(M, T, full); P_ess_dis sdpvar(M, T, full); SOC sdpvar(M, T1, full); u_ch binvar(M, T, full); u_dis binvar(M, T, full); P_line sdpvar(2, T, full); % 两条联络线正方向定义清楚 P_buy sdpvar(M, T, full); P_sell sdpvar(M, T, full); %% 约束集合 Constraints []; for i 1:M for t 1:T % 功率平衡 Constraints [Constraints, P_pv(i,t) P_wt(i,t) P_mt(i,t) P_ess_dis(i,t) P_buy(i,t) ... P_line_in(i,t) P_load(i,t) P_ess_ch(i,t) P_sell(i,t)]; % 燃机出力上下限和爬坡 Constraints [Constraints, mt_min(i) P_mt(i,t) mt_max(i)]; if t 1 Constraints [Constraints, -ramp_down(i) P_mt(i,t) - P_mt(i,t-1) ramp_up(i)]; end % 储能互补 Constraints [Constraints, 0 P_ess_ch(i,t) ess_ch_max(i) * u_ch(i,t)]; Constraints [Constraints, 0 P_ess_dis(i,t) ess_dis_max(i) * u_dis(i,t)]; Constraints [Constraints, u_ch(i,t) u_dis(i,t) 1]; end % SOC递推 for t 1:T Constraints [Constraints, SOC(i,t1) SOC(i,t) ... eta_ch * P_ess_ch(i,t) * dt - P_ess_dis(i,t) * dt / eta_dis]; end % 初始和终值SOC Constraints [Constraints, SOC(i,1) soc_init(i)]; Constraints [Constraints, SOC(i,T1) soc_end(i)]; end % 联络线约束 for t 1:T Constraints [Constraints, -line_max(1) P_line(1,t) line_max(1)]; Constraints [Constraints, -line_max(2) P_line(2,t) line_max(2)]; end %% 目标函数 Objective sum(sum(fuel_coef(1) * P_mt.^2 fuel_coef(2) * P_mt)) ... sum(sum(om_cost .* P_mt)) ... sum(price_buy .* sum(P_buy,1) * dt) - sum(price_sell .* sum(P_sell,1) * dt) ... carbon_price * (total_emission - total_quota); %% 求解 ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); result optimize(Constraints, Objective, ops); %% 结果提取与绘图这里有几个细节值得说明。燃料成本如果保留了二次项模型就变成二次规划QP求解器处理起来也还行但如果想稳妥地用MILP就把二次项分段线性化。代码里的fuel_coef(1) * P_mt.^2只是一个示意实际项目中我更推荐预先算好分段线性化的系数再用辅助变量表示。联络线功率的方向和归属需要特别小心。三条微网两条联络线我建议在定义变量之前先画一张系统拓扑图把联络线编号和正方向标清楚写代码的时候就不容易乱。不少人在这个环节弄反了方向结果互联功率变成“互坑”功率仿真结果完全失真。目标函数里sum(P_buy,1)是对三个微网同一时段的购电功率求和这样恰好对应分时电价向量。这个地方如果维度不对Matlab会直接报错或者给出错误结果每次运行后都要检查结果里购电费用的数量级是否合理。储能SOC的终值约束也很重要。研究表明如果允许SOC环境随意变化调度模型往往会倾向于在调度周期结束时把电量耗尽因为这种“吃干榨净”的策略看起来成本最低但对长周期运行是不可持续的。所以我会加上SOC(i,T1) soc_end(i)让储能系统在一个调度周期后回到初始电量保证调度的可重复性。3.3 算例设计与结果分析独立vs互联对比为了验证互联调度的价值我设计了两个场景做对比。场景一三个微网独立运行各自只能从主网购售电微网之间不交换功率。场景二三个微网通过联络线互联可以互相交换功率。其他参数完全一致。基础数据方面我选取某地典型日的24小时负荷、光伏出力和风电出力数据。三个微网的参数设计如下微网编号负荷峰值/kW光伏装机/kW风电装机/kW燃机容量/kW储能容量/kWh特征说明微网132010060200100负荷峰谷差大晚高峰明显微网22002604010080光伏充裕午间出力溢出微网325080160150100风电充足夜间出力大跑完两个场景结果差异非常明显对比指标独立运行互联运行变化幅度总运行成本/元2687.52421.3-9.9%碳排放总量/kg1842.61597.4-13.3%弃风弃光率/%7.82.1-73.1%主网购电总量/kWh862604-29.9%数字背后的逻辑很清楚互联之后微网2午间富余的光伏不再被白白丢弃而是通过联络线送到微网1和微网3微网3夜间的富余风电也能在微网1晚高峰时发挥作用。新能源利用率提高主网购电和燃机出力双降碳排放自然就下来了。调度的细节规律也很典型。联络线功率在午间和傍晚呈现出明显的方向性午间总体上从微网2流向微网1和微网3傍晚则部分反转。这是分时电价引导下的合理行为电价高峰时段互联系统整体减少主网购电电价低谷时段则适当购电配合储能充电。如果画一张24小时联络线功率曲线图能看到非常规律的“潮汐式”流动这就是多微网互联的经济性体现。4. 实操中的坑和排查经验这部分是我最想分享的。模型本身照书可以搭出来但真正让模型“跑得通、解得对、算得快”需要大量的实测调试经验。4.1 模型无解或结果离谱先查这五个地方遇到Infeasible problem或者结果数值异常先别怀疑求解器大概率是模型本身的问题。按照优先级排查这五个位置功率平衡等式两边的项是否写齐。我调试时经常发现少写了一项比如储能充电功率没有加到等号右边导致系统凭空少了一部分负荷求解器为了满足等式会强行压低燃机出力结果看起来“很绿色”但物理上根本不成立。SOC递推关系是否用了上一时段的值。Yalmip里的变量是符号性的SOC(i,t1)和SOC(i,t)都必须在约束中定义清楚。很多人复制代码时把SOC索引写错递推关系循环没对上导致整个储能模块失效。联络线功率方向是否一致。三条微网两条联络线每个变量都必须明确谁是送端、谁是受端。如果正方向定义不一致约束内部会相互矛盾模型无解。爬坡约束的第一时段。爬坡约束需要用到 t-1 的燃机出力t1 时没有前值代码里必须用if t 1跳过否则索引越界。有些人不做判断直接把所有时段的爬坡约束一次性加上Matlab直接报维度错误。免费碳排放配额是否设置得过低。配额设得过低会导致约束无解特别是新能源渗透率较低的场景。可以先跑一个不带碳交易约束的版本看看实际碳排放水平再根据这个数值合理设置配额。4.2 求解太慢从模型简化和参数调整入手对于24时段的三个微网模型Cplex求解时间通常是秒级几秒钟就出结果。如果你发现求解时间动辄几分钟甚至上小时问题一定出在模型本身。最需要警惕的是二次项。如果目标函数里包含了 $P_{\text{mt}}^2$ 这样的二次项模型就变成MIQP求解难度比MILP高一个档次。改用分段线性近似之后求解速度会提升几个数量级。然后是二进制变量的规模。每个时段每个储能设备都有两个二进制变量24时段三个储能一共144个0-1变量规模不大。但如果你的模型扩展到了96时段二进制变量数量变成576个求解时间会急剧上升。这时候可以检查一下是否真的需要0-1互补约束。如果充放电同时还伴随损耗惩罚求解器在大多数情况下不会出现同时充放电的情况我实践过可以先不加二进制变量跑一版看看如果结果没有出现充放电重叠就没必要加。另外Yalmip里给变量设置边界约束也能加速收敛。Cplex对变量边界是非常敏感的边界越紧分支定界树的剪枝效果越好。我给每个决策变量都设置了上下界特别是功率变量精确到设备容量这对求解速度有显著帮助。4.3 结果合理性检验的三板斧模型跑出结果后千万不要直接拿去写论文先做几个基本的合理性检验。第一板斧功率平衡是否严格成立。把最优解代回功率平衡等式检查每时段每个微网的残差是否为零。由于Cplex求解精度很高残差应该在1e-6量级。残差大说明约束书写有误大概率是某个变量的索引错了。第二板斧储能SOC曲线是否在一个合理范围。正常调度下SOC应该在[20%, 100%]区间内平滑变化充放电切换不会过于频繁。如果SOC曲线出现锯齿状剧烈波动说明储能的投资回报分析可能不准确或者目标函数里对储能启停的惩罚力度不够。第三板斧碳排放和用量的计算结果是否可复现。手动选一个时段按公式算一遍购电产生的排放和燃机产生的排放和代码输出的数值对比。这一步能揪出量纲错误和排放因子设置错误。最后再分享一点个人经验跑完整套流程之后我最大的体会是多微网能量互联优化调度这个题目真正有价值的不是模型本身而是分布式能源之间“互通有无”这个理念的量化体现。通过联络线功率的时空调配把富余的可再生能源送到最需要的地方经济效益和低碳效益就这样同时实现了。后续如果想扩展可以从冷热电联供、需求响应或者不确定优化这几个方向入手这套Matlab代码框架完全支持在这些方向上继续改造。对正在做类似课题的朋友我建议先把你手头的算例数据整理干净然后把本文的模型和代码框架消化透再针对你的具体问题做调整会比从零开始自己摸索高效得多。