企业官网建设流程全解析

简介：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，覆盖从数据准备到结果落盘的完整链路。直流潮流的本质是用一组线性方程P=Bθ替代交流潮流非线性迭代：固定电压幅值为1，忽略电阻和支路对地电纳，保留节点相角作为唯一状态量。因此它特别适合做断面有功校核、机组组合内嵌潮流约束、静态安全分析这类需要反复求解的场景。读者至少要有Matlab矩阵运算基础，并知道节点有功注入等于发电机出力减负荷，否则建议先补这两块。

2. 直流潮流数据与导纳矩阵:data.m和FormY.m到底做了什么

2.1 直流潮流的输入数据比交流潮流少得多

交流潮流需要完整的母线参数、变压器变比、并联补偿、无功上下限，直流潮流只需要基准容量、节点有功注入和支路电抗。下面表格可参考用于数据准备。

这个程序包里的data.m是一段脚本而不是函数，它直接在工作区里写入bus和branch矩阵。常见写法如下：

% data.m - 三节点直流潮流算例 % bus 列定义：编号 节点类型 有功出力 有功负荷 % 节点类型: 1=参考节点, 2=PV, 3=PQ 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，只要知道每个不是参考节点的节点有功注入即可。第三列是有功出力，第四列是有功负荷，两者的差值才是净注入。单位推荐用标幺值，baseMVA=100时，100MW对应注入1.0。如果从工程数据直接抄MW，需要除以baseMVA，否则求解结果会整体偏大。最简单的错误是把负荷填成正数，导致节点注入方向反了；正确做法是负荷在第四列，而FormP.m里会用第三列减第四列。

2.2 FormY.m组装的是电纳阵而不是导纳阵

FormY.m从分支数据建立节点电纳矩阵。严格按照直流潮流假设，应该忽略电阻和并联电纳，因此这个函数通常有两种实现：先形成复数节点导纳矩阵再取虚部，或者直接按B=1/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和第四列负荷Pd，FormP.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。

对于三节点示例，节点注入如下表：

参考节点本身P值任意，程序会忽略它。这里有一个常见误区：有人会在FormP.m中把参考节点强制为0，认为参考节点不参与计算。实际上参考节点有功注入求解后会自动等于其余节点注入之和的相反数，这一步不用手写。若在直流潮流中强行令参考节点P=0，整个方程的功率平衡就会被破坏。

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

逻辑说明：直流潮流中P=Bθ是线性关系，因此这是一个纯粹的矩阵乘法。参数说明：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求解线性方程组Ax=b的标准写法，内部会做LU分解，数值稳定性比inv(Bp)*Pn好得多，也不容易受到矩阵病态影响。Delta理论上不会超过1e-10，如果达到1e-4量级，优先检查Bp和Pn是否用了不同基准容量。参数说明：ref必须是对应完整角度向量的位置，假设ref=1时，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做二次分析。输出结构对应关系如下：

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列是QD；gen第2列是PG。这里把PG除以baseMVA转为标幺值，PD同样处理。BUS_TYPE列保留下来只是为了维持data.m的列结构。参数说明：如果原始case里多台机组挂在同一母线，需要用累加而不是直接赋值，上面用的是bus(id,3)+，不会覆盖已有值。

5.2 不能用直流潮流硬算的情况

直流潮流的实用性有边界，换算例时要逐条检查：

参数上，当线路电抗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.csv，P是之前计算出的节点注入向量，若没有保留，可在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/x12+1/x13, -1/x12, -1/x13; -1/x12, 1/x12+1/x23, -1/x23; -1/x13, -1/x23, 1/x13+1/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.5273，Pflow13约为0.3273，Pflow23约为0.4727。与程序包output.dat对应栏目逐项比对，若偏差超过1e-6，说明FormY或FormP里出现了单位错误。

6.2 用角度差判断是否该换算法

直流潮流省时间，但也省精度。在验证阶段，我一般会计算所有支路的(θi-θj)弧度绝对值，取其最大值。如果最大角度差超过0.3弧度，约合17度，直流潮流实际上已不适用。这时的“正确结果”不是调程序的参数，而是回到交流潮流。另一种快速评估是把所得角度代入交流有功方程，观察不考虑电压幅值时的误差量级；误差主要出现在低电压母线和重载线路上，这为后续进一步优化提供了对象。

6.3 稀疏化，让同一套程序能跑大网络

FormY.m用zeros生成全矩阵，节点数到5000时内存约200MB，还能接受，到5万节点就非常吃紧。常见做法是把B初始化为sparse：

B = 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里把节点和支路数组定义为稀疏数据，通常能把内存占用降一个量级。

本文还有配套的精品资源，点击获取