1. 项目概述与核心需求解析1.1 这个模型到底在解决什么问题电力系统里面有个老生常谈但又绕不开的话题叫“无功优化”。你要是刚接触这个方向可能被各种论文里的名词绕晕有功调度、无功优化、日前计划、两阶段、分布式电源...其实把这些词拆开看就是一件事配电网里接入了越来越多的小型发电设备光伏、风机这些电网运行人员需要在“明天”到来之前提前安排好这些设备和传统配电网里各种调压手段的工作状态让电网在满足电压、潮流等安全约束的前提下运行成本最低、线损最小、电压质量最好。这里的关键词有两个“日前”和“两阶段”。日前指的是提前一天做计划把24小时有些模型取96个时段每15分钟一个点的调度策略先定下来两阶段则意味着这个优化不是一步到位的而是先做第一阶段的有功调度决定分布式电源出力、储能充放电再做第二阶段的无功优化调整无功补偿装置、逆变器无功出力、有载调压变压器分接头等两个阶段之间通过边界条件相互耦合。为什么要拆成两步因为如果把所有决策变量都塞进一个大模型里求解难度会指数级上升而且有功和无功的时间尺度、控制手段本身就不同分开处理更符合工程实际。这种模型目前在学术圈和工程圈都属于很常见但含金量不低的研究方向。如果你是电力系统方向的研究生正在找毕设题目或者小论文方向或者你在配电网规划、调度运行相关的岗位上想系统梳理一下分布式电源接入后的优化调度方法这个题目都值得好好拆解一遍。1.2 为什么选择Matlab做载体Matlab在电力系统领域的使用率极高Power System Toolbox、Matpower这些工具箱几乎是标配。就这个题目而言用Matlab搭建优化调度模型有两个天然优势一是矩阵运算和约束条件表达非常直观尤其是配合YALMIP这类建模语言写约束条件几乎和数学公式一一对应二是Matlab生态里有现成的优化求解器接口无论是开源的还是商业的都能直接调。不过我要先泼一盆冷水Matlab本身不提供“配电网优化调度”的现成工具箱你要做的是用Matlab作为平台把数学模型代码化再调用求解器去算。所以这篇文章里所有讲的“代码”本质上是数学模型在Matlab中的实现方式核心是在讲模型而不是某个现成工具。2. 两阶段优化调度模型的整体设计思路2.1 阶段划分的工程逻辑与数学逻辑先讲清楚两阶段划分的动机。配电网的调度问题其实是一个混合整数非线性规划问题——里面有连续变量各节点电压、功率有整数变量变压器分接头档位、电容器投切组数还有非线性约束潮流方程本质是非线性的。如果把所有东西全部一起优化问题规模大、求解慢而且整数变量和非线性约束混在一起很容易陷入局部最优或者干脆解不出来。所以工程上就采取“分层决策”的思路。第一阶段解决“明天每时段分布式电源出多少有功、储能充多少放多少”的问题这个阶段可以先把无功侧的整数变量松弛掉重点考虑有功功率平衡、线路传输容量、储能SOC等约束目标是最小化运行成本或者最大化可再生能源消纳。第二阶段固定第一阶段的有功决策结果在这一前提下去优化无功补偿装置投入、变压器分接头位置、逆变器无功出力目标是降低网络损耗、改善电压分布。两个阶段之间怎么衔接关键在第一阶段的边界条件要传递到第二阶段作为固定参数。比如第一阶段算出来某时段光伏有功出力是500kW第二阶段做无功优化时光伏逆变器的有功输出就固定为500kW在这个基础上再决定它的无功出力。这种串行结构的好处是每个阶段的问题规模都比原问题小得多第二阶段甚至可以独立求解计算效率大幅提升。2.2 目标函数怎么选——成本和损耗的双重考量目标函数的设计要匹配阶段目标。第一阶段目标函数最常见的选择是系统总运行成本最小化包括从上级电网购电的成本、分布式电源的发电成本如果有、储能充放电的折旧成本等。这里要特别留意分布式电源的发电成本在不同文献里处理方式差异很大有的把它设成0光伏、风机燃料成本确实为0有的会加一个运维成本系数比如每千瓦时0.02元还有的会用二次函数表示——但对一个配电网层面的日前调度模型来说用线性成本就已经足够二次函数反而会让模型更难解。第二阶段的目标函数首选网络损耗最小化。配电网的损耗计算是基于潮流结果的Ploss sum(I_ij^2 * R_ij)即所有支路电流平方乘以电阻之和。优化无功出力的本质就是改变系统无功分布从而改变电流分布和损耗。这里有个很实用的视角网络损耗是关于节点电压幅值的函数而无功补偿直接影响电压。你把损耗表达式展开会发现它和电压幅值的平方相关这是典型的凸函数在正常运行范围内所以第二阶段在保持线性约束的结构下目标函数是一个二次规划或者二阶锥规划的形式求解起来相对友好。2.3 为什么不是“一把梭”单阶段优化可能有人要问直接一个混合整数非线性规划把所有变量一次性优化掉理论上结果不是更全局最优吗这话没错但实际中不可行。在配电网规模较大比如IEEE 33节点、IEEE 123节点甚至更大、调度时段又长达24小时的情况下单阶段MINLP的求解规模会达到数百上千个整数变量加数万个连续变量现有的求解器在有限时间内很难收敛而且很容易陷入局部最优出不来。两阶段策略本质上是一种“分解协调”思想虽然牺牲了理论上的全局最优性但换来了计算可行性和工程可实施性。在实际调度场景中运行人员更需要的是一个在可接受时间内给出的、满足所有安全约束的可行解而不是一个理论上最优但算不出来的解。IEEE 33节点系统上的大量论文数据也表明两阶段模型和单阶段全局模型之间的目标函数差距通常在5%以内这个损失是可以接受的交换。3. 关键技术环节拆解配电网潮流与分布式电源建模3.1 配电网潮流计算——DistFlow模型的妙处要做优化调度首先得能描述电网的运行状态这就离不开潮流计算。传统输电网的牛顿-拉夫逊法在配电网里也能用但配电网的结构绝大多数是辐射状的闭环设计、开环运行更适合用DistFlow模型——一种专门针对辐射状配电网的潮流方程形式。DistFlow方程长这样P_ij - sum(P_jk) P_j_load - P_j_gen有功功率平衡Q_ij - sum(Q_jk) Q_j_load - Q_j_gen无功功率平衡V_j^2 V_i^2 - 2(R_ijP_ij X_ijQ_ij) (R_ij^2 X_ij^2)*I_ij^2电压降方程这套方程的精髓在于它不用去求解复杂的雅可比矩阵迭代直接用支路功率和节点电压作为变量非常适合嵌入优化模型。不过在优化模型里第三个方程里有个I_ij^2项而且电压是平方形式这会让约束变成非线性的。怎么处理两步走。第一步做变量替换用U_i替代V_i^2用L_ij替代I_ij^2这样一来约束形式就清爽很多。但还有一个问题L_ij (P_ij^2 Q_ij^2) / V_i^2这个表达式在U变量下会变成L_ij * U_i P_ij^2 Q_ij^2这仍然是一个非凸的等式约束。第二步就来了把这个等式松弛成不等式L_ij (P_ij^2 Q_ij^2) / U_i再变形为L_ij * U_i P_ij^2 Q_ij^2。这个不等式的形式恰好是一个二阶锥约束SOCP而二阶锥规划在现代求解器里非常成熟求解速度快、收敛稳定性好。这种“DistFlow 二阶锥松弛”的组合是当前配电网优化调度文献里的绝对主流做法我在实操中实测下来在IEEE 33节点系统上收敛非常稳定基本没有碰到过数值问题。3.2 分布式电源的建模细节有功出力和无功容量怎么设分布式电源里光伏和风机是最常见的。光伏出力建模相对简单可以用典型日光照强度曲线乘以装机容量再乘一个效率系数来得到有功出力的预测值。无功部分才是重头戏。现在配电网里的光伏逆变器和风机变流器通常都具备无功调节能力调节范围取决于逆变器的视在功率容量S和当前有功出力P即Q_max sqrt(S^2 - P^2)。这个约束在代码里看起来简单但它是一个非线性函数直接写进模型会让整个问题更难解。好在有工程近似手段。最常用的做法是把逆变器无功范围简化为一个固定上下限比如Q ∈ [-0.4S, 0.4S]这是基于大量实际运行数据得出的保守范围在优化模型里处理起来非常方便。如果非要精确建模就得引入分段线性化或者用旋转锥约束逼近但在我看来对于配电网优化调度这个尺度固定比例法完全够用还会让模型计算效率提升明显。另有少量文献用功率因数来约束比如要求功率因数不低于0.95这个可以换算成无功范围的限制——如果你的审稿人或者导师对这块有要求你要能灵活转换。储能系统的建模相对独立一些核心在于SOC递推方程SOC(t1) SOC(t) η_ch*P_ch(t)*Δt/E_cap - P_dis(t)Δt/(η_disE_cap)。这个约束让储能的时间耦合性很强——如果你在第一阶段不加这个约束储能就可能在同一时段同时充放电或者出现“今天放完明天充”这种不合理行为这都是初做模型时特别容易踩的坑。3.3 场景生成与不确定性考虑进阶内容再往深做一步的话光伏出力和负荷都有很强的不确定性。如果只按单场景的预测值去做日前调度到了当天实际值偏差很大决策就可能失效。更学术一点的处理方式是用场景法通过拉丁超立方采样或者蒙特卡洛模拟生成多个光伏出力和负荷场景每个场景赋一个概率把确定性问题扩展成随机优化问题。这样目标函数就变成所有场景下期望成本的最优值每个场景下都有潮流约束但决策变量只在第一阶段统一决定非预期性约束。这种做法代码里就是多跑几套数据加一个期望算子实际工作量不大但写进论文里档次明显不一样。如果你的题目里导师追求“创新点”这绝对是个投入产出比很高的方向。不过第一次做的话我的建议是先把确定性模型跑通跑透再考虑做随机规划别一口吃成胖子。4. 实操过程Matlab代码结构、YALMIP建模与求解器配置4.1 代码整体结构与数据准备整个Matlab代码的架构我按下面的模块来组织——单看层级结构就是一个标准的“数据加载→模型构建→求解→结果分析”流。%% 主程序入口run_scheduling.m clc; clear; close all; addpath(genpath(./functions)); % 添加函数库路径 %% 1. 加载配电网基础数据IEEE 33节点标准算例 [bus, branch] load_case_ieee33(); %% 2. 加载分布式电源和负荷预测曲线24时段数据 [load_profile, pv_profile] load_prediction_data(typical_day.csv); %% 3. 构建并求解第一阶段日前有功调度模型 [x_phase1, info_phase1] solve_phase1(bus, branch, load_profile, pv_profile); disp([第一阶段求解完成总运行成本 , num2str(info_phase1.obj), 元]); %% 4. 构建并求解第二阶段无功优化模型固定第一阶段有功 [x_phase2, info_phase2] solve_phase2(bus, branch, load_profile, pv_profile, x_phase1); disp([第二阶段求解完成网络损耗 , num2str(info_phase2.obj), kW]); %% 5. 结果可视化 plot_results(bus, branch, x_phase1, x_phase2);数据准备是每天最容易忽略但最容易出错的地方。IEEE 33节点系统的数据网上到处都有但你要注意单位统一阻抗单位是欧姆功率基准值取多少电压基准值取10kV还是12.66kV这些都要提前定清楚不然算出来的结果五花八门你根本没法判断谁对谁错。我的习惯是统一用标幺值pu功率基准值取1MVA电压基准值取12.66kV这样数值都在0.9到1.1之间求解器的数值稳定性最好。4.2 核心建模代码YALMIP让数学公式直接可读YALMIP的使用我觉得是让Matlab代码质量产生质的飞跃的关键工具。以前用纯Matlab写优化模型得把所有约束展开成矩阵形式那简直是一场灾难。现在用YALMIP做建模语言写出来的代码跟论文里的数学表达式几乎一一对应。下面给一个完整的第一阶段核心建模片段。function [x, info] solve_phase1(bus, branch, load_profile, pv_profile) %% 读网络参数 nb size(bus, 1); % 节点数 33 nl size(branch, 1); % 支路数 32 nt 24; % 时段数 24 S_base 1e6; % 功率基准值 1MVA V_base 12.66e3; % 电压基准值 12.66kV %% 定义决策变量 P_buy sdpvar(1, nt, full); % 上级电网购入有功 P_pv sdpvar(nb, nt, full); % 各节点光伏有功注入 P_ch sdpvar(nb, nt, full); % 储能充电功率 P_dis sdpvar(nb, nt, full); % 储能放电功率 SOC sdpvar(nb, nt, full); % 储能荷电状态 P_br sdpvar(nl, nt, full); % 支路有功功率 U sdpvar(nb, nt, full); % 节点电压幅值平方标幺值 %% 添加约束条件 constraints []; for t 1:nt for k 1:nl i branch(k, 1); j branch(k, 2); R_ij branch(k, 3) / (V_base^2 / S_base); % 电阻标幺值 X_ij branch(k, 4) / (V_base^2 / S_base); % 电抗标幺值 constraints [constraints, ... P_br(k, t) P_buy(t) * (i 1), ... % 根节点功率注入 P_br(k, t) -5/S_base, ... P_br(k, t) 5/S_base]; % 支路传输容量限值 end constraints [constraints, ... sum(P_br(:, t)) sum(load_profile(:, t)) / S_base ... % 有功平衡 - sum(P_pv(:, t)) - sum(P_dis(:, t)) sum(P_ch(:, t))]; constraints [constraints, ... 0 P_pv(:, t) pv_profile(:, t) / S_base, ... % 光伏出力上限 0 P_ch(:, t) 0.2/S_base, ... 0 P_dis(:, t) 0.2/S_base, ... SOC(:, t) 0.2, SOC(:, t) 0.9]; end % SOC递推约束储能时间耦合 for k 1:nb for t 2:nt constraints [constraints, ... SOC(k, t) SOC(k, t-1) 0.9*P_ch(k, t)*1/0.5 - P_dis(k, t)*1/(0.9*0.5)]; end end %% 目标函数购电成本 储能损耗成本 弃光惩罚 c_buy 0.5; % 购电电价 0.5元/kWh c_pv 0.02; % 光伏运维成本 c_pen 0.8; % 弃光惩罚系数 objective sum(c_buy * P_buy * 1e3) sum(c_pv * P_pv * 1e3) ... sum(c_pen * (pv_profile/S_base - P_pv) * 1e3); %% 配置求解器并求解 ops sdpsettings(solver, cplex, verbose, 0, showprogress, 0); optimize(constraints, objective, ops); %% 结果回传 x.P_buy value(P_buy) * S_base / 1e3; % 换算成kW x.P_pv value(P_pv) * S_base / 1e3; x.P_ch value(P_ch) * S_base / 1e3; x.P_dis value(P_dis) * S_base / 1e3; x.SOC value(SOC); info.obj value(objective) * S_base / 1e3; info.status yalmiperror(optimize(constraints, objective, ops)); end注意看SOC递推约束我加了储能效率充电效率0.9、放电效率0.9这个系数会直接影响储能一天下来“充进去多少、放出来多少”的能量守恒关系。不加效率的话储能就变成了永动机这是很多入门级代码会犯的明显错误。第二阶段的YALMIP建模思路跟第一阶段几乎一样区别在于一是有功变量P_pv要固定为第一阶段算出来的值二是决策变量换成无功补偿装置的无功出力Q_c、逆变器无功Q_inv、变压器分接头Tap三是目标函数从“成本最小”换成“网损最小”。核心约束里要加上潮流方程DistFlow在YALMIP里可以直接用二阶锥约束表达。% 第二阶段核心片段DistFlow二阶锥约束 constraints [constraints, ... U(j,t) U(i,t) - 2*(R_ij*P_br(k,t) X_ij*Q_br(k,t)) (R_ij^2 X_ij^2)*L_br(k,t)]; constraints [constraints, ... L_br(k,t) (P_br(k,t)^2 Q_br(k,t)^2) / U(i,t)]; % SOCP松弛这里的U、L、P_br、Q_br都是二阶锥变量YALMIP会自动识别并传给求解器代码层面的改动成本很低但模型性质完全变了——从MINLP变成了SOCP求解难度下降了一个量级。4.3 求解器选型Cplex、Gurobi还是开源方案做了这么多年优化模型求解器算是我最看重的环节。商业求解器方面Cplex和Gurobi是主流选择两个都支持大规模线性规划、混合整数规划、二次约束规划含二阶锥性能在伯仲之间。我个人实际用下来Gurobi在大规模二阶锥问题上经常比Cplex快10%到20%而且Gurobi的学术授权申请特别方便用学校邮箱就能申请到一年期的免费license对在校学生非常友好。如果暂时没有商业求解器license也不建议放弃这个方向。开源的CBCCOIN-OR分支切制求解器能处理LP和MILP配合YALMIP可以直接用。SCS和ECOS这两个开源求解器则专门针对锥规划在二阶锥问题上表现不错虽然跟商业求解器比在大型算例上还有差距但如果你的算例只是IEEE 33节点这个规模SCS完全够用求解时间通常在几秒到十几秒之间。安装求解器的建议先在Matlab里运行yalmiptest命令它会自动检测当前环境里已经安装了哪些求解器然后你再根据需要去下载对应的Matlab接口包。这里有个常见的坑——下载完Gurobi之后光有安装包不行还要在Matlab里运行gurobi_setup进行路径配置否则YALMIP会报“solver not found”。5. 实操经验总结参数整定、结果分析与避坑指南5.1 参数怎么调才合理——从一次失败案例说起我印象很深的一次经历是刚开始搭这个模型的时候求解器怎么都不收敛要么迭代半天卡在某个地方要么直接报infeasible。折腾了大半天最后发现问题是节点电压约束范围设得太死——我把电压下限设成了0.95标幺值但第一阶段的模型里根本没加任何电压约束导致光伏大发时段部分节点电压被算得离0.95很远第一阶段决策出来的有功值到了第二阶段在电压约束下根本无解。后来我把电压约束的处理方式改成两阶段都加第一阶段加上松弛版的电压范围约束比如0.93到1.07让有功调度“知道”有电压安全底线这回事第二阶段再收紧到0.95到1.05做精细的无功调节。这样两阶段之间设备的边界条件就一致了模型再也没有出现过不可行的情况。这个案例说明一个道理两阶段调度的“通信”非常关键。第一阶段不能只盯着经济性完全忽略电网安全约束第二阶段也不能在固定了有功之后发现电压稳不住就束手无策。合理的做法是第一阶段就给第二阶段的可行域留足裕量。5.2 结果分析怎么看——网损、电压分布与DG消纳求解完成后很多人习惯直接看总目标函数值就完事了这不行。优化调度模型算出来的海量结果信息才是评估方案好坏的真正钥匙。我建议至少画三张图第一张是24时段的有功调度曲线图横轴是时间纵轴是各类电源出力和负荷曲线。从这张图能直观看到光伏出力高峰期是否被合理消纳储能是否在低谷时段充电、高峰时段放电。如果发现光伏出力大中午的被大量削减而储能却明明还有容量没用说明目标函数里的弃光惩罚系数设得不够要么调高惩罚要么检查储能容量约束。第二张是各节点电压幅值分布图用semilogy或者plot画出来时序曲线。正常情况下所有节点电压应该在0.95到1.05之间波动分布式电源接入点的电压在正午时段会抬升——这是高渗透率配电网的典型特征。如果某几个节点电压持续偏高但无功补偿装置没有任何动作说明第二阶段的目标函数权重设置有问题或者无功补偿装置的容量配置不足。第三张是网损的时段分布图。无功优化后的网损曲线应该明显低于未优化前的方案尤其在重负荷时段。用这两条曲线的差就能算出一天节省的电量换算成电费这是你论文里写“节能减排效益”的数据来源。5.3 代码调试的常见坑与排查技巧要我说这个模型能踩的坑基本都在下面这几个地方。遇到问题了先别急着改代码照着清单逐项排一遍。不可行错误Infeasible problem——这个最让人头大一出现就说明约束之间打架了。优先检查三个地方一是两阶段之间传递的边界变量第一阶段的有功出力是不是超出了第二阶段允许的范围二是储能SOC的初始值和终止值是否设置得一致我习惯设SOC(1)0.5、SOC(24)也限制在0.5左右避免储能“白嫖”能量三是光伏出力上限是否跟实际的容量数据一致。收敛慢或者求解时间异常长——大概率是模型中隐藏的非线性约束把问题的凸性破坏了。YALMIP里可以用check(constraints)命令快速检查每个约束的类型看看有没有意外的非凸等式约束。还有一种可能是约束里有变量和变量相乘比如P*Q这种双线性项会让求解器直接卡死一定要用变量替换或者松弛把它消掉。Numerical issue或者数值溢出告警——这一般是单位没有统一造成的。我见过有人把电阻值直接填成欧姆量级比如0.5欧姆但功率却是千瓦量级在标幺化之后数值差了百倍量级求解器内部的容差设置可能就直接崩溃。解决办法是严格按照4.1节说的做标幺化处理并且给所有约束设置合理的显式上下界比如U的范围设在0.8到1.2L的范围设在0到1帮求解器把搜索空间压缩到合理范围内。YALMIP版本兼容性问题——如果你用的是老版本的YALMIP2020年之前的版本某些二阶锥约束的传递可能有问题报一些莫名其妙的错。此外如果你用到的是R2023a以上的较新Matlab版本我建议把YALMIP更新到GitHub上的最新版本业界更新速度比Matlab新版本发布的速度还活跃兼容性更好。6. 扩展方向从两阶段到多阶段、从确定到鲁棒两阶段日前调度做熟练之后往哪走我这里推荐三个方向难度递增、发表潜力也递增。第一个是日内滚动修正。日前计划做得再精细到了当天总会有偏差所以工程上更常见的是“日前计划日内调整”的模式。日内阶段每15分钟滚动执行一次基于最新的实时量测数据和超短期预测结果只修正未来1到4小时的调度指令——这个在代码层面做得更省力一点的话就需要将模型从离线转成循环调用每轮都重新求解但建模内核跟两阶段模型是一致的。第二个是考虑不确定性的随机优化或者鲁棒优化。随机优化需要生成大量场景代码量适中但场景数多了求解时间会直线上涨鲁棒优化则用不确定集合来描述预测误差不需要场景概率模型复杂度略高但解出来的方案保守性更强。两阶段模型本身就是一个天然的鲁棒优化结构把第二阶段写成“在最坏场景下调整无功手段”就成了经典的两阶段鲁棒优化这个方向近年在顶级期刊上的热度一直很高如果你手里已经有两阶段的确定性代码往鲁棒方向扩展的路径相当顺畅。第三个是考虑网络重构的联合优化。配电网里除了无功补偿还有联络开关可以进行网络拓扑调整这在故障恢复和降损场景下效果显著。但一旦引入开关状态变量问题就变成混合整数二阶锥规划MISOCP求解难度又上一个台阶——不过这恰恰是很多论文能发得出来的原因算力不够、方法来凑分解算法、Benders分解、拉格朗日松弛这些高级解法都能派上用场。我个人实际做下来比较推荐的进阶路线是先把确定性两阶段模型做扎实跑通IEEE 33节点和IEEE 123节点两个标准算例然后加一个日内滚动模块模拟全天运行最后再上一套场景法做不确定性分析。这三步走完别说硕士毕业论文就是一篇还不错的期刊小论文的主体内容也都能拿得出手了。 基于两阶段优化的配电网日前调度与无功优化模型详解（Matlab实现）
1. 项目概述与核心需求解析1.1 这个模型到底在解决什么问题电力系统里面有个老生常谈但又绕不开的话题叫“无功优化”。你要是刚接触这个方向可能被各种论文里的名词绕晕有功调度、无功优化、日前计划、两阶段、分布式电源...其实把这些词拆开看就是一件事配电网里接入了越来越多的小型发电设备光伏、风机这些电网运行人员需要在“明天”到来之前提前安排好这些设备和传统配电网里各种调压手段的工作状态让电网在满足电压、潮流等安全约束的前提下运行成本最低、线损最小、电压质量最好。这里的关键词有两个“日前”和“两阶段”。日前指的是提前一天做计划把24小时有些模型取96个时段每15分钟一个点的调度策略先定下来两阶段则意味着这个优化不是一步到位的而是先做第一阶段的有功调度决定分布式电源出力、储能充放电再做第二阶段的无功优化调整无功补偿装置、逆变器无功出力、有载调压变压器分接头等两个阶段之间通过边界条件相互耦合。为什么要拆成两步因为如果把所有决策变量都塞进一个大模型里求解难度会指数级上升而且有功和无功的时间尺度、控制手段本身就不同分开处理更符合工程实际。这种模型目前在学术圈和工程圈都属于很常见但含金量不低的研究方向。如果你是电力系统方向的研究生正在找毕设题目或者小论文方向或者你在配电网规划、调度运行相关的岗位上想系统梳理一下分布式电源接入后的优化调度方法这个题目都值得好好拆解一遍。1.2 为什么选择Matlab做载体Matlab在电力系统领域的使用率极高Power System Toolbox、Matpower这些工具箱几乎是标配。就这个题目而言用Matlab搭建优化调度模型有两个天然优势一是矩阵运算和约束条件表达非常直观尤其是配合YALMIP这类建模语言写约束条件几乎和数学公式一一对应二是Matlab生态里有现成的优化求解器接口无论是开源的还是商业的都能直接调。不过我要先泼一盆冷水Matlab本身不提供“配电网优化调度”的现成工具箱你要做的是用Matlab作为平台把数学模型代码化再调用求解器去算。所以这篇文章里所有讲的“代码”本质上是数学模型在Matlab中的实现方式核心是在讲模型而不是某个现成工具。2. 两阶段优化调度模型的整体设计思路2.1 阶段划分的工程逻辑与数学逻辑先讲清楚两阶段划分的动机。配电网的调度问题其实是一个混合整数非线性规划问题——里面有连续变量各节点电压、功率有整数变量变压器分接头档位、电容器投切组数还有非线性约束潮流方程本质是非线性的。如果把所有东西全部一起优化问题规模大、求解慢而且整数变量和非线性约束混在一起很容易陷入局部最优或者干脆解不出来。所以工程上就采取“分层决策”的思路。第一阶段解决“明天每时段分布式电源出多少有功、储能充多少放多少”的问题这个阶段可以先把无功侧的整数变量松弛掉重点考虑有功功率平衡、线路传输容量、储能SOC等约束目标是最小化运行成本或者最大化可再生能源消纳。第二阶段固定第一阶段的有功决策结果在这一前提下去优化无功补偿装置投入、变压器分接头位置、逆变器无功出力目标是降低网络损耗、改善电压分布。两个阶段之间怎么衔接关键在第一阶段的边界条件要传递到第二阶段作为固定参数。比如第一阶段算出来某时段光伏有功出力是500kW第二阶段做无功优化时光伏逆变器的有功输出就固定为500kW在这个基础上再决定它的无功出力。这种串行结构的好处是每个阶段的问题规模都比原问题小得多第二阶段甚至可以独立求解计算效率大幅提升。2.2 目标函数怎么选——成本和损耗的双重考量目标函数的设计要匹配阶段目标。第一阶段目标函数最常见的选择是系统总运行成本最小化包括从上级电网购电的成本、分布式电源的发电成本如果有、储能充放电的折旧成本等。这里要特别留意分布式电源的发电成本在不同文献里处理方式差异很大有的把它设成0光伏、风机燃料成本确实为0有的会加一个运维成本系数比如每千瓦时0.02元还有的会用二次函数表示——但对一个配电网层面的日前调度模型来说用线性成本就已经足够二次函数反而会让模型更难解。第二阶段的目标函数首选网络损耗最小化。配电网的损耗计算是基于潮流结果的Ploss sum(I_ij^2 * R_ij)即所有支路电流平方乘以电阻之和。优化无功出力的本质就是改变系统无功分布从而改变电流分布和损耗。这里有个很实用的视角网络损耗是关于节点电压幅值的函数而无功补偿直接影响电压。你把损耗表达式展开会发现它和电压幅值的平方相关这是典型的凸函数在正常运行范围内所以第二阶段在保持线性约束的结构下目标函数是一个二次规划或者二阶锥规划的形式求解起来相对友好。2.3 为什么不是“一把梭”单阶段优化可能有人要问直接一个混合整数非线性规划把所有变量一次性优化掉理论上结果不是更全局最优吗这话没错但实际中不可行。在配电网规模较大比如IEEE 33节点、IEEE 123节点甚至更大、调度时段又长达24小时的情况下单阶段MINLP的求解规模会达到数百上千个整数变量加数万个连续变量现有的求解器在有限时间内很难收敛而且很容易陷入局部最优出不来。两阶段策略本质上是一种“分解协调”思想虽然牺牲了理论上的全局最优性但换来了计算可行性和工程可实施性。在实际调度场景中运行人员更需要的是一个在可接受时间内给出的、满足所有安全约束的可行解而不是一个理论上最优但算不出来的解。IEEE 33节点系统上的大量论文数据也表明两阶段模型和单阶段全局模型之间的目标函数差距通常在5%以内这个损失是可以接受的交换。3. 关键技术环节拆解配电网潮流与分布式电源建模3.1 配电网潮流计算——DistFlow模型的妙处要做优化调度首先得能描述电网的运行状态这就离不开潮流计算。传统输电网的牛顿-拉夫逊法在配电网里也能用但配电网的结构绝大多数是辐射状的闭环设计、开环运行更适合用DistFlow模型——一种专门针对辐射状配电网的潮流方程形式。DistFlow方程长这样P_ij - sum(P_jk) P_j_load - P_j_gen有功功率平衡Q_ij - sum(Q_jk) Q_j_load - Q_j_gen无功功率平衡V_j^2 V_i^2 - 2(R_ijP_ij X_ijQ_ij) (R_ij^2 X_ij^2)*I_ij^2电压降方程这套方程的精髓在于它不用去求解复杂的雅可比矩阵迭代直接用支路功率和节点电压作为变量非常适合嵌入优化模型。不过在优化模型里第三个方程里有个I_ij^2项而且电压是平方形式这会让约束变成非线性的。怎么处理两步走。第一步做变量替换用U_i替代V_i^2用L_ij替代I_ij^2这样一来约束形式就清爽很多。但还有一个问题L_ij (P_ij^2 Q_ij^2) / V_i^2这个表达式在U变量下会变成L_ij * U_i P_ij^2 Q_ij^2这仍然是一个非凸的等式约束。第二步就来了把这个等式松弛成不等式L_ij (P_ij^2 Q_ij^2) / U_i再变形为L_ij * U_i P_ij^2 Q_ij^2。这个不等式的形式恰好是一个二阶锥约束SOCP而二阶锥规划在现代求解器里非常成熟求解速度快、收敛稳定性好。这种“DistFlow 二阶锥松弛”的组合是当前配电网优化调度文献里的绝对主流做法我在实操中实测下来在IEEE 33节点系统上收敛非常稳定基本没有碰到过数值问题。3.2 分布式电源的建模细节有功出力和无功容量怎么设分布式电源里光伏和风机是最常见的。光伏出力建模相对简单可以用典型日光照强度曲线乘以装机容量再乘一个效率系数来得到有功出力的预测值。无功部分才是重头戏。现在配电网里的光伏逆变器和风机变流器通常都具备无功调节能力调节范围取决于逆变器的视在功率容量S和当前有功出力P即Q_max sqrt(S^2 - P^2)。这个约束在代码里看起来简单但它是一个非线性函数直接写进模型会让整个问题更难解。好在有工程近似手段。最常用的做法是把逆变器无功范围简化为一个固定上下限比如Q ∈ [-0.4S, 0.4S]这是基于大量实际运行数据得出的保守范围在优化模型里处理起来非常方便。如果非要精确建模就得引入分段线性化或者用旋转锥约束逼近但在我看来对于配电网优化调度这个尺度固定比例法完全够用还会让模型计算效率提升明显。另有少量文献用功率因数来约束比如要求功率因数不低于0.95这个可以换算成无功范围的限制——如果你的审稿人或者导师对这块有要求你要能灵活转换。储能系统的建模相对独立一些核心在于SOC递推方程SOC(t1) SOC(t) η_ch*P_ch(t)*Δt/E_cap - P_dis(t)Δt/(η_disE_cap)。这个约束让储能的时间耦合性很强——如果你在第一阶段不加这个约束储能就可能在同一时段同时充放电或者出现“今天放完明天充”这种不合理行为这都是初做模型时特别容易踩的坑。3.3 场景生成与不确定性考虑进阶内容再往深做一步的话光伏出力和负荷都有很强的不确定性。如果只按单场景的预测值去做日前调度到了当天实际值偏差很大决策就可能失效。更学术一点的处理方式是用场景法通过拉丁超立方采样或者蒙特卡洛模拟生成多个光伏出力和负荷场景每个场景赋一个概率把确定性问题扩展成随机优化问题。这样目标函数就变成所有场景下期望成本的最优值每个场景下都有潮流约束但决策变量只在第一阶段统一决定非预期性约束。这种做法代码里就是多跑几套数据加一个期望算子实际工作量不大但写进论文里档次明显不一样。如果你的题目里导师追求“创新点”这绝对是个投入产出比很高的方向。不过第一次做的话我的建议是先把确定性模型跑通跑透再考虑做随机规划别一口吃成胖子。4. 实操过程Matlab代码结构、YALMIP建模与求解器配置4.1 代码整体结构与数据准备整个Matlab代码的架构我按下面的模块来组织——单看层级结构就是一个标准的“数据加载→模型构建→求解→结果分析”流。%% 主程序入口run_scheduling.m clc; clear; close all; addpath(genpath(./functions)); % 添加函数库路径 %% 1. 加载配电网基础数据IEEE 33节点标准算例 [bus, branch] load_case_ieee33(); %% 2. 加载分布式电源和负荷预测曲线24时段数据 [load_profile, pv_profile] load_prediction_data(typical_day.csv); %% 3. 构建并求解第一阶段日前有功调度模型 [x_phase1, info_phase1] solve_phase1(bus, branch, load_profile, pv_profile); disp([第一阶段求解完成总运行成本 , num2str(info_phase1.obj), 元]); %% 4. 构建并求解第二阶段无功优化模型固定第一阶段有功 [x_phase2, info_phase2] solve_phase2(bus, branch, load_profile, pv_profile, x_phase1); disp([第二阶段求解完成网络损耗 , num2str(info_phase2.obj), kW]); %% 5. 结果可视化 plot_results(bus, branch, x_phase1, x_phase2);数据准备是每天最容易忽略但最容易出错的地方。IEEE 33节点系统的数据网上到处都有但你要注意单位统一阻抗单位是欧姆功率基准值取多少电压基准值取10kV还是12.66kV这些都要提前定清楚不然算出来的结果五花八门你根本没法判断谁对谁错。我的习惯是统一用标幺值pu功率基准值取1MVA电压基准值取12.66kV这样数值都在0.9到1.1之间求解器的数值稳定性最好。4.2 核心建模代码YALMIP让数学公式直接可读YALMIP的使用我觉得是让Matlab代码质量产生质的飞跃的关键工具。以前用纯Matlab写优化模型得把所有约束展开成矩阵形式那简直是一场灾难。现在用YALMIP做建模语言写出来的代码跟论文里的数学表达式几乎一一对应。下面给一个完整的第一阶段核心建模片段。function [x, info] solve_phase1(bus, branch, load_profile, pv_profile) %% 读网络参数 nb size(bus, 1); % 节点数 33 nl size(branch, 1); % 支路数 32 nt 24; % 时段数 24 S_base 1e6; % 功率基准值 1MVA V_base 12.66e3; % 电压基准值 12.66kV %% 定义决策变量 P_buy sdpvar(1, nt, full); % 上级电网购入有功 P_pv sdpvar(nb, nt, full); % 各节点光伏有功注入 P_ch sdpvar(nb, nt, full); % 储能充电功率 P_dis sdpvar(nb, nt, full); % 储能放电功率 SOC sdpvar(nb, nt, full); % 储能荷电状态 P_br sdpvar(nl, nt, full); % 支路有功功率 U sdpvar(nb, nt, full); % 节点电压幅值平方标幺值 %% 添加约束条件 constraints []; for t 1:nt for k 1:nl i branch(k, 1); j branch(k, 2); R_ij branch(k, 3) / (V_base^2 / S_base); % 电阻标幺值 X_ij branch(k, 4) / (V_base^2 / S_base); % 电抗标幺值 constraints [constraints, ... P_br(k, t) P_buy(t) * (i 1), ... % 根节点功率注入 P_br(k, t) -5/S_base, ... P_br(k, t) 5/S_base]; % 支路传输容量限值 end constraints [constraints, ... sum(P_br(:, t)) sum(load_profile(:, t)) / S_base ... % 有功平衡 - sum(P_pv(:, t)) - sum(P_dis(:, t)) sum(P_ch(:, t))]; constraints [constraints, ... 0 P_pv(:, t) pv_profile(:, t) / S_base, ... % 光伏出力上限 0 P_ch(:, t) 0.2/S_base, ... 0 P_dis(:, t) 0.2/S_base, ... SOC(:, t) 0.2, SOC(:, t) 0.9]; end % SOC递推约束储能时间耦合 for k 1:nb for t 2:nt constraints [constraints, ... SOC(k, t) SOC(k, t-1) 0.9*P_ch(k, t)*1/0.5 - P_dis(k, t)*1/(0.9*0.5)]; end end %% 目标函数购电成本 储能损耗成本 弃光惩罚 c_buy 0.5; % 购电电价 0.5元/kWh c_pv 0.02; % 光伏运维成本 c_pen 0.8; % 弃光惩罚系数 objective sum(c_buy * P_buy * 1e3) sum(c_pv * P_pv * 1e3) ... sum(c_pen * (pv_profile/S_base - P_pv) * 1e3); %% 配置求解器并求解 ops sdpsettings(solver, cplex, verbose, 0, showprogress, 0); optimize(constraints, objective, ops); %% 结果回传 x.P_buy value(P_buy) * S_base / 1e3; % 换算成kW x.P_pv value(P_pv) * S_base / 1e3; x.P_ch value(P_ch) * S_base / 1e3; x.P_dis value(P_dis) * S_base / 1e3; x.SOC value(SOC); info.obj value(objective) * S_base / 1e3; info.status yalmiperror(optimize(constraints, objective, ops)); end注意看SOC递推约束我加了储能效率充电效率0.9、放电效率0.9这个系数会直接影响储能一天下来“充进去多少、放出来多少”的能量守恒关系。不加效率的话储能就变成了永动机这是很多入门级代码会犯的明显错误。第二阶段的YALMIP建模思路跟第一阶段几乎一样区别在于一是有功变量P_pv要固定为第一阶段算出来的值二是决策变量换成无功补偿装置的无功出力Q_c、逆变器无功Q_inv、变压器分接头Tap三是目标函数从“成本最小”换成“网损最小”。核心约束里要加上潮流方程DistFlow在YALMIP里可以直接用二阶锥约束表达。% 第二阶段核心片段DistFlow二阶锥约束 constraints [constraints, ... U(j,t) U(i,t) - 2*(R_ij*P_br(k,t) X_ij*Q_br(k,t)) (R_ij^2 X_ij^2)*L_br(k,t)]; constraints [constraints, ... L_br(k,t) (P_br(k,t)^2 Q_br(k,t)^2) / U(i,t)]; % SOCP松弛这里的U、L、P_br、Q_br都是二阶锥变量YALMIP会自动识别并传给求解器代码层面的改动成本很低但模型性质完全变了——从MINLP变成了SOCP求解难度下降了一个量级。4.3 求解器选型Cplex、Gurobi还是开源方案做了这么多年优化模型求解器算是我最看重的环节。商业求解器方面Cplex和Gurobi是主流选择两个都支持大规模线性规划、混合整数规划、二次约束规划含二阶锥性能在伯仲之间。我个人实际用下来Gurobi在大规模二阶锥问题上经常比Cplex快10%到20%而且Gurobi的学术授权申请特别方便用学校邮箱就能申请到一年期的免费license对在校学生非常友好。如果暂时没有商业求解器license也不建议放弃这个方向。开源的CBCCOIN-OR分支切制求解器能处理LP和MILP配合YALMIP可以直接用。SCS和ECOS这两个开源求解器则专门针对锥规划在二阶锥问题上表现不错虽然跟商业求解器比在大型算例上还有差距但如果你的算例只是IEEE 33节点这个规模SCS完全够用求解时间通常在几秒到十几秒之间。安装求解器的建议先在Matlab里运行yalmiptest命令它会自动检测当前环境里已经安装了哪些求解器然后你再根据需要去下载对应的Matlab接口包。这里有个常见的坑——下载完Gurobi之后光有安装包不行还要在Matlab里运行gurobi_setup进行路径配置否则YALMIP会报“solver not found”。5. 实操经验总结参数整定、结果分析与避坑指南5.1 参数怎么调才合理——从一次失败案例说起我印象很深的一次经历是刚开始搭这个模型的时候求解器怎么都不收敛要么迭代半天卡在某个地方要么直接报infeasible。折腾了大半天最后发现问题是节点电压约束范围设得太死——我把电压下限设成了0.95标幺值但第一阶段的模型里根本没加任何电压约束导致光伏大发时段部分节点电压被算得离0.95很远第一阶段决策出来的有功值到了第二阶段在电压约束下根本无解。后来我把电压约束的处理方式改成两阶段都加第一阶段加上松弛版的电压范围约束比如0.93到1.07让有功调度“知道”有电压安全底线这回事第二阶段再收紧到0.95到1.05做精细的无功调节。这样两阶段之间设备的边界条件就一致了模型再也没有出现过不可行的情况。这个案例说明一个道理两阶段调度的“通信”非常关键。第一阶段不能只盯着经济性完全忽略电网安全约束第二阶段也不能在固定了有功之后发现电压稳不住就束手无策。合理的做法是第一阶段就给第二阶段的可行域留足裕量。5.2 结果分析怎么看——网损、电压分布与DG消纳求解完成后很多人习惯直接看总目标函数值就完事了这不行。优化调度模型算出来的海量结果信息才是评估方案好坏的真正钥匙。我建议至少画三张图第一张是24时段的有功调度曲线图横轴是时间纵轴是各类电源出力和负荷曲线。从这张图能直观看到光伏出力高峰期是否被合理消纳储能是否在低谷时段充电、高峰时段放电。如果发现光伏出力大中午的被大量削减而储能却明明还有容量没用说明目标函数里的弃光惩罚系数设得不够要么调高惩罚要么检查储能容量约束。第二张是各节点电压幅值分布图用semilogy或者plot画出来时序曲线。正常情况下所有节点电压应该在0.95到1.05之间波动分布式电源接入点的电压在正午时段会抬升——这是高渗透率配电网的典型特征。如果某几个节点电压持续偏高但无功补偿装置没有任何动作说明第二阶段的目标函数权重设置有问题或者无功补偿装置的容量配置不足。第三张是网损的时段分布图。无功优化后的网损曲线应该明显低于未优化前的方案尤其在重负荷时段。用这两条曲线的差就能算出一天节省的电量换算成电费这是你论文里写“节能减排效益”的数据来源。5.3 代码调试的常见坑与排查技巧要我说这个模型能踩的坑基本都在下面这几个地方。遇到问题了先别急着改代码照着清单逐项排一遍。不可行错误Infeasible problem——这个最让人头大一出现就说明约束之间打架了。优先检查三个地方一是两阶段之间传递的边界变量第一阶段的有功出力是不是超出了第二阶段允许的范围二是储能SOC的初始值和终止值是否设置得一致我习惯设SOC(1)0.5、SOC(24)也限制在0.5左右避免储能“白嫖”能量三是光伏出力上限是否跟实际的容量数据一致。收敛慢或者求解时间异常长——大概率是模型中隐藏的非线性约束把问题的凸性破坏了。YALMIP里可以用check(constraints)命令快速检查每个约束的类型看看有没有意外的非凸等式约束。还有一种可能是约束里有变量和变量相乘比如P*Q这种双线性项会让求解器直接卡死一定要用变量替换或者松弛把它消掉。Numerical issue或者数值溢出告警——这一般是单位没有统一造成的。我见过有人把电阻值直接填成欧姆量级比如0.5欧姆但功率却是千瓦量级在标幺化之后数值差了百倍量级求解器内部的容差设置可能就直接崩溃。解决办法是严格按照4.1节说的做标幺化处理并且给所有约束设置合理的显式上下界比如U的范围设在0.8到1.2L的范围设在0到1帮求解器把搜索空间压缩到合理范围内。YALMIP版本兼容性问题——如果你用的是老版本的YALMIP2020年之前的版本某些二阶锥约束的传递可能有问题报一些莫名其妙的错。此外如果你用到的是R2023a以上的较新Matlab版本我建议把YALMIP更新到GitHub上的最新版本业界更新速度比Matlab新版本发布的速度还活跃兼容性更好。6. 扩展方向从两阶段到多阶段、从确定到鲁棒两阶段日前调度做熟练之后往哪走我这里推荐三个方向难度递增、发表潜力也递增。第一个是日内滚动修正。日前计划做得再精细到了当天总会有偏差所以工程上更常见的是“日前计划日内调整”的模式。日内阶段每15分钟滚动执行一次基于最新的实时量测数据和超短期预测结果只修正未来1到4小时的调度指令——这个在代码层面做得更省力一点的话就需要将模型从离线转成循环调用每轮都重新求解但建模内核跟两阶段模型是一致的。第二个是考虑不确定性的随机优化或者鲁棒优化。随机优化需要生成大量场景代码量适中但场景数多了求解时间会直线上涨鲁棒优化则用不确定集合来描述预测误差不需要场景概率模型复杂度略高但解出来的方案保守性更强。两阶段模型本身就是一个天然的鲁棒优化结构把第二阶段写成“在最坏场景下调整无功手段”就成了经典的两阶段鲁棒优化这个方向近年在顶级期刊上的热度一直很高如果你手里已经有两阶段的确定性代码往鲁棒方向扩展的路径相当顺畅。第三个是考虑网络重构的联合优化。配电网里除了无功补偿还有联络开关可以进行网络拓扑调整这在故障恢复和降损场景下效果显著。但一旦引入开关状态变量问题就变成混合整数二阶锥规划MISOCP求解难度又上一个台阶——不过这恰恰是很多论文能发得出来的原因算力不够、方法来凑分解算法、Benders分解、拉格朗日松弛这些高级解法都能派上用场。我个人实际做下来比较推荐的进阶路线是先把确定性两阶段模型做扎实跑通IEEE 33节点和IEEE 123节点两个标准算例然后加一个日内滚动模块模拟全天运行最后再上一套场景法做不确定性分析。这三步走完别说硕士毕业论文就是一篇还不错的期刊小论文的主体内容也都能拿得出手了。