简介面向电力系统谐波分析场景的MATLAB程序包专注谐波潮流计算与谐波解耦算法适合电气工程专业学生、科研人员以及从事电能质量治理的工程师使用。当前电网中开关电源、整流器等非线性负载大量接入谐波畸变已成为影响电能质量的主要因素这套代码即可用于求解各节点谐波电压与电流分布辅助开展谐波评估与治理方案设计。压缩包共9个文件核心是main.m算法脚本同时配有8个txt格式的节点、负荷、线路参数等数据文件整体体积仅15KB结构清晰、易于上手。目前已有586人学习/下载。通过运行代码既能理解谐波解耦将复杂多变量问题拆分为单变量子问题的矩阵对角化思想也能结合不同规模数据完成基波与谐波联合求解在此基础上用户可以针对仿真结果安排滤波器配置、调整负荷分配从而提升系统稳定性和电能质量是一份兼顾理论学习与工程验证的实用资料。1. 谐波潮流计算到底在算什么从一次“变压器超温”说起我刚工作那会儿遇到一件怪事一台 10kV/0.4kV 干式变压器负荷率才 65%外壳却烫得能煎鸡蛋。拿电能质量分析仪一测谐波电流畸变率 THD 干到 32%5 次和 7 次谐波占大头。那时候我才意识到传统基波潮流算出来的“设备不超载”在谐波面前跟笑话差不多。所谓谐波潮流计算就是回答一个问题电网里的各次谐波功率从哪儿来、流到哪儿去、在哪个节点把电压畸变成什么样。而谐波解耦是这个计算里最关键的一步——把混在一起的 31 种谐波拆成 31 个互不干扰的独立网络逐次求解。这篇文章我打算把谐波潮流计算的原理、解耦方法、可复现的最小计算流程以及我踩过的几个坑一次讲完。2. 为什么谐波能解耦先理解频域隔离2.1 线性假设下的叠加原理谐波潮流计算能成立靠的是电力系统里绝大多数元件在稳态下满足线性或准线性假设。变压器、输电线路、电抗器、电容器这些元件的阻抗虽然是频率的函数但在同一个频率下它们就是线性阻抗。谐波源如整流器、变频器可以等效成一个恒流源或者恒功率源注入对应频率的电流。有了这个前提叠加原理就能用了含多个谐波源的系统某节点电压的 5 次谐波分量等于所有谐波源在 5 次频率下分别作用再叠加。我在实际项目里最深的体会是谐振分析比谐波源建模容易得多。谐振是线性问题阻抗扫频一下就出来了谐波源建模才是真正的黑匣子不同负载类型、不同控制策略谐波电流的幅值和相角都不一样。解耦计算本身是成熟技术真正的变量在上游的谐波源测量与拟合上。2.2 不同频率之间的弱耦合理论上谐波之间并非完全不耦合。比如变压器铁芯饱和会产生谐波且各次谐波之间有一定的交互电弧炉这类非线性源更是极端。但工程实践里除了专门研究铁磁谐振和电弧炉的特殊工况绝大多数电能质量评估项目都接受一个简化各次谐波之间解耦同一次谐波在不同节点之间的相互作用是主导项不同次谐波之间的交互忽略。这就把问题从“一个巨大的非线性方程组”降维成了“31 个或 50 个互不干扰的线性方程组”。基波潮流用牛顿拉夫逊或潮流计算软件求解然后对每一个整数次谐波单独构建导纳矩阵、单独求解节点电压。最后把基波和各次谐波电压做向量合成得到总的电压畸变。这个过程就是我们常说的谐波解耦。2.3 网格谐波变形与频率扫描说到“网格谐波变形”这个词我理解的是系统阻抗随频率变化呈现的响应形态。电力系统在不同频率下感抗随频率线性增加容抗随频率反比减小L-C 组合会在特定频率产生并联谐振或串联谐振。这个阻抗随频率的曲线形状就叫阻抗频率特性。5 次谐波时系统阻抗是感性还是容性决定了同样大的谐波电流注入节点电压畸变会被放大还是吸收。根据我的经验解耦后每个频率的网络就是一组线性方程 Z(f)·I(f)V(f)其中 Z 随频率变化。所以谐波潮流要做的就是把基波潮流计算网架结构拿过来修改每个元件的阻抗模型重新生成各次频率下的 Z 矩阵。解耦在这里的意义是不需要同时求解所有频率的耦合方程一次算一个频率极大降低计算规模。3. 从基波潮流到谐波潮流扩展与修改3.1 用确定性方法而不是概率方法谐波潮流计算有两类路径确定性方法和概率方法。确定性方法输入固定的谐波源发射水平输出确定性的各节点谐波电压概率方法考虑谐波源的波动性输入一组统计分布输出谐波电压的概率分布。我做项目时按以下原则选型如果做电能质量评估、滤波器设计校核、或者看某个新接入的整流器对周边节点的影响——用确定性方法就够了计算简单、结果可解释、工程师惯于用国标 GB/T 14549-93 限值去对比。如果是研究规模化电动汽车充电桩、光伏逆变器群接入后配电网谐波的整体分布——概率方法更贴合实际但建模工作量大、对实测数据的依赖高。对大多数落地场景确定性方法足够撑起评估结论。3.2 基波潮流网络参数的频率修正写谐波潮流程序时我习惯把基波潮流当成“第 0 次迭代”来复用。先在基波频率下完成潮流求解取得各节点基波电压幅值和相角。然后切换到谐波计算模式线路/变压器阻抗做频率修正电阻要计及趋肤效应感抗按频率线性放大容抗按频率反比缩小。发电机和电动机在谐波网络中简化成次暂态电抗对应的阻抗。负荷在谐波网络中常用静态阻抗模型近似取基波负荷功率 P、Q 反算一个等效阻抗再按频率修正。频率修正公式我放在下面说明里。这是最容易写出“看起来对、实际错”的一步因为修正系数取错了谐振点位置会偏好几百赫兹滤波器设计直接翻车。3.3 节点导纳矩阵的构建对每个谐波次数 h构建节点导纳矩阵对角线元素 Y(h,ii) 是节点 i 上所有支路导纳和所有接地导纳之和非对角线元素 Y(h,ij) 是连接节点 i 和 j 的支路导纳的负值支路导纳用频率修正后的阻抗求倒数矩阵规模等于节点数稀疏矩阵存储。然后谐波源在相应节点注入电流求解方程Y(h) · V(h) I(h)V(h) 是该次谐波各节点电压相量。构建导纳矩阵的核心代码用 Python 写出来大概是下面这个样子的import numpy as np from scipy.sparse import lil_matrix, csc_matrix from scipy.sparse.linalg import spsolve def build_y_matrix(buses, branches, h): 构建 h 次谐波的节点导纳矩阵。 buses: 节点数组, 每个元素为 (idx, base_kv) branches: 支路数组, 每个元素为 (from_bus, to_bus, r, x, b, type) type: line 线路, xfmr 变压器, reactor 电抗器 h: 谐波次数 返回值: Y 矩阵 (csc 稀疏格式), 以及节点编号到数组下标的映射 n len(buses) bus_index {b[0]: i for i, b in enumerate(buses)} Y lil_matrix((n, n), dtypenp.complex128) for frm, to, r, x, b in branches: i bus_index[frm] j bus_index[to] # 频率修正: 感抗按 h 倍放大, 电阻计及趋肤效应α 一般取 0.5~0.8 r_h r * (h ** 0.6) x_h x * h if h 1: # 谐波下线路对地电容基本不修正或按频率反比微调 b_h b / h else: b_h b z_series r_h 1j * x_h y_series 1.0 / z_series y_shunt 1j * b_h / 2 # 线路对地导纳均分到两端 Y[i, i] y_series y_shunt Y[j, j] y_series y_shunt Y[i, j] - y_series Y[j, i] - y_series return csc_matrix(Y)这段代码是我实际项目里简化后的版本。逻辑上主要做了三件事第一建立稀疏矩阵并用lil_matrix存放第二对每条支路的阻抗做频率修正第三把支路导纳填入相应位置注意对地导纳是均分到两端。csc_matrix是为了后续用spsolve求解。我把r_h的修正系数取 0.6这是钢芯铝绞线的常见取值如果你处理的是纯电缆线路趋肤效应更显著这个指数可以取到 0.8。在换成一个新系统时先取一个保守值再与实测结果对比微调。3.4 谐波源注入的等效模型谐波源建模直接决定谐波潮流计算的准确性这部分不能简单套用一个固定公式。我实际处理时把谐波源分成三类电流源型整流器、变频器、UPS 等输出特性近似于恒流源。工程上按实测特征频谱给出各次谐波电流的幅值和相角。6 脉动整流器的特征谐波是 6k±1 次12 脉动是 12k±1 次。电压源型部分 PWM 逆变器在谐波频率下更接近电压源。此时需要构建谐波 Norton 或 Thevenin 等效。不可控型电弧炉这类强非线性、非平稳负载谐波电流随机波动大实测数据只是统计代表值。我自己最常用的做法阶段性实测 查标准谱。实测可以获得准确的谐波电流频谱查标准谱如 IEEE 519 推荐的发射限值做规划计算。建模时把电流源注入放到对应节点注入幅值给实测均方根值相角给实测统计值。注意相角是谐波潮流计算里“玄学”的重灾区——取错了相角各源之间的叠加结果会完全不同甚至方向反转。4. 用 Python 在本地跑通最小谐波潮流计算4.1 一个最小可复现的算例与脚本结构我设计了一个最简单的 4 节点系统来演示谐波潮流计算流程一个 110kV 电源节点节点 1、一个 10kV 主变节点节点 2、一条配电线路带一个整流负荷节点节点 3、还有一个无功补偿电容器节点节点 4。电容器节点加进来的目的就是制造一个潜在的谐振条件——配电网里并联电容器是 5 次、7 次谐波放大的常见来源。最小脚本由三个文件构成bus_data.py放节点数据branch_data.py放支路数据harmonic_flow.py放主计算流程。脚本不需要做任何基波潮流迭代直接把基波电压设为 1.0pu 附近因为作为演示核心是展示谐波网络的构建、解耦求解和结果合成。# bus_data.py # (节点编号, 节点名称, 基准电压kV, 基波电压幅值pu, 基波电压相角度) buses [ (1, source, 110.0, 1.00, 0.0), (2, main_xfmr_10kv, 10.0, 1.02, -2.0), (3, rectifier_load, 10.0, 0.99, -5.5), (4, cap_bank, 10.0, 1.01, -3.2), ]branch_data.py里变压器阻抗按短路电压百分数折算线路按每公里参数折算电容器组按容量折算。# branch_data.py # (起点, 终点, 电阻Ohm, 感抗Ohm, 对地导纳S, 类型) branches [ (1, 2, 0.40, 12.0, 0.0, xfmr), (2, 3, 0.85, 1.20, 5e-6, line), (2, 4, 0.10, 5.00, 0.0, reactor), ]注意实际搭建时参数要在统一基准下折算到同一电压等级否则计算出来的导纳矩阵比例失衡。我这里写的参数已经折算到 10kV 侧。4.2 逐次谐波扫描与结果合成主计算流程harmonic_flow.py做了这么几件事先读两个数据文件然后循环扫描 3 到 19 次奇次谐波工程评估最常关注的次数对每次谐波调用build_y_matrix构建导纳矩阵再按照谐波源参数设置注入电流求解线性方程组获得各节点该次谐波电压。# harmonic_flow.py import numpy as np from scipy.sparse.linalg import spsolve from bus_data import buses from branch_data import branches from build_y import build_y_matrix # 整流负荷的谐波电流特征频谱(幅值相对于基波电流的标幺值, 相角/度) # h次谐波相角按典型6脉动整流器经验取值 harmonic_source { 3: (0.055, -108.0), 5: (0.230, 55.0), 7: (0.110, 62.0), 9: (0.025, -130.0), 11: (0.045, 27.0), 13: (0.033, 40.0), 17: (0.015, 15.0), 19: (0.012, 20.0), } base_current 100.0 # 整流器基波电流幅值, 单位: A def solve_one_harmonic(h, Iinject_pu): Y build_y_matrix(buses, branches, h) n len(buses) I np.zeros(n, dtypenp.complex128) src_idx 2 I[src_idx] Iinject_pu * base_current V spsolve(Y, I) return V # 主循环: 逐次谐波解耦求解, 并计算总畸变率 V_sum np.zeros((len(buses), len(harmonic_source)), dtypenp.complex128) for k, (h, (mag_per_unit, angle_deg)) in enumerate(harmonic_source.items()): I_mag mag_per_unit * base_current I_phase np.deg2rad(angle_deg) I_ph I_mag * np.exp(1j * I_phase) V solve_one_harmonic(h, I_ph) V_sum[:, k] V print(fh{h:2d}: V3{abs(V[2]):.4f}V, ang{np.angle(V[2], degTrue):.1f}deg) # 计算 THD: 各次谐波电压有效值之和与基波有效值的比值 V_base np.array([b[3] for b in buses]) * 10000.0 # 10kV节点基波相电压近似值 V_base[0] 110.0 * 10000.0 # 源节点 THD np.sqrt(np.sum(abs(V_sum) ** 2, axis1)) / V_base * 100.0 for i, b in enumerate(buses): print(fBus {b[1]}: THD {THD[i]:.2f}%)代码的逻辑不复杂但有几个点需要特别注意。第一个是spsolve直接求解了一次谐波方程如果你要算的节点数上万建议改用factorized对相同节点结构的矩阵做分解复用可以节省大量重复分解时间。第二个是注入电流相角的单位换算np.deg2rad不能漏。我在实际项目里踩过一次相角漏转弧度的坑算出来的 5 次谐波结果完全失真既偏大又偏小毫无规律排查了整整一个下午。第三个是V_base取了基波相电压的近似值严格做法是直接用基波潮流结果里各节点的实际基波电压相量但演示算例里把幅值近似为 1pu 足够说明流程。4.3 谐波解耦的计算域选择在前面这段代码里我用了相量域解耦每个频率单独建矩阵、单独求解频率之间不设交互。这是工程上最通用的做法速度快、内存小、适合大规模系统。另一个做法是模态域解耦也叫模态分析——通过特征值分解把谐波网络分解成独立模态用于识别谐振放大最严重的频率和位置。这个我在滤波器设计和谐振抑制时用得更多模态分析能把“哪个节点、哪个频率危险”直接告诉你比扫频曲线的眼力判断更可靠。但模态分析计算量比相量域大不适合频繁测试。表格对比一下这个差别方便选型解耦方式计算量输出形式适用场景相量域逐次扫描低各节点各次谐波电压常规谐波潮流、电能质量评估模态分解中高模态阻抗、关键参与节点谐振风险识别、滤波器设计时域波形高波形级谐波交互控制器交互研究、非特征谐波分析我做谐波潮流的时候80% 的场景用相量域就够剩下 20% 才是模态分析补上谐振风险排查。5. 谐波潮流计算翻车记录五条值得收藏的避坑经验5.1 频率修正系数不是对所有设备一视同仁现象同样的系统用不同软件算出来的谐波电压5 次谐波差 40%。原因变压器泄漏电感在谐波下的频率修正不同算法取不同系数线路电阻的趋肤效应修正有的软件取 0.5有的取 0.8。这些系数直接影响谐振点位置。解决在建模时统一约定修正规则并写明在报告里做滤波器设计时至少用两组系数做敏感性分析看谐振点偏移范围。5.2 电容器组的谐波放大不是算出来的现象一台并联电容器组投运后5 次谐波电流从 12A 涨到 71A电容器保险熔断而仿真结果却显示谐波电流只涨了 20%。原因仿真用的系统阻抗是 50Hz 基波阻抗外推的但实际系统在 250Hz 时存在一个未建模的并联谐振点。解决谐波潮流计算之前一定要先做 50Hz 到 2500Hz 的阻抗频率扫描看清楚系统在哪几个频率存在谐振峰再决定建模深度。这个步骤能省下来后面的返工。这是我从项目失败里收获的习惯之后每个项目都强制先扫描阻抗曲线。5.3 电容量的基准折算不一致导致结果量纲错乱现象某分区网谐波报告里节点电压畸变率 0.5%后来复核发现全算错了。原因变压器支路的数据用 110kV 侧欧姆值线路数据用 10kV 侧欧姆值没有统一折算到同一基准电压。解决程序里加一个断言检查所有支路阻抗的电压等级标记不匹配时直接报错中断。我用了一个笨办法——每条支路在数据文件里标记基准电压程序启动时逐个比对比值超出 1.05 倍就提示。5.4 相角问题现象谐波电压畸变率计算结果和实测差 3 倍。原因谐波电流相角给得不对。6 脉动整流器的 5 次谐波相角不是固定值它随触发角和直流侧电感变化不同工况可能差 120 度。解决从设计院拿到谐波源数据时第一件事确认相角的计量参考点如果厂家只给幅值不给相角那就用经验谱比如 5 次取 55 度、7 次取 62 度。千万别自己“猜”一个看起来舒服的相角。相角是谐波计算里最玄学、但偏偏又分量最重的输入参数。5.5 谐波潮流计算结果无法直接用于滤波器设计现象谐波潮流显示 5 次谐波电压超标PLC 直接按这个结果选了一个 5 次单调谐滤波器参数装上后 5 次谐波反而更大了。原因滤波器设计必须考虑滤波器本身接入后改变了系统阻抗阻抗频率特性从“原来无滤波器”变成了“有滤波器”谐振峰重新分布。显然不能拿没装滤波器的潮流结果来设计滤波器本身。解决选好滤波器参数后必须把滤波器作为固定支路重新构建导纳矩阵重跑一遍谐波潮流对比加装前后的谐波电压分布和滤波器自身电流。理想情况下这一步应该和滤波器优化迭代联合求解而不是单向串行。6. 往回收一步滤波器参数校验的三个关键量谐波潮流的最终目的往往是治理。滤波器设计好之后我习惯回到谐波潮流做三件事。第一滤波器支路电流核算——看流入滤波器的谐波电流是否超过其热额定第二滤波后剩余谐波电压是否满足国标限值同时对邻近节点也做一次快速复算第三偏移工况敏感性检查——系统阻抗在最大和最小运行方式下谐振点会跑谐波潮流结果必须包住这两个极端。这一套流程走下来谐波潮流计算的核心就是三句话逐次频率解耦、准确的谐波源注入、有效的边界校核。只要有可靠的基波潮流数据和一个能算稀疏线性方程组的脚本基本上任何配电网或输电网都能照这个框架跑起来。我在自己电脑上搭过一套纯 Python 的谐波评估工具读 Data、建矩阵、扫频半天能出结果。希望这个方法也能帮你把谐波评估从黑匣子变成手里可复算的脚本。如果你跑通了最小算例下一步可以试试往里面接入实测谐波频谱替换掉经验谱做出来的结果会更贴近现场。希望帮到你。本文还有配套的精品资源点击获取
