提到含分布式电源的配电网日前两阶段优化调度模型很多刚接触电力系统方向的同学第一反应是先找个粒子群算法套上去跑个曲线出来。但如果你真正在Matlab里从头搭过一个完整的、可复现的调度模型就会知道粒子群只是最后一步的花架子真正花时间的地方是建模逻辑、约束处理、求解器接口和两阶段之间的数据衔接。这个项目标题看上去只是一个模型加一段代码实际上它把配电网潮流、分布式电源出力特性、日前预测不确定性、调度策略分层这几个大问题全部揉在了一起既适合做毕业设计、论文仿真也适合刚入门的工程师用它理解日前计划日内修正这套调度体系到底是怎么运转的。这篇文章我会从模型设计思路、关键数学方法、Matlab具体实现、常见调参排错四个方面展开尽量把每一步为什么这么做讲清楚给出一套可以直接复用的框架。你不需要提前掌握所有配电网细节但懂一点Matlab矩阵操作和基础的优化概念会更顺。如果你正在为怎么写好调度代码发愁或者论文里差一个拿得出手的仿真结果这篇文章应该能帮你少走不少弯路。1. 项目概述与问题拆解1.1 分布式电源接入后配电网调度的核心矛盾先聊一个最基本的场景。过去配电网是单电源、单向潮流调度员只需要关心从变电站往下送多少电。但光伏、风电、储能、电动汽车充电桩大量接入后配电网从被动末端变成了有源网络潮流方向可能反转节点电压可能越限变压器可能倒送过载。更麻烦的是分布式电源出力受天气影响具有明显的间歇性和随机性如果调度策略只是简单按预测值安排出力到了实际运行时刻难免出现偏差要么弃光弃风要么电压越限。所以配电网调度问题的本质是在有限的可控资源联络线功率、储能充放电、可调分布式电源、可中断负荷和不确定的扰动源光伏波动、负荷波动之间找到一组最优的决策序列使得运行成本最低、网络损耗最小、电压质量最好同时满足所有安全和物理约束。这个多目标、多约束、含不确定性的问题单靠一个全局优化模型硬解很容易算出理论上最优、实际上无法执行的愚蠢结果。这就是两阶段调度被广泛使用的根本原因——把问题按决策时域和不确定性来源拆成两半。1.2 日前两阶段到底解决什么问题日前意思是提前一天做计划两阶段指的是决策过程分成前后两个层次。第一阶段通常是指前一日基于预测信息典型日曲线决定那些需要提前确定、响应慢或启停成本高的决策量比如火电机组或上级电网购电的日前计划出力、储能日前充放电计划、联络线交换功率等。第二阶段则是日内或实时运行阶段在日前计划的基础上根据实际的光伏出力、负荷实测值对可控资源做再调整修正第一阶段计划带来的偏差。注意这里两阶段在不同文献里含义有差异有的指日前日内有的指日前机组组合实时经济调度有的指同一优化问题引入场景的两阶段随机规划。做Matlab代码之前建议先把你的场景定义清楚不然很容易把程序结构写乱。我下面讲的方案是日前预测调度日内偏差修正的确定性框架适合作为项目基础后续再扩展成多场景随机规划也方便。1.3 模型整体技术路线整个模型的实现路径可以拆成下面几个环节算例系统构建选定一个配电网拓扑IEEE 33节点、IEEE 69节点或自己搭的典型馈线配好线路参数、负荷曲线、分布式电源参数。不确定性建模为光伏、风电、负荷设置预测误差场景或典型日曲线确定日前预测值和日内修正所需的扰动范围。第一阶段建模构建日前调度优化模型目标函数取运行成本最低或网损最小约束包括潮流约束、节点电压约束、支路容量约束、DG出力约束、储能约束、联络线功率约束。第二阶段建模在日前计划固定部分变量的前提下针对预测误差场景做再调度目标通常是偏差修正成本最小或调整量最小。求解与衔接用Matlab调用YalmipCplex/Gurobi或Cplex for MATLAB求解把第一阶段的部分结果作为第二阶段的输入参数。结果后处理绘制电压分布、DG出力曲线、储能SOC曲线、各时段成本构成等图表。这套路线最大的优势是模块化。你在处理第二阶段时不必把第一阶段的全部逻辑重写一遍只需在已有模型上固定变量或加偏差项Matlab里用几分钟就能完成。下面每一章都对应路线里的一个环节我逐个展开。2. 核心建模思路与关键技术解析2.1 两阶段优化框架的第一阶段日前计划先定义第一阶段要决策的对象。对配电网而言最核心的决策变量通常包括上级电网联络线有功无功交换功率、储能各时段充放电功率和充放电状态、可调DG比如柴油机或燃气轮机的有功无功出力、可中断负荷的削减量。光伏和风电因为不可控在日前阶段一般按预测值作为负的负荷处理不参与决策。目标函数比较常见的有两种写法。一种是运行成本最低 [ \min \sum_{t1}^{T} \left[ c_{grid,t} P_{grid,t} \sum_{g \in \mathcal{G}{DG}} c{DG,g} P_{DG,g,t} c_{ess}^{dis} P_{ess,t}^{dis} c_{curtail} P_{curtail,t} \right] ] 另一种是网损最小化 [ \min \sum_{t1}^{T} \sum_{(i,j) \in \Omega_{br}} I_{ij,t}^2 r_{ij} ] 两种目标也可以加权合并。论文里经常会开多目标的噱头但实际用Matlab实现时建议还是先用单目标把模型跑通再扩展成加权多目标调试成本会低很多。第一阶段的关键在于计划值要留有余地。比如储能日前计划不必把SOC充满到100%再放完因为日内出现偏差时你手里已经没有可调节的储能容量了。这类专业性细节在论文里不太会有人强调但实际仿真时你会发现如果日前计划不留备用第二阶段的修正空间会非常小甚至出现无解。2.2 第二阶段日内修正与实际执行第二阶段解决的问题是到了当天某个时段光伏和负荷的实际值偏离了预测值此时在电网安全约束下以最小调整成本将各可控资源从日前计划值调整到新的运行点。实现方式有两种常见做法。第一种做法是把第二阶段建模为一个重新优化问题。将日前计划中已确定的联络线功率、储能SOC轨迹等作为已知参数重新求解一次优化允许DG功率、储能功率、切负荷量等变量在日前计划基础上做增量调整目标是对调整量施加惩罚成本 [ \min \sum_{t} w_P|\Delta P_{t}| w_Q|\Delta Q_{t}| c_{penalty} P_{curtail,t} ] 第二种做法是机会约束或鲁棒修正在约束里预留不确定区间。这种方法更高级但求解难度也大前几版代码不建议直接上鲁棒容易把自己绕晕。从Matlab代码结构来看第二阶段最省事的方法是复用第一阶段的约束函数。写一个构建约束的函数输入参数是一个结构体包含哪些变量是固定值第二阶段生成约束时把日前计划值传入对相应变量设等式约束其余变量保留自由度求解器自动给出修正解。这个函数复用的设计在做两类模型时极其好用。2.3 配电网潮流约束与凸松弛处理配电网通常呈辐射状所以潮流计算常用DistFlow方程。对每条支路(i,j)[ P_{ij,t} - r_{ij} I_{ij,t}^2 \sum_{k: j \to k} P_{jk,t} P_{j,t}^{load} - P_{j,t}^{dg} ] [ Q_{ij,t} - x_{ij} I_{ij,t}^2 \sum_{k: j \to k} Q_{jk,t} Q_{j,t}^{load} - Q_{j,t}^{dg} ] [ V_{i,t}^2 - V_{j,t}^2 2(r_{ij} P_{ij,t} x_{ij} Q_{ij,t}) - (r_{ij}^2 x_{ij}^2) I_{ij,t}^2 ]这三个方程里面存在 (I_{ij,t}^2) 和 (P_{ij,t}^2 Q_{ij,t}^2 I_{ij,t}^2 V_{i,t}^2) 这样的非线性项直接求解是非凸的Matlab里的求解器一般解不了非凸问题。工程上常用的做法是二阶锥松弛把 (I_{ij,t}^2) 和 (V_{i,t}^2) 当作新变量然后用凸的锥约束近似原方程。松弛之后原问题变成一个二阶锥规划SOCPCplex和Gurobi都能高效求解。二阶锥松弛偶尔会带来误差原因在于它把潮流方程放松了。如果网络重载或者参数不合理松弛后的解可能在物理上不可行。此时需要检查是否满足松弛间隙即 ( \lVert P_{ij}^2 Q_{ij}^2 - l_{ij}v_i \rVert_\infty \leq \epsilon )如果间隙太大说明这个算例不适合直接SOCP或者参数设置有误。2.4 不确定性处理场景还是鲁棒关于分布式电源出力的不确定性论文里主要有三类处理方式确定性预测、多场景随机、鲁棒区间。你项目标题里只写了日前两阶段没说具体用哪种所以我建议第一阶段先用场景法做也就是用几个典型场景比如晴天、多云、阴雨天分别做日前计划日内再随机采样若干误差场景做修正。这种做法的优点是可以直接看出调度策略对不同天气的适应能力仿真图也好看。如果你想更进一步把问题做成两阶段随机规划那么第一阶段变量是here and now决策第二阶段变量是wait and see决策目标函数变成 [ \min c^T x \sum_{s\in\mathcal{S}} \pi_s Q(x, \xi_s) ] 其中 (\pi_s) 是场景概率(Q(x,\xi_s)) 是场景 (s) 下第二阶段的最优成本。这在Matlab里也能实现只需把场景循环先展开再将所有约束拼进一个大的优化模型里或者写成Benders分解。对大多数毕业设计而言先写确定性两阶段再扩展场景循环是最稳的路线。3. 基于Matlab的代码实现过程3.1 环境准备与工具箱选型先交代我的运行环境方便你对照Matlab R2022b及以上2026b也没问题核心逻辑一致Yalmip工具箱Cplex 12.10 或 Gurobi 9.5以上二选一学生可用学术版可选Matpower用于潮流对比验证有几个热搜里提到的点值得单独说。很多同学问Matlab装了Gurobi还能装Cplex吗答案是完全可以。Yalmip只是建模层后端求解器可以同时注册多个用ops sdpsettings(solver,cplex)或gurobi指定即可互不冲突。另外如果你只是学习建模不想折腾商业求解器可以先用Matlab自带的linprog或quadprog跑通线性化版本甚至fmincon也能跑小算例但涉及SOCP时强烈建议商业求解器。Yalmip安装很简单把文件夹解压放到任意目录进入该目录后运行addpath(genpath(pwd)); savepath;即可。Cplex for MATLAB安装时需要注意版本匹配建议直接到IBM官网下载与Matlab版本对应的安装包安装完成后在Matlab里执行cplex.getVersion验证。3.2 算例系统构建与参数设计我推荐用IEEE 33节点配电系统做算例原因有三个节点规模适中学生电脑几秒就能出解线路和负荷参数公开文献里到处都能找到对照它是一条辐射状馈线非常方便验证DistFlow约束。构建系统时建议把所有参数写成一个结构体mpc struct(); mpc.bus busdata(:, 1:3); % 列节点编号、有功负荷、无功负荷 mpc.branch branchdata(:, 1:4); % 列首端节点、末端节点、电阻r、电抗x mpc.baseMVA 10; mpc.Vmax 1.05; % 电压上限标幺值 mpc.Vmin 0.95; mpc.Sbase 12.66; % 电压基准/kV然后在指定节点接入分布式电源。以33节点系统为例常见做法是在节点18接入光伏因为它在馈线末端电压问题最严重、节点22接入风电、节点25接入储能配一台小型燃气轮机在节点8作为可调DG。光伏和风电的日出力曲线可以自己构造也可以用实测数据。如果做论文建议用典型日的标幺曲线再乘上装机容量例如tpv 0:1:24; pv_profile 0.5 * max(0, sin(pi * (tpv - 6) / 12)).^1.2; wind_profile 0.35 0.25 * sin(2*pi*tpv/24 - 0.5) 0.1*randn(1,25);这里加了一点随机项是为了体现功率波动。注意randn每次运行结果不同如果你要求结果可复现记得在代码最前面写rng(2024);。这个细节做实验时极重要审稿人和导师都会看你的结果能不能复现。3.3 第一阶段模型代码实现先说决策变量的定义方式。用Yalmip定义24时段、n个节点的优化变量最常见的是二维变量sdpvar(n, T)其中行是节点或机组编号列是时段。下面给出一个简化但完整的日前调度模型核心代码T 24; nb 33; % 节点数 ng 3; % 燃气轮机台数可调DG nw 2; % 风机数量 npv 2; % 光伏数量 ness 1; % 储能数量 % 决策变量 P_grid sdpvar(1, T, full); % 联络线有功 Q_grid sdpvar(1, T, full); % 联络线无功 P_dg sdpvar(ng, T, full); % 燃气轮机有功 Q_dg sdpvar(ng, T, full); P_ch sdpvar(ness, T, full); % 储能充电功率 P_dis sdpvar(ness, T, full); % 储能放电功率 SOC sdpvar(ness, T1, full); % 储能SOC z_ch binvar(ness, T, full); % 充电状态指示 z_dis binvar(ness, T, full); % 放电状态指示目标函数定义为购电成本、燃气轮机燃料成本、储能充放电老化成本之和并加入弃光和弃风惩罚c_grid 0.5 * (1 0.1*sin((1:T)/T*pi)); % 分时电价 cost sum(c_grid .* P_grid) ... sum(sum(c_dg .* P_dg)) ... sum(sum(0.02 * P_ch 0.03 * P_dis)) ... sum(sum(50 * P_curtail_pv)) sum(sum(40 * P_curtail_wind));注意到这里P_curtail_pv和P_curtail_wind一般定义为非负连续变量作为光伏逆变器弃功率和风机弃功量它们同时要满足实际出力弃功率预测出力这样的约束。约束条件按块写。潮流约束我用DistFlow并做了二阶锥松弛其中l_ij代表(I_{ij,t}^2)v_i代表(V_{i,t}^2)for t 1:T for k 1:nb % 节点功率平衡DistFlow形式 % 父支路流进的功率 - 子支路流出的功率 注入功率 end % 二阶锥约束P_ij^2 Q_ij^2 l_ij * v_i for branch 1:nb-1 Pij sdpvar(1,1); Qij sdpvar(1,1); lij sdpvar(1,1); vi sdpvar(1,1); Constraints [Constraints, Pij^2 Qij^2 lij * vi]; end end写Yalmip代码时有个很实用的技巧不要在循环里直接追加Constraints [Constraints, expr]太多次那样会拖慢Matlab的解析过程。正确做法是先用cell数组收集约束最后Constraints [Constraints{:}]一次性合并速度能快很多。这个优化在你把T改成96时段15分钟一个点时尤其明显。储能约束需要特别注意时间耦合因为SOC从第1个时段传到第t1个时段Constraints [Constraints, SOC(:, 1) 0.3]; % 初始SOC for t 1:T Constraints [Constraints, SOC(:, t1) SOC(:, t) eta_ch * P_ch(:, t) - P_dis(:, t) / eta_dis]; Constraints [Constraints, 0 SOC(:, t1) 1]; Constraints [Constraints, P_ch(:, t) P_ch_max * z_ch(:, t)]; Constraints [Constraints, P_dis(:, t) P_dis_max * z_dis(:, t)]; Constraints [Constraints, z_ch(:, t) z_dis(:, t) 1]; % 不能同时充放电 end严格来说充放电同时性约束在SOC不为零的情况下会让问题变成一个MILP因为有二进制变量求解速度会下降但这是物理约束不能省。3.4 第二阶段模型与两阶段衔接实现第二阶段不是从零开始。假设你已经把第一阶段的优化结果保存到了结构体result_day里第二阶段的模型只需要做三件事把光伏、风电、负荷不再是确定的预测曲线而是当前时段的实测曲线或误差场景。例如设P_pv_real P_pv_pred delta_pv其中delta_pv按正态分布抽样。把第一阶段已经确定的变量作为已知参数传进来P_grid result_day.P_grid; % 日前计划第二阶段固定 SOC_fixed result_day.SOC; % 日前计划SOC轨迹固定把可调DG、储能功率、切负荷量作为修正决策变量允许他们偏离日前计划但偏离要付出惩罚成本。第二阶段目标函数可以写成cost_pen sum(P_price_pen .* abs(P_dg - P_dg_scheduled)) ... sum(sum(P_curtail_pen .* P_curtail_real));如果你嫌abs()在Yalmip里不好处理两个解决办法一是引入辅助变量delta要求delta P_dg - P_dg_scheduled且delta -(P_dg - P_dg_scheduled)目标里写sum(delta)二是如果你装了Gurobi/Cplex这类支持MILP的求解器abs()本身也能处理只是中间会引入额外变量速度稍慢。第二阶段模型里潮流约束、电压约束、线路容量约束必须重新生成因为变量变了运行点也变了。这里最笨但最不容易出错的办法是写一个函数function [cons, vars] build_power_flow_constraints(mpc, t, var_struct) % 输入当前时段的变量结构体输出该时段潮流约束这样第一阶段和第二阶段都可以调它传入不同的变量结构体减少大量重复代码。这个方法我建议所有做配电网优化的同学都学会它本质上是把场景和模型解耦后续加新场景只是多调一次函数而不是复制粘贴一大段代码。当你跑完第二阶段后还会想验证一下结果在真实网络的可行性。做法是把求得的各节点注入功率退回潮流计算器比如用Matpower做一次潮流验证runpf看电压和线路是否越限。这一步属于闭环验证论文里加一张表列出优化结果电压vs潮流验证电压专家一眼就觉得你工程功底扎实。3.5 求解结果与可视化分析仿真跑完后至少要画四张图电压分布图、各电源出力时序图、储能SOC曲线、成本箱线图多场景时用。画出这些图不是单纯为了好看而是用于验证三个事情电压是否在0.95到1.05之间有没有节点电压越下限。光伏出力高峰时段是否出现倒送联络线功率是否为负向上一级电网送电。储能SOC曲线是否在安全范围充放电切换是否频繁频繁切换在真实系统里是不可接受的。电压分布图用二维热力图展示最直观横轴是时间纵轴是节点号颜色代表节点电压色条贴在图右侧。Matlab代码如下V_plot reshape(V_result, nb, T); figure; imagesc(1:T, 1:nb, V_plot); colorbar; xlabel(时间/h); ylabel(节点编号); title(节点电压时序分布); caxis([0.95 1.05]);如果发现电压越限不要慌先检查是哪个节点在哪个时段越限。通常解决办法包括增加储能充电功率、提高DG无功吸收能力储能或逆变器提供无功支撑、调整有载调压变压器的挡位。4. 常见问题与调试经验4.1 求解器报错的典型原因与处理Yalmip报错最常见的几类我列在下面你遇到了可以直接对照报错信息可能原因解决方法No suitable solver found模型里有二进制变量或锥约束但未安装对应求解器安装Cplex/Gurobi或把二进制变量松驰化不推荐/改用线性化方案NaN in constraints变量定义维数不一致或某参数为NaN检查所有load数据里有没有空值打印nnz(isnan(data))Infeasible problem约束冲突常见于储能初始SOC与目标SOC矛盾逐条注释约束排查或先用optimize返回的info查看Yalmip给出的诊断Solver timed out规模太大、二进制变量太多减少时段数96改24、减少场景数、或给求解器加TimeLimit参数Cplex error 5002模型包含非凸二次约束检查是否有P_ij*V_i这类双线性项确认是否已做变量替换和松弛一个非常实用的技巧第一次求解先用ops sdpsettings(verbose, 2, solver, cplex)把完整日志打开求解器会告诉你卡在哪个约束上。别嫌日志刷屏它是你调试最好的朋友。4.2 二阶锥松弛不收敛怎么办SOCP模型最常见的玄学问题是求解器提示求解成功但你把解代回原潮流方程发现节点电压对不上或者存在较大的松弛间隙。出现这个问题的原因通常是你把潮流方程简化过头了比如把V_i^2和I_ij^2当成完全独立的变量却没有严格执行锥约束[P_ij; Q_ij; (V_i/sqrt(2)); ...] 属于旋转二阶锥Yalmip写旋转锥约束时用cone([P_ij; Q_ij; (v_i - l_ij)/sqrt(2)], (v_i l_ij)/sqrt(2))这种形式或者[P_ij^2 Q_ij^2 l_ij * v_i]都行但后者如果v_i接近0数值稳定性会变差。建议给节点电压平方加个下限比如v_i 0.8^2让求解器不用在病态区域探索。如果松弛间隙还是很大还有一种工程上较稳的兜底方案把最优解作为初值重新用交流潮流方程做一次还原也就是通过逐次逼近或内点法把SOCP解投影到真实潮流可行域上。这就是所谓的恢复算法在你做论文时可以写一小节松弛解的可行性校验来体现严谨性。4.3 运行时间过长如何提速配电网两阶段模型如果规模太大Matlab运行时间可能从几秒飙到几分钟甚至更久。提速的核心是减少维度和变量数。把一天96个时段压缩到24个时段。如果你只是研究算法15分钟粒度和1小时粒度在策略层面的差别并不大先把模型跑通再加密。选择稀疏参数。Cplex和Gurobi对稀疏模型很友好Yalmip会自动检测稀疏度但如果你在Matlab脚本里用了大量full()转换会破坏稀疏结构。除非必要不要轻易把sdpvar转成full矩阵。删除对结果影响极小的冗余变量。比如一些节点上并没有接入任何DG或负荷那它的注入功率就是零不需要定义为决策变量直接在潮流方程里把对应的P_load设成0即可。求解器参数设置。给Cplex设置ops sdpsettings(solver,cplex,cplex.mip.tolerances.mipgap,0.001)允许MIP目标在0.1%的间隙内停止速度提升非常明显。4.4 结果不合理时的排查清单有时候模型有解、求解器也不报错但结果一看就离谱比如DG全部停机、储能一次都不充电、弃光率100%而联络线却满载送电。这时候不急着改代码先按下面清单逐项排查目标函数量纲是否一致购电成本单位是元/kWhDG成本单位是不是元/MWh如果差1000倍优化结果会疯狂偏向某一项。负荷曲线和DG出力曲线是否有单位坐标不一致的问题。我见过很多同学把负荷功率当成kW、光伏功率当成MW来计算结果配网潮流严重失衡。储能SOC的初值和终值设定是否合理。如果要求SOC一天结束后回到初始值那储能一天的总充电量必须等于总放电量再考虑效率这个约束会让储能调度空间变小。如果你不想加终值约束也可以不加但需要明确储能只是日内的调节工具。电价时段是否在目标函数里正确使用。有些同学把分时电价向量和P_grid维度搞反导致目标函数出现size mismatchYalmip不报错但会隐式广播结果完全不可解释。最后检查一下神经网络常用的归一化问题。在做算例分析时如果你对数据做了归一化别忘了在结果展示阶段反归一化回来。电压本来就是标幺值别再归一化一遍否则图上的电压范围完全不符合物理意义。排查完以上项目绝大多数诡异结果都能找到原因。如果还找不到就从头打印每一时段各变量的数值用一个小时时段看注入功率的平衡关系该节点注入功率下级支路流出本地负荷。手算一遍通常能立刻发现问题。我在做这个项目时最大的体会是写Matlab代码本身不难难的是把每个环节的物理含义翻译成数学约束再翻译成代码约束最后用合理的可视化检查结果。尤其是第一和第二阶段之间哪些变量要固定、哪些变量要放开这个边界如果没想清楚代码很容易写出一个第一第二阶段的优化结果完全没关系的缝合怪。最后分享一个小技巧把整个模型的输入参数和求解配置全部抽到一个config.m脚本里运行主程序之前先运行config.m。这样你实验时只需要改一个文件里的数值不需要去主程序里面翻变量别人拿到你的代码复现论文结果也会非常顺手。分布式电源配电网调度这个方向真正的门槛不在算法本身而在你是否能把一套模型用工程化思维实现得足够干净、可控、可复用。希望这篇文章能帮你把那个门槛跨过去。 含分布式电源的配电网日前两阶段优化调度Matlab实现
提到含分布式电源的配电网日前两阶段优化调度模型很多刚接触电力系统方向的同学第一反应是先找个粒子群算法套上去跑个曲线出来。但如果你真正在Matlab里从头搭过一个完整的、可复现的调度模型就会知道粒子群只是最后一步的花架子真正花时间的地方是建模逻辑、约束处理、求解器接口和两阶段之间的数据衔接。这个项目标题看上去只是一个模型加一段代码实际上它把配电网潮流、分布式电源出力特性、日前预测不确定性、调度策略分层这几个大问题全部揉在了一起既适合做毕业设计、论文仿真也适合刚入门的工程师用它理解日前计划日内修正这套调度体系到底是怎么运转的。这篇文章我会从模型设计思路、关键数学方法、Matlab具体实现、常见调参排错四个方面展开尽量把每一步为什么这么做讲清楚给出一套可以直接复用的框架。你不需要提前掌握所有配电网细节但懂一点Matlab矩阵操作和基础的优化概念会更顺。如果你正在为怎么写好调度代码发愁或者论文里差一个拿得出手的仿真结果这篇文章应该能帮你少走不少弯路。1. 项目概述与问题拆解1.1 分布式电源接入后配电网调度的核心矛盾先聊一个最基本的场景。过去配电网是单电源、单向潮流调度员只需要关心从变电站往下送多少电。但光伏、风电、储能、电动汽车充电桩大量接入后配电网从被动末端变成了有源网络潮流方向可能反转节点电压可能越限变压器可能倒送过载。更麻烦的是分布式电源出力受天气影响具有明显的间歇性和随机性如果调度策略只是简单按预测值安排出力到了实际运行时刻难免出现偏差要么弃光弃风要么电压越限。所以配电网调度问题的本质是在有限的可控资源联络线功率、储能充放电、可调分布式电源、可中断负荷和不确定的扰动源光伏波动、负荷波动之间找到一组最优的决策序列使得运行成本最低、网络损耗最小、电压质量最好同时满足所有安全和物理约束。这个多目标、多约束、含不确定性的问题单靠一个全局优化模型硬解很容易算出理论上最优、实际上无法执行的愚蠢结果。这就是两阶段调度被广泛使用的根本原因——把问题按决策时域和不确定性来源拆成两半。1.2 日前两阶段到底解决什么问题日前意思是提前一天做计划两阶段指的是决策过程分成前后两个层次。第一阶段通常是指前一日基于预测信息典型日曲线决定那些需要提前确定、响应慢或启停成本高的决策量比如火电机组或上级电网购电的日前计划出力、储能日前充放电计划、联络线交换功率等。第二阶段则是日内或实时运行阶段在日前计划的基础上根据实际的光伏出力、负荷实测值对可控资源做再调整修正第一阶段计划带来的偏差。注意这里两阶段在不同文献里含义有差异有的指日前日内有的指日前机组组合实时经济调度有的指同一优化问题引入场景的两阶段随机规划。做Matlab代码之前建议先把你的场景定义清楚不然很容易把程序结构写乱。我下面讲的方案是日前预测调度日内偏差修正的确定性框架适合作为项目基础后续再扩展成多场景随机规划也方便。1.3 模型整体技术路线整个模型的实现路径可以拆成下面几个环节算例系统构建选定一个配电网拓扑IEEE 33节点、IEEE 69节点或自己搭的典型馈线配好线路参数、负荷曲线、分布式电源参数。不确定性建模为光伏、风电、负荷设置预测误差场景或典型日曲线确定日前预测值和日内修正所需的扰动范围。第一阶段建模构建日前调度优化模型目标函数取运行成本最低或网损最小约束包括潮流约束、节点电压约束、支路容量约束、DG出力约束、储能约束、联络线功率约束。第二阶段建模在日前计划固定部分变量的前提下针对预测误差场景做再调度目标通常是偏差修正成本最小或调整量最小。求解与衔接用Matlab调用YalmipCplex/Gurobi或Cplex for MATLAB求解把第一阶段的部分结果作为第二阶段的输入参数。结果后处理绘制电压分布、DG出力曲线、储能SOC曲线、各时段成本构成等图表。这套路线最大的优势是模块化。你在处理第二阶段时不必把第一阶段的全部逻辑重写一遍只需在已有模型上固定变量或加偏差项Matlab里用几分钟就能完成。下面每一章都对应路线里的一个环节我逐个展开。2. 核心建模思路与关键技术解析2.1 两阶段优化框架的第一阶段日前计划先定义第一阶段要决策的对象。对配电网而言最核心的决策变量通常包括上级电网联络线有功无功交换功率、储能各时段充放电功率和充放电状态、可调DG比如柴油机或燃气轮机的有功无功出力、可中断负荷的削减量。光伏和风电因为不可控在日前阶段一般按预测值作为负的负荷处理不参与决策。目标函数比较常见的有两种写法。一种是运行成本最低 [ \min \sum_{t1}^{T} \left[ c_{grid,t} P_{grid,t} \sum_{g \in \mathcal{G}{DG}} c{DG,g} P_{DG,g,t} c_{ess}^{dis} P_{ess,t}^{dis} c_{curtail} P_{curtail,t} \right] ] 另一种是网损最小化 [ \min \sum_{t1}^{T} \sum_{(i,j) \in \Omega_{br}} I_{ij,t}^2 r_{ij} ] 两种目标也可以加权合并。论文里经常会开多目标的噱头但实际用Matlab实现时建议还是先用单目标把模型跑通再扩展成加权多目标调试成本会低很多。第一阶段的关键在于计划值要留有余地。比如储能日前计划不必把SOC充满到100%再放完因为日内出现偏差时你手里已经没有可调节的储能容量了。这类专业性细节在论文里不太会有人强调但实际仿真时你会发现如果日前计划不留备用第二阶段的修正空间会非常小甚至出现无解。2.2 第二阶段日内修正与实际执行第二阶段解决的问题是到了当天某个时段光伏和负荷的实际值偏离了预测值此时在电网安全约束下以最小调整成本将各可控资源从日前计划值调整到新的运行点。实现方式有两种常见做法。第一种做法是把第二阶段建模为一个重新优化问题。将日前计划中已确定的联络线功率、储能SOC轨迹等作为已知参数重新求解一次优化允许DG功率、储能功率、切负荷量等变量在日前计划基础上做增量调整目标是对调整量施加惩罚成本 [ \min \sum_{t} w_P|\Delta P_{t}| w_Q|\Delta Q_{t}| c_{penalty} P_{curtail,t} ] 第二种做法是机会约束或鲁棒修正在约束里预留不确定区间。这种方法更高级但求解难度也大前几版代码不建议直接上鲁棒容易把自己绕晕。从Matlab代码结构来看第二阶段最省事的方法是复用第一阶段的约束函数。写一个构建约束的函数输入参数是一个结构体包含哪些变量是固定值第二阶段生成约束时把日前计划值传入对相应变量设等式约束其余变量保留自由度求解器自动给出修正解。这个函数复用的设计在做两类模型时极其好用。2.3 配电网潮流约束与凸松弛处理配电网通常呈辐射状所以潮流计算常用DistFlow方程。对每条支路(i,j)[ P_{ij,t} - r_{ij} I_{ij,t}^2 \sum_{k: j \to k} P_{jk,t} P_{j,t}^{load} - P_{j,t}^{dg} ] [ Q_{ij,t} - x_{ij} I_{ij,t}^2 \sum_{k: j \to k} Q_{jk,t} Q_{j,t}^{load} - Q_{j,t}^{dg} ] [ V_{i,t}^2 - V_{j,t}^2 2(r_{ij} P_{ij,t} x_{ij} Q_{ij,t}) - (r_{ij}^2 x_{ij}^2) I_{ij,t}^2 ]这三个方程里面存在 (I_{ij,t}^2) 和 (P_{ij,t}^2 Q_{ij,t}^2 I_{ij,t}^2 V_{i,t}^2) 这样的非线性项直接求解是非凸的Matlab里的求解器一般解不了非凸问题。工程上常用的做法是二阶锥松弛把 (I_{ij,t}^2) 和 (V_{i,t}^2) 当作新变量然后用凸的锥约束近似原方程。松弛之后原问题变成一个二阶锥规划SOCPCplex和Gurobi都能高效求解。二阶锥松弛偶尔会带来误差原因在于它把潮流方程放松了。如果网络重载或者参数不合理松弛后的解可能在物理上不可行。此时需要检查是否满足松弛间隙即 ( \lVert P_{ij}^2 Q_{ij}^2 - l_{ij}v_i \rVert_\infty \leq \epsilon )如果间隙太大说明这个算例不适合直接SOCP或者参数设置有误。2.4 不确定性处理场景还是鲁棒关于分布式电源出力的不确定性论文里主要有三类处理方式确定性预测、多场景随机、鲁棒区间。你项目标题里只写了日前两阶段没说具体用哪种所以我建议第一阶段先用场景法做也就是用几个典型场景比如晴天、多云、阴雨天分别做日前计划日内再随机采样若干误差场景做修正。这种做法的优点是可以直接看出调度策略对不同天气的适应能力仿真图也好看。如果你想更进一步把问题做成两阶段随机规划那么第一阶段变量是here and now决策第二阶段变量是wait and see决策目标函数变成 [ \min c^T x \sum_{s\in\mathcal{S}} \pi_s Q(x, \xi_s) ] 其中 (\pi_s) 是场景概率(Q(x,\xi_s)) 是场景 (s) 下第二阶段的最优成本。这在Matlab里也能实现只需把场景循环先展开再将所有约束拼进一个大的优化模型里或者写成Benders分解。对大多数毕业设计而言先写确定性两阶段再扩展场景循环是最稳的路线。3. 基于Matlab的代码实现过程3.1 环境准备与工具箱选型先交代我的运行环境方便你对照Matlab R2022b及以上2026b也没问题核心逻辑一致Yalmip工具箱Cplex 12.10 或 Gurobi 9.5以上二选一学生可用学术版可选Matpower用于潮流对比验证有几个热搜里提到的点值得单独说。很多同学问Matlab装了Gurobi还能装Cplex吗答案是完全可以。Yalmip只是建模层后端求解器可以同时注册多个用ops sdpsettings(solver,cplex)或gurobi指定即可互不冲突。另外如果你只是学习建模不想折腾商业求解器可以先用Matlab自带的linprog或quadprog跑通线性化版本甚至fmincon也能跑小算例但涉及SOCP时强烈建议商业求解器。Yalmip安装很简单把文件夹解压放到任意目录进入该目录后运行addpath(genpath(pwd)); savepath;即可。Cplex for MATLAB安装时需要注意版本匹配建议直接到IBM官网下载与Matlab版本对应的安装包安装完成后在Matlab里执行cplex.getVersion验证。3.2 算例系统构建与参数设计我推荐用IEEE 33节点配电系统做算例原因有三个节点规模适中学生电脑几秒就能出解线路和负荷参数公开文献里到处都能找到对照它是一条辐射状馈线非常方便验证DistFlow约束。构建系统时建议把所有参数写成一个结构体mpc struct(); mpc.bus busdata(:, 1:3); % 列节点编号、有功负荷、无功负荷 mpc.branch branchdata(:, 1:4); % 列首端节点、末端节点、电阻r、电抗x mpc.baseMVA 10; mpc.Vmax 1.05; % 电压上限标幺值 mpc.Vmin 0.95; mpc.Sbase 12.66; % 电压基准/kV然后在指定节点接入分布式电源。以33节点系统为例常见做法是在节点18接入光伏因为它在馈线末端电压问题最严重、节点22接入风电、节点25接入储能配一台小型燃气轮机在节点8作为可调DG。光伏和风电的日出力曲线可以自己构造也可以用实测数据。如果做论文建议用典型日的标幺曲线再乘上装机容量例如tpv 0:1:24; pv_profile 0.5 * max(0, sin(pi * (tpv - 6) / 12)).^1.2; wind_profile 0.35 0.25 * sin(2*pi*tpv/24 - 0.5) 0.1*randn(1,25);这里加了一点随机项是为了体现功率波动。注意randn每次运行结果不同如果你要求结果可复现记得在代码最前面写rng(2024);。这个细节做实验时极重要审稿人和导师都会看你的结果能不能复现。3.3 第一阶段模型代码实现先说决策变量的定义方式。用Yalmip定义24时段、n个节点的优化变量最常见的是二维变量sdpvar(n, T)其中行是节点或机组编号列是时段。下面给出一个简化但完整的日前调度模型核心代码T 24; nb 33; % 节点数 ng 3; % 燃气轮机台数可调DG nw 2; % 风机数量 npv 2; % 光伏数量 ness 1; % 储能数量 % 决策变量 P_grid sdpvar(1, T, full); % 联络线有功 Q_grid sdpvar(1, T, full); % 联络线无功 P_dg sdpvar(ng, T, full); % 燃气轮机有功 Q_dg sdpvar(ng, T, full); P_ch sdpvar(ness, T, full); % 储能充电功率 P_dis sdpvar(ness, T, full); % 储能放电功率 SOC sdpvar(ness, T1, full); % 储能SOC z_ch binvar(ness, T, full); % 充电状态指示 z_dis binvar(ness, T, full); % 放电状态指示目标函数定义为购电成本、燃气轮机燃料成本、储能充放电老化成本之和并加入弃光和弃风惩罚c_grid 0.5 * (1 0.1*sin((1:T)/T*pi)); % 分时电价 cost sum(c_grid .* P_grid) ... sum(sum(c_dg .* P_dg)) ... sum(sum(0.02 * P_ch 0.03 * P_dis)) ... sum(sum(50 * P_curtail_pv)) sum(sum(40 * P_curtail_wind));注意到这里P_curtail_pv和P_curtail_wind一般定义为非负连续变量作为光伏逆变器弃功率和风机弃功量它们同时要满足实际出力弃功率预测出力这样的约束。约束条件按块写。潮流约束我用DistFlow并做了二阶锥松弛其中l_ij代表(I_{ij,t}^2)v_i代表(V_{i,t}^2)for t 1:T for k 1:nb % 节点功率平衡DistFlow形式 % 父支路流进的功率 - 子支路流出的功率 注入功率 end % 二阶锥约束P_ij^2 Q_ij^2 l_ij * v_i for branch 1:nb-1 Pij sdpvar(1,1); Qij sdpvar(1,1); lij sdpvar(1,1); vi sdpvar(1,1); Constraints [Constraints, Pij^2 Qij^2 lij * vi]; end end写Yalmip代码时有个很实用的技巧不要在循环里直接追加Constraints [Constraints, expr]太多次那样会拖慢Matlab的解析过程。正确做法是先用cell数组收集约束最后Constraints [Constraints{:}]一次性合并速度能快很多。这个优化在你把T改成96时段15分钟一个点时尤其明显。储能约束需要特别注意时间耦合因为SOC从第1个时段传到第t1个时段Constraints [Constraints, SOC(:, 1) 0.3]; % 初始SOC for t 1:T Constraints [Constraints, SOC(:, t1) SOC(:, t) eta_ch * P_ch(:, t) - P_dis(:, t) / eta_dis]; Constraints [Constraints, 0 SOC(:, t1) 1]; Constraints [Constraints, P_ch(:, t) P_ch_max * z_ch(:, t)]; Constraints [Constraints, P_dis(:, t) P_dis_max * z_dis(:, t)]; Constraints [Constraints, z_ch(:, t) z_dis(:, t) 1]; % 不能同时充放电 end严格来说充放电同时性约束在SOC不为零的情况下会让问题变成一个MILP因为有二进制变量求解速度会下降但这是物理约束不能省。3.4 第二阶段模型与两阶段衔接实现第二阶段不是从零开始。假设你已经把第一阶段的优化结果保存到了结构体result_day里第二阶段的模型只需要做三件事把光伏、风电、负荷不再是确定的预测曲线而是当前时段的实测曲线或误差场景。例如设P_pv_real P_pv_pred delta_pv其中delta_pv按正态分布抽样。把第一阶段已经确定的变量作为已知参数传进来P_grid result_day.P_grid; % 日前计划第二阶段固定 SOC_fixed result_day.SOC; % 日前计划SOC轨迹固定把可调DG、储能功率、切负荷量作为修正决策变量允许他们偏离日前计划但偏离要付出惩罚成本。第二阶段目标函数可以写成cost_pen sum(P_price_pen .* abs(P_dg - P_dg_scheduled)) ... sum(sum(P_curtail_pen .* P_curtail_real));如果你嫌abs()在Yalmip里不好处理两个解决办法一是引入辅助变量delta要求delta P_dg - P_dg_scheduled且delta -(P_dg - P_dg_scheduled)目标里写sum(delta)二是如果你装了Gurobi/Cplex这类支持MILP的求解器abs()本身也能处理只是中间会引入额外变量速度稍慢。第二阶段模型里潮流约束、电压约束、线路容量约束必须重新生成因为变量变了运行点也变了。这里最笨但最不容易出错的办法是写一个函数function [cons, vars] build_power_flow_constraints(mpc, t, var_struct) % 输入当前时段的变量结构体输出该时段潮流约束这样第一阶段和第二阶段都可以调它传入不同的变量结构体减少大量重复代码。这个方法我建议所有做配电网优化的同学都学会它本质上是把场景和模型解耦后续加新场景只是多调一次函数而不是复制粘贴一大段代码。当你跑完第二阶段后还会想验证一下结果在真实网络的可行性。做法是把求得的各节点注入功率退回潮流计算器比如用Matpower做一次潮流验证runpf看电压和线路是否越限。这一步属于闭环验证论文里加一张表列出优化结果电压vs潮流验证电压专家一眼就觉得你工程功底扎实。3.5 求解结果与可视化分析仿真跑完后至少要画四张图电压分布图、各电源出力时序图、储能SOC曲线、成本箱线图多场景时用。画出这些图不是单纯为了好看而是用于验证三个事情电压是否在0.95到1.05之间有没有节点电压越下限。光伏出力高峰时段是否出现倒送联络线功率是否为负向上一级电网送电。储能SOC曲线是否在安全范围充放电切换是否频繁频繁切换在真实系统里是不可接受的。电压分布图用二维热力图展示最直观横轴是时间纵轴是节点号颜色代表节点电压色条贴在图右侧。Matlab代码如下V_plot reshape(V_result, nb, T); figure; imagesc(1:T, 1:nb, V_plot); colorbar; xlabel(时间/h); ylabel(节点编号); title(节点电压时序分布); caxis([0.95 1.05]);如果发现电压越限不要慌先检查是哪个节点在哪个时段越限。通常解决办法包括增加储能充电功率、提高DG无功吸收能力储能或逆变器提供无功支撑、调整有载调压变压器的挡位。4. 常见问题与调试经验4.1 求解器报错的典型原因与处理Yalmip报错最常见的几类我列在下面你遇到了可以直接对照报错信息可能原因解决方法No suitable solver found模型里有二进制变量或锥约束但未安装对应求解器安装Cplex/Gurobi或把二进制变量松驰化不推荐/改用线性化方案NaN in constraints变量定义维数不一致或某参数为NaN检查所有load数据里有没有空值打印nnz(isnan(data))Infeasible problem约束冲突常见于储能初始SOC与目标SOC矛盾逐条注释约束排查或先用optimize返回的info查看Yalmip给出的诊断Solver timed out规模太大、二进制变量太多减少时段数96改24、减少场景数、或给求解器加TimeLimit参数Cplex error 5002模型包含非凸二次约束检查是否有P_ij*V_i这类双线性项确认是否已做变量替换和松弛一个非常实用的技巧第一次求解先用ops sdpsettings(verbose, 2, solver, cplex)把完整日志打开求解器会告诉你卡在哪个约束上。别嫌日志刷屏它是你调试最好的朋友。4.2 二阶锥松弛不收敛怎么办SOCP模型最常见的玄学问题是求解器提示求解成功但你把解代回原潮流方程发现节点电压对不上或者存在较大的松弛间隙。出现这个问题的原因通常是你把潮流方程简化过头了比如把V_i^2和I_ij^2当成完全独立的变量却没有严格执行锥约束[P_ij; Q_ij; (V_i/sqrt(2)); ...] 属于旋转二阶锥Yalmip写旋转锥约束时用cone([P_ij; Q_ij; (v_i - l_ij)/sqrt(2)], (v_i l_ij)/sqrt(2))这种形式或者[P_ij^2 Q_ij^2 l_ij * v_i]都行但后者如果v_i接近0数值稳定性会变差。建议给节点电压平方加个下限比如v_i 0.8^2让求解器不用在病态区域探索。如果松弛间隙还是很大还有一种工程上较稳的兜底方案把最优解作为初值重新用交流潮流方程做一次还原也就是通过逐次逼近或内点法把SOCP解投影到真实潮流可行域上。这就是所谓的恢复算法在你做论文时可以写一小节松弛解的可行性校验来体现严谨性。4.3 运行时间过长如何提速配电网两阶段模型如果规模太大Matlab运行时间可能从几秒飙到几分钟甚至更久。提速的核心是减少维度和变量数。把一天96个时段压缩到24个时段。如果你只是研究算法15分钟粒度和1小时粒度在策略层面的差别并不大先把模型跑通再加密。选择稀疏参数。Cplex和Gurobi对稀疏模型很友好Yalmip会自动检测稀疏度但如果你在Matlab脚本里用了大量full()转换会破坏稀疏结构。除非必要不要轻易把sdpvar转成full矩阵。删除对结果影响极小的冗余变量。比如一些节点上并没有接入任何DG或负荷那它的注入功率就是零不需要定义为决策变量直接在潮流方程里把对应的P_load设成0即可。求解器参数设置。给Cplex设置ops sdpsettings(solver,cplex,cplex.mip.tolerances.mipgap,0.001)允许MIP目标在0.1%的间隙内停止速度提升非常明显。4.4 结果不合理时的排查清单有时候模型有解、求解器也不报错但结果一看就离谱比如DG全部停机、储能一次都不充电、弃光率100%而联络线却满载送电。这时候不急着改代码先按下面清单逐项排查目标函数量纲是否一致购电成本单位是元/kWhDG成本单位是不是元/MWh如果差1000倍优化结果会疯狂偏向某一项。负荷曲线和DG出力曲线是否有单位坐标不一致的问题。我见过很多同学把负荷功率当成kW、光伏功率当成MW来计算结果配网潮流严重失衡。储能SOC的初值和终值设定是否合理。如果要求SOC一天结束后回到初始值那储能一天的总充电量必须等于总放电量再考虑效率这个约束会让储能调度空间变小。如果你不想加终值约束也可以不加但需要明确储能只是日内的调节工具。电价时段是否在目标函数里正确使用。有些同学把分时电价向量和P_grid维度搞反导致目标函数出现size mismatchYalmip不报错但会隐式广播结果完全不可解释。最后检查一下神经网络常用的归一化问题。在做算例分析时如果你对数据做了归一化别忘了在结果展示阶段反归一化回来。电压本来就是标幺值别再归一化一遍否则图上的电压范围完全不符合物理意义。排查完以上项目绝大多数诡异结果都能找到原因。如果还找不到就从头打印每一时段各变量的数值用一个小时时段看注入功率的平衡关系该节点注入功率下级支路流出本地负荷。手算一遍通常能立刻发现问题。我在做这个项目时最大的体会是写Matlab代码本身不难难的是把每个环节的物理含义翻译成数学约束再翻译成代码约束最后用合理的可视化检查结果。尤其是第一和第二阶段之间哪些变量要固定、哪些变量要放开这个边界如果没想清楚代码很容易写出一个第一第二阶段的优化结果完全没关系的缝合怪。最后分享一个小技巧把整个模型的输入参数和求解配置全部抽到一个config.m脚本里运行主程序之前先运行config.m。这样你实验时只需要改一个文件里的数值不需要去主程序里面翻变量别人拿到你的代码复现论文结果也会非常顺手。分布式电源配电网调度这个方向真正的门槛不在算法本身而在你是否能把一套模型用工程化思维实现得足够干净、可控、可复用。希望这篇文章能帮你把那个门槛跨过去。