1、算法的原理
谐波潮流计算作为分析配电网谐波的主要手段。其中，基于解耦法的谐波潮流计算是配电网谐波潮流主要的计算方法之一。基波潮流和谐波潮流在计算过程中存在复杂的耦合关系，可采用较为复杂的交替迭代谐波分析方法进行求解，但考虑到基波潮流对谐波潮流影响较大、谐波潮流对基波潮流影响较小的特点。 在实际工程应用方面，进行基波潮流计算时可以忽略谐波潮流的影响，在计算过程中将基波潮流和谐波潮流解耦分开处理，先运用牛顿-拉夫逊法迭代求解基波潮流，再运用高斯消元法求解谐波潮流。两者通过节点有功、无功功率方程建立关系，谐波潮流与基波潮流功率平衡方程表达式如下式所示：
基于解耦法谐波潮流算法，在实现基波潮流和谐波潮流完全解耦的同时，具 有线性法和非线性法两种算法的优点。在计算过程中不需反复对基波潮流和谐波 潮流进行更新迭代操作，而是在基波潮流结果的基础上求解谐波潮流。基波潮流和谐波潮流分别运用牛顿-拉夫逊法和高斯消元法方式求解，这样可以避免线性 计算过程中精度不高和非线性计算过程中运算量大的问题，该算法广泛应用于复 杂配电系统的谐波潮流分析。
2、求解步骤
步骤 1：将分布式电源类以及非线性负荷类谐波源作为 PQ 节点，运用牛顿-拉夫逊法迭代求解基波潮流。 根据配电系统的网络结构、线路参数、节点参数以及分布式电源和非线性负荷的接入位置。对于分布式电源，建立系统阻抗电路模型，在基波潮流计算时相当于在分布式电源接入的 PQ 节点并联一个等效导纳 Yz_0；对于非线性负荷，根据其功率大小更新其接入位置的 PQ 节点参数。综合考虑分布式电源类和非线性负荷类谐波源建立配电系统基波下的导纳矩 阵 Y1。潮流计算中牛顿-拉夫逊法的修正方程可以表示为：
式中，P、Q 分别表示配电系统节点的有功和无功功率；U、θ分别表示配电系统节点的电压幅值和相角。在配电系统中的任意第 z 节点，其有功和无功功率不平衡量 ΔPz 和 ΔQz 可以表示为：
对于配电系统中的 PV 节点，节点有功功率 Pz 和电压幅值 Uz 是已知的，需要求解的未知量是节点电压相角θz 和无功功率 Qz。利用 ΔPz 的修正方程求得节点电压幅值修正量 ΔUz 和电压相角修正量 Δθz；然后组成牛顿法计算潮流的迭代格式；循环修正、迭代过程以求解，直至收敛为止，即可获取配电系统的基波潮流结果。
步骤 2：根据配电系统分布式电源类和非线性负荷类谐波源接入位置，建立包含谐波源的系统谐波导纳矩阵 Yh。对配电网进行谐波分析时，可将分布式电源和非线性负荷均看作是一个向电网注入特定次谐波电流的整体，不考虑其复杂的内部结构与运行特性，用诺顿等效模型表示。分布式光伏等效阻抗模型，可以得到不同阶次谐波下的等效导纳 ，用以建立配电系统的谐波导纳矩阵 Yh 如下式：
式中：h 表示谐波阶次；下标表示节点编号；x 和 y 节点表示谐波源接入位 置。在实际工程应用中，考虑到谐波源接入节点的基波电压和谐波电压均会对谐 波电流造成影响，无论是分布式电源类谐波源还是非线性负荷类谐波源，可以用 诺顿等效模型进行等效处理，即用各类谐波源的谐波等效阻抗来模拟其谐波电流 的线性部分，用一个谐波电流源去模拟谐波电流的非线性部分。实际计算中可根 据各类谐波源的谐波特性，利用基波潮流结果中的基波电流计算各阶次谐波电流。 h 次谐波电流、h 次谐波电压以及 h 次谐波电流非线性部分的诺顿等效模型如下所示：
式中，Ih 表示谐波源注入电网 h 次谐波电流， I0_h表示谐波源 h 次谐波电流的非线性部分大小， Uh表示谐波源接入节点 h 次谐波电压,ZhZ_0 表示谐波源 h 次谐波阻抗。
步骤 3：利用高斯消元法求解谐波潮流，获取配电系统各节点的谐波电压。 将谐波源等效成谐波电流源，不同阶次谐波电流作用于系统谐波阻抗在配电网各节点产生谐波电压。节点谐波电流、谐波电压以及系统导纳关系如下式所示：
步骤 4：计算获取的各节点谐波电压得到谐波源在配电系统各节点的谐波电流修正量；
步骤 5：根据谐波电流的修正量和节点电压方程，得到更新后的系统各节点谐波电压，反复迭代计算直至满足收敛精度为止。
详细流程图如下：
解耦谐波潮流算法流程图
3、代码实现
部分代码：
%% 清空环境变量
clear; clc;
%% 基础设置
% 基准电压（kV）
UB = 12.66;
% 基准视在功率（MVA）
S_base = 100;
% 基准阻抗（Ohm）
Z_base = (UB^2) / S_base; % 1.6035 Ohm
%% 节点数据
% bus 矩阵
% 列说明：
% 1 - bus_i: 节点编号
% 2 - type: 节点类型 (3: 平衡节点, 1: PQ 节点)
% 3 - Pd: 有功负荷 (p.u.)
% 4 - Qd: 无功负荷 (p.u.)
% 5 - Gs: 导纳 (p.u.)
% 6 - Bs: 电纳 (p.u.)
% 7 - area: 区域号
% 8 - Vm: 电压幅值 (p.u.)
% 9 - Va: 电压相角 (度)
% 10 - baseKV: 基准电压 (kV)
% 11 - zone: 区域号
% 12 - Vmax: 电压最大值 (p.u.)
% 13 - Vmin: 电压最小值 (p.u.)
bus = [
1 3 0 0 0 0 1 1.0 0 12.66 1 1.1 0.9;
2 1 100 60 0 0 1 1.0 0 12.66 1 1.1 0.9;
3 1 90 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
4 1 120 80 0 0 1 1.0 0 12.66 1 1.1 0.9;
5 1 60 30 0 0 1 1.0 0 12.66 1 1.1 0.9;
6 1 60 20 0 0 1 1.0 0 12.66 1 1.1 0.9;
7 1 200 100 0 0 1 1.0 0 12.66 1 1.1 0.9;
8 1 200 100 0 0 1 1.0 0 12.66 1 1.1 0.9;
9 1 60 20 0 0 1 1.0 0 12.66 1 1.1 0.9;
10 1 60 20 0 0 1 1.0 0 12.66 1 1.1 0.9;
11 1 45 30 0 0 1 1.0 0 12.66 1 1.1 0.9;
12 1 60 35 0 0 1 1.0 0 12.66 1 1.1 0.9;
13 1 60 35 0 0 1 1.0 0 12.66 1 1.1 0.9;
14 1 120 80 0 0 1 1.0 0 12.66 1 1.1 0.9;
15 1 60 10 0 0 1 1.0 0 12.66 1 1.1 0.9;
16 1 60 20 0 0 1 1.0 0 12.66 1 1.1 0.9;
17 1 60 20 0 0 1 1.0 0 12.66 1 1.1 0.9;
18 1 90 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
19 1 90 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
20 1 90 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
21 1 90 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
22 1 90 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
23 1 90 50 0 0 1 1.0 0 12.66 1 1.1 0.9;
24 1 420 200 0 0 1 1.0 0 12.66 1 1.1 0.9;
25 1 420 200 0 0 1 1.0 0 12.66 1 1.1 0.9;
26 1 60 25 0 0 1 1.0 0 12.66 1 1.1 0.9;
27 1 60 25 0 0 1 1.0 0 12.66 1 1.1 0.9;
28 1 60 20 0 0 1 1.0 0 12.66 1 1.1 0.9;
29 1 120 70 0 0 1 1.0 0 12.66 1 1.1 0.9;
30 1 200 600 0 0 1 1.0 0 12.66 1 1.1 0.9;
31 1 150 70 0 0 1 1.0 0 12.66 1 1.1 0.9;
32 1 210 100 0 0 1 1.0 0 12.66 1 1.1 0.9;
33 1 60 40 0 0 1 1.0 0 12.66 1 1.1 0.9;
];
% 将 Pd 和 Qd 从 kW 和 kVAr 转换为 MW 和 MVAr，然后标幺化
bus(:,3) = bus(:,3) / 1000 / S_base; % Pd in p.u.
bus(:,4) = bus(:,4) / 1000 / S_base; % Qd in p.u.
%% 支路数据
% branch 矩阵
% 列说明：
% 1 - fbus: 起始节点
% 2 - tbus: 终止节点
% 3 - r: 电阻 (p.u.)
% 4 - x: 电抗 (p.u.)
% 5 - b: 电纳 (p.u.)
% 6 - rateA: 线路容量 (MVA)
% 7 - rateB: 线路容量 (MVA)
% 8 - rateC: 线路容量 (MVA)
% 9 - ratio: 变压器变比
% 10 - angle: 变压器角度
% 11 - status: 线路状态 (1: 投入, 0: 退出)
% 12 - angmin: 最小相角差
% 13 - angmax: 最大相角差
branch = [
1 2 0.0922 0.0470 0 0 0 0 0 0 1 -360 360;
2 3 0.4930 0.2511 0 0 0 0 0 0 1 -360 360;
3 4 0.3660 0.1864 0 0 0 0 0 0 1 -360 360;
4 5 0.3811 0.1941 0 0 0 0 0 0 1 -360 360;
5 6 0.8190 0.7070 0 0 0 0 0 0 1 -360 360;
6 7 0.1872 0.6188 0 0 0 0 0 0 1 -360 360;
7 8 0.7114 0.2351 0 0 0 0 0 0 1 -360 360;
8 9 1.0300 0.7400 0 0 0 0 0 0 1 -360 360;
9 10 1.0440 0.7400 0 0 0 0 0 0 1 -360 360;
10 11 0.1966 0.0650 0 0 0 0 0 0 1 -360 360;
11 12 0.3744 0.1238 0 0 0 0 0 0 1 -360 360;
12 13 1.4680 1.1550 0 0 0 0 0 0 1 -360 360;
13 14 0.5416 0.7129 0 0 0 0 0 0 1 -360 360;
14 15 0.5910 0.5260 0 0 0 0 0 0 1 -360 360;
15 16 0.7463 0.5450 0 0 0 0 0 0 1 -360 360;
16 17 1.2890 1.7210 0 0 0 0 0 0 1 -360 360;
17 18 0.7320 0.5740 0 0 0 0 0 0 1 -360 360;
2 19 0.1640 0.1565 0 0 0 0 0 0 1 -360 360;
19 20 1.5042 1.3554 0 0 0 0 0 0 1 -360 360;
20 21 0.4095 0.4784 0 0 0 0 0 0 1 -360 360;
21 22 0.7089 0.9373 0 0 0 0 0 0 1 -360 360;
3 23 0.4512 0.3083 0 0 0 0 0 0 1 -360 360;
23 24 0.8980 0.7091 0 0 0 0 0 0 1 -360 360;
24 25 0.8960 0.7011 0 0 0 0 0 0 1 -360 360;
7 26 0.2030 0.1034 0 0 0 0 0 0 1 -360 360;
26 27 0.2842 0.1447 0 0 0 0 0 0 1 -360 360;
27 28 1.0590 0.9337 0 0 0 0 0 0 1 -360 360;
28 29 0.8042 0.7006 0 0 0 0 0 0 1 -360 360;
29 30 0.5075 0.2585 0 0 0 0 0 0 1 -360 360;
30 31 0.9744 0.9630 0 0 0 0 0 0 1 -360 360;
31 32 0.3105 0.3619 0 0 0 0 0 0 1 -360 360;
32 33 0.3410 0.5302 0 0 0 0 0 0 1 -360 360;
];
% 将电阻和电抗从 Ohm 转换为标幺值
branch(:,3) = branch(:,3) / Z_base; % r in p.u.
branch(:,4) = branch(:,4) / Z_base; % x in p.u.
%% 接入新能源
% 指定接入新能源的节点
gen_buses = [1, 7, 8, 9, 14, 22, 24, 25, 30]; % 包括平衡节点 1
% 新能源的有功和无功出力（单位：p.u.）
Pg_gen = [0; 0.0010; 0.0010; 0.0008; 0.0012; 0.0009; 0.0020; 0.0020; 0.0010]; % 有功出力
Qg_gen = [0; 0.0002; 0.0002; 0.00016; 0.00024; 0.00018; 0.0004; 0.0004; 0.0002]; % 无功出力
4、仿真图像
节点谐波电压水平图
5、代码获取
代码获取地址：IEEE33节点谐波潮流计算
扫一扫
2216
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
