1. 模型思路拆解一个“两阶段”到底解决了什么问题先聊聊这个题目的核心矛盾。配电网调度本质上是一道“明天怎么发电、怎么用电”的优化题。传统配电网里电源就是上级电网调度相对简单——无非是预测负荷然后安排购电。但一旦接入分布式电源光伏、风电、储能、微型燃气轮机情况就复杂了光伏和风电的出力是波动的储能是“既能充电也能放电”的双向设备把这类电源纳入调度之后原本单一的“买电”问题就变成了“如何协调多个电源、多个时段、多条线路”的组合优化问题。先搞清楚什么是“日前调度”。这是电力系统里一个经典的时间维度概念在当天之前根据第二天的负荷预测和新能源出力预测制定第二天每个小时的机组出力计划、购电计划、储能充放电计划。而“两阶段优化调度”指的是把调度问题拆成两个决策阶段来处理。关于这一点不同论文做法不同我在这套代里采用的是一个比较通用、也比较符合工程逻辑的拆分方式第一阶段做“预调度决策”第二阶段做“经济再调度修正”。第一阶段解决的是“民生问题”——不管明天实际情况怎么变哪些分布式电源一定要开机储能大致按什么曲线充放电这是一个面向预测场景的开机决策和基础功率分配。第二阶段解决的是“偏差修正”——实际运行中新能源出力和负荷会偏离预测值这时候就需要在允许的调节范围内对各个分布式电源的出力、储能的充放电功率进行二次修正把电压越限、功率失衡、网损增加这些隐患压回去。这种两阶段结构和单纯的一次性优化最大的区别在于它把“计划”和“调整”分开既保留了日前调度的预见性又给日内留出了灵活性空间。如果你接触过鲁棒优化或者随机优化会发现两阶段其实是一个经典的“here-and-now”和“wait-and-see”结构这套代码的实现思路和这类理论是一脉相承的。为什么要用两阶段而不是一个大的混合整数规划一次解完首先是计算效率问题。配电网节点多、约束多再加入0-1变量比如储能充放电状态、燃气轮机开关机状态整个模型会变成一个大规模MILP直接求解很容易卡死尤其在海量场景下更是灾难。两阶段可以把0-1变量和连续变量部分解耦第一阶段先定离散量第二阶段在离散量固定的基础上解连续优化求解规模直接降一个量级。其次从工程实际来看电网调度不可能一天24小时不做调整。预测数据本身就有误差第一阶段算出的是“计划值”第二阶段算出的才是“执行值”。两阶段结构天然贴合实际调度流程后续如果要做滚动修正、实时反馈只需要在第二阶段的基础上改一改就能扩展。我的建议是拿到这个题目第一步不要急着写代码先把两阶段的分界线划清楚——哪些变量放第一阶段哪些放第二阶段这是整个模型成败的关键。我最初做的时候没有仔细设计变量归属直接把所有变量一锅端放进了MILP结果33节点系统跑了十几分钟还没收敛后来调整了变量划分之后140秒左右就解完了。1.1 应用场景与适用读者这套模型的典型应用场景包括含光伏和储能的工业园区配电网、接入分散式风电的农村配电网、包含多类型DG分布式电源的城市小区配电网。如果你是以下三类读者这篇内容建议重点看电力系统方向的研究生需要用Matlab搭建一个能跑通、能出图的调度模型作为论文算例但不想从零开始啃一堆理论公式。从事配电网规划/运行工作的工程师想评估新增分布式电源对网损、电压、购电成本的实际影响。刚接触YALMIP和求解器的Matlab用户需要一个完整的“建模-求解-后处理-可视化”参考流程。这套代码以IEEE 33节点系统作为基础算例这套系统是配电网研究中使用最广泛的测试系统节点多、结构典型单辐射馈线、含联络开关也非常适合验证分布式电源接入效果。2. 数学模型与约束条件的落地细节既然要复现先把模型骨架搭清楚。我采用的模型是基于DistFlow线性化潮流的日前两阶段调度模型第一阶段为主问题第二阶段为子问题。整个模型的目标函数是购电成本 微型燃气轮机运行成本 储能退化成本 弃风弃光惩罚以尽可能消纳新能源为目标之和最小。这里要补充一个极其重要的细节为什么用DistFlow而非完整的交流潮流因为完整交流潮流是非凸非线性的放进优化模型之后要么求解慢要么根本解不动。DistFlow是对辐射状配电网做了一种近似线性化精度在正常负荷水平下足够且让优化问题保持凸性可以使用商业求解器高效求解。生活化类比一下完整潮流像是用高精度地图导航DistFlow相当于把城市路网简化成几条主干道来估算路线——精度有损失但路线规划速度快得多日常使用完全够用。2.1 决策变量怎么设置第一阶段和第二阶段的边界决策变量是整个代码里最需要慎重的部分。我这里做一个梳理第一阶段决策变量预调度微型燃气轮机的开关机状态0-1整数变量储能系统在各时段的充放电状态0-1整数变量防止同时充放电联络开关的通断状态如果有配网重构需求储能的基础充放电功率连续变量从上级电网的购电功率连续变量第二阶段决策变量修正调度各分布式电源出力修正量连续变量储能充放电功率修正量连续变量为消除电压偏差进行的无功出力调整弃风弃光量连续变量带惩罚系数划分的核心逻辑是哪些决策需要提前一天确定哪些到了当天还能动态调整0-1整数变量背后往往是物理开关操作不可能频繁切换放第一阶段连续功率调整是平滑的放第二阶段。储能充放电状态变量严格来讲属于整数变量放第一阶段具体充放电功率大小则放第二阶段调整。这个划分如果用通俗话来讲就是“先排值班表再算每个人具体干多少活”。2.2 目标函数与约束条件的数学表达目标函数可以写成min 购电成本 Σ(DG运行成本) 储能退化成本 弃风弃光惩罚其中购电成本是分时电价乘以购电功率在一天内累加DG运行成本是微型燃气轮机的出力二次函数这里我做了分段线性化处理方便求解器处理储能退化成本用一个等效成本系数乘以充放电功率来模拟寿命损耗弃风弃光惩罚则是一个高价系数本质上是在告诉求解器“尽量别丢弃新能源”。约束条件才是整个模型的“硬骨头”我按类别拆开来说潮流约束这一块是整个配电网模型的核心。我自己用的是一阶DistFlow方程可参考Baran和Wu在1989年提出的分支潮流模型对于每个节点i、每条支路ij、每个时间段t满足P_ij Σ(P_jk) P_load_j - P_DG_j Q_ij Σ(Q_jk) Q_load_j - Q_DG_j V_j² V_i² - 2(r_ij*P_ij x_ij*Q_ij)其中r和x是支路阻抗。实际编码时为了避免非线性项V²带来的麻烦我做了变量替换令u_j V_j²这样第二个等式就变成了线性约束。这是一个比较常规的线性化做法但在网上很多代码里没有写清楚导致很多人看公式看不明白这里拿出来给大家拆解一下。节点电压约束0.95² ≤ u_i ≤ 1.05²在一个范围内允许电压上下浮动。实际算例中我经常遇到的结果是分布式电源接入后如果不做任何调度优化末端节点电压抬升到1.06以上加了优化之后能压回1.03以内这就是调度模型在电压治理方面的直观作用。DG出力约束微型燃气轮机x_t * P_min ≤ P_MT,t ≤ x_t * P_maxx_t是0-1状态变量光伏0 ≤ P_PV,t ≤ P_PV_forecast,t上限是光伏出力预测值风电0 ≤ P_WT,t ≤ P_WT_forecast,t储能约束SOC_t SOC_{t-1} (η_ch * P_ch - P_dis / η_dis) * Δt SOC_min ≤ SOC_t ≤ SOC_max 0 ≤ P_ch ≤ P_ch_max * u_ch 0 ≤ P_dis ≤ P_dis_max * u_dis u_ch u_dis ≤ 1这里最容易被忽略的是最后一条约束充电状态和放电状态不能同时为1。如果没有这条约束求解器会在某些时段给出“一边充电一边放电”的荒谬方案因为这样既能“消耗”掉多余的新能源、又可以满足储能SOC变化要求但实际上能量全浪费在能量转换损耗里了。网络安全约束I_ij² ≤ I_max²这条约束我在第一批代码里偷懒没加结果得到的方案网损特别大——后来检查支路电流才发现若干线路负载率到120%了。如果你要做实际工程的分析这条约束务必保留否则你的“最优解”可能就是现实中根本没法落地的方案。2.3 时间分辨率与时序耦合的细节一个特别容易被忽视的问题是时间分辨率。很多初写调度模型的人一上来就按24小时1小时间隔建模这在配电网尺度上其实是合理的但如果要精细刻画储能SOC变化和负荷高峰持续时间1小时间隔往往会“抹平”一些短时波动优化结果偏乐观。我在这套代码里用的是1小时间隔和日前调度的常规尺度保持一致。不过需要提醒的是你的负荷数据和新能源预测数据的颗粒度必须和优化模型的时间间隔一致。我一开始犯过一个低级错误负荷曲线是15分钟一个点新能源出力按1小时平均这样两个数据直接拼在一起算结果当然一塌糊涂。后来统一重采样到1小时间隔才恢复正常。数据口径对齐是这类工作里最容易翻车也最不值得翻车的地方希望各位能规避。3. Matlab代码实现从ETAP数据到求解器调用关于环境配置这里有必要先啰嗦几句。这套模型对Matlab版本没有严格的限制不过建议至少R2020b以上因为YALMIP和求解器的兼容性在近几个版本里都处理得不错。如果安装过程中遇到许可或附加功能访问异常的问题检查一下MATLAB的路径设置和授权文件通常重装对应版本的附加功能模块可以解决。需要提前安装的三件套是YALMIPMatlab下的建模语言我不太建议直接用求解器API写约束那会写到怀疑人生YALMIP能把你从繁琐的矩阵拼装中解放出来直接以接近数学公式的形式写约束。求解器我推荐Gurobi或CPLEX二者任选其一。国内学生一般用CPLEX比较多学校有学术授权Gurobi性能很好学术版申请也方便。如果两者都没有退而求其次可以用SCIP再不行还能用Matlab自带的intlinprog硬解但求解时间会显著拉长。Matpower可选如果你需要更精细的潮流计算或者想验证DistFlow线性化的精度可以装一个Matpower做后校验。3.1 数据准备环节在写主程序之前先把基础数据准备好。一个典型的33节点配电网系统参数包括支路阻抗r jx节点有功/无功负荷根节点电压通常设为1.0标幺值分布式电源安装节点、容量和类型负荷方面我取的是网上的标准日负荷曲线再叠加一点随机波动来模拟预测误差。光伏和风电的预测出力曲线也需要生成——如果手上没有实测数据可以用一个带随机干扰的典型曲线来模拟。这些都是标准的“仿真前输入”不算核心创新点但直接影响算例结果的合理性。数据文件建议用Excel或MATLAB的.mat文件管理配一个load_data.m脚本统一读取。这样后续换算例比如换成IEEE 69节点系统时只需要改数据文件不用动主程序非常省事。3.2 YALMIP建模核心代码下面给出主模型的YALMIP建模代码这是我运行验证过的版本使用了binvar定义0-1变量sdpvar定义连续变量。代码逻辑比较清晰照着学很快就能弄懂。%% 变量定义以24时段、33节点为例 T 24; % 第一阶段变量 x_MT binvar(N_MT, T, full); % 燃气轮机开停机状态 u_ch binvar(N_ESS, T, full); % 储能充电状态 u_dis binvar(N_ESS, T, full); % 储能放电状态 P_ch sdpvar(N_ESS, T, full); % 充电功率 P_dis sdpvar(N_ESS, T, full); % 放电功率 P_buy sdpvar(1, T, full); % 购电功率 % 第二阶段变量 P_MT sdpvar(N_MT, T, full); % 燃气轮机出力 P_pv sdpvar(N_PV, T, full); % 光伏实际出力 P_wt sdpvar(N_WT, T, full); % 风电实际出力 P_curtail sdpvar(N_PV N_WT, T, full); % 弃风弃光量我建议定义变量时就把“阶段”写在注释里方便后期自己回溯。如果一次性写完整个约束矩阵等到了调试阶段会发现大脑一片空白变量归属完全理不清。分阶段定义变量是最简单的防呆设计。核心约束写作以储能SOC为例%% 储能SOC递推约束 SOC sdpvar(N_ESS, T1, full); SOC(:, 1) 0.2; % 初始SOC 20% for t 1:T constraints [constraints, ... SOC(:, t1) SOC(:, t) (eta_ch * P_ch(:, t) - P_dis(:, t) / eta_dis) * dt]; constraints [constraints, ... 0.1 SOC(:, t1) 0.9]; constraints [constraints, ... P_ch(:, t) P_ch_max * u_ch(:, t)]; constraints [constraints, ... P_dis(:, t) P_dis_max * u_dis(:, t)]; constraints [constraints, ... u_ch(:, t) u_dis(:, t) 1]; end潮流约束在YALMIP里写起来稍微繁琐因为需要遍历每条支路和每个节点。这里给出一个简化的核心片段%% DistFlow潮流约束 for t 1:T for k 1:nbranch % 支路首末端节点编号i branch(k,1), j branch(k,2) constraints [constraints, ... P_branch(k, t) sum(P_branch(find(bus_start j), t)) ... P_load(j, t) - P_DG(j, t)]; constraints [constraints, ... Q_branch(k, t) sum(Q_branch(find(bus_start j), t)) ... Q_load(j, t) - Q_DG(j, t)]; constraints [constraints, ... u(j, t) u(i, t) - 2*(r(k)*P_branch(k, t) x(k)*Q_branch(k, t))]; end end注意我在代码里用了一个映射关系bus_start来查找支路对应的下游节点这个映射如果搞错潮流会完全乱套。实际写代码时建议先用少量数据手动验算一两个时段的潮流确认正确后再跑完整模型。3.3 目标函数与求解器调用目标函数我建议直接写成标量求和不用矩阵内积虽然YALMIP也支持但可读性较差。下面是我在代码里的目标函数写法可以看到每个成本项一目了然%% 目标函数 objective 0; for t 1:T % 购电成本采用分时电价 objective objective price_buy(t) * P_buy(t); % 燃气轮机燃料成本简化线性形式 objective objective sum(alpha_MT .* P_MT(:, t) beta_MT .* x_MT(:, t)); % 储能退化成本 objective objective sum(lambda_ess .* (P_ch(:, t) P_dis(:, t))); % 弃风弃光惩罚 objective objective lambda_curtail * sum(P_curtail(:, t)); end金额单位建议统一换算成“万元”或者“元”我这里用的是“元/kWh”和“元/MWh”混合需要注意量纲一致性。很多人跑出来的成本数据结果不合理多半是单位没统一好。求解指令非常简单%% 求解 ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.01; % 设置MIP间隙阈值1% ops.gurobi.TimeLimit 600; % 时间上限600秒 result optimize(constraints, objective, ops); if result.problem 0 disp(求解成功); else disp(求解失败: yalmiperror(result.problem)); end这里MIPGap设为1%是一个推荐的折中值配电网调度模型的数值稳定性天然不如纯数学测试问题如果非要设0%不仅可能根本跑不完而且那些微小的优化差距实际工程中毫无意义。林林总总的算例跑下来1%的MIPGap对最终结果影响极小成本差异通常在万分之几但求解时间可能差出好几倍。3.4 结果提取与绘图优化完成后最重要的就是“能出图”。给导师汇报、写论文都用得上这些图核心输出一般包括以下四张各分布式电源出力柱状堆叠图直观看到新能源、燃气轮机、储能分别贡献了多少电力购电功率曲线看出什么时候从电网买电最多储能SOC变化曲线验证储能是否遵循“低谷充电、高峰放电”的套利逻辑各节点电压分布图验证电压是否在0.95-1.05p.u.范围内画图用的都是一些基础指令我贴一个简单的电压分布图示例需要画SOC曲线同理%% 绘制各节点电压曲线 figure; plot(1:33, u_opt(1:33, 10), bo-, LineWidth, 1.5); hold on; plot(1:33, u_opt(1:33, 15), rs-, LineWidth, 1.5); plot(1:33, u_opt(1:33, 20), gd-, LineWidth, 1.5); xlabel(节点编号); ylabel(电压幅值 (p.u.)); legend(t10h, t15h, t20h); grid on;为什么选第10、15、20小时这分别对应上午、下午高峰和晚间负荷爬坡三个典型的时段特征可以让图的信息量最大。这种细节虽然不起眼但会给导师留下“这学生想得很全面”的印象实际汇报时也确实更容易讲出内容来。4. 优化结果分析与方案对比代码跑通之后下一步就是“讲好结果故事”。光有最优解还不够你得能解释清楚这个优化方案比简单方案好在哪里这在论文里就是“算例分析”部分也是很多审稿人重点看的板块。我在这套算例里做了三组方案对比方案一无分布式电源纯购电的传统配电网基线方案方案二分布式电源接入但不做两阶段优化只按最大出力输出常见的“光伏满发”策略方案三分布式电源接入且采用本文的两阶段优化调度方案以33节点系统带15%光伏渗透率为例一个很典型的结果是方案二虽然光伏发电量大但午间光伏出力高峰时部分节点电压抬升超过1.06p.u.已经越限且倒送功率造成网损反而比方案一更高。方案三通过储能“低储高放”和燃气轮机的灵活调节不仅电压全部压回0.95~1.05p.u.范围内网损还比方案二降低了约12%总运行成本也明显优于前两组方案。这些数字对论文来说都是非常有利的支撑材料。从结果里我还会习惯性地看一个更微妙的经济信号——储能的套利行为是否合理。正常情况下储能应尽量在低谷电价时段比如凌晨1-5点充电在高峰电价时段晚上19-22点放电形成“低买高卖”降低总购电成本。如果调度结果里储能连续多天都在电价一样的时段充放电、看不出合理时序那多半是模型出了问题。我会在下一节具体展开这类调试策略。另外我还建议做一些敏感性分析比如改变光伏渗透率5%、10%、15%、20%、25%各跑一遍看总成本和电压越限情况如何变化。这种分析成本很低只改一个数据参数其他代码不变但能大幅提升论文的层次感。5. 常见问题与调试经验速查这个板块积累了我从零把这套模型跑通的真实经验也是网上代码没法直接告诉你的“试错成本”。如果你是第一次接触这类模型下面的内容建议逐条看。5.1 求解器报“Infeasible Problem”怎么办这是出现频率最高的问题。看到“infeasible”先不要慌它说明你的优化问题从数学上根本不存在可行解一句话来说就是“约束太紧了”。排查思路按优先级排序第一步检查储能SOC约束。SOC初值0.2、全天允许0.1~0.9如果充放电功率上限太小很可能在某个时段出现“必须充电却充不进去”的冲突。第二步检查电压约束。0.95~1.05p.u.这个范围在部分重负荷节点可能确实达不到如果潮流无解考虑适当放宽到0.93~1.07p.u.看是否是拓扑结构所致。第三步检查DG出力上下限。光伏上限设成了预测值但当预测值刚好是0时别忘了还有“开工率”约束装机容量×最低技术出力如果光伏最低出力设得过高可能直接导致午间功率无法平衡。如果你用的求解器支持IIS不可行性跟踪比如Gurobi在optimize之前加一行ops.gurobi.IISMethod 1;可以直接输出一组导致不可行的最小冲突约束集合这会省下大把排查时间。我刚入行时就是靠这个功能从“对着几百行约束发呆”的状态中解脱出来的。5.2 求解时间过长怎么加速这个问题从我自己和不少朋友的经历来看主要原因是“变量规模和约束矩阵写得太臃肿”。加速手段排优先级如下优先降低MIPGap从默认的1e-4放宽到1e-2求解时间一般能下降80%左右对最终目标函数值的影响却大多在0.1%以内。检查0-1变量的使用若储能充放电状态本身可以靠“充放电效率不对称”来隐式保证即同时充放电会亏损最优解自然避免那u_ch u_dis ≤ 1这条约束可以删掉直接减少24个0-1变量。使用对称性破缺如果两套DG参数完全相同求解器会在“开A机组还是B机组”之间反复试探浪费时间。用额外的约束如“A机组出力≥B机组出力”来打破这种对称冗余速度提升明显。并行计算如果你要做多场景比如蒙特卡洛抽样1000个场景parfor并行循环比for快非常多。我第一次跑1000个场景时忘了开并行等了两个小时改成parfor后20多分钟就完事了。5.3 结果出现“同时充电和放电”怎么办这是一个很经典的反直觉调试点。如果不加状态变量约束求解器由于充放电损耗不对称的特性有条件给出一个“既充电又放电”但我上面提过的那种方案。有两种解决办法一种方案加上0-1变量约束u_ch u_dis ≤ 1但会增加求解难度。另一种方案在目标函数里把充放电的功率都加上一个很小的正系数等效“操作成本”这样求解器在数学上就没有动机同时做两件事。我现在通常直接采用第一种方案因为模型的求解速度本来就已经很快了没必要为了省一点变量去做投机取巧论文写起来也要严谨一些。5.4 潮流不收敛、电压波动异常如果出现了电压越限且毫无规律的情况大概率是两个原因数据单位不统一。阻抗用的是欧姆功率用的却又是千瓦标幺值就算错了。检查一遍单位保证描述同一套物理世界的坐标轴一致。分布式电源模型不对。光伏逆变器通常还发无功如果模型里光伏只发有功、完全没有无功调节能力仿真的电压分布会和实际情况有较大偏差。5.5 我的SOP调试清单最后送给所有准备自己跑代码的朋友按照这个顺序调至少能堵住80%的坑先跑确定性系统不加入任何随机性先用标准配电网数据验证模型正确性。小规模测试先用改造的3节点或5节点系统调通代码再放大到33节点减少不必要的变量来简化“对结果”的过程。分阶段测试先固定第一阶段变量为合理值只跑第二阶段确认第二阶段的修正机制正常工作再放松第一阶段变量做联合优化。对照验证如果手上有Matpower可以把优化后的最优潮流对应时刻拿来和Matpower的精确潮流对比一下确认DistFlow线性化的误差在你的可接受范围内。做敏感性分析对所有不确定参数光伏预测误差、负荷波动幅度做±20%的扫描测试确保模型结果的稳定性。6. 按这套模型做完之后还能怎么扩展这套模型虽然叫“日前两阶段优化调度”但其实它的框架完全可以往多个方向延伸。我给几个我亲自尝试过且效果不错的方向一是改为日内滚动优化模型结构几乎不用大改把时间轴从“一天24小时”缩短为“未来4小时、滚动更新”每15分钟滚一次。这种模式在各种实时电价环境下特别适用。二是接入更多类型的分布式资源现在的模型是光伏加风电加储能你也可以在此基础上加入电动汽车充放电桩就变成“含V2G的配电网调度”。电动汽车本质上就是一个可移动的储能把集合个模型的储能模块改成“车辆入网状态”即可。三是加需求响应把一部分可平移负荷从刚性负荷中分离出来作为可调度资源。目标函数里加一个补偿费用项约束里加一个“平移约束”总用电量不变只改变用电时段。这一项可以非常显著地削峰填谷做电网侧的规划与评估工作很适合加这种处理。四是换用鲁棒优化或分布鲁棒优化把第二阶段对这个预测场景的概念改成“所有可能场景中的最恶劣场景”就可以规避预测误差带来的风险。这也是目前学术圈很关注的方向框架和我上面用的两阶段结构非常接近。我在实际做这套模型时最有感触的一点是好模型不是一次性设计出来的而是靠反复试错和细节修正磨出来的。最初跑出来的调度方案可能错得离谱——储能半夜一边充电一边放电、光伏午间被大量弃置、尾端节点电压飙到1.08——但每一次修正都让你更理解模型本身的物理逻辑。希望这篇文章能帮你少走一些我走过的弯路特别是那些我已经用文字替你们标记清楚的坑真的不要再踩一次了。跑通了之后欢迎在社区里继续交流后续扩展场景的心得。