月初帮一位师弟把一套“含分布式电源的配电网日前两阶段优化调度模型”的完整版Matlab代码从头到尾梳理了一遍顺手还给他补了环境配置、数据替换和结果分析这几个最容易卡壳的环节。这套代码典型的用途就是在已知分布式光伏、风电预测出力曲线的前提下把未来24小时配电网里可控机组、储能、向上级电网购电的计划一次性算出来同时用两阶段结构应对新能源出力不确定性让计划既经济又不至于在实际运行中被偏差打穿。无论你是做配电网调度、微电网能量管理、新能源并网研究的硕博生还是刚接配电网规划的工程人员把这套代码吃透基本就等于把日前调度这个方向最主流的一套研究范式拿下了。下面我把拆解过程和实操经验完整写出来照着走基本都能跑通。1. 项目整体拆解这个模型到底在优化什么1.1 为什么要做“日前两阶段”调度很多刚开始接触配电网调度的人会问直接把每个时刻的负荷、光伏、风电给定做一个24小时的最优潮流不就行了理论上确实可以但实际工程中光伏和风电的预测误差是绕不开的。早上8点的出力预测到了当天正午可能差出去20%甚至更多如果你把预测值当成确定值去安排储能充放电和机组启停那当天实际运行时电压越限、功率倒送、弃风弃光等问题会接二连三地冒出来。“两阶段”的思路就是为这个问题设计的。第一阶段叫“日前决策”在新能源出力未知量还没完全实现之前先决定那些一旦定了就很难快速调整的变量比如机组启停0-1整数变量、储能是否充放电的开关逻辑、向上级电网的购电渠道分配等。第二阶段叫“运行校正”或“场景校核”在第一阶段决策固定的前提下用若干个具有代表性的不确定性场景去检验系统运行是否安全、经济并根据场景反馈对每小时的可调变量做再优化。用一个生活化类比你每天上班前会先定一个大致路线第一阶段知道早高峰某条路常常堵车就提前把项目地点、会议时间这些“硬约束”定好真正开车上路后再根据实时路况临时变道第二阶段。两阶段模型比单阶段多出来的正是这种“先定大方向、再应对偏差”的冗余能力在数学上表现为你不是优化一条出力曲线而是优化一个“在若干场景下都可行”的运行策略。1.2 代码的整体架构与数据流这套代码从文件结构上看大致分为四个模块数据输入模块、模型构建模块、求解模块、结果输出模块。正常跑通一次完整流程数据流的走向是下面这样的第一步读入配电网拓扑参数节点数、支路阻抗、负荷基准值第二步读入预测数据24小时负荷曲线、光伏出力预测、风电出力预测第三步读入储能和机组参数容量、功率上下限、爬坡率、效率、SOC初值第四步进入YALMIP建模环节定义决策变量、目标函数、等式/不等式约束第五步调用求解器Cplex或Gurobi等求出最优解第六步从求解结果中提取数据绘制机组出力堆叠图、储能SOC曲线、电压分布图。这套架构最值得借鉴的地方在于建模层和求解层完全分离。模型用YALMIP描述求解器只是后端接口这意味着你写好一套模型后想换求解器只需改一行命令而不是把整个模型推倒重来。对于做研究的人来说这个设计能省下大量重复工作时间。1.3 为什么选Matlab而不是Python这两年Python在优化领域势头很猛Pyomo、pulp用户越来越多但配电网日前调度这个方向Matlab的生态仍然是很多研究者的首选原因很实际一是Matlab的矩阵计算和稀疏矩阵处理对大规模配电网问题先天友好写节点系数矩阵、变量索引比较顺手二是YALMIP这套建模工具箱在电力领域深耕多年网上能找到大量可直接借鉴的代码模板其中包括大量基于YALMIP的配电网/微电网调度开源代码三是Solution Pool、可行性容差、MIP gap等求解器控制参数的设置在Matlab环境下调试起来更直观。另外如果你没有商业求解器Matlab自带intlinprog也能应付中小规模算例我后面会专门说这个零成本替代方案。2. 数学建模与代码实现的关键细节2.1 目标函数到底在惩罚什么这套代码的目标函数从物理意义上看是在最小化整个配电网次日运行的总成本。最常见的写法包含四项向上级电网购电成本各时段购电功率乘以分时电价按时段累加常规机组燃料成本通常用机组出力的二次函数或线性函数近似机组启停成本为了抑制频繁启停目标函数中要惩罚启停动作弃风弃光惩罚成本当消纳不完时给一个较大的惩罚系数或者反过来在约束里限制弃风弃光率上限。用数学符号表达就是Objective sum(price(t) * P_buy(t) * dt) ... sum(a_g * P_g(t)^2 b_g * P_g(t) c_g) * dt ... sum(StartCost * u_start(t) StopCost * u_stop(t)) ... sum(penalty * P_curtail(t) * dt);这里要注意一个容易出错但非常重要的细节目标函数各项的单位必须统一。如果你的电价值单位是元/MWh功率单位是MW时间步长是1小时那计算结果是元没问题但如果你把时间步长改成15分钟dt0.25小时电价值和成本项都没变乘积就会差4倍结果直接跑偏。我见过不少人踩这个坑检查半天都找不到问题最后发现是dt没跟着改。2.2 约束条件潮流、储能、爬坡是怎么写进代码的这部分是整个模型的核心也是代码正文里占比最大的部分。先说潮流约束。为配合优化求解代码通常采用线性化的DistFlow分支潮流模型。完整DistFlow方程是非线性的但配电网辐射状结构下可以做合理简化把电压降方程中的二次项忽略就得到线性化的DistFlow% 支路潮流线性化 for t 1:T for k 1:BranchNum % 功率平衡支路流入功率 下游节点净负荷 P_branch_out(k,t) P_branch_in(k,t) ... - r(k) * (P_branch_in(k,t)^2 Q_branch_in(k,t)^2) / V0^2; % 电压降落线性近似 V_bus_end(k,t) V_bus_start(k,t) ... - (r(k) * P_branch_in(k,t) x(k) * Q_branch_in(k,t)) / V0; end end节点功率平衡约束也很直观注入功率等于负荷消耗加流出功率。在含分布式电源的情况下PV和风机出力作为一种“负的负荷”接入对应节点储能充电时是负荷、放电时是电源。储能约束是另一块重点。标准的储能模型包含SOC状态转移方程以及充放电上下限约束% 储能SOC递推 for t 1:T-1 SOC(t1) SOC(t) (eta_c * P_ch(t) * dt - P_dis(t) * dt / eta_d) / Energy_cap; end % 充放电功率上下限 for t 1:T 0 P_ch(t) u_ch(t) * P_ch_max; 0 P_dis(t) u_dis(t) * P_dis_max; end % 同一时刻不能同时充电和放电 for t 1:T u_ch(t) u_dis(t) 1; end这里有个工程上的小心得SOC初值不能随便设置为0或1因为日前优化是一个滚动过程今天的结束点就是明天的起点。代码里一般考虑日终SOC回到初值附近或者设置一个SOC终值范围避免储能建成“一次性电池”。机组爬坡约束属于典型的时序耦合约束表达的是机组在两个相邻时段之间的出力变化限制for t 1:T-1 P_g(t1) - P_g(t) R_up * u_g(t) P_g_max * (1 - u_g(t1)); P_g(t) - P_g(t1) R_down * u_g(t1) P_g_max * (1 - u_g(t)); end这两条约束里各带一个“大M项”目的是让机组在关闭状态下允许出力从0启动时不完全受爬坡率限制。这个写法很典型也是很多开源代码里最容易被抄错的地方——如果把大M项漏了模型会把“关停机组的启停过渡”和“在线机组的正常爬坡”混为一谈导致无解或次优解。2.3 分布式电源的不确定性处理场景法与多面体不确定集“两阶段”的不同实现方式直接决定了代码里不确定性建模的写法。目前主流的范式有两类范式核心思想求解复杂度适合场景随机规划构造若干典型出力场景两阶段求期望成本最优较高场景数决定规模有历史数据能生成概率场景鲁棒优化用不确定集描述出力范围求最恶劣场景下的可行解较低需要线性对偶或列约束生成数据不足侧重运行安全性完整版代码里常见的是“盒式不确定集 线性对偶”的鲁棒优化写法因为配电网日前调度强调的是“我在任何可能的出力偏差下都能安全运行”而不是“平均来看成本最低”。盒式不确定集的表达式很简单各时段光伏/风电出力在预测值附近一个给定偏差范围内波动。第一阶段决策完成后第二阶段针对不确定集里的最恶劣场景进行校核与变量调整。还有一部分完整版代码用的是随机规划写法先生成光伏出力的多个典型场景通常用K-means聚类历史出力数据得到再对每个场景分别求解第二阶段问题目标函数是各场景成本的期望。这种做法代码结构上更直观但对场景数量比较敏感——场景太少结果太乐观场景太多求解速度拖到怀疑人生。3. 实操过程从环境准备到跑通全套代码3.1 软件环境配置清单先把环境踩平再谈跑代码。我建议的配置方案如下Matlab版本R2020b及以上都可以我自己实测过R2020b、R2022a、R2023b都能正常跑最新版本也不用担心兼容性YALMIP工具箱这是一个免费开源的工具箱从这里下载后把文件夹加到Matlab路径即可求解器优先选Cplex或者Gurobi。Cplex安装相对简单但老的12.10版本和Matlab高版本比如R2024b以上存在兼容性问题Gurobi需要申请免费学术许可证但性能通常更好零成本替代如果装不上商业求解器Matlab自带的intlinprog也能用来解MILP。唯一的问题是一旦算例规模超过几百个整数变量求解时间会明显拉长小算例验证代码逻辑是够用的。配置完成后在Matlab命令行输入yalmiptest看到类似“Found installed solver: GUROBI”的提示就说明环境OK了。3.2 数据准备与参数设置详解整套代码里最有“工程含量”的不是建模语句而是数据文件的整理。以最常见的IEEE33节点配电网为例你需要准备、核对下面这些参数的取值与单位缺一个都会导致结果失真支路阻抗参数各条支路的电阻和电抗r、x值单位是Ω基准值通常取基准功率100MVA、基准电压12.66kV负荷曲线各节点24小时的有功/无功负荷基准值有些代码用标幺值有些用有名值要看清说明光伏/风电预测出力以小时为单位单位通常为MW或标幺值储能参数额定容量、最大充放电功率、充放电效率、SOC上下限、初始SOC常规机组参数出力上下限、爬坡速率、燃料成本系数、启停成本分时电价峰平谷时段划分和对应电价。建议不要跳过数据准备环节直接跑代码。拿到完整版代码后第一件事是把数据文件从示例状态改为你自己算例的参数用数据表的格式一行一行核对。很多“代码跑不通”的问题根源不是模型错了而是数据单位错位、节点号不连续、表里多了空格等低级问题。3.3 核心代码段逐块讲解下面给出一段浓缩版的核心代码展示YALMIP建模的整体框架这是整套代码的骨架%% 决策变量定义 P_buy sdpvar(1, T, full); % 向上级电网购电功率 P_g sdpvar(N_g, T, full); % 常规机组出力 u_g binvar(N_g, T); % 机组启停状态0-1变量 P_ch sdpvar(N_ess, T, full); % 储能充电功率 P_dis sdpvar(N_ess, T, full); % 储能放电功率 SOC sdpvar(N_ess, T, full); % 储能荷电状态 P_curt sdpvar(N_dg, T, full); % 弃风弃光功率 %% 目标函数 Objective sum(price .* P_buy) * dt ... sum(sum(fuel_cost * P_g)) * dt ... sum(sum(start_cost * u_g)) ... penalty * sum(sum(P_curt)) * dt; %% 约束 Constraints []; % 功率平衡约束 for t 1:T Constraints [Constraints, sum(P_buy(t)) sum(P_g(:,t)) ... sum(P_dis(:,t)) - sum(P_ch(:,t)) sum(P_dg(:,t)) - sum(P_curt(:,t)) ... sum(P_load(:,t))]; end % 储能与机组约束 % 潮流、电压约束 % 不确定集与两阶段耦合约束 %% 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.01); optimize(Constraints, Objective, ops);如果你在代码里看到的模型比这个复杂千万别慌——核心仍然逃不出上面这些变量和约束的组合。增加的细节一般体现在分布式电源的容量上限约束、变压器容量约束、更细粒度的电压约束比如用节点电压上下限V_min V_bus V_max这些都是在基本骨架上做加法。3.4 从代码跑到出图的结果解读代码跑完后记住先看求解器输出信息再画图。求解器输出的mipgap和solverTime是判断结果可信度的第一道关口如果gap一直是100%或者求解时间异常短多半是模型在求解一开始就被证明无解或者被求解器提前截断成不可行解。出图环节我习惯五张图一起看机组出力堆叠图看出力的时序合理性储能SOC曲线看SOC是否在上下限之间是否有不合理的突跳弃风弃光功率图看哪个时段消纳压力最大节点电压分布图重点检查有没有电压越限购电功率曲线结合分时电价判断购电计划是否“削峰填谷”。如果SOC曲线在凌晨充电时段直接顶到上限又瞬间跌到下限大概率不是真实调度需求而是约束写软化或初值设置有误。4. 常见问题与排查技巧实录4.1 建模报错速查表下面这份速查表里的问题是我和师弟们在实际跑这套代码时遇到频率最高的几种报错信息背后原因解决办法Undefined function or variable yalmipYALMIP没加入Matlab路径添加工具箱路径并保存pathdefThe solver is not installed求解器接口没配置好运行yalmiptest确认求解器是否被发现Warning: Solver not applicable模型是MILP但你调了非线性求解器确认变量含binvar时求解器选Cplex/Gurobi/intlinprogNaN or Inf in objective变量未初始化或数据存在空值检查数据文件添加变量初始化Infeasible problem约束过紧或大M系数选择不当先用0-1变量松弛测试逐步排查核心约束求解时间过长整数变量过多或场景数太大增大mipgap、压缩场景数、检查冗余约束4.2 求解器与电脑环境问题这部分属于纯环境问题但卡人的概率非常高。先说最常见的很多人装完Matlab后一运行就弹“由于找不到msvcp140.dll无法继续执行代码”这通常是Visual C运行库缺失去装对应版本的Microsoft Visual C Redistributable就能解决跟代码本身没关系。再就是Gurobi许可证过期。Gurobi的学术License一般12个月有效过期后Matlab里调用会报许可错误。解决办法是登录官网重新申请申请完刷新环境变量即可。Cplex在Matlab新版本上偶尔出现动态库加载失败优先选用Cplex官方安装包配套推荐的Matlab版本而不是随便拖一个旧版cplexlink去Matlab里硬加载。另外新版Matlab对代码诊断做得越来越强推荐打开编辑器自带的代码诊断功能“代码分析器”很多括号不匹配、变量未定义的问题在运行前就会用红线标出来这比运行后看报错效率高得多。4.3 数值稳定性与模型调试技巧做日前优化的人绕不开数值病态问题。最典型的是大M系数设置不当如果某个约束里M取得太大比如比所有变量量级大五个数量级求解器的数值稳定性会显著下降表现为同样的代码换个参数就得到完全不同的解。我的经验是M取该约束中可能出现的最大物理量的1.5到2倍即可不要盲目取一个巨大的数。如果模型报无解调试顺序非常重要。第一步把目标函数去掉只求解可行性问题确认约束本身是否自洽第二步把0-1整数变量全部固定为1或0转成一个纯LP问题排除整数组合爆炸导致的不可行第三步逐个删除约束找矛盾源。这套“三步定位法”每次都很管用比在几百行代码里瞎找要快得多。4.4 如何把你的算例替换到代码里完整版代码通常都带一个示例算例33节点或69节点替换成自己的系统时我建议按以下步骤操作把系统节点数、支路数、基准容量改成你的系统数据把线路阻抗矩阵换成你的拓扑注意辐射状配电网支路方向要和拓扑一致重新设置负荷曲线和DG接入位置别把DG功率加到纯负荷节点储能参数按你的设备实际参数填写初值务必设置合理约束里的电压上下限依据你的系统额定电压和允许偏差来定。改动过程中最隐蔽的错误是节点编号错位。很多开源代码的节点编号从1开始而且支路的首末端节点编号顺序决定了DistFlow方程里功率流方向的正负。你要是把支路首末端写反了模型也能解出结果但潮流方向会乱画出来的电压分布图会出现大量异常。还有一个小技巧把基准功率换成100MVA后欧姆值转标幺值要除以基准阻抗这个算式一定要自己核对一遍不同代码里基值定义不完全一致。算错基值会导致全部结果出现系统性的比例偏差。5. 个人实操的一点体会最后聊几句掏心窝的经验。这套代码我通读下来最有价值的不是最终的求解结果而是里面标准的建模范式——线性化DistFlow怎么处理、储能时序耦合怎么表达、两阶段耦合怎么通过大M项衔接、以及不确定性集合如何影响模型保守性。这些几乎是配电网日前调度方向的公共知识不管你以后是自己搭新模型还是看别人论文去复现这套基本功都用得上。跑通代码之后我强烈建议你做两个实验一是把分时电价从24时段改成96时段15分钟粒度观察求解时间变化和SOC曲线平滑度二是把光伏出力的不确定集范围从正负5%逐步调整到正负20%看系统总成本如何上升。做完这两个实验你对“两阶段”模型的理解会比只看代码深刻得多。毕竟代码是静态的模型对参数的反应才是活的。