简介Matlab直流潮流计算程序包面向电力系统初学者及需要快速开展稳态分析的工程人员通过线性化直流模型解决大规模网络功率分布的近似计算问题。资源共9个文件以8个Matlab脚本为核心辅以1个dat数据文件压缩包仅4KB轻量易部署。脚本覆盖原始数据读取、节点导纳矩阵生成、功率残差计算、电压相角修正与直流潮流主迭代等环节输出模块可查看节点注入功率与线路潮流结果关键函数划分清晰便于二次开发与调试。已有399人学习适合配合教材、课程设计或论文仿真进行验证。通过修改数据文件和计算参数学习者能深入理解直流潮流的建模假设与求解流程进一步将程序扩展到标准节点系统或交直流混联场景为后续电力系统分析打下扎实基础。1. 直流潮流不是直流:交流电网快速有功分析的Matlab实现做直流潮流计算时Matlab程序包最常见的问题不是算法不会而是不知道每个.m文件该干什么。这套Matlab直流潮流程序把直流潮流计算拆成了data.m、FormY.m、FormP.m、CalP.m、CalDelta.m、DCPF.m、openfile.m、Output.m和output.dat覆盖从数据准备到结果落盘的完整链路。直流潮流的本质是用一组线性方程PBθ替代交流潮流非线性迭代固定电压幅值为1忽略电阻和支路对地电纳保留节点相角作为唯一状态量。因此它特别适合做断面有功校核、机组组合内嵌潮流约束、静态安全分析这类需要反复求解的场景。读者至少要有Matlab矩阵运算基础并知道节点有功注入等于发电机出力减负荷否则建议先补这两块。2. 直流潮流数据与导纳矩阵:data.m和FormY.m到底做了什么2.1 直流潮流的输入数据比交流潮流少得多交流潮流需要完整的母线参数、变压器变比、并联补偿、无功上下限直流潮流只需要基准容量、节点有功注入和支路电抗。下面表格可参考用于数据准备。参数交流潮流直流潮流说明电压幅值求解变量固定为1.0 p.u.直流潮流不考虑电压调节无功功率求解变量不参与节点P方程中无Q电阻用于Ybus可忽略线路损耗不计并联电纳用于Ybus忽略对地充电电容不计节点相角求解变量求解变量唯一状态量这个程序包里的data.m是一段脚本而不是函数它直接在工作区里写入bus和branch矩阵。常见写法如下% data.m - 三节点直流潮流算例 % bus 列定义编号 节点类型 有功出力 有功负荷 % 节点类型: 1参考节点, 2PV, 3PQ baseMVA 100; bus [ 1 1 0 0; 2 2 1.0 0; 3 3 0 0.8 ]; % branch 列定义送端 受端 电阻 电抗 branch [ 1 2 0 0.10; 1 3 0 0.20; 2 3 0 0.25 ];这里bus第二列的类型字段主要用于兼容交流潮流习惯直流潮流不需要区分PV和PQ只要知道每个不是参考节点的节点有功注入即可。第三列是有功出力第四列是有功负荷两者的差值才是净注入。单位推荐用标幺值baseMVA100时100MW对应注入1.0。如果从工程数据直接抄MW需要除以baseMVA否则求解结果会整体偏大。最简单的错误是把负荷填成正数导致节点注入方向反了正确做法是负荷在第四列而FormP.m里会用第三列减第四列。2.2 FormY.m组装的是电纳阵而不是导纳阵FormY.m从分支数据建立节点电纳矩阵。严格按照直流潮流假设应该忽略电阻和并联电纳因此这个函数通常有两种实现先形成复数节点导纳矩阵再取虚部或者直接按B1/X组装。后者更直观代码如下function B FormY(branch, nb) % 仅考虑支路电抗的直流潮流节点电纳矩阵 B zeros(nb, nb); for k 1:size(branch, 1) f branch(k, 1); t branch(k, 2); x branch(k, 4); if x 0 error(直流潮流不处理零阻抗支路请检查数据); end b 1 / x; B(f, f) B(f, f) b; B(t, t) B(t, t) b; B(f, t) B(f, t) - b; B(t, f) B(t, f) - b; end end逻辑说明每个循环只处理一条支路支路电抗x以标幺值给出倒数就是该支路的串联电纳b。对角元累加与该节点相连所有支路的电纳非对角元乘以负一。这个B矩阵在直流潮流中相当于“导纳”的角色但因为忽略了实部它是对称的实矩阵。参数上需要特别注意x必须是标幺值如果输入的是有名值欧姆还要先除以基准阻抗才能参与组装。为什么不用复数Ybus直流潮流本就不关心电阻损耗和并联充电功率用复数Ybus取imag虽然也能得到同样结果但会带来两个小风险一是如果Ybus中包含了非零电阻某些实现会错误地把实部保留进B矩阵二是复数运算多花一倍内存。直接按电抗组装更干净。2.3 参考节点的行和列必须删掉节点电纳矩阵B是奇异的因为所有行之和为0直接B\P会得到NaN或Inf。直流潮流的惯例是删去参考节点所在行和列保留N-1阶方程求解后再补参考节点的0角度。以2.1的三节点算例为例B矩阵如下 B FormY(branch, 3) B 15 -10 -5 -10 14 -4 -5 -4 9删去第一行第一列后得到Bp [14 -4; -4 9]行列式为110非奇异。Bp的物理含义是去掉参考点后的等效电纳映射任何节点角度叠加一个常数都不会改变支路有功因此必须让参考节点角度固定为0来消除自由度。程序里通常会这样取子矩阵ref find(bus(:, 2) 1); % 找到参考节点所在行号 Bp B; Bp(ref, :) []; Bp(:, ref) [];这五行代码本身不是函数放在DCPF.m中即可。注意find返回的是参考节点在bus矩阵中的行号不是节点编号本身如果节点编号不从1开始连续排列建议先重排节点编号再建矩阵否则删行删列时会错位。2.4 组装前的数据校验在更大的算例上我一般会在FormY后加两行校验assert(norm(B - B., inf) 1e-10, B矩阵不对称检查支路编号); assert(max(abs(sum(B, 2))) 1e-10, B矩阵行和应接近零);第一行检查对称性支路两端写反不会造成不对称但data.m里若出现重复支路或漏写支路对角元会异常。第二行检查B矩阵的零和特性它来自基尔霍夫定律每一行元素和表示节点注入对单位角度偏移的响应理论上应为0。如果这两项不通过问题几乎都在data.m的branch矩阵典型情况是节点编号越界或电抗填成了电阻。3. 直流潮流注入功率与相角求解:FormP.m、CalP.m和CalDelta.m的分工3.1 FormP.m只做一件事净注入向量直流潮流方程右侧是节点有功注入向量列向量维度为节点总数。data.m中的bus矩阵包含第三列出力Pg和第四列负荷PdFormP.m把它们按节点汇总成P向量function P FormP(bus, nb) % 返回节点净注入有功单位与data.m一致 P zeros(nb, 1); for i 1:nb P(i) bus(i, 3) - bus(i, 4); end end逻辑说明bus的第3、4列顺序必须与data.m注释一致。如果程序包实际版本中第3列是负荷、第4列是出力这里就是反向减法。参数说明nb由DCPF.m用size(bus,1)取得不传容易出现维度不一致。循环写法比向量化慢但对于几百节点完全可接受而且便于在循环内检查每条数据是否存在NaN。对于三节点示例节点注入如下表节点PgPdPnet100021.001.0300.8-0.8参考节点本身P值任意程序会忽略它。这里有一个常见误区有人会在FormP.m中把参考节点强制为0认为参考节点不参与计算。实际上参考节点有功注入求解后会自动等于其余节点注入之和的相反数这一步不用手写。若在直流潮流中强行令参考节点P0整个方程的功率平衡就会被破坏。3.2 CalP.m可能是“用当前角度算功率”的回代函数在这个压缩包中CalP.m和CalDelta.m容易混淆。建议把CalP.m理解为计算“由当前角度导出的功率注入”P_calc而不是形成外部注入P。因为形成外部P的工作已经由FormP.m完成CalP.m更合理的定位是function Pcalc CalP(B, theta) % 根据当前节点角度计算各节点注入有功 % 输入B为完整节点电纳矩阵theta为完整角度列向量 Pcalc B * theta; end逻辑说明直流潮流中PBθ是线性关系因此这是一个纯粹的矩阵乘法。参数说明B是n×n维的奇异矩阵theta必须包含参考节点的0元素否则结果没意义。这个函数在求解后的残差校验里很有用比如在DCPF.m里比较CalP(B,theta)与P的差别。由于我们忽略了电阻和电压幅值影响这种残差通常很小但如果网络中含大量重载支路直流近似误差会集中体现在有电阻的线路上。3.3 CalDelta.m求解角度而不是“修正量”程序命名为Delta容易让人误以为要做牛顿-拉夫逊式迭代修正。实际上直流潮流是线性方程一次求解直接得到角度不需要迭代。CalDelta.m的常见实现function [theta, Delta] CalDelta(Bp, Pn, ref) % Bp删去参考节点的电纳矩阵 % Pn删去参考节点的注入有功 % ref参考节点在完整角度向量中的位置 n length(Pn) 1; theta zeros(n, 1); idx true(n, 1); idx(ref) false; theta(idx) Bp \ Pn; Delta max(abs(Pn - Bp * theta(idx))); end逻辑说明第一步构建逻辑索引idx参考节点对应的位置被排除第二步theta(idx)就是非参考节点角度用左除求解第三步Delta是残差最大值用于判断求解质量。左除(Bp \ Pn)是Matlab求解线性方程组Axb的标准写法内部会做LU分解数值稳定性比inv(Bp)*Pn好得多也不容易受到矩阵病态影响。Delta理论上不会超过1e-10如果达到1e-4量级优先检查Bp和Pn是否用了不同基准容量。参数说明ref必须是对应完整角度向量的位置假设ref1时Bp删除第一行第一列Pn删除第一个元素。在三节点算例中Bp [14 -4; -4 9]; Pn [1; -0.8]; theta_noref Bp \ Pn;得到角度为0、0.0527、-0.0655弧度。角度很小说明直流潮流的线性化假设成立。如果某条支路两侧角度差超过0.3弧度结果已经严重偏离交流潮流此时不要硬用直流潮流而应退回交流潮流。3.4 为什么不需要牛顿迭代交流潮流之所以需要迭代是因为P和V、θ之间存在正弦和乘积关系。直流潮流把V固定为1把sinθ近似为θ方程退化为线性。即使程序里某个版本把CalDelta写成循环修正也只是重复“Bp\Pn”若干次属于冗余。真正需要迭代的场景是直流二次规划或考虑线路损耗修正那已经超出本程序范围。因此主流程里DCPF.m会直接把CalDelta的结果当作最终角度不会再做外层迭代。4. 把散文件串成主流程:DCPF.m、openfile.m和Output.m的运行逻辑4.1 DCPF.m按顺序调用所有函数理解了各文件后DCPF.m的工作就是把它们串成一条流水线装载数据→建B→算P→解角度→算支路潮流→写文件。一个可运行的DCPF.m骨架如下% DCPF.m clear; clc; data; % 在工作区生成 bus、branch、baseMVA nb size(bus, 1); B FormY(branch, nb); P FormP(bus, nb); ref find(bus(:, 2) 1); Bp B; Bp(ref, :) []; Bp(:, ref) []; Pn P; Pn(ref, :) []; [theta, residual] CalDelta(Bp, Pn, ref); % 支路有功: Pij (theta_i - theta_j) / x Pflow zeros(size(branch, 1), 1); for k 1:size(branch, 1) f branch(k, 1); t branch(k, 2); Pflow(k) (theta(f) - theta(t)) / branch(k, 4); end fid openfile(output.dat); Output(fid, bus, branch, theta, Pflow, baseMVA, residual); fclose(fid);逻辑说明data是脚本所以它定义的所有变量会留在当前工作区后面的函数调用可直接使用。Bp和Pn的删行操作必须放在求解前因为CalDelta内部不再处理参考节点。支路有功计算放在角度求解后公式里的电抗必须和FormY.m用同一个branch(:,4)。代码最后的fclose不能省略否则文件内容可能没有完全写入磁盘。4.2 openfile.m负责检查文件句柄openfile.m名字容易让人以为它读取文件实际在这套流程里它负责创建并打开输出文件function fid openfile(fname) fid fopen(fname, w); if fid -1 error(无法创建 %s请检查当前目录是否可写, fname); end end使用方式是在主脚本里传一个名字给它返回值fid是文件标识符。如果当前目录没有写权限fopen返回-1程序会立即报错而不是在后面的fprintf里弹一堆难懂的异常。参数说明w模式会覆盖已有output.dat如果需要保留历史结果可把模式改成a追加。该程序包默认输出output.dat所以openfile的实参通常是output.dat。4.3 Output.m写什么格式Output.m需要把节点角度和支路有功按人类可读的格式写入文件。一种常见实现function Output(fid, bus, branch, theta, Pflow, baseMVA, residual) fprintf(fid, 节点角度(弧度, 度)\n); for i 1:size(bus, 1) fprintf(fid, %4d %10.6f %10.4f\n, ... i, theta(i), theta(i) * 180 / pi); end fprintf(fid, 支路有功(MW)\n); for k 1:size(Pflow, 1) fprintf(fid, %4d %4d %12.4f\n, ... branch(k, 1), branch(k, 2), Pflow(k) * baseMVA); end fprintf(fid, 最大残差: %.2e\n, residual); end逻辑说明角度先输出弧度再输出度支路有功是标幺值乘以baseMVA得到MW。参数说明branch和baseMVA通过函数参数传入避免在函数体内引用工作区变量导致“未定义变量”。output.dat作为纯文本可直接被Python的pandas.read_csv用分隔符解析也可以直接把第一段粘贴到Excel做二次分析。输出结构对应关系如下输出段字段实际单位节点角度节点编号、弧度、度rad、deg支路有功送端、受端、PflowMW残差max residualp.u.4.4 跑不动时的三个检查点运行时最常见的报错是“未定义函数或变量FormY”。这通常不是程序写错而是Matlab当前路径没有包含这些文件需要在DCPF.m所在目录执行addpath(genpath(pwd))或者把当前文件夹切换正确。第二个常见问题是“矩阵维度不一致”出现在data.m中bus和branch行数不匹配比如bus只有3个节点但branch引用了节点4。此时先检查branch中max(max(branch(:,1:2)))是否等于size(bus,1)。第三个问题是“矩阵接近奇异”除了没用删除参考节点外还可能因为两条并联支路电抗相近导致Bp条件数变大。并联线路的Bp仍可计算但残差会从1e-15升到1e-12依然可用如果残差到1e-4多半是参考节点没删干净。5. 直流潮流算例替换:把data.m改成你自己的电网拓扑5.1 从交流潮流算例里取数而不是手填如果你自己没有一个现成的直流潮流数据文件最快的方式是从MATPOWER的case9中抽取。MATPOWER不是本压缩包的一部分但它输出的矩阵格式稳定常见做法是这样mpc loadcase(case9); baseMVA mpc.baseMVA; % 节点: 编号、类型、有功出力、有功负荷 bus mpc.bus(:, [1, 2]); bus [bus, zeros(size(bus, 1), 2)]; % 先占列 for g 1:size(mpc.gen, 1) id mpc.gen(g, 1); bus(id, 3) bus(id, 3) mpc.gen(g, 2) / baseMVA; % PG累加 end bus(:, 4) mpc.bus(:, 3) / baseMVA; % PD % 支路: 送端、受端、电阻(填0)、电抗 branch [mpc.branch(:, 1), mpc.branch(:, 2), ... zeros(size(mpc.branch, 1), 1), mpc.branch(:, 4)];逻辑说明MATPOWER的bus第3列是PD第4列是QDgen第2列是PG。这里把PG除以baseMVA转为标幺值PD同样处理。BUS_TYPE列保留下来只是为了维持data.m的列结构。参数说明如果原始case里多台机组挂在同一母线需要用累加而不是直接赋值上面用的是bus(id,3)不会覆盖已有值。5.2 不能用直流潮流硬算的情况直流潮流的实用性有边界换算例时要逐条检查场景直流潮流行为建议线路电阻不可忽略结果偏乐观忽略损耗导致断面功率偏低用交流潮流或只用于初筛电压过低或过高电压幅值偏离1造成误差检查潮流后电压若低于0.9p.u.慎用重载线路相角差大sinθ约等于θ失效相角差超过30度时换交流变压器变比非1需折算到统一基准将变比折算到电抗零阻抗支路1/x无穷合并节点或给一个很小x参数上当线路电抗x很小时直流潮流会算出很大的相角差结果不稳定。处理方式把两条并联支路用等值电抗x/2合并不要保留两个端点完全相同的支路。如果只是做灵敏度分析可以在当前算例上把负荷都按比例放大到120%观察角度是否线性放大若放大后角度分布与交流潮流明显偏离说明已到边界。5.3 在DCPF.m里追加CSV导出很多时候output.dat的格式不适合直接画图与其改Output.m的文件格式不如在DCPF.m里再补一段CSV导出fid2 fopen(result.csv, w); fprintf(fid2, bus_id,theta_rad,theta_deg,Pbus_mw\n); for i 1:length(theta) fprintf(fid2, %d,%.6f,%.4f,%.2f\n, ... i, theta(i), theta(i)*180/pi, P(i)*baseMVA); end fclose(fid2);代码说明这个循环把theta和P打包到result.csvP是之前计算出的节点注入向量若没有保留可在DCPF.m里复制一份。fprintf用逗号作为分隔符Excel直接双击可打开。请注意这里的P(i)如果是第i个节点的净注入在参考节点上会是待求的实际注入不应写死为0。5.4 替换算例后的前后一致性检查换完数据后我会先跑一遍原始三节点算例记录output.dat里的角度和Pflow再替换自己的数据。如果新结果里支路潮流量级与手动估算差10倍以上多半是baseMVA没乘回去或电抗用了欧姆值。标幺值电抗的折算方式是实际电抗欧姆乘以baseMVA除以基准电压平方这块是最容易出错的建议单独写一行注释放到data.m顶部。6. 验证与排错:用几行matlab命令复现直流潮流结果6.1 绕开整个程序包独立复算为了确认程序包结果不是“自己算自己”我会在Matlab命令行里用最原始的方式重算一遍三节点x12 0.10; x13 0.20; x23 0.25; B [ 1/x121/x13, -1/x12, -1/x13; -1/x12, 1/x121/x23, -1/x23; -1/x13, -1/x23, 1/x131/x23]; P [0; 1; -0.8]; theta B(2:end, 2:end) \ P(2:end); theta [0; theta]; Pflow12 (theta(1) - theta(2)) / x12; Pflow13 (theta(1) - theta(3)) / x13; Pflow23 (theta(2) - theta(3)) / x23;这段代码没有任何函数调用完全从定义出发。theta计算结果为0、0.0527、-0.0655弧度Pflow12约为-0.5273Pflow13约为0.3273Pflow23约为0.4727。与程序包output.dat对应栏目逐项比对若偏差超过1e-6说明FormY或FormP里出现了单位错误。6.2 用角度差判断是否该换算法直流潮流省时间但也省精度。在验证阶段我一般会计算所有支路的(θi-θj)弧度绝对值取其最大值。如果最大角度差超过0.3弧度约合17度直流潮流实际上已不适用。这时的“正确结果”不是调程序的参数而是回到交流潮流。另一种快速评估是把所得角度代入交流有功方程观察不考虑电压幅值时的误差量级误差主要出现在低电压母线和重载线路上这为后续进一步优化提供了对象。6.3 稀疏化让同一套程序能跑大网络FormY.m用zeros生成全矩阵节点数到5000时内存约200MB还能接受到5万节点就非常吃紧。常见做法是把B初始化为sparseB sparse(nb, nb); for k 1:size(branch, 1) f branch(k, 1); t branch(k, 2); b 1 / branch(k, 4); B(f, f) B(f, f) b; B(t, t) B(t, t) b; B(f, t) B(f, t) - b; B(t, f) B(t, f) - b; end后续赋值逻辑不变左除求解时Matlab会自动用稀疏LU分解。另一个节省内存的点是在DCPF.m里不再显式构造Bp而是用索引方式取子阵idx true(nb, 1); idx(ref) false; theta(idx, 1) B(idx, idx) \ P(idx);这样既删去了参考节点又避免了先复制再删行的临时变量。在data.m里把节点和支路数组定义为稀疏数据通常能把内存占用降一个量级。本文还有配套的精品资源点击获取