CN116467556B - 考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法 - Google Patents
考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法 Download PDFInfo
- Publication number
- CN116467556B CN116467556B CN202310316761.7A CN202310316761A CN116467556B CN 116467556 B CN116467556 B CN 116467556B CN 202310316761 A CN202310316761 A CN 202310316761A CN 116467556 B CN116467556 B CN 116467556B
- Authority
- CN
- China
- Prior art keywords
- harmonic
- impedance
- value
- mutation
- pcc
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 47
- 238000000034 method Methods 0.000 claims abstract description 120
- 230000008859 change Effects 0.000 claims abstract description 26
- 230000011218 segmentation Effects 0.000 claims abstract description 25
- 238000012216 screening Methods 0.000 claims abstract description 12
- 230000035772 mutation Effects 0.000 claims description 77
- 239000011159 matrix material Substances 0.000 claims description 23
- 230000008929 regeneration Effects 0.000 claims description 18
- 238000011069 regeneration method Methods 0.000 claims description 18
- 238000001514 detection method Methods 0.000 claims description 14
- 238000012360 testing method Methods 0.000 claims description 9
- 238000012417 linear regression Methods 0.000 claims description 7
- 238000006073 displacement reaction Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 5
- 238000005070 sampling Methods 0.000 claims description 5
- 230000008878 coupling Effects 0.000 claims description 4
- 238000010168 coupling process Methods 0.000 claims description 4
- 238000005859 coupling reaction Methods 0.000 claims description 4
- 238000007689 inspection Methods 0.000 claims description 4
- 238000013519 translation Methods 0.000 claims description 4
- 238000005315 distribution function Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 238000013459 approach Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 6
- 230000001172 regenerating effect Effects 0.000 abstract description 6
- 230000002411 adverse Effects 0.000 abstract description 5
- 238000004088 simulation Methods 0.000 abstract description 4
- 238000004458 analytical method Methods 0.000 abstract description 3
- 230000002159 abnormal effect Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000010998 test method Methods 0.000 description 4
- 239000003990 capacitor Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 239000000243 solution Substances 0.000 description 3
- 238000002347 injection Methods 0.000 description 2
- 239000007924 injection Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000010248 power generation Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000009774 resonance method Methods 0.000 description 1
- 238000000528 statistical test Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
- Y02E40/40—Arrangements for reducing harmonics
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Measurement Of Resistance Or Impedance (AREA)
Abstract
本发明公开了一种考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，首先利用时序聚类中的互相关系数原理使谐波电压和谐波电流曲线平移对齐并根据相似度筛选出背景谐波波动较小的数据段，再利用基于二元分割的Pettitt方法检验出谐波阻抗发生突变的时间点进行分段，最后使用再生权最小二乘方法对分段后的数据进行谐波阻抗计算。通过仿真分析验证了本发明方法能够较好减小背景谐波电压波动和阻抗变化给***谐波阻抗计算带来的不利影响，具有更广的适用范围。
Description
技术领域
本发明涉及电力谐波阻抗计算技术领域，具体涉及一种考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法。
背景技术
建设“以新能源为主体的新型电力***”成为了当前的研究热点，而新型电力***由于高比例新能源、高比例电力电子设备的接入使***谐波含量增加、频次增高，从而引起谐波分量放大、谐波分量超限的问题。准确的***谐波阻抗计算是解决谐波问题的前提。但是在电力***实际运行下，由于风力发电容易受到风速影响、光伏电源容易受到光照强度影响以及电力电子器件开断等原因，背景谐波电压会发生波动，而背景谐波电压波动会对谐波阻抗的精确计算带来困难。不仅如此，电力***运行方式发生变化、投切电容器组或无功补偿方式发生改变等原因都有可能导致***谐波阻抗发生变化，而目前传统计算方法对***谐波阻抗变化的研究相对较少。
目前传统的谐波阻抗计算方法主要分为干预式和非干预式两类。干预式方法利用人为产生的附加扰动量注入***或者改变***网络拓扑结构，向***注入谐波/间谐波电流、通断支路或启停负荷来估算***谐波阻抗。目前干预式***谐波阻抗计算的主要方法有：投切电容器法、晶闸管支路投切法、谐波电流源注入法、开关元件法、用户侧并联阻抗法和LCL谐振法等等。此类干预式法所获得的结果较为精确，但注入式方法注入的扰动有可能受到***的影响或者对***运行产生不利影响，而且此类方法的成本较高，因此没有被推广而得到广泛使用。
非干预式谐波阻抗计算方法指在不干扰***正常运行情况下，利用负荷或***本身的自然扰动及可测量参数等来计算谐波阻抗，这类方法对***正常运行不产生影响而被广泛采用。现有方法分为：波动量法、线性回归法、独立分量法和协方差法。波动法通过测量出一段时间内的谐波电压与谐波电流变化量，使用其比值近似替代谐波阻抗，方法使用简单，但该法仅适用于背景谐波相对稳定且谐波电压和谐波电流的测量精度较高的场合；回归法在***背景谐波基本稳定的情况下，以PCC处实际测量的谐波电压与谐波电流数据为基础，根据电力***谐波分析的等效电路构建回归方程，进而推算出***侧谐波阻抗，但这类方法在背景谐波波动较大或冲击性负荷较多时难以保持稳健性；独立分量法利用公共耦合点的谐波电压和谐波电流与***和用户侧谐波电流的关系方程来计算***谐波阻抗，但是源信号需要具有非高斯性；协方差法根据公共连接点处谐波电流与***背景谐波弱相关的特点，运用随机矢量协方差特性来抵消偏差量方程中背景谐波变动项，解得***侧谐波阻抗，适用于***侧谐波阻抗远小于用户侧的场合。
由于线性回归法在单谐波源和多谐波源的阻抗计算上同样适用而被广泛使用，但是当背景谐波大幅波动时必须先进行数据筛选。目前普遍认为***谐波阻抗是阶段性变化的，即只在某一时间点发生突变而在时间段内基本保持不变，目前的研究主要是通过检测突变点进行分段计算，但是并没有关注背景谐波大幅波动会给检测带来干扰，从而影响突变点的检测结果。
鉴于此，本发明提出了一种考虑背景谐波电压波动和阻抗变化的非干预式谐波阻抗计算方法，利用时序聚类中的互相关系数和基于二元分割的Pettitt方法筛选出背景谐波电压波动较小和***谐波阻抗相对恒定的数据段，再采用再生权最小二乘法稳健回归，减小异常值对计算结果的影响，该方法具有更广的应用范围。
发明内容
为解决现有技术中存在的问题，本发明提供了一种考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，能够较好减小背景谐波电压波动和阻抗变化给***谐波阻抗计算带来的不利影响，能够在背景谐波电压波动和***阻抗变化的情况下精确的计算***谐波阻抗，具有更广的适用范围，解决了上述背景技术中提到的问题。
为实现上述目的，本发明提供如下技术方案：一种考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，包括如下步骤：
S1、基于时序聚类中的互相关系数筛选出背景谐波电压波动较小的时段；
S2、基于二元分割的Pettitt方法检测出***谐波阻抗发生突变的时间点并分段；
S3、采用再生权最小二乘法进行***谐波阻抗计算。
优选的，在步骤S1中，所述基于时序聚类中的互相关系数筛选出背景谐波电压波动较小的时段具体是：将谐波电压U和谐波电流I分为n个子序列，对每个子序列采用时序聚类K-Shape方法中互相关系数进行筛选，计算两个子序列的距离测度SBD来表示其相似度的强弱。
优选的，所述步骤S1具体包括：
根据两序列互相关关系，固定公共耦合点谐波电压曲线时间窗，平移另一条谐波电流曲线/>时间窗，对谐波数据曲线/>和/>进行全局对齐：
其中，s为谐波电流曲线时间窗的平移量；im为谐波电流曲线/>上的数据点m；m为谐波电流曲线/>的数据长度。
根据谐波电流曲线的所有平移量s，s∈[-m,m]，得到谐波电流曲线/>和/>互相关序列：
其中，il为谐波电流曲线的数据点l；ul为谐波电压曲线/>的数据点l；w为电流曲线/>的平移量，w由-m到m依次递增；k＝w-m；k∈[-m,m]。
当达到最大值时，w的位置相对于/>的最佳位移量为s'＝w-m，将谐波数据曲线/>和/>互相关序列系数归一化后进行标幺值处理：
其中，为/>的最大值；/>为完全相似的谐波电流曲线/>不发生相对位移的互相关值；/>为完全相似的谐波电压曲线/>不发生相对位移的互相关值；的波动范围为[-1,1]。
根据得到SBD距离的计算公式：
其中，SBD值范围为[0,2]，SBD越小两曲线相似程度越高，筛选出SBD趋近于0的数据时段。
优选的，所述步骤S2具体包括：
S21、根据步骤S1已经筛选出背景谐波电压波动较小的时段，通过公式Upcc(n)＝ZsIpcc(n)+Us可知，当背景谐波电压波动较小时，***谐波阻抗由式(6)近似估计：
其中，n为采样点序号；Upcc(n)为第n个采样点的PCC处谐波电压；Ipcc(n)为第n个采样点的PCC处谐波电流；Zs为***谐波阻抗；Us为背景谐波电压。
S22、为减小计算量，将近似估计的谐波阻抗值均分为n小段，运用基于二元分割的Pettitt方法进行突变点检验；
S23、若没有检测出突变点或检测显著性概率p＞0.01，则认为在该时间段内，***谐波阻抗没有发生突变；
S24、若检测出突变点或检测显著性概率p＜0.01时，则认为该时间段内***谐波阻抗发生突变；
S25、步骤S24中是有一个突变点还是多个突变点需要再判断，从突变点将该时间序列分段；
S26、对分段后的数据序列再进行Pettitt突变检验，如果没有发现突变点，则步骤S24检测的突变点为真突变点；
如果检测出新突变点，则步骤S24中序列存在多个突变点，且检测出的突变点发生偏移，则回到步骤S25进行循环，直到检测出所有真突变点；
S27、将所有检测到的真突变点汇总排序，作为检测最终结果。
优选的，所述的Pettitt突变检验具体包括：将样本数据任意分为两个子序列，两个子序列均值相等，即表明当随机变量序列x1,x2,…,xn被任意分割成x1,x2,…,xm和xm+1,xm+2,…,xn两部分以后，如果两部分中随机变量的分布函数F1(x)≠F2(x)，则认为突变发生在m处。
优选的，所述步骤S3具体包括：通过诺顿电路等效模型可知***谐波阻抗计算公式为：
Upcc(n)＝ZsIpcc(n)+Us (7)
利用时序聚类中的互相关系数和基于二元分割的Pettitt方法筛选出的数据段谐波电压Upcc和Ipcc通过再生权最小二乘法计算出系数Zs和Us。
优选的，所述再生权最小二乘法的步骤包括如下：
1)把PCC处谐波电流作为自变量，谐波电压作为因变量可列出线性回归观测方程：
其中，表示观测值矩阵；β为观测方程的系数矩阵；/>为未知数矩阵；d为观测方程的常数项矩阵；X0为未知数的近似值；/>为未知数的值；/>
误差方程为：
其中，V表示观测值Y的残差；l为误差方程常数项；l＝-(d+βX0-Y)。
2)由最小二乘法得到未知数的解、观测值的残差和单位权方差的估值为：
V＝β(βTPβ)-1βTPl-l (11)
其中，P表示观测值Y的权阵，其主对角线元素为pj，初始值为用n表示观测值的数量，t表示未知数的数量，r表示多余观测值的数量(自由度)，即r＝n-t,r＞0。
3)将式(9)分为以下两部分，其中βt为t×t阶的满秩矩阵，通过对系数矩阵β进行线性变换来确定：
4)通过对式(13)、(14)变形可得：
Vr＝βrtVt-Wrt (16)
式中，
5)偶然误差的绝对值不会超过一定限值，使用来限制单位权真误差的范围，η称为单位权真误差取值范围的阈值系数，由式(16)可解得满足限制条件(18)的m(m→∞)组真误差的估值V(1),V(2),…,V(m)为：
其中，以/>为初值、/>为终值、/>为步长取值，/>则由式(16)确定，2θ+1是/>在区间/>中取值的节点数，θ则是/>在区间或/>中取值的节点数，即半截点数；
6)再生方差是用同一个观测值真误差的多个估值计算得到的方差，再生权是用观测值的再生方差计算得到的该观测值的权；不同观测值的再生方差、再生方差的平均值和再生权分别按式(19)、式(20)和式(21)计算：
本发明的有益效果是：本发明首先利用时序聚类中的互相关系数原理使谐波电压和谐波电流曲线平移对齐并根据相似度筛选出背景谐波波动较小的数据段，再利用基于二元分割的Pettitt方法检验出谐波阻抗发生突变的时间点进行分段，最后使用再生权最小二乘方法对分段后的数据进行谐波阻抗计算。能够较好减小背景谐波电压波动和阻抗变化给***谐波阻抗计算带来的不利影响，具有更广的适用范围。
附图说明
图1是本发明方法中的时序聚类中互相关系数方法滑动过程图；
图2是本发明方法中的诺顿等效电路图；
图3是本发明方法中的***谐波阻抗变化示意图；
图4是本发明方法中Pettitt检验突变点发生偏移示意图；
图5是本发明方法中基于二元分割的Pettitt方法分段检验出真突变点示意图；
图6是本发明提出的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法流程图；
图7是本发明实施例中PCC处5次谐波电压和电流幅值曲线图，图7(a)为谐波电压幅值曲线，图7(b)为谐波电流幅值曲线；
图8是本发明实施例中基于二元分割的Pettitt方法检验结果图，图8(a)为突变点在4秒，图8(b)为突变点在8秒，图8(c)为突变点在12秒；
图9是本发明实施例方法与其他各方法计算结果对比图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
本发明提出的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，首先利用时序聚类中互相关系数对公共连接点(Point of Common Connection,PCC)处谐波电压与谐波电流曲线进行平移对齐并检验其相似性，筛选出背景谐波电压比较稳定的时段，再对筛选后的数据利用基于二元分割的Pettitt方法进行谐波阻抗突变检验，根据检验出的真突变点进行分段，最后采用再生权最小二乘法分别计算各段的***谐波阻抗，本发明提出的方法具有更广的适用范围。
本发明提出的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，如图6所示，包括以下步骤：
Step1基于互相关系数筛选出背景谐波电压波动较小的时段；
Step2基于二元分割的Pettitt方法检测出***谐波阻抗发生突变的时间点并分段；
Step3采用再生权最小二乘法进行***谐波阻抗计算。
通过本发明的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，能够在背景谐波电压波动和阻抗变化的情况下比较精确的计算***谐波阻抗。
进一步的，Step1包括：将谐波电压U和谐波电流I分为n个子序列，对每个子序列采用时序聚类(K-Shape)方法中互相关系数进行筛选，计算两个子序列的距离测度(ShapeBased Distance,SBD)来表示其相似度的强弱。
在一个具体实施方式中，本发明的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法Step1中，运用时序聚类中互相关系数筛选背景谐波电压波动较小的数据段包括：
一般电力***比较复杂，可以将其***侧和用户侧等效为如图2所示的诺顿电路进行分析。在图2中，UPCC和IPCC分别代表PCC处的谐波电压和谐波电流，Is和Ic分别代表***侧和用户侧的等效谐波电流源，Zs和Zc分别代表***侧和用户侧等效谐波阻抗。
由下面关系式可知
UPCC＝ZsIPCC+Us (22)
当***侧谐波源波动较小时，PCC处的谐波电流IPCC大部分由用户侧Ic提供，此时UPCC和IPCC呈现很强的正相关性，由于Is和Ic在一段时间内是相互独立的，上述结论的逆结论依然成立，即如果UPCC和IPCC呈现很强的正相关性，背景谐波相对稳定，因此可以利用这种方法筛选出满足条件的数据段进行计算。
由于PCC处的谐波电压UPCC和谐波电流IPCC存在相位漂移，直接计算两者的相关性存在较大误差，故选择使用时序聚类中互相关系数进行曲线相似度检验来表示两者的相关性强弱较为准确。
根据两序列互相关关系，固定公共耦合点谐波电压曲线时间窗，平移另一条谐波电流曲线/>时间窗，如图1所示，对谐波数据曲线/>和/>进行全局对齐：
其中，il为谐波电流曲线的数据点l；ul为谐波电压曲线/>的数据点l；w为电流曲线/>的平移量，w由-m到m依次递增；k＝w-m；k∈[-m,m]。
根据谐波电流曲线的所有平移量s，s∈[-m,m]，得到谐波电流曲线/>和/>互相关序列：
其中，il为谐波电流曲线的数据点l；ul为谐波电压曲线/>的数据点l；k∈[-m,m]；w由-m到m依次递增。
当达到最大值时，w的位置相对于/>的最佳位移量为s'＝w-m，将谐波数据曲线/>和/>互相关序列系数归一化后进行标幺值处理：
其中，为/>的最大值；/>为完全相似的谐波电流曲线/>不发生相对位移的互相关值；/>为完全相似的谐波电压曲线/>不发生相对位移的互相关值；的波动范围为[-1,1]。
根据得到SBD距离的计算公式：
其中，SBD值范围为[0,2]，SBD越小两曲线相似程度越高，筛选出SBD趋近于0的数据段，即可减小背景谐波电压波动的影响而可以较为准确的计算出***谐波阻抗值。
在一个具体实施方式中，本发明的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法Step2中，运用基于二元分割的Pettitt方法检验谐波阻抗突变点包括：
实际上，电力***运行方式发生变化、投切电容器组或无功补偿方式发生改变等原因都可能导致***谐波阻抗发生变化，目前普遍认为***谐波阻抗是阶段性变化的，即只在某一时间点发生突变而在时间段内基本保持不变，如图3所示，即只需要检验发生突变的时间点并分段计算***谐波阻抗值。
Pettitt突变检验法是一种基于秩的非参数统计检验方法，计算简便，可以明确突变时间点，并且能够较好的识别序列的突变点。Pettitt突变检验法的基本原理是假设将样本数据任意分为两个子序列，两个子序列均值相等，即表明当随机变量序列x1,x2,…,xn被任意分割成x1,x2,…,xm和xm+1,xm+2,…,xn两部分以后，如果两部分中随机变量的分布函数F1(x)≠F2(x)，则认为突变发生在m处，其定义如下：
式中：
统计量Ut,N计算了第一个样本序列大于第二个样本序列的次数，则最显著突变点t处对应|Ut,n|的最大值，这里用Kn表示：
Kn＝Max1≤t≤n|Ut,n| (25)
对应的显著性概率可以通过下式近似计算得到：
由于Pettitt突变检验方法只能识别长序列中一个突变点，但谐波阻抗序列中可能存在多个突变点，而且如果谐波阻抗序列较大时，该方法的计算量较大。本发明基于二元分割法对Pettitt突变检验法进行改进，使其可以识别全部的突变点且减少计算量。二元分割法(Binary Segmentation,BS),常用于数据分段与图像分割领域，是一种迭代检测的处理方法。
进一步的，基于二元分割的Pettitt方法进行突变检验的具体步骤包括：
(1)由于Step1过程已经筛选出背景谐波电压波动较小的时段，通过公式Upcc(n)＝ZsIpcc(n)+Us可知，当背景谐波电压较小时，***谐波阻抗可由下式近似估计：
(2)为了减小计算量，将近似估计的谐波阻抗值均分为n小段，运用基于二元分割的Pettitt方法进行突变点检验。
(3)若没有检测出突变点或检测显著性概率p＞0.01，则认为在该时间段内，***谐波阻抗没有发生突变。
(4)若检测出突变点或检测显著性概率p＜0.01时，则认为该时间段内***谐波阻抗发生突变。
(5)步骤(4)中是有一个突变点还是多个突变点需要再判断，从突变点将该时间序列分段。
(6)对分段后的数据序列再进行Pettitt突变检验，如果没有发现突变点，则步骤(4)检测的突变点为真突变点，如果检测出新突变点，则步骤(4)中序列存在多个突变点，且检测出的突变点发生偏移，如图4中m点所示，则回到步骤(5)进行循环，直到检测出所有真突变点，如图5中i、j点所示。
(7)将所有检测到的真突变点汇总排序，作为检测最终结果。
在一个具体实施方式中，本发明的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法Step3中，采用再生权最小二乘法计算***谐波阻抗包括：
再生权最小二乘法(self-born weighted least squares method,SBWLS)是一种新的稳健回归估计方法，它利用观测值的残差条件方程提供的有效信息构造观测值的权，而不是像一般传统方法用最小二乘法得到的观测值残差构造观测值的权，当数据存在多个异常值时，相比一般稳健回归方法，再生权最小二乘法能够更有效地消除或减弱异常值对参数估计带来的影响。
把PCC处谐波电流作为自变量，谐波电压作为因变量可列出线性回归观测方程：
其中，表示观测值矩阵；β为观测方程的系数矩阵；/>为未知数矩阵；d为观测方程的常数项矩阵；X0为未知数的近似值；/>为未知数的值；/>
误差方程为：
其中，V表示观测值Y的残差；l为误差方程常数项；l＝-(d+βX0-Y)。
由最小二乘法可得未知数的解、观测值的残差和单位权方差的估值为：
V＝β(βTPβ)-1βTPl-l (11)
其中，P表示观测值Y的权阵，其主对角线元素为pj，初始值为 用n表示观测值的数量，t表示未知数的数量，r表示多余观测值的数量(自由度)，即r＝n-t,r＞0
将式(9)分为以下两部分，其中βt为t×t阶的满秩矩阵，可通过对系数矩阵β进行线性变换来确定。
将式(13)变形可得
将式(15)代入到式(14)中可得
Vr＝βrtVt-Wrt (16)
式中，
由于偶然误差的绝对值不会超过一定限值，所以用限制单位权真误差的范围，η称为单位权真误差取值范围的阈值系数。由式(16)可解得满足限制条件(18)的m(m→∞)组真误差的估值V(1),V(2),…,V(m)为：
在实际计算中，以/>为初值、/>为终值、/>为步长取值，/>则由式(16)确定。2θ+1是/>在区间/>中取值的节点数，θ则是/>在区间/>或/>中取值的节点数，称为半截点数。
再生方差是用同一个观测值真误差的多个(m个)估值计算得到的方差。再生权是用观测值的再生方差计算得到的该观测值的权。
不同观测值的再生方差、再生方差的平均值和再生权分别按式(19)、(20)和(21)计算。
式(16)称为再生权函数式，式(18)称为再生权函数限制条件。不同的参数估计问题具有不同的再生权函数式和相同的再生权函数限制条件。两个基本参数值阈系数η和半节点数θ通过仿真实验的方法确定。由于半节点数θ的限制，m不会是无穷大。βt和βr具有不唯一性，但是对参数估计结果没有显著影响。
由***谐波阻抗线性关系式可知
其中，为第i个谐波电压的估计值；ii为第i个谐波电流观测值；us为待求背景谐波电压；zs为待求***谐波阻抗。
由可得
ui+vi＝1·us+ii·zs (28)
其中，ui为第i个谐波电压真实值；vi为第i个谐波电压残差。
将上式矩阵化
在此线性回归模型中，n为观测值个数，其中待求量个数t＝2。一元线性回归模型误差方程如下：
对误差方程(9)中的系数矩阵β作线性变换，确定组成最大线性无关组的t＝2个误差方程式。组成最大线性无关组的2个误差方程不是唯一的，任取一组即可。
则该模型误差方程式的分块矩阵如下
由式(30)可得再生权函数式为：
Vr＝AtVt-W (32)
式中，Vr＝[vu3 … vun]T＝[v3 … vn]T,W＝-(Atlt-lr)，Vt＝[vu1 vu2]T＝[v1 v2]T。
以式(32)得到的观测值的再生权作为下一次迭代的观测值的权，用最小二乘法进行迭代计算，迭代的终止条件为相邻两次观测值残差小于0.1，则此时根据式(10)解得的背景谐波电压us和***谐波阻抗zs逼近真实值。
实施例2
以下通过一个具体实施例来说明本发明的技术效果。在MATLAB软件平台搭建如图2所示的等效电路模型，以五次谐波为例，仿真模型中谐波源与谐波阻抗参数设置如下所示，每0.02s测量一次谐波电压与谐波电流，仿真共16s。
谐波源：的幅值为120A，/>的幅值为12A，/>的相角为-30°，/>的相角为70°。在整个时段内给/>加上相对于其幅值大小±10％的随机扰动；其相角加上相对于相角大小±10％的随机扰动。给/>幅值加上±5％的随机扰动；给/>相角加上±5％的随机扰动。
谐波阻抗：为了模拟谐波阻抗发生突变，每4s大幅改变谐波阻抗值，0～4s***侧谐波阻抗Zs设为(15+j20)Ω，4～8s***侧谐波阻抗Zs设为(30+j40)Ω，8～12s***侧谐波阻抗Zs设为(45+j60)Ω，12～16s***侧谐波阻抗Zs设为(15+10)Ω。
测量PCC处的谐波电压和电流值，幅值曲线如图7(a)、图7(b)所示，将其分为n个子序列，运用时序聚类的互相关系数分别计算每个子序列的SBD值，并筛选出SBD小于0.2的数据段，然后用谐波电压比谐波电流近似代替谐波阻抗值，使用基于二元分割的Pettitt方法检验突变点。
检验结果如图8(a)、图8(b)、图8(c)所示。从图中可以清楚的看出谐波阻抗突变点发生在4s、8s、12s处，与上文设置的阻抗突变点相符，故本发明方法可以有效的识别谐波阻抗发生突变的时间点，然后根据突变点进行分段，使用再生权最小二乘法稳健回归计算各段的谐波阻抗值，计算结果和误差百分比如表1所示，由表1可知，本发明方法能够较为准确的计算出***谐波阻抗值。
表1
当数据存在多个异常值时，本发明采用的再生权最小二乘法相对于其他方法更具稳健性，利用波动量法(方法1)，稳健回归法(方法2)，独立分量法(方法3)，协方差法(方法4)，与本发明所提考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法(方法5)作对比，选择使用超过3个异常值的数据段，进行50次重复实验，计算结果如下图9。
从图9可知，本发明通过时序聚类互相关系数筛选的数据然后使用再生权最小二乘法计算***谐波阻抗的方法优于其他的四种方法，且对于存在多个异常值时，更具稳健性。
本发明提出考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，适用于单谐波源和多谐波源的阻抗计算。通过仿真结果可知，在背景谐波电压波动时，所提出的时序聚类的互相关系数能够将谐波电压和电流平移对齐筛选出满足条件的数据段，相对于传统的皮尔逊相关系数方法，能够更好的减小误差。在***阻抗变化时，所提出的基于二元分割的Pettitt方法能够有效筛选出谐波阻抗发生突变的时间点，效果优于一般检验方法。当存在多个异常值时，再生权最小二乘法相比其他方法更具稳健性。因此，本发明所提出的方法能够减小背景谐波电压波动和***谐波阻抗变化给阻抗计算带来的不利影响，具有更广的适用范围。
尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换,凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (4)
1.一种考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，其特征在于，包括如下步骤：
S1、基于时序聚类中的互相关系数筛选出背景谐波电压波动较小的时段；
S2、基于二元分割的Pettitt方法检测出***谐波阻抗发生突变的时间点并分段；
S3、采用再生权最小二乘法进行***谐波阻抗计算；
所述步骤S1具体包括：
根据两序列互相关关系，固定公共耦合点谐波电压曲线时间窗，平移另一条谐波电流曲线/>时间窗，对谐波数据曲线/>和/>进行全局对齐：
其中，s为谐波电流曲线时间窗的平移量；im为谐波电流曲线/>上的数据点m；m为谐波电流曲线/>的数据长度；
根据谐波电流曲线的所有平移量s，s∈[-m,m]，得到谐波电流曲线/>和/>互相关序列：
其中，il为谐波电流曲线的数据点l；ul为谐波电压曲线/>的数据点l；w为电流曲线/>的平移量，w由-m到m依次递增；k＝w-m；k∈[-m,m]；
当达到最大值时，w的位置相对于/>的最佳位移量为s'＝w-m，将谐波数据曲线/>和/>互相关序列系数归一化后进行标幺值处理：
其中，为/>的最大值；/>为完全相似的谐波电流曲线/>不发生相对位移的互相关值；/>为完全相似的谐波电压曲线/>不发生相对位移的互相关值；的波动范围为[-1，1]；
根据得到SBD距离的计算公式：
其中，SBD值范围为[0,2]，SBD越小两曲线相似程度越高，筛选出SBD趋近于0的数据时段；
所述步骤S2具体包括：
S21、根据步骤S1已经筛选出背景谐波电压波动较小的时段，通过公式Upcc(n)＝ZsIpcc(n)+Us可知，当背景谐波电压较小时，***谐波阻抗由式(6)近似估计；
其中，n为采样点序号；Upcc(n)为第n个采样点的PCC处谐波电压；Ipcc(n)为第n个采样点的PCC处谐波电流；Zs为***谐波阻抗；Us为背景谐波电压；
S22、为减小计算量，将近似估计的谐波阻抗值均分为n小段，运用基于二元分割的Pettitt方法进行突变点检验；
S23、若没有检测出突变点或检测显著性概率p＞0.01，则认为在该时间段内，***谐波阻抗没有发生突变；
S24、若检测出突变点或检测显著性概率p＜0.01时，则认为该时间段内***谐波阻抗发生突变；
S25、步骤S24中是有一个突变点还是多个突变点需要再判断，从突变点将该时间序列分段；
S26、对分段后的数据序列再进行Pettitt突变检验，如果没有发现突变点，则步骤S24检测的突变点为真突变点；
如果检测出新突变点，则步骤S24中序列存在多个突变点，且检测出的突变点发生偏移，则回到步骤S25进行循环，直到检测出所有真突变点；
S27、将所有检测到的真突变点汇总排序，作为检测最终结果；
所述步骤S3具体包括：通过诺顿电路等效模型可知***谐波阻抗计算公式为：
Upcc(n)＝ZsIpcc(n)+Us (7)
利用时序聚类中的互相关系数和基于二元分割的Pettitt方法筛选出的数据段谐波电压Upcc和Ipcc通过再生权最小二乘法计算出系数Zs和Us。
2.根据权利要求1所述的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，其特征在于：在步骤S1中，所述基于时序聚类中的互相关系数筛选出背景谐波电压波动较小的时段具体是：将谐波电压U和谐波电流I分为n个子序列，对每个子序列采用时序聚类K-Shape方法中互相关系数进行筛选，计算两个子序列的距离测度SBD来表示其相似度的强弱。
3.根据权利要求1所述的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，其特征在于：所述的Pettitt突变检验具体包括：将样本数据任意分为两个子序列，两个子序列均值相等，即表明当随机变量序列x1,x2,…,xn被任意分割成x1,x2,…,xm和xm+1,xm+2,…,xn两部分以后，如果两部分中随机变量的分布函数F1(x)≠F2(x)，则认为突变发生在m处。
4.根据权利要求1所述的考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法，其特征在于：所述再生权最小二乘法的步骤包括如下：
1)把PCC处谐波电流作为自变量，谐波电压作为因变量列出线性回归观测方程为：
其中，表示观测值矩阵；β为观测方程的系数矩阵；/>为未知数矩阵；d为观测方程的常数项矩阵；X0为未知数的近似值；/>为未知数的值；/>
误差方程为：
其中，V表示观测值Y的残差；l为误差方程常数项；l＝-(d+βX0-Y)；
2)由最小二乘法得到未知数的解、观测值的残差和单位权方差的估值为：
V＝β(βTPβ)-1βTPl-l (11)
其中，P表示观测值Y的权阵，其主对角线元素为pj，初始值为用n表示观测值的数量，t表示未知数的数量，r表示多余观测值的数量，即r＝n-t,r＞0；
3)将式(9)分为以下两部分，其中βt为t×t阶的满秩矩阵，通过对系数矩阵β进行线性变换来确定：
4)通过对式(13)、(14)变形可得：
Vr＝βrtVt-Wrt (16)
式中，
5)偶然误差的绝对值不会超过一定限值，使用来限制单位权真误差的范围，η称为单位权真误差取值范围的阈值系数，由式(16)可解得满足限制条件(18)的m(m→∞)组真误差的估值V(1),V(2),…,V(m)为：
其中，以/>为初值、/>为终值、/>为步长取值，/>则由式(16)确定，2θ+1是/>在区间/>中取值的节点数，θ则是/>在区间/>或中取值的节点数，即半截点数；
6)再生方差是用同一个观测值真误差的多个估值计算得到的方差，再生权是用观测值的再生方差计算得到的该观测值的权；不同观测值的再生方差、再生方差的平均值和再生权分别按式(19)、式(20)和式(21)计算：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310316761.7A CN116467556B (zh) | 2023-03-28 | 2023-03-28 | 考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310316761.7A CN116467556B (zh) | 2023-03-28 | 2023-03-28 | 考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116467556A CN116467556A (zh) | 2023-07-21 |
| CN116467556B true CN116467556B (zh) | 2023-12-15 |
Family
ID=87183417
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310316761.7A Active CN116467556B (zh) | 2023-03-28 | 2023-03-28 | 考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116467556B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116663431B (zh) * | 2023-07-28 | 2023-11-14 | 山东大学 | 基于分段仿射的风电场阻抗确定方法、***、设备及介质 |
| CN117630487B (zh) * | 2023-10-25 | 2024-06-07 | 国网甘肃省电力公司电力科学研究院 | 基于lof筛选-启发式分割算法的***谐波阻抗估计方法 |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106443285A (zh) * | 2016-11-30 | 2017-02-22 | 济南大学 | 基于总体最小二乘法的多谐波源谐波责任量化分析方法 |
| CN107565559A (zh) * | 2017-09-29 | 2018-01-09 | 西南交通大学 | 一种配电网负荷谐波贡献评估的分段有界约束优化方法 |
| EP3361272A1 (en) * | 2017-02-09 | 2018-08-15 | ABB Schweiz AG | Grid impedance estimation based on total least squares |
| CN109409762A (zh) * | 2018-11-07 | 2019-03-01 | 国网四川省电力公司成都供电公司 | 基于改进数据关联分析的谐波责任划分方法 |
| CN109709395A (zh) * | 2017-10-26 | 2019-05-03 | 华北电力大学（保定） | 基于双重筛选三点法的背景谐波源等效方法 |
| CN110059295A (zh) * | 2019-04-08 | 2019-07-26 | 国网江苏省电力有限公司电力科学研究院 | 一种基于三点筛选法的***谐波阻抗估计方法 |
| CN110456159A (zh) * | 2019-08-21 | 2019-11-15 | 四川大学 | 基于修正独立随机矢量的***侧谐波阻抗估计方法及*** |
| CN110850164A (zh) * | 2019-10-25 | 2020-02-28 | 华北电力大学 | 考虑背景谐波电压波动和阻抗变化的多谐波责任划分方法 |
| CN111693773A (zh) * | 2020-04-29 | 2020-09-22 | 国网江苏省电力有限公司电力科学研究院 | 一种基于互近似熵数据筛选的谐波源责任划分方法 |
| CN111898499A (zh) * | 2020-07-17 | 2020-11-06 | 东南大学 | 基于互近似熵与聚类的多谐波源谐波责任划分方法 |
| CN115343535A (zh) * | 2022-08-05 | 2022-11-15 | 国网冀北电力有限公司唐山供电公司 | 一种基于卡尔曼滤波的***谐波阻抗估计方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11221353B2 (en) * | 2018-07-06 | 2022-01-11 | Schneider Electric USA, Inc. | Systems and methods for analyzing power quality events in an electrical system |
| WO2021020985A1 (en) * | 2019-07-31 | 2021-02-04 | Schlumberger Canada Limited | A method and system for monitoring a wellbore object using a reflected pressure signal |
-
2023
- 2023-03-28 CN CN202310316761.7A patent/CN116467556B/zh active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106443285A (zh) * | 2016-11-30 | 2017-02-22 | 济南大学 | 基于总体最小二乘法的多谐波源谐波责任量化分析方法 |
| EP3361272A1 (en) * | 2017-02-09 | 2018-08-15 | ABB Schweiz AG | Grid impedance estimation based on total least squares |
| CN107565559A (zh) * | 2017-09-29 | 2018-01-09 | 西南交通大学 | 一种配电网负荷谐波贡献评估的分段有界约束优化方法 |
| CN109709395A (zh) * | 2017-10-26 | 2019-05-03 | 华北电力大学（保定） | 基于双重筛选三点法的背景谐波源等效方法 |
| CN109409762A (zh) * | 2018-11-07 | 2019-03-01 | 国网四川省电力公司成都供电公司 | 基于改进数据关联分析的谐波责任划分方法 |
| CN110059295A (zh) * | 2019-04-08 | 2019-07-26 | 国网江苏省电力有限公司电力科学研究院 | 一种基于三点筛选法的***谐波阻抗估计方法 |
| CN110456159A (zh) * | 2019-08-21 | 2019-11-15 | 四川大学 | 基于修正独立随机矢量的***侧谐波阻抗估计方法及*** |
| CN110850164A (zh) * | 2019-10-25 | 2020-02-28 | 华北电力大学 | 考虑背景谐波电压波动和阻抗变化的多谐波责任划分方法 |
| CN111693773A (zh) * | 2020-04-29 | 2020-09-22 | 国网江苏省电力有限公司电力科学研究院 | 一种基于互近似熵数据筛选的谐波源责任划分方法 |
| CN111898499A (zh) * | 2020-07-17 | 2020-11-06 | 东南大学 | 基于互近似熵与聚类的多谐波源谐波责任划分方法 |
| CN115343535A (zh) * | 2022-08-05 | 2022-11-15 | 国网冀北电力有限公司唐山供电公司 | 一种基于卡尔曼滤波的***谐波阻抗估计方法 |
Non-Patent Citations (5)
| Title |
|---|
| Study on Harmonic Impedance Estimation and Harmonic Contribution Evaluation Index;Ziteng Liu等;IEEE;第8卷;第59114页-59125页 * |
| 基于时序聚类及改进Pettitt的***谐波阻抗估计;李文涛等;电力***及其自动化学报;第1页-10页 * |
| 基于波形匹配的谐波责任划分方法;顾伟;邱海峰;尹香;陈兵;袁晓冬;王旭冲;;电力***自动化;第41卷(第02期);第129页-134页 * |
| 基于混合寻优算法及分离一致性判断的***谐波阻抗估计;林顺富等;电力***保护与控制;第50卷(第20期);第106页-116页 * |
| 改进变系数回归在***谐波阻抗估计中的应用;程卫健等;电力电容器与无功补偿;第44卷(第1期);第34页-48页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116467556A (zh) | 2023-07-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116467556B (zh) | 考虑背景谐波电压波动和阻抗变化的谐波阻抗计算方法 | |
| Su et al. | Power harmonic and interharmonic detection method in renewable power based on Nuttall double‐window all‐phase FFT algorithm | |
| CN103326358A (zh) | 基于同步相角测量装置的电力***动态状态估计方法 | |
| CN111596232B (zh) | 一种基于电压序列变点检测的户变关系校验方法 | |
| CN110333389A (zh) | 基于插值dft的正弦信号频率估计方法 | |
| CN106443285B (zh) | 基于总体最小二乘法的多谐波源谐波责任量化分析方法 | |
| CN110676940B (zh) | 一种提高参数辨识精度的配电网pmu配置方法及*** | |
| CN107453351B (zh) | 一种基于节点注入功率的配电网运行拓扑辨识方法 | |
| CN110086168A (zh) | 一种基于多源数据融合的配电网运行拓扑辨识方法 | |
| CN108199375A (zh) | 基于同步相量量测的智能配电网拓扑辨识方法 | |
| CN107658881A (zh) | 基于戴维南等值方法的电压稳定临界点判断方法 | |
| CN106980044B (zh) | 一种适应风电接入的电力***谐波电流估计方法 | |
| CN108054757A (zh) | 一种内嵌无功和电压的n-1闭环安全校核方法 | |
| CN111541255B (zh) | 基于动力学***的低频振荡模态识别方法及*** | |
| Li et al. | Grid-constrained data cleansing method for enhanced bus load forecasting | |
| Shu et al. | A new algorithm for calculating utility harmonic impedance | |
| CN104393588B (zh) | 一种基于最小均方误差的最优pmu布置方法 | |
| CN110210690B (zh) | 一种配电***微型同步相量测量单元优化配置方法 | |
| CN117786446A (zh) | 一种基于多判据融合的单相接地故障免阈值选段方法 | |
| Zheng et al. | Non-intrusive load monitoring based on the graph least squares reconstruction method | |
| CN102156783B (zh) | 电力***仿真精度综合评估方法 | |
| CN107069710B (zh) | 计及新能源时空相关性的电力***状态估计方法 | |
| CN117630487B (zh) | 基于lof筛选-启发式分割算法的***谐波阻抗估计方法 | |
| CN113553538A (zh) | 一种递推修正混合线性状态估计方法 | |
| Li et al. | Two-stage fault section location for distribution networks based on compressed sensing with estimated voltage measurements |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |