Joint probability distribution model construction method and harmonic sampling method

Abstract

The invention relates to a joint probability distribution construction method and a harmonic sampling method, and the method comprises the steps: obtaining load current data of different working conditions, enabling fundamental wave effective values corresponding to the load current data to be different, carrying out the classification statistics of the load current data, and classifying the load current data, so as to improve the accuracy of joint probability distribution. Each similar class is obtained, and joint probability distribution is established for each similar class. And decomposing the load current of each similar class to obtain a fundamental component and a harmonic component of each order corresponding to each similar class. According to the invention, density distribution estimation is carried out on each order harmonic component through the preset kernel density function to obtain the probability density distribution of each order harmonic component, so that the probability density distribution of each order harmonic component is prevented from being generally identified as obeying normal distribution according to experience, subjective identification of each order harmonic component is reduced, and the probability density distribution of each order harmonic component is improved. And the probability density distribution of the harmonic components of each order is closer to the distribution of a real situation.

Description

Technical Field

[0001] This disclosure relates to the field of rail transit traction system technology, and in particular to a method for constructing a joint probability distribution model and a method for sampling harmonics. Background Technology

[0002] Currently, electrified EMUs both domestically and internationally primarily employ AC-DC-AC traction systems. The current drawn from the traction power supply network by the EMU varies within a certain range under different operating conditions. Compared to AC-DC traction systems, AC-DC-AC traction systems have a higher power factor, thus reducing the low-order harmonic current content injected into the system. However, high-order harmonics still exist on the AC side. Simultaneously, as harmonic currents flow with the fundamental current through various power conversion stages, problems such as harmonic resonance and electromagnetic compatibility (EMC) issues arise in the traction system. Therefore, establishing a reasonable joint probability distribution model of harmonics in the EMU traction load provides rich reference information for harmonic identification, prediction, and suppression control algorithms in the traction system, possessing significant theoretical research value and practical engineering implications.

[0003] Existing technologies for addressing the harmonic distribution of traction system grid-side current primarily rely on statistical analysis based on measured data. However, these statistical analysis methods only analyze the probability distribution of specific current harmonics, failing to consider the interactions between different current harmonics. This leads to biases in providing reference information for subsequent harmonic identification, prediction, and suppression control in the traction system. Summary of the Invention

[0004] This disclosure provides a method for constructing a joint probability distribution model and a method for sampling harmonics, providing more accurate reference information for subsequent harmonic identification in traction systems.

[0005] Firstly, this disclosure provides a method for constructing a joint probability distribution model, including:

[0006] The load current data under different operating conditions are classified and statistically analyzed to classify the load currents under the same operating condition and whose fundamental effective values ​​differ within a specified range into similar categories.

[0007] For each similar class of load current, it is decomposed into fundamental component and harmonic components of each order;

[0008] Based on the preset kernel density function, the probability density distribution of each order harmonic component is determined;

[0009] The probability density distributions of each order harmonic component are sorted and classified, and the sorted marginal objective functions under each category are used as the first layer of the vine structure.

[0010] The copula function is used to correlate the two adjacent marginal objective functions in each category of the first layer to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, the copula function is used to correlate the two adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2 and n is a positive integer greater than m.

[0011] By using the copula function to correlate the joint probability distribution corresponding to each category, the joint probability distribution of each order harmonic component is obtained.

[0012] In some embodiments, the operating conditions include traction, braking, or coasting conditions.

[0013] In some embodiments, the step of sorting and classifying the probability density distributions of the harmonic components of each order to obtain the sorted marginal objective functions under each category as the first layer of the vine structure includes:

[0014] Calculate the correlation between the probability density distributions of the harmonic components of each order to obtain the correlation coefficient between the probability density distributions of the harmonic components of each order.

[0015] Based on the correlation coefficient, the probability density distribution of each order harmonic component is classified to obtain the marginal objective function under each category.

[0016] Based on the correlation coefficient, the marginal objective functions under the same category are sorted, and the sorted marginal objective functions under each category are used as the first layer of the vine structure.

[0017] In some embodiments, the step of determining the probability density distribution of each order harmonic component based on a preset kernel density function includes:

[0018] With the goal of minimizing the kernel density estimation error, the density distribution of each order harmonic component is calculated based on the preset kernel density function to obtain the probability density distribution of each order harmonic component; wherein, the kernel density estimation error is calculated based on the distribution error between the estimated density distribution and the actual density distribution.

[0019] In some embodiments, the kernel density estimation error is calculated based on a comparison coefficient between the integral of the square of the estimated density distribution and the integral of the estimated density distribution; wherein, when the comparison coefficient is greater than a first threshold, the kernel density estimation error is calculated based on the integral of the square of the estimated density distribution; when the comparison coefficient is less than or equal to the first threshold and greater than or equal to a second threshold, the kernel density estimation error is calculated based on the integral of the square of the estimated density distribution and the integral of the estimated density distribution; when the comparison coefficient is less than the second threshold, the kernel density estimation error is calculated based on the integral of the estimated density distribution.

[0020] Secondly, this disclosure provides a method for sampling harmonic currents, including:

[0021] The joint probability distribution of each order harmonic component obtained using the construction method described in any of the above embodiments;

[0022] Obtain the current object to be studied;

[0023] The correlation between time and current of the current object under study is analyzed to determine the periodic variation law of the current object under study.

[0024] Random sampling is performed on the joint probability distribution to obtain random current arrays corresponding to each preset current component;

[0025] Based on the periodic variation relationship, the sample current values ​​in each of the random current arrays are sorted to obtain harmonic random samples;

[0026] Based on the inverse transformation of the preset kernel density function, sample extraction is performed on the edge function constructed from the random harmonic samples to obtain the harmonic components contained in the current object under study.

[0027] Thirdly, this disclosure provides an apparatus for constructing the joint probability distribution of harmonic components of various orders, comprising:

[0028] The current classification module is used to classify and statistically analyze load current data under different operating conditions, so as to divide load currents under the same operating condition and whose fundamental frequencies differ within a specified range into similar categories.

[0029] The current decomposition module is used to decompose the load current of each similar class into the fundamental component and harmonic components of each order.

[0030] The harmonic density module is used to determine the probability density distribution of each order of harmonic components based on a preset kernel density function.

[0031] The sorting and classification module is used to sort and classify the probability density distribution of each order harmonic component, and obtain the sorted marginal objective function under each category as the first layer of the vine structure.

[0032] The density correlation module is used to correlate adjacent marginal objective functions in each category of the first layer using the copula function to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, the copula function is used to correlate adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2 and n is a positive integer greater than m.

[0033] The category association module is used to associate a joint probability distribution corresponding to each category using the copula function to obtain the joint probability distribution of each order harmonic component.

[0034] Fourthly, this disclosure provides a computer device including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method described in the preceding aspects.

[0035] Fifthly, this disclosure provides a computer-readable storage medium having a computer program stored thereon that, when executed by a processor, implements the steps of the methods described in the preceding aspects.

[0036] In a sixth aspect, this disclosure provides a computer program product, including a computer program / instructions that, when executed by a processor, implement the steps of the methods described in the foregoing aspects.

[0037] This disclosure provides a method for constructing a joint probability distribution and a method for sampling harmonics. It acquires load current data under different operating conditions. Since the fundamental RMS values ​​corresponding to the load current data may differ, to improve the accuracy of the joint probability distribution, the load current data is classified and statistically analyzed to obtain similarity classes. A joint probability distribution is then established for each similarity class. The load current of each similarity class is decomposed to obtain the fundamental component and harmonic components of each order corresponding to each similarity class. The density distribution of each harmonic component is estimated using a preset kernel density function to obtain the probability density distribution of each harmonic component. This avoids the assumption that the probability density distribution of each harmonic component follows a normal distribution based on experience, reducing subjective assumptions about the harmonic components and making the probability density distribution of each harmonic component closer to the actual distribution.

[0038] The probability density distributions of each harmonic component are sorted and classified. Harmonic components with high correlation are grouped into the same category, with highly correlated components grouped close together during sorting. The sorted marginal objective functions under each category form the first layer of the vine structure. Copula functions are used to correlate adjacent marginal objective functions within each category. Each pair of adjacent marginal objective functions correlated yields a joint probability distribution in the second layer. After all the marginal objective functions in the first layer of the vine structure are correlated, the joint probability distributions under each category in the second layer are obtained. Similarly, the copula function is used to correlate adjacent joint probability distributions in each category of the m-th layer, resulting in joint probability distributions for each category in the (m+1)-th layer of the vine structure. The number of joint probability distributions in the (m+1)-th layer is less than that in the m-th layer. During this process, the number of joint probability distributions for each category decreases layer by layer until a unique joint probability distribution is obtained for each category. Finally, the copula function is used to correlate the joint probability distributions for each category to obtain the joint probability distributions for each order harmonic component. Because the joint probability distributions for each order harmonic component can simultaneously consider the distribution of each order harmonic component and the mutual influence between them, the description of each order harmonic component is more accurate, providing more accurate reference information for subsequent harmonic identification, harmonic prediction, and harmonic suppression control of the traction system. Attached Figure Description

[0039] The present disclosure will be described in more detail below based on embodiments and with reference to the accompanying drawings:

[0040] Figure 1 A flowchart illustrating a method for constructing a joint probability distribution model of harmonic components of various orders, provided in an embodiment of this disclosure;

[0041] Figure 2 A schematic diagram of the structure of a device for constructing a joint probability distribution model of harmonic components of various orders provided in an embodiment of this disclosure;

[0042] Figure 3 The present invention provides an embodiment of obtaining a fundamental frequency and total current harmonic diagram from the original load current through a low-pass filter;

[0043] Figure 4 This disclosure provides a current harmonic diagram for each frequency range after wavelet transform of the original load current in an embodiment of the present disclosure.

[0044] Figure 5 This is a schematic diagram illustrating the frequency and time domain correspondence of current harmonics obtained by passing the load current through a filter, according to an embodiment of this disclosure.

[0045] Figure 6This is a schematic diagram of wavelet transform corresponding to the current frequency domain interval division provided in an embodiment of the present disclosure;

[0046] Figure 7 A time-domain signal diagram of obtaining specific sub-current harmonics by passing a load current through a filter, provided in an embodiment of this disclosure;

[0047] Figure 8 This embodiment of the present disclosure provides a time-domain diagram of the current harmonics corresponding to the frequency domain interval obtained by wavelet transform of the load current;

[0048] Figure 9 A schematic diagram of a cannonical vine structure for an improved Pair-Copula method provided in this disclosure embodiment;

[0049] Figure 10-A A schematic diagram illustrating an overall implementation framework provided in this disclosure.

[0050] Figure 10-B This is another schematic diagram illustrating an overall implementation framework provided by an embodiment of this disclosure.

[0051] In the accompanying drawings, the same parts are referred to by the same reference numerals, and the drawings are not drawn to scale. Detailed Implementation

[0052] To enable those skilled in the art to better understand the technical solutions of this disclosure, and to fully understand and implement the process of how this disclosure applies technical means to solve technical problems and achieve corresponding technical effects, the technical solutions in the embodiments of this disclosure will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this disclosure, not all embodiments. The embodiments of this disclosure and the various features within them can be combined with each other without conflict, and the resulting technical solutions are all within the protection scope of this disclosure. All other embodiments obtained by those skilled in the art based on the embodiments of this disclosure without creative effort should fall within the protection scope of this disclosure.

[0053] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this disclosure are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this disclosure described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0054] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0055] Existing patents and papers on the distribution of harmonic currents on the grid side of traction systems primarily rely on statistical analysis based on measured data. Firstly, they fail to describe the probability distribution of each current harmonic within a certain fundamental current range. Furthermore, most literature assumes an empirical normal distribution for each harmonic, introducing subjectivity and leading to significant discrepancies between the harmonic distribution function and actual conditions. Secondly, accurate extraction of harmonic currents of each order increases the difficulty of designing ideal low-pass and band-pass filters. Thirdly, they fail to establish a joint probability distribution of current harmonics on the inverter side of the traction system, considering only the probability distribution of specific harmonics. Fourthly, they lack analysis of the cross-correlation between different harmonics, such as the correlation between lower-order harmonics and between lower and higher-order harmonics, resulting in discrepancies between the obtained harmonic current probability distribution and actual conditions. Therefore, researching a joint probability distribution model for train current harmonics, reasonably assessing the impact of current harmonics generated by traction loads on the grid-side power quality of the traction system, and verifying the effectiveness of harmonic suppression control algorithms have practical value.

[0056] The paper "Research on Harmonic Distribution Characteristics and Probabilistic Model of High-Speed ​​EMU Based on Measured Data" (Journal of Railway Engineering, Vol. 32, No. 3, 2010) establishes a current harmonic model for the CRH2-200 high-speed EMU based on a large amount of measured data, using curve fitting and probabilistic statistical methods to simulate the harmonic probability distribution characteristics of the EMU affected by numerous random factors. This model, based on measured data, considers not only the variation of harmonic content and phase angle with the fundamental current, but also the probability distribution characteristics caused by the influence of many random factors. Although this paper establishes a current harmonic probability density model based on measured EMU data, it still assumes that the current harmonics follow a normal distribution for parameter estimation, resulting in differences between the obtained probability density functions of each current harmonic and the actual situation. Furthermore, the extracted harmonic currents are highly dependent on high-precision low-pass and band-pass filters, which cannot guarantee the purity and accuracy of the extracted harmonics. Although this paper mainly establishes low-order and high-order current harmonic probability models, the modeling process does not consider the cross-correlation between current harmonics of different orders, nor does it construct a joint distribution function for current harmonics. Therefore, it still differs from a true reflection of the load current harmonic distribution. Consequently, current harmonic probability models established using this method will have a certain impact on actual current harmonic identification, harmonic prediction, or harmonic suppression control algorithms and parameter adjustments. This will cause simulation software and real-time control programs to fail to achieve the expected results in practical applications, and may even fail to meet EMC and other standard requirements.

[0057] Example 1

[0058] Figure 1 This is a flowchart illustrating a method for constructing a joint probability distribution model of harmonic components of various orders, as provided in an embodiment of this disclosure. Figure 1 As shown, a method for constructing a joint probability distribution model of harmonic components of various orders includes:

[0059] Step 110: Classify and statistically analyze the load current data under different operating conditions to classify the load currents under the same operating condition and whose fundamental effective values ​​differ within a specified range into similar categories.

[0060] In this embodiment, the operating conditions and fundamental RMS values ​​of the collected load current data may differ. When the operating conditions are different or the fundamental RMS values ​​differ significantly, the intrinsic correlation between the load current data is weak, making it difficult to uncover the inherent relationships between the load current data and establish a joint probability distribution. Therefore, it is necessary to classify and statistically analyze the load current data, dividing load currents under the same operating condition and whose fundamental RMS values ​​differ within a specified range into similar categories, and establishing joint probability distributions for different similar categories. The specified range is the range of differences in the fundamental RMS values, used to filter out load current data whose fundamental RMS values ​​differ within the specified range. When these load current data belong to the same operating condition, they can be divided into similar categories.

[0061] Step 120: For each similar class of load current, decompose it into fundamental component and harmonic components of each order.

[0062] In this embodiment, in order to provide more accurate reference information for harmonic identification, harmonic prediction, and harmonic suppression control of the subsequent traction system, the load current in the similar class is decomposed to obtain the fundamental component and harmonic components of each order, so as to analyze the harmonic components of each order and the joint probability distribution of the harmonic components of each order.

[0063] Step 130: Determine the probability density distribution of each order harmonic component based on the preset kernel density function.

[0064] In this embodiment, based on a preset kernel density function, kernel density estimation is performed on each harmonic component to obtain its probability density distribution. In this step, it is not assumed that all harmonic components follow a normal distribution; rather, kernel density estimation is performed specifically for each harmonic component to ensure that the resulting probability density distribution accurately reflects the distribution of each harmonic component. The preset kernel density function is used to perform kernel density estimation on each harmonic component, calculating the kernel density estimate for each harmonic component at each current value, thus obtaining the probability density distribution of each harmonic component over the current value range.

[0065] Step 140: Sort and classify the probability density distributions of each order harmonic component, and obtain the sorted marginal objective functions under each category as the first layer of the vine structure.

[0066] In this embodiment, to improve the utilization of the probability density distribution of each order harmonic component, the probability density distributions of each order harmonic component are divided into the same category if they meet the preset conditions. For example, the probability density distributions of each order harmonic component with a correlation coefficient greater than a certain value are divided into one category, and the probability density distributions of each order harmonic component with a correlation coefficient less than a certain value are divided into different categories. This allows the probability density distributions of each order harmonic component with a large correlation coefficient under the same category to be combined in the vine structure, effectively utilizing the information of the inherent correlation between the probability densities of each order harmonic component, and constructing a more accurate joint probability distribution of each order harmonic component. In a vine structure, the marginal objective functions or joint probability distributions of the current layer are correlated with each other. The resulting joint probability distribution becomes the next layer of the vine structure. This next layer is then treated as the current layer, and the joint probability distributions of that layer are correlated with each other, resulting in the next layer of the vine structure. This process continues, and as the number of layers increases, the number of joint probability distributions in each layer gradually decreases, eventually leading to a joint probability distribution that combines all marginal objective functions. The vine structure is a method for building multivariate statistical models. It decomposes the correlation between multiple variables into two parts: uncorrelated and univariate correlated, to better describe the dependency structure between variables.

[0067] Step 150: Use the copula function to correlate the two adjacent marginal objective functions in each category of the first layer to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, use the copula function to correlate the two adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2 and n is a positive integer greater than m.

[0068] In this embodiment, marginal objective functions under different categories are correlated. In the first layer, the correlation between marginal objective functions within the same category is relatively high. The copula function is used to correlate adjacent marginal objective functions within the same category, and the resulting joint probability distribution more accurately reflects the correlation between the two marginal objective functions. However, the correlation between marginal objective functions from different categories is relatively low, and the joint probability distribution obtained from these correlations does not accurately describe the correlation between the two marginal objective functions. Therefore, in this step, marginal objective functions from different categories are not correlated. This improves the extraction rate of information from the previous layer of the vine structure by the joint probability distribution obtained after correlation, which is beneficial for improving the extraction of effective information from the first to the second layer and from the m to the (m+1)th layer of the vine structure.

[0069] It should be understood that before or during the association process, there are multiple joint probability distributions corresponding to each category. Each time the joint probability distribution corresponding to each category is associated by the copula function, the number of joint probability distributions obtained in the next layer will decrease until a unique joint probability distribution corresponding to each category is obtained.

[0070] For example, after obtaining the joint probability distributions for each category in the second layer, the copula function is used to correlate adjacent joint probability distributions in each category of the second layer to obtain the joint probability distributions for each category in the third layer of the vine structure. Subsequently, the copula function is used to correlate adjacent joint probability distributions in each category of the third layer to obtain the joint probability distributions for each category in the fourth layer of the vine structure... Following the above process, the correlation of joint probability distributions in each layer is performed sequentially until a unique joint probability distribution corresponding to each category is obtained.

[0071] Step 160: Use the copula function to correlate the joint probability distribution corresponding to each category to obtain the joint probability distribution of each order harmonic component.

[0072] In this embodiment, as each category is interconnected layer by layer, the number of joint probability distributions at each layer gradually decreases until a unique joint probability distribution for that category is obtained. To comprehensively consider the interaction between the joint probability distributions corresponding to different categories, a copula function is used to associate the joint probability distribution corresponding to each category, resulting in the joint probability distribution of each order harmonic component. The joint probability distribution of each order harmonic component can reflect the correlation between the harmonic components of each order in the corresponding similar classes, and can more accurately reflect the differences between the harmonic components of each order. The joint probability distribution of each order harmonic component can be used for subsequent research on harmonic identification, harmonic prediction, and harmonic suppression control of the traction system.

[0073] This disclosure provides a method for constructing a joint probability distribution and a method for sampling harmonics. It acquires load current data under different operating conditions. Since the fundamental RMS values ​​corresponding to the load current data may differ, to improve the accuracy of the joint probability distribution, the load current data is classified and statistically analyzed to obtain similarity classes. A joint probability distribution is then established for each similarity class. The load current of each similarity class is decomposed to obtain the fundamental component and harmonic components of each order corresponding to each similarity class. The density distribution of each harmonic component is estimated using a preset kernel density function to obtain the probability density distribution of each harmonic component. This avoids the assumption that the probability density distribution of each harmonic component follows a normal distribution based on experience, reducing subjective assumptions about the harmonic components and making the probability density distribution of each harmonic component closer to the actual distribution.

[0074] The probability density distributions of each harmonic component are sorted and classified. Harmonic components with high correlation are grouped into the same category, with highly correlated components grouped close together during sorting. The sorted marginal objective functions under each category form the first layer of the vine structure. Copula functions are used to correlate adjacent marginal objective functions within each category. Each pair of adjacent marginal objective functions correlated yields a joint probability distribution in the second layer. After all the marginal objective functions in the first layer of the vine structure are correlated, the joint probability distributions under each category in the second layer are obtained. Similarly, the copula function is used to correlate adjacent joint probability distributions in each category of the m-th layer, resulting in joint probability distributions for each category in the (m+1)-th layer of the vine structure. The number of joint probability distributions in the (m+1)-th layer is less than that in the m-th layer. During this process, the number of joint probability distributions for each category decreases layer by layer until a unique joint probability distribution is obtained for each category. Finally, the copula function is used to correlate the joint probability distributions for each category to obtain the joint probability distributions for each order harmonic component. Because the joint probability distributions for each order harmonic component can simultaneously consider the distribution of each order harmonic component and the mutual influence between them, the description of each order harmonic component is more accurate, providing more accurate reference information for subsequent harmonic identification, harmonic prediction, and harmonic suppression control of the traction system.

[0075] Example 2

[0076] Based on the above embodiments, a method for constructing the joint probability distribution of harmonic components of each order is provided, which includes:

[0077] Step 1: Classify and statistically analyze the load current data under different operating conditions to classify load currents under the same operating condition that differ in fundamental effective values ​​within a specified range into similar categories.

[0078] In one embodiment, the operating conditions include traction, braking, or coasting conditions. In this embodiment, in order to more accurately construct the joint probability distribution of each order harmonic component and obtain load current data under traction, braking, or coasting conditions, the data distribution of each order harmonic component in the load current data differs under different operating conditions. Obtaining load current data under different operating conditions is beneficial for providing richer statistical data. By constructing the joint probability distribution of each order harmonic component for different operating conditions, the joint probability distribution of each order harmonic component corresponding to different operating conditions is obtained, providing a more targeted joint probability distribution of each order harmonic component for subsequent research on harmonic identification, etc.

[0079] Step two: For each similar load current, decompose it into the fundamental component and harmonic components of each order.

[0080] In one embodiment, a filter or wavelet transform is used to process the load current of each similar class, decomposing it into a fundamental component and harmonic components of each order. In this embodiment, the load current includes a fundamental component and harmonic components of each order. To extract information about the harmonic components of each order in the load current, a filter or wavelet transform is used to decompose the load current of each similar class, obtaining the fundamental component and harmonic components of each order. The harmonic components of different orders are separated so that the probability density distribution of each harmonic component can be calculated, more accurately reflecting the distribution of each harmonic component.

[0081] Step 3: Determine the probability density distribution of each order harmonic component based on the preset kernel density function.

[0082] In one embodiment, the step of determining the probability density distribution of each order harmonic component based on a preset kernel density function includes:

[0083] With the goal of minimizing the kernel density estimation error, the density distribution of each harmonic component is calculated based on the preset kernel density function, resulting in the probability density distribution of each harmonic component. The kernel density estimation error is calculated based on the distribution error between the estimated and actual density distributions. In this embodiment, to accurately quantify the density distribution of each harmonic component, different window widths are traversed during density distribution calculation. The density distributions calculated for each harmonic component are based on the preset kernel density function. Different window widths may yield different density distributions. The smoothness of the preset kernel density function is strongly correlated with the window width; therefore, selecting a suitable window width is crucial for accurately calculating the density distribution of each harmonic component. To obtain accurate probability density distributions of each harmonic component, different initial window widths are traversed, and the kernel density estimation error corresponding to each initial window width is calculated. When the kernel density estimation error is minimized, the corresponding initial window width is determined as the optimized window width. The preset kernel density function corresponding to the optimized window width is then determined, and the probability density distribution of each harmonic component is calculated based on the preset kernel density function corresponding to the optimized window width.

[0084] Furthermore, the preset kernel density function is shown in equation (1):

[0085]

[0086] In equation (1), Let represent the kernel density estimate of each order harmonic component at current x, b represent the window width, m represent the number of samples, W() represent the kernel function, Xi represent the samples from the total distribution X, and x represent the current value under each specific order current harmonic or each frequency domain interval current harmonic.

[0087] The kernel function W() satisfies the condition shown in equation (2):

[0088]

[0089] In equation (2), x represents the current value under each specific current harmonic or the current harmonic in each frequency range. That is to say, the function value of the kernel function W() is greater than or equal to 0, and the area of ​​the kernel function W() on the integration interval from negative infinity to positive infinity is 1.

[0090] Furthermore, the kernel function W() is subject to the conditions shown in equation (3):

[0091]

[0092] In equation (3), u represents the current value under each specific sub-current harmonic or each frequency range current harmonic.

[0093] In one embodiment, the kernel density estimation error is calculated based on a comparison coefficient between the expected value of the integral of the square of the predicted density distribution and the expected value of the integral of the predicted density distribution. Specifically, when the comparison coefficient is greater than a first threshold, the kernel density estimation error is calculated based on the expected value of the integral of the square of the predicted density distribution; when the comparison coefficient is less than or equal to the first threshold and greater than or equal to a second threshold, the kernel density estimation error is calculated based on the expected value of the integral of the square of the predicted density distribution and the expected value of the integral of the predicted density distribution; when the comparison coefficient is less than the second threshold, the kernel density estimation error is calculated based on the expected value of the integral of the predicted density distribution. In this embodiment, a preset kernel density function is used. The smoothness of the kernel density function is strongly correlated with the window width *b*. The appropriate selection of the window width *b* is a crucial technical aspect of kernel density estimation. Furthermore, the window width *b* significantly affects the performance of the preset kernel density function. The accuracy is greatly affected. To minimize the kernel density estimation error, the calculation formula for the kernel density estimation error is set for different situations when iterating through the window width b. First, the method to solve for the optimal window width b can be described as minimizing the initial error, which is shown in equation (4):

[0094]

[0095] In equation (4), Indicates the initial error, used to describe The mean square error of the integral between f(x) and f(x) is the actual population distribution density function, and f(x) is the predicted population distribution density function. From equation (4), we can see that the third term E{∫f 2 The value of dx (x) is independent of the window width b, so equation (4) can be converted into a formula that includes the kernel density estimation error of at least one of the first two terms. The optimization problem of the window width b can be converted into finding the minimum of the first two terms in equation (4) to determine the optimal window width b. That is, minimizing the following equation (5) to determine the corresponding window width b as the optimal solution. Equation (5) is shown below:

[0096]

[0097] In equation (5), the difference between the expected value of the integral of the square of the predicted density distribution and the expected value of the integral of the predicted density distribution is used as the error. Both terms are related to the window width b. Since the ratio of these two terms varies greatly under different working conditions, in order to flexibly calculate the kernel density estimation error for different working conditions, a corresponding kernel density estimation error calculation formula is set according to the comparison coefficient between the expected value of the integral of the square of the predicted density distribution and the expected value of the integral of the predicted density distribution. A first threshold and a second threshold are preset, and the first threshold is greater than the second threshold. When the comparison coefficient is greater than the first threshold, it indicates that the expected value of the integral of the square of the predicted density distribution is much greater than the expected value of the integral of the predicted density distribution. The error calculated in equation (5) mainly comes from the first term, namely the expected value of the integral of the square of the predicted density distribution. The kernel density estimation error should be calculated based on the expected value of the integral of the square of the predicted density distribution. When the comparison coefficient is less than or equal to the first threshold and greater than or equal to the second threshold, it indicates that the expected value of the integral of the square of the predicted density distribution is close to the expected value of the integral of the predicted density distribution. The error calculated in equation (5) originates from the first and second terms, namely the expected value of the integral of the square of the predicted density distribution and the expected value of the integral of the predicted density distribution. The kernel density estimation error should be calculated based on these two terms. When the comparison coefficient is less than the second threshold, it indicates that the expected value of the integral of the square of the predicted density distribution is much smaller than the expected value of the integral of the predicted density distribution. The error calculated in equation (5) mainly originates from the second term, namely the expected value of the integral of the predicted density distribution. The kernel density estimation error should be calculated based on the expected value of the integral of the predicted density distribution.

[0098] Furthermore, the formula for calculating the kernel density estimation error is shown in equation (6):

[0099]

[0100] In equation (6), The values ​​represent the kernel density estimation error, where t1 represents the first threshold and t2 represents the second threshold. This represents the expectation of the integral of the square of the predicted density distribution. This represents the expectation of the integral of the predicted density distribution.

[0101] Step four: Sort and classify the probability density distributions of each order harmonic component, and obtain the sorted marginal objective functions under each category as the first layer of the vine structure.

[0102] In one embodiment, the step of sorting and classifying the probability density distributions of the harmonic components of each order to obtain the sorted marginal objective functions under each category as the first layer of the vine structure includes:

[0103] Calculate the correlation between the probability density distributions of the harmonic components of each order to obtain the correlation coefficient between the probability density distributions of the harmonic components of each order.

[0104] Based on the correlation coefficient, the probability density distribution of each order harmonic component is classified to obtain the marginal objective function under each category.

[0105] Based on the correlation coefficients, the marginal objective functions within the same category are sorted, and the sorted marginal objective functions within each category are used as the first layer of the vine structure. In this embodiment, to effectively extract the relevant information between the probability density distributions of each order harmonic component, the correlation coefficients between the probability density distributions of each order harmonic component are obtained. The larger the correlation coefficient, the greater the correlation between the probability density distributions of the corresponding order harmonic components, and they are classified into the same category. For example, the probability density distributions of each order harmonic component with a correlation coefficient greater than or equal to a preset correlation threshold can be classified into the same category. Furthermore, the probability density distributions of each order harmonic component within the same category are sorted according to their correlation from largest to smallest, thereby obtaining the sorted marginal objective functions within each category as the first layer of the vine structure. Therefore, when subsequently associating adjacent marginal objective functions or adjacent joint probability distributions in the vine structure, two marginal objective functions or two joint probability distributions with a high degree of correlation can be associated, more effectively extracting association information and improving information utilization.

[0106] Step 5: Use the copula function to correlate the two adjacent marginal objective functions in each category of the first layer to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, use the copula function to correlate the two adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2, and n is a positive integer greater than m.

[0107] Step 6: Use the copula function to correlate the joint probability distribution corresponding to each category to obtain the joint probability distribution of each order harmonic component.

[0108] Example 3

[0109] Based on the above embodiments, a method for constructing a joint probability distribution model of harmonic components of various orders is provided, including:

[0110] Step (1), extraction of specific subcurrent harmonics / division of current harmonic frequency domain intervals and extraction of frequency domain harmonics.

[0111] Step (2), probability distribution of specific subcurrent harmonics or frequency domain current harmonics.

[0112] Step (3): Correlation analysis of specific subcurrent harmonics or frequency domain current harmonics.

[0113] Step (4), the joint probability distribution of specific subcurrent harmonics or frequency domain interval current harmonics.

[0114] In step (1), specific subcurrent harmonics are extracted / current harmonic frequency domain intervals are divided and frequency domain harmonics are extracted.

[0115] Because the current harmonics in rail transit traction systems exhibit characteristics such as fluctuation, periodicity, superposition, and wide frequency domain bandwidth distribution, frequency domain analysis methods are required to achieve in-depth analysis of current harmonics. Based on measured data from trains operating on actual lines, load current data under different load fundamental currents, such as traction, braking, and coasting, are classified and statistically analyzed. Load currents under the same operating conditions with fundamental currents differing within a certain range are grouped into similar categories.

[0116] The load current sample data from similar classes are sequentially passed through low-pass, band-pass, and other filters to obtain specific harmonic current signals, such as... Figure 3 As shown. Figure 3 The curves from top to bottom represent: the original load current signal, the fundamental current obtained after the load current passes through a low-pass filter, and the total harmonic current remaining after the fundamental current is filtered out.

[0117] Alternatively, wavelet transform can be applied sequentially to the load current sample data in the similar classes, and frequency domain partitioning analysis can be performed on the load current frequency domain, such as... Figure 4 As shown. Figure 4 The curves from top to bottom represent: the original load current signal graph, the fundamental current graph obtained by wavelet transform of the load current at the corresponding fundamental frequency, the current harmonic graph corresponding to the fifth frequency domain interval of the load current outside the fundamental frequency domain, the current harmonic graph corresponding to the fourth frequency domain interval of the load current outside the fundamental frequency domain, the current harmonic graph corresponding to the third frequency domain interval of the load current outside the fundamental frequency domain, the current harmonic graph corresponding to the second frequency domain interval of the load current outside the fundamental frequency domain, and the current harmonic graph corresponding to the first frequency domain interval of the load current outside the fundamental frequency domain.

[0118] Therefore, see Figure 3 and Figure 4 It can be seen that the inverter-side current harmonics of the traction system do exhibit certain characteristics such as volatility, periodicity, superposition, and wide frequency distribution, which are significantly different from other intermittent random variables (such as wind power and photovoltaic power). This invention patent will simultaneously consider statistical analysis and probabilistic modeling of the original current signal to obtain each harmonic current through a filter and the corresponding current harmonic signals in each frequency range through wavelet transform. Figure 5 A schematic diagram showing the one-to-one correspondence between the frequency domain and time domain of each current harmonic obtained by passing the load current through a filter is given.

[0119] At the same time, such as Figure 6 To divide the load current frequency domain into n+1 frequency intervals using wavelet transform, the first interval is the fundamental frequency component of the load current, denoted by 'an'. The second interval is assumed to include the third and fifth harmonics of the load current, denoted by 'd1'. Subsequent harmonic frequency interval divisions follow the same assumption. The advantages of using wavelet transform for frequency interval division are: firstly, it eliminates the need to design low-pass and band-pass filters for each order of current harmonics, significantly reducing design workload and difficulty; secondly, it considers the interaction between different current harmonics within the same frequency interval, avoiding the impact of abrupt changes in specific harmonics on the overall harmonic periodicity and stability of the frequency interval.

[0120] Therefore, this invention proposes for the first time to use wavelet transform to perform frequency domain partitioning analysis on the load current of the traction system, which can obtain the fundamental current of the corresponding frequency domain interval and the current harmonic components in each frequency domain interval.

[0121] In step (2), the probability distribution of a specific subcurrent harmonic or a frequency domain current harmonic is obtained.

[0122] By using filters to obtain specific sub-current harmonics corresponding to similar load currents, or by obtaining current harmonics in various frequency domain intervals based on wavelet transform, empirical density functions can be obtained, denoted as follows, for specific sub-current harmonics or current harmonics in various frequency domain intervals corresponding to the fundamental current within a certain range, based on statistical principles. Based on empirical density function An improved kernel density estimation method is proposed to estimate the kernel density of each specific current harmonic or each frequency domain current harmonic f(x) at current x as follows:

[0123]

[0124] In equation (1) above, b represents the window width, m represents the number of samples, and W() represents the kernel function, where the kernel function W() needs to satisfy the requirements of equation (2). Xi represents the sample from the total distribution X, and x represents the current value under each specific current harmonic or current harmonic in each frequency domain interval.

[0125]

[0126] Referring to relevant literature on kernel density estimation methods, most studies use the Gaussian kernel function to estimate the marginal distribution of current harmonics at specific sub-current harmonics or in different frequency ranges.

[0127]

[0128] Overall analysis shows that the kernel density estimation function... The smoothness of the kernel density estimation function is strongly correlated with the window width *b*. Therefore, appropriately selecting the window width *b* is a crucial technical aspect of kernel density estimation methods, and also affects the smoothness of the kernel density estimation function. The accuracy is very important. The method for solving the optimal window width b can be described as follows:

[0129]

[0130] The above formula describes The mean square error of the integral between f(x) and f(x) is the actual population distribution density function. From equation (4), it can be seen that the third term is independent of the window width b. Therefore, equation (4) can be transformed into an optimization problem where only the first two terms contain the window width b, i.e.

[0131]

[0132] Therefore, for solving the window width b in the optimization model of equation (5), this application proposes an adaptive optimization method for solving the kernel function window width b. That is,

[0133]

[0134] Where t represents the adaptive selection factor. The above equation shows that the adaptive optimization method for solving the kernel function window width b will... By rationally selecting the window width b closest to the optimization for the evaluation index, the computational load of the optimization model can be greatly reduced, making it a feasible optimization solution method for engineering applications with high real-time computation requirements. By solving with this adaptive window width b, more accurate probability densities of specific current harmonics or current harmonics in various frequency domain intervals can be obtained.

[0135] (3) Correlation analysis of specific sub-current harmonics or current harmonics in the frequency domain.

[0136] according to Figure 5 The load current is filtered to obtain the time-domain signals of each specific current harmonic, such as... Figure 7 As shown.

[0137] In this embodiment, Figure 7 I1 represents the time-domain signal with a frequency of 3. This is because, after passing through a low-pass filter, this method can obtain a current harmonic time-domain signal with only the third fundamental frequency.

[0138] according to Figure 6 The load current is transformed using wavelet transform to obtain the time-domain signal of the corresponding current harmonics in the frequency domain interval, as shown below. Figure 8 As shown.

[0139] In this embodiment, Figure 8 This diagram can be interpreted as follows: After wavelet transform, the time-domain signals of current harmonics corresponding to different frequency ranges are obtained. Specifically, `an` represents the fundamental frequency. `d1` represents the harmonics within this frequency range that include 3fc and 5fc; therefore, `I1(t)` represents the time-domain signal of the current containing these harmonics. `I0(t)` represents the time-domain signal of the current corresponding to the fundamental frequency `fc`. `f1`, `f2`, ... only represent the endpoints of the frequency ranges; that is, the frequency range (f1, f2) only contains the fundamental current `fc`, and the corresponding time-domain signal is `I0(t)`. Similarly, (f2, f3) represents the frequency range containing only the frequencies of 3fc and 5fc, and the corresponding time-domain signal is `I1(t)`. The correlation coefficient refers to the correlation between time-domain signals in different frequency ranges.

[0140] In this embodiment, Figure 7 and Figure 8 The expressions are different. Figure 7 Using a low-pass filter, only the time-domain signals corresponding to specific sub-current harmonics can be obtained, such as 3fc corresponding to I1(t), 5fc corresponding to I2(t), ... Figure 8 This indicates that the original current signal is divided into frequency domain intervals using the wavelet transform method. For example, the frequency domain interval d1 may include different low-frequency current harmonics such as 3fc and 5fc. That is, I1(t) represents the time domain signal of the current harmonics corresponding to the two current harmonics of 3fc and 5fc. The frequency domain interval d2 may include the time domain signal of the current harmonics corresponding to the three current harmonics of 7fc, 9fc and 11fc.

[0141] Therefore, the harmonic signal corresponding to the frequency domain can be understood as the frequency range containing harmonic currents of different harmonic frequencies. Figure 5 It represents harmonic currents of only a single frequency. In the frequency domain, it is not a simple addition, but rather it is understood as dividing the entire frequency domain into different segments, and each segment or frequency domain interval contains current harmonics of different frequencies.

[0142] Meanwhile, a correlation evaluation criterion is defined. To account for the differences in current harmonic amplitudes within different frequency ranges, the rank correlation coefficient is used to perform correlation analysis on the time-domain signals of current harmonics in different frequency ranges. Taking I0(t) and I1(t) as examples, the following explanation is provided:

[0143] If the rank correlation coefficient ρ between I0(t) and I1(t) is I0(t)、I1(t) If I0(t) ≥ a, then I1(t) is considered to be strongly correlated.

[0144] If the rank correlation coefficient b between I0(t) and I1(t) is ≤ ρ I0(t)、I1(t)If <a, it is considered that there is a weak correlation between I0(t) and I1(t);

[0145] If the rank correlation coefficient ρ between I0(t) and I1(t) I0(t)、I1(t) <b, it is considered that there is no correlation between I0(t) and I1(t).

[0146] Therefore, by using the rank correlation coefficient to characterize the cross-correlation degree of current harmonics in different frequency domain intervals, time-domain signals with similar correlation degree ranges can be divided into the same class, and then the joint probability distribution of such time-domain signals can be described. The different classes of time-domain signals obtained here can be used to obtain the corresponding probability distribution function by using the improved kernel density estimation method. The magnitude of the correlation coefficient here is for the division of a certain class in the first layer of the subsequent improved pair-copula function. For example, if the rank correlation coefficient between I0, I1, and I3 is large, the time-domain signal of this current harmonic can be used as the first class in the first layer of the pair-copula; if the rank correlation coefficient between I2 and I4 is large, the time-domain signal of this current harmonic can be used as the second class in the first layer of the Pair-copula,....

[0147] (4) Joint probability distribution of specific order current harmonics or current harmonics in frequency domain intervals

[0148] Selecting the same type of copula function can characterize the joint distribution function between specific order current harmonics or current harmonics in different frequency domain intervals. However, the copula function does not consider the cross-correlation between current harmonics in any two frequency domain intervals. Therefore, the present invention specifically proposes an improved vine structure pair-copula method considering spatio-temporal cross-correlation to construct the joint distribution function of multivariate copulas. On the one hand, this improved method uses the vine structure pair-copula method to merge different specific order current harmonics or current harmonics in different frequency domain intervals pairwise and can introduce different types of bivariate copula functions. On the other hand, it performs hierarchical merging of different classes according to the rank correlation coefficients between variables, which can save a large amount of computational work for constructing the high-dimensional correlation of variables in the conventional pair-copula method. The improved pair-copula method is as follows Figure 9 as shown.

[0149] above Figure 9 Among them, F i is the marginal distribution function of the i-th specific order current harmonic or current harmonic in the frequency domain interval of the current, which can be obtained by the kernel density estimation method. C ij is the bivariate copula function with the marginal distribution function F i of the i-th specific order current harmonic or current harmonic in the frequency domain interval of the current and the marginal distribution function F jThe binary joint distribution function is called the copula function. The cross-correlation between different specific current harmonics or current harmonics in different frequency ranges can be determined by selecting different types of copula functions based on the actual degree of correlation. Therefore, compared to the previous method, the improved pair-copula method can still consider using different copula functions for pairs of variables with high correlation, while avoiding merging variables with correlation below a certain threshold with the root node. This is equivalent to having multiple parallel, nested pair-copula structures in the entire pair-copula vine structure diagram. This satisfies the need to consider the cross-correlation between variables while saving a significant amount of computational work.

[0150] use Figure 9 The improved pair-copula method shown can establish the joint probability distribution function of specific sub-current harmonics or current harmonics in various frequency domains. The main calculation steps are as follows:

[0151] Step 1: Calculate the rank correlation coefficient of specific current harmonics or current harmonics in various frequency ranges. Based on the correlation degree evaluation criteria, classify the random variables with high correlation into one category. Then, calculate and classify the correlation coefficients of the remaining random variables, ultimately obtaining the following result: Figure 9 The variable partitioning shown is similar to... Figure 7 U1-U4 are classified into one class, and Uj-1-Un are classified into another class.

[0152] Step 2: Based on the correlation coefficients of each type of random variable in Step 1, determine the order of the specific sub-current harmonics or frequency range current harmonics in the first layer of the canonical vine structure, and denote each specific sub-current harmonic or frequency range current harmonic in the sorted order as {(I0,I1,...,I...}}. k ),(I k+1 ,I k+2 ,...,I j ),...,(I m ,I m+ 2,...,I n ))};

[0153] Step 3: Based on the number of first-type random variables, select a binary copula function to combine the reordered first-type specific current harmonics or frequency domain harmonics with the specific current harmonics or frequency domain harmonics corresponding to the first root node of the first layer, thus establishing the first-type first-layer pair-copula sequence {C}. 01 (F 1 (I0),F2(I 1 )),...,C 0k (F 1(I0),F k+1 (I k ))};

[0154] Step 4: Following the method in Step 3, establish the first-level pair-copula sequence for the random variables of the second and subsequent categories respectively;

[0155] Step 5: Use the distribution functions of each category obtained in steps 3 and 4 as new variable sequences, and repeat steps 3 and 4 until all conditional copula sequences {C} of the second level of random variables for each category are obtained. 12 |0(F(I1|I0),F(I2|I0)),C 1,3 |0(F(I1|I0),F(I3|I0))},...{C j,n|j-1 (F(I j |I j-1 ),F(I n |I j-1 The calculation ends here. The correlation between variables C01, C02, and C03 of the first layer, first class, has already been considered. Subsequent calculations only require pairwise merging; further correlation checks are unnecessary. This step already assesses correlation at the first layer, first class level, ensuring a certain correlation exists between variables. Subsequent calculations only require pairwise merging. When each class is merged into a joint function, the final distribution functions of each class are combined to obtain the joint probability density function f(I0,...,I...). n ).

[0156] Step 6: Repeat step 5 until each category has only one binary conditional copula distribution remaining.

[0157] {C k-1,k|0,1,...,k-2 (F(I k-1 |I0,I1,...,I k-2 ),F(I k |I0,I1,...,I k-2 ))},...,...function{C n-1,n|j-1,j-2,...,n-2 (F(I n-1 |I j-1 ,I j-2 ,...,I n-2 ),F(I n |I j-1 ,I j-2 ,...,I n-2 ))}until.

[0158] Step 7: Using the last binary conditional copula function calculated for each category as the distribution function of the category variable, and based on the cross-correlation between the category variables, the joint probability density function f(I0,...,I) of all specific sub-current harmonics or current harmonics in each frequency domain interval is finally obtained using a multivariate copula function. n ).

[0159] Taking the calculation of the joint probability distribution function of current harmonics in the frequency domain as an example, the overall implementation framework is illustrated in Figures 10A and 10B. To calculate the joint probability distribution function of a specific current harmonic, it is only necessary to replace the current signal in each frequency domain interval with the corresponding current harmonic signal of that specific frequency.

[0160] Because the inverter current harmonics in the traction system exhibit periodic characteristics, and the implementation frameworks shown in Figures 10A and 10B only consider the cross-correlation between random variables, the corresponding periodic characteristics cannot be reflected. Therefore, sampling only within the joint probability distribution function of the current harmonics in the frequency domain corresponding to Figures 10A and 10B can only guarantee the cross-correlation of the current harmonics in the frequency domain at a single moment, but cannot reflect its periodicity. This patent, based on this, will also consider the temporal correlation of the current harmonics, specifically as follows:

[0161] Step 1: Analyze the sample data of current harmonics in each frequency domain by period. Taking a sinusoidal current waveform as an example, record the current harmonic values ​​in ascending order within the ranges (0, π / 2), (π / 2, 3π / 2), and (3π / 2, 2π). The harmonic value refers to the magnitude of the harmonic current. Theoretically, current harmonics should also be sinusoidal signals, only the instantaneous value of the harmonic amplitude changes with time.

[0162] Step 2: Random sampling is performed from the joint probability distribution function of the frequency domain current harmonics corresponding to Figures 10A and 10B. Before performing the inverse kernel density estimation transform to obtain the harmonic values ​​of each frequency domain current interval, the random sample numbers are sequentially assigned to the sample permutation index values. The random sample numbers are then rearranged according to the permutation index values ​​to give the random sample numbers a periodic characteristic. The previous wavelet transform only obtained the harmonic currents of each frequency domain interval, while this step refers to considering the correlation between the current harmonics after considering the joint probability distribution. Sampling is then performed from the joint probability distribution to mainly consider their cross-correlation. The sample size is set according to actual needs when sampling from the constructed joint probability distribution. Figures 10A and 10B are based on the modeling needs of harmonic current samples and can consider the construction of a joint probability distribution with cross-correlation. The sampling is performed in subsequent practical applications to simulate the real harmonic current situation, considering sampling from the constructed joint probability distribution to ensure that its spatiotemporal correlation is considered simultaneously. If the overall probability distribution of harmonic current is obtained through linearization, then sampling can be performed according to different categories, and sampling should be conducted separately for different operating conditions. The sample sequence number is arranged in ascending order of period according to the time domain signal corresponding to the current harmonics. Therefore, after random sampling, the randomly sampled samples need to be reordered according to the sequence number before performing an inverse kernel density transformation to obtain sample values ​​that take into account spatial and temporal correlations. If the values ​​of current harmonics in a certain period are [2,4,7,9,10,14,16,20,18,15], based on the order of these sample values, assuming the random numbers for random sampling are [0.2,0.4,0.1,0.05,0.6,0.75,0.23,0.55,0.67,0.70], then matching the original order of sample values, the random numbers for random sampling are reordered to [0.05,0.1,0.2,0.23,0.4,0.55,0.6,0.75,0.70,0.67].

[0163] Step 3: The rearranged random sample numbers from Step 2 are inversely transformed using kernel density estimation to obtain sampled current harmonics for each frequency domain interval. Random sampling is performed within the constructed joint probability distribution, using random numbers. The samples need to be rearranged according to their sorting values ​​before undergoing inverse kernel density transformation to ensure that the final sampled data takes into account both spatial and temporal correlations. The final sampled data needs to be transformed to the time domain signal corresponding to a specific frequency domain, or the time domain signal corresponding to a specific current harmonic.

[0164] Therefore, by following the above approach, the cross-correlation between current harmonics in each frequency domain is fully considered, and the pair-copula model is improved. This saves a lot of computational workload while ensuring the required computational accuracy. Furthermore, the periodicity characteristics of current harmonics are considered simultaneously through the sample sequence sorting method. As a result, the joint probability distribution of current harmonics in each frequency domain can be well applied to practical engineering projects.

[0165] This embodiment relates to the field of rail transit traction systems, focusing on modeling the joint probability distribution of harmonic currents on the inverter side of the traction converter. This provides reference information for subsequent harmonic identification, prediction, and verification of harmonic suppression control algorithms on the inverter side. This method is also applicable to modeling the joint probability distribution of current harmonics under different operating conditions, such as the power supply grid side of rail transit traction systems, power transmission networks, distribution networks, and microgrids.

[0166] In this embodiment, to ensure the accuracy of the traction load current harmonic probability density model and the effectiveness of the harmonic suppression control algorithm used in engineering projects, this invention proposes a joint probability distribution model of inverter current harmonics in a train traction system. First, based on different operating conditions of the train and corresponding measured data, statistics and analysis are performed according to the range of the fundamental current. Low-pass and band-pass filters that meet the computational accuracy requirements can be designed to extract the fundamental current and each harmonic current from the current data; alternatively, wavelet transform can be used to perform frequency domain partitioning analysis on the measured current data, dividing the frequency domain corresponding to the current into specific frequency domain intervals. Each frequency domain interval contains not only a single frequency harmonic of a specific current, but may also contain multiple harmonics of specific frequencies. Then, an improved nonparametric kernel density estimation method is used to obtain the probability distribution of specific current harmonics or the probability distribution of current harmonics in each frequency domain interval. Finally, the cross-correlation between current harmonics in each frequency domain interval is statistically analyzed, and an improved Pair-Copula method is proposed to characterize the cross-correlation of each specific sub-current harmonic or the current harmonics in each frequency domain interval. Based on the current linearization segmentation method, the joint probability distribution function of the specific sub-current harmonics or the current harmonics in each frequency domain interval of the line is obtained.

[0167] This embodiment provides a joint distribution function of current harmonics that can describe the traction load of a train, offering reference information for more accurate assessment of the current harmonic characteristics of the traction load, the impact of current harmonics on grid-side power quality, and the effectiveness verification of current harmonic suppression control algorithms. The overall approach involves designing low-pass and band-pass filters that meet computational accuracy requirements based on measured train operation data, extracting the fundamental current and various current harmonics from the current data; or using wavelet transform for frequency domain partitioning analysis, i.e., statistically analyzing the current harmonics in each frequency domain interval. An improved nonparametric kernel density estimation method is used to obtain the probability distribution of specific current harmonics or current harmonics in each frequency domain interval. An improved pair-copula method is used to consider the cross-correlation of specific current harmonics or current harmonics in different frequency domain intervals. Finally, the joint probability distribution function of specific current harmonics or current harmonics in each frequency domain interval is obtained based on the current linearization piecewise method. This scheme is mainly divided into several parts, including specific harmonic extraction or frequency domain division and frequency domain harmonic extraction, probability distribution of specific harmonic current or frequency domain interval current harmonics, correlation analysis of specific harmonic current or frequency domain interval current harmonics, and joint probability distribution of specific harmonic current or frequency domain interval current harmonics.

[0168] To reduce the impact of train traction load current harmonics on train operation and the resulting electromagnetic compatibility issues, and to improve the accuracy of train traction load modeling and verify the effectiveness of traction load current harmonic suppression control algorithms, this invention proposes a joint probability distribution model for current harmonics. The main innovations of this application are as follows:

[0169] (1) Based on the test data of actual train operation conditions, an innovative scheme is proposed to divide the train traction load current according to different operating conditions and different fundamental current ranges. Load currents with the same operating conditions and fundamental currents within a certain range are divided into similar categories, so as to realize the current harmonic statistics and analysis under the same operating conditions and the same fundamental current range.

[0170] (2) Based on the load current sample data divided into similar classes, specific current harmonics are obtained through filters or wavelet transforms are used to obtain current harmonics in each frequency range. The cross-correlation coefficients of specific current harmonics or current harmonics in each frequency range are analyzed using rank correlation coefficients. The probability density distribution of current harmonics is obtained by applying kernel density estimation methods. Compared with existing parameter estimation methods, it is not necessary to assume a definite distribution function of current harmonics in each frequency range, and a more accurate probability density can be obtained.

[0171] (3) An innovative improved pair-copula method is proposed to describe the cross-correlation between specific sub-current harmonics of the load current or between current harmonics in different frequency domains. Simultaneously, the joint probability density function of each specific sub-current harmonic or current harmonic in different frequency domains can be obtained. Compared with existing literature reports using linear correlation and copula methods, this application can consider the high-dimensional correlation of current harmonics in different frequency domains. Furthermore, compared with existing pair-copula methods, it can save a significant amount of computational work, breaking the limitation of existing methods that only consider the pairwise cross-correlation between variables, and accurately reflecting the intrinsic relationship between current harmonics in different frequency domains.

[0172] (4) Simultaneously consider the periodicity of current harmonics. Based on the sequence number of the current harmonic sample data changing from small to large within the period, reorder the random number sequence obtained from the improved pair-copula method to ensure that the sampled data still maintains the periodicity.

[0173] (5) The joint probability distribution of current harmonics in each frequency domain interval of the line is obtained by the current linearization segmentation method, which provides reference information for subsequent current harmonic prediction, harmonic identification and verification of specific current harmonic suppression control algorithms.

[0174] Through application retrieval and literature review, current modeling of traction load current harmonics in rail transit primarily relies on statistical analysis. This involves performing FFT transformations on measured load current data and assuming deterministic distributions for each order of current harmonics based on empirical considerations. Probability densities for each order of current harmonics are then obtained through parameter fitting. This invention, building upon existing current harmonic probability distributions, innovatively proposes a classification scheme based on operating conditions and the range of the fundamental current. It also employs kernel density estimation to avoid the limitations of empirical distribution considerations and utilizes an improved pair-copula method to achieve high-dimensional correlation of current harmonics across different frequency domains. Ultimately, this yields a joint probability density function for current harmonics that more closely approximates actual conditions. Furthermore, the joint probability distribution function exhibits periodicity based on the variation sequence of the sample data. This invention provides valuable reference information for predicting harmonic currents on train lines. Simultaneously, it provides a more accurate load model for verifying specific current harmonic suppression control algorithms, enabling more accurate assessment of the effectiveness and applicability of these algorithms. This application does not require additional hardware facilities. It only reserves sample storage space in the underlying control software of the equipment for storing load current under actual operating conditions and calculating the joint probability density distribution model. The method is simple, safe, and reliable.

[0175] In this embodiment, a current harmonic joint probability distribution model is proposed. The object of application is the traction load of a train. However, it is also applicable to rail transit subways, locomotives, power system transmission networks, distribution networks, etc., and is also within the scope of protection of this invention.

[0176] In this embodiment, a joint probability distribution model considering specific sub-current harmonics or current harmonics in various frequency domain intervals of the traction load current is proposed. Currently, relevant patents or literature reports all take the probability density of specific sub-current harmonics as an example for research, while this patent can also quickly obtain the probability density by taking the partial derivative of the joint probability distribution, which also falls within the scope of protection of this invention.

[0177] Example 4

[0178] Based on the above embodiments, an apparatus for constructing the joint probability distribution of harmonic components of each order is provided, comprising:

[0179] The current classification module is used to classify and statistically analyze load current data under different operating conditions, so as to divide load currents under the same operating condition and whose fundamental frequencies differ within a specified range into similar categories.

[0180] The current decomposition module is used to decompose the load current of each similar class into the fundamental component and harmonic components of each order.

[0181] The harmonic density module is used to determine the probability density distribution of each order of harmonic components based on a preset kernel density function.

[0182] The sorting and classification module is used to sort and classify the probability density distribution of each order harmonic component, and obtain the sorted marginal objective function under each category as the first layer of the vine structure.

[0183] The density correlation module is used to correlate adjacent marginal objective functions in each category of the first layer using the copula function to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, the copula function is used to correlate adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2 and n is a positive integer greater than m.

[0184] The category association module is used to associate a joint probability distribution corresponding to each category using the copula function to obtain the joint probability distribution of each order harmonic component.

[0185] Example 5

[0186] Based on the above embodiments, a method for sampling harmonic current is provided, comprising:

[0187] Step 510: Obtain the joint probability distribution of each order harmonic component using the construction method of the joint probability distribution model of each order harmonic component obtained by the construction method described in any of the above embodiments.

[0188] Step 520: Obtain the current object to be studied.

[0189] In this embodiment, the current object to be studied is the current object that needs to be studied in the application scenario when the joint probability distribution of harmonic components of each order is applied. The current object to be studied contains harmonic components of each order, and the harmonic components contained therein are the reference information required in subsequent research problems such as harmonic identification.

[0190] Step 530: Analyze the correlation between time and current for the current object under study, and determine the periodic variation law of the current object under study.

[0191] In this embodiment, the frequencies and periods of the harmonic components corresponding to different current objects under study are different. Within the same time period, the current magnitude of different current objects under study shows the periodic variation pattern over time.

[0192] Step 540: Randomly sample the joint probability distribution to obtain random current arrays corresponding to each preset current component.

[0193] In this embodiment, in order to obtain harmonic components more accurately, the joint probability distribution f(I0,I1,...,I) is used. n Random sampling is performed to obtain random current arrays {i1}, {i2}, {i3}, ..., {i...} corresponding to each preset current component. n These random current arrays follow the interrelationships between harmonic components of each order in the joint probability distribution, and can more accurately measure each harmonic component.

[0194] Step 550: Based on the periodic change relationship, sort the sample current values ​​in each of the random current arrays to obtain harmonic random samples.

[0195] In this embodiment, in order to make the variation pattern of the random current array sampled in step 540 conform to the periodic variation pattern of the current object under study, the sample current values ​​in the random current array are sorted so that the sorted harmonic random samples exhibit the same periodic variation pattern as the current object under study.

[0196] Step 560: Based on the inverse transform of the preset kernel density function, sample extraction is performed on the edge function constructed by the random harmonic samples to obtain the harmonic components contained in the current object under study.

[0197] In this embodiment, the edge function of the harmonic can be calculated based on random samples of the wave. By using the inverse transformation of the preset kernel density function, samples can be extracted from the edge function to obtain the harmonic components contained in the current object under study, so that the harmonic components of the current object under study can be identified, predicted, and studied.

[0198] In this embodiment, the application process of the joint probability distribution of each order harmonic component is as follows:

[0199] (1) The joint distribution function f(I0,I1,...,I) obtained based on the constructed improved pair-copula model n Generate n sets of uniformly distributed random numbers: {i1}, {i2}, {i3}, ..., {i... n Based on the conditional probability distribution and the method of partial derivatives, n sets of specific sub-current harmonic random samples or n sets of frequency domain interval current harmonic random samples {I0}, {I1}, ..., {I...} are calculated. n}

[0200] (2) In order to ensure the periodicity of each current harmonic or each frequency range current harmonic, the n random samples of specific current harmonics or the n random samples of frequency range current harmonics are reordered according to the sequence number of the current harmonic sampling sample increasing (or decreasing) in the monotonic interval, so that the n random samples of specific current harmonics or the n random samples of frequency range current harmonics also have the characteristics of periodic change.

[0201] (3) Marginal distribution function: The n sets of random samples of specific sub-current harmonics or n sets of random samples of current harmonics in the frequency domain are subjected to the inverse transformation of the corresponding kernel density estimation method to obtain the relevant and periodic sampling samples of each specific sub-current harmonic or each frequency domain current harmonic: {I'1}, {I'2}, ..., {I' n}

[0202] (4) Based on the obtained current harmonic sampling samples, subsequent harmonic prediction, development of specific current harmonic suppression control algorithms, and research on technical issues related to current harmonics can be carried out.

[0203] Example 6

[0204] Based on the above embodiments, a harmonic current sampling device is provided, comprising:

[0205] The joint probability module is used to obtain the joint probability distribution of each order harmonic component using the construction method described in any of the above embodiments;

[0206] The research object acquisition module is used to acquire the current object to be studied.

[0207] The periodicity module is used to analyze the correlation between time and current of the current object under study and to determine the periodic variation pattern of the current object under study.

[0208] The distribution sampling module is used to randomly sample the joint probability distribution to obtain a random current array corresponding to each preset current component.

[0209] The current sorting module is used to sort the sample current values ​​in each of the random current arrays based on the periodic change relationship to obtain harmonic random samples;

[0210] The edge extraction module is used to extract the edge function constructed from the random harmonic samples based on the inverse transformation of the preset kernel density function, so as to obtain the harmonic components contained in the current object under study.

[0211] Example 7

[0212] Based on the above embodiments, this embodiment provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method described in the above embodiments.

[0213] In some embodiments of this example, a computer-readable storage medium is provided, on which a computer program is stored, which, when executed by a processor, implements the steps of the method described in the above embodiments.

[0214] In some embodiments of this example, a computer program product is provided, including a computer program / instructions, which, when executed by a processor, implements the steps of the method described in the above embodiments.

[0215] The processor may include, but is not limited to, one or more processors or microprocessors. Each processor may be implemented as an Application Specific Integrated Circuit (ASIC), Digital Signal Processor (DSP), Digital Signal Processing Device (DSPD), Programmable Logic Device (PLD), Field Programmable Gate Array (FPGA), controller, microcontroller, microprocessor, or other electronic component, for executing the methods described in the above embodiments.

[0216] Computer-readable storage media can be implemented by any type of volatile or non-volatile storage device or a combination thereof. Computer-readable storage media may include, but are not limited to, random access memory (RAM), read-only memory (ROM), flash memory, EPROM memory, EEPROM memory, registers, and computer storage media (e.g., hard disks, floppy disks, solid-state drives, removable disks, CD-ROMs, DVD-ROMs, Blu-ray discs, etc.).

[0217] Computer-readable storage media may also store at least one computer-executable program / instruction, such as computer-readable instructions. Computer-readable storage media include, but are not limited to, volatile memory and / or non-volatile memory. Volatile memory may include, for example, random access memory (RAM) and / or cache memory. Computer-readable storage media may include, for example, read-only memory (ROM), hard disk, flash memory, etc. For example, a non-transitory computer-readable storage medium may be connected to a computing device such as a computer, and then, when the computing device executes the computer-readable instructions stored on the computer-readable storage medium, the various methods described above can be performed.

[0218] In addition, the computer device may include (but is not limited to) a data bus, an input / output (I / O) bus, a display, and input / output devices (e.g., keyboard, mouse, speakers, etc.).

[0219] The processor can communicate with external devices via the I / O bus through wired or wireless networks.

[0220] In one embodiment, the at least one computer-executable instruction may also be compiled into or comprise a software product / computer program product, wherein one or more computer-executable instructions are executed by a processor to perform the steps of the various functions and / or methods in the embodiments described herein.

[0221] In the embodiments provided in this disclosure, it should be understood that the disclosed apparatus and methods can also be implemented in other ways. The apparatus embodiments described above are merely illustrative; for example, the flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of apparatus, methods, and computer program products according to various embodiments of this disclosure. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of code containing one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions marked in the blocks may occur in a different order than those marked in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in a block diagram and / or flowchart, and combinations of blocks in block diagrams and / or flowcharts, can be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.

[0222] It should be noted that, in this disclosure, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element limited by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0223] While the embodiments disclosed herein are as described above, the foregoing content is merely for the purpose of facilitating understanding of this disclosure and is not intended to limit this disclosure. Any person skilled in the art to which this disclosure pertains may make any modifications and changes in form and detail of the implementation without departing from the spirit and scope of this disclosure; however, the scope of patent protection of this disclosure shall still be determined by the scope defined in the appended claims.

Claims

1. A method for constructing a joint probability distribution model of harmonic components of various orders, characterized in that, include: The load current data under different operating conditions are classified and statistically analyzed to classify the load currents under the same operating condition and whose fundamental effective values ​​differ within a specified range into similar categories. For each similar class of load current, it is decomposed into fundamental component and harmonic components of each order; Based on the preset kernel density function, the probability density distribution of each order harmonic component is determined; The probability density distributions of each order harmonic component are sorted and classified, and the sorted marginal objective functions under each category are used as the first layer of the vine structure. The copula function is used to correlate the two adjacent marginal objective functions in each category of the first layer to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, the copula function is used to correlate the two adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2 and n is a positive integer greater than m. By using the copula function to correlate the joint probability distribution corresponding to each category, the joint probability distribution of each order harmonic component is obtained.

2. The method according to claim 1, characterized in that, The operating conditions include traction, braking, or coasting conditions.

3. The method according to claim 1, characterized in that, The step of sorting and classifying the probability density distributions of the harmonic components of each order, and obtaining the sorted marginal objective functions under each category as the first layer of the vine structure, includes: Calculate the correlation between the probability density distributions of the harmonic components of each order to obtain the correlation coefficient between the probability density distributions of the harmonic components of each order. Based on the correlation coefficient, the probability density distribution of each order harmonic component is classified to obtain the marginal objective function under each category. Based on the correlation coefficient, the marginal objective functions under the same category are sorted, and the sorted marginal objective functions under each category are used as the first layer of the vine structure.

4. The method according to claim 1, characterized in that, The step of determining the probability density distribution of each order harmonic component based on a preset kernel density function includes: With the goal of minimizing the kernel density estimation error, the density distribution of each order harmonic component is calculated based on the preset kernel density function to obtain the probability density distribution of each order harmonic component; wherein, the kernel density estimation error is calculated based on the distribution error between the estimated density distribution and the actual density distribution.

5. The method according to claim 4, characterized in that, The kernel density estimation error is calculated based on a comparison coefficient between the integral of the square of the estimated density distribution and the integral of the estimated density distribution. Specifically, when the comparison coefficient is greater than a first threshold, the kernel density estimation error is calculated based on the integral of the square of the estimated density distribution; when the comparison coefficient is less than or equal to the first threshold and greater than or equal to a second threshold, the kernel density estimation error is calculated based on the integral of the square of the estimated density distribution and the integral of the estimated density distribution; and when the comparison coefficient is less than the second threshold, the kernel density estimation error is calculated based on the integral of the estimated density distribution.

6. A method for sampling harmonic current, characterized in that, include: The joint probability distribution of each order harmonic component obtained by the construction method described in any one of claims 1-5; Obtain the current object to be studied; The correlation between time and current of the current object under study is analyzed to determine the periodic variation law of the current object under study. Random sampling is performed on the joint probability distribution to obtain random current arrays corresponding to each preset current component; Based on the periodic variation relationship, the sample current values ​​in each of the random current arrays are sorted to obtain harmonic random samples; Based on the inverse transformation of the preset kernel density function, sample extraction is performed on the edge function constructed from the random harmonic samples to obtain the harmonic components contained in the current object under study.

7. A device for constructing the joint probability distribution of harmonic components of various orders, characterized in that, include: The current classification module is used to classify and statistically analyze load current data under different operating conditions, so as to divide load currents under the same operating condition and whose fundamental frequencies differ within a specified range into similar categories. The current decomposition module is used to decompose the load current of each similar class into the fundamental component and harmonic components of each order. The harmonic density module is used to determine the probability density distribution of each order of harmonic components based on a preset kernel density function. The sorting and classification module is used to sort and classify the probability density distribution of each order harmonic component, and obtain the sorted marginal objective function under each category as the first layer of the vine structure. The density correlation module is used to correlate adjacent marginal objective functions in each category of the first layer using the copula function to obtain the joint probability distributions under each category of the second layer. From the second layer to the nth layer of the vine structure, the copula function is used to correlate adjacent joint probability distributions in each category of the mth layer to obtain the joint probability distributions under each category of the (m+1)th layer of the vine structure, until a unique joint probability distribution is obtained for each category. The total number of layers in the vine structure is n, where m is a positive integer greater than or equal to 2 and n is a positive integer greater than m. The category association module is used to associate a joint probability distribution corresponding to each category using the copula function to obtain the joint probability distribution of each order harmonic component.

8. A computer device, comprising a memory, a processor, and a computer program stored in the memory, characterized in that, The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 5.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 5.

10. A computer program product comprising a computer program / instructions, characterized in that, When executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 5.

Images