CC BY
33Nonlinear Load Modeling for Analysis of Non-Sinusoidal Conditions in Electrical Networks Based on Measurements of Harmonic Parameters
L.I. Kovernikova1*, VC. Luong2
1 Melentiev Eneigy Systems Institute of the Siberian Branch of the Russian Academy of Sciences, Irkutsk, Russian Federation
2 Institute of Ship Design, Hanoi, Vietnam
Abstract — Non-sinusoidal conditions in electrical networks need to be calculated for their control and development of technical measures to maintain harmonic parameters according to the requirements of regulatory documents. These calculations are impossible without electrical network and nonlinear load models that adequately reflect them in computational programs. Nonlinear load models have been developed for a long time. Some studies present general modeling principles and models of various nonlinear devices. Others consider some nonlinear devices as equivalent nonlinear loads connected to low and medium voltage networks. A whole host of highpower nonlinear electrical equipment is connected to high voltage networks. Modeling nonlinear loads connected to these networks is a problem. Research of measured parameters of harmonic conditions in electrical networks has shown that they are random values. The probabilistic nature is determined by the network configuration, a range of network components, the number of nonlinear loads, wave and frequency properties of the network, harmonic source phase currents, voltage at terminals of nonlinear electrical equipment, changes in operating conditions and load power, and many other factors. Nonlinear loads can only be modeled based on the measurements of parameters of harmonic conditions due to their unpredictability. The paper presents an overview of existing methods for modeling nonlinear loads, a methodological approach to modeling nonlinear loads based on measured parameters, an algorithm for modeling harmonics of active and reactive currents, a computational program
* Corresponding author. E-mail: [email protected]
http://dx.doi.org/10.38028/esr.2021.03.0001
Received September 23, 2021. Revised October 07, 2021.
Accepted October 28, 2021. Available online November 28, 2021.
This is an open access article under a Creative Commons Attribution-NonCommercial 4.0 International License.
© 2021 ESI SB RAS and authors. All rights reserved.
algorithm designed to identify distribution functions of measured current harmonics, and modeling results for current harmonics of railway transformers supplying power to electric locomotives.
Index Terms. Harmonics, harmonic parameter analysis, measurements, nonlinear load modeling, power quality.
I. Introduction
Operating parameters of electric power systems should ensure economic feasibility, reliability of electrical networks, and quality of power supplied to consumers. When the latter is compromised, electrical power losses in electrical networks increase, the service life of electrical equipment of electrical networks and consumers shortens, the performance of process equipment of industrial enterprises decreases, which causes economic damage. Currently, the most relevant and most acute power quality issue is associated with the non-sinusoidal nature of voltage. This problem is typical not only of Russia but also of many countries around the world. The number of electrical equipment units with nonlinear current-voltage characteristics (nonlinear loads) that are the source of current harmonics and the cause of the voltage waveform distortion at network nodes is rapidly and continuously growing. In addition to the conventional nonlinear electrical equipment of aluminum smelters, railways, pulp and paper mills, and others, the amount of electronic equipment used in the digitalization of all areas of human life is drastically increasing. Nonlinear loads are present in networks of all voltages. Harmonics penetrate from high voltage networks into medium and low voltage networks and vice versa.
Currently, three legislative documents of the Russian Federation (the Civil Code [1], the Law "On Protection of Consumers' Rights" [2], the Law "On the Electric Power Industry" [3]) contain provisions on ensuring proper power quality and the responsibility for it. In real life, the parameters of electrical network operation often differ from those required for electrical equipment [4]. In particular, they do not meet the requirements outlined in [5]. Control of non-sinusoidal conditions in electrical networks, their analysis, and development of engineering measures to maintain the harmonic parameters in line
with the requirements of regulatory documents require calculations of harmonic parameters, which cannot be done without models that faithfully represent electrical networks and nonlinear loads in computational programs. At present, increasingly favorable conditions are being created for this, as the objective has been set to establish electric power systems equipped with numerous measuring instruments. Individual devices and systems are installed in electrical networks for long-term continuous monitoring of power quality indices and operating parameters, which allows obtaining information for developing nonlinear load models based on measurement results.
This paper deals with an analysis of existing methods designed to model the nonlinear loads, describes a methodological approach and an algorithm for modeling nonlinear loads based on measured harmonic parameters, a computational program algorithm for identification of cumulative distribution functions of a measured series of random harmonics of nonlinear load currents and calculation of their values with a given probability, as well as the results of modeling backed by the methodological approach, algorithms, and a computational program developed by us for currents of railway traction transformers that supply electric power to nonlinear loads of high-power, i.e., electric locomotives.
II. Analysis of nonlinear load modeling methods
AVAILABLE IN THE LITERATURE
Nonlinear load models have been developed for a long time. In some papers, the general principles of modeling are given and models of various individual nonlinear devices are presented on their basis, in others, individual nonlinear devices are considered as equivalent nonlinear loads connected to low and medium voltage networks. A large amount of nonlinear electrical equipment is connected to high voltage networks. Modeling of nonlinear loads connected to these networks remains an unsolved issue.
Models of nonlinear electrical equipment can be divided into three groups according to the principle of its operation: models of ferromagnetic devices, models of arc devices, and models of electronic devices. Ferromagnetic electrical equipment includes power transformers, synchronous machines, and induction motors. To model power transformers, the authors in [6, 7] propose a T-type transformer equivalent circuit. The magnetizing branch is represented as a source of no-load current harmonics, i.e., l0n = Icn + jl , where l0n is the complex of the no-load current for the nth harmonic, Icn and I are the active and reactive components of harmonics of the no-load current for the nth harmonic. Their definition is presented in detail in [7-9]. In [10] the transformer is represented as a source of harmonic currents of the two windings:
I - = Z 2 n /(kT Z ln + Z 2 n )0 n - for the primary winding,
12„ = Zln /(kT Z ln + Z 2n )i0n
- for the secondary winding, where Z 1n = R1 + jnXl , Z2n = R2 + jnX2 are complex impedances of the transformer windings; R1, R2, X1, X2 are the resistance and reactance values of the transformer windings; kT is the transformer transformation ratio.
The values of the winding resistance are calculated as R1 = 0.1026KnX1(J + n)
and
R2 = 0.1026KnX2(J + n), where the coefficients J and K are determined using the expressions in [11].
When developing models of synchronous machines, the authors of [12, 13] propose assuming that synchronous generators do not create voltage harmonics. Synchronous generators can be represented by equivalent impedance with respect to the busbars. The equivalent impedance is determined by the expression in [12] as Zn = nX, where X = 1/2(X" + X"), X"d is the direct-axis sub transient
reactance, and X" is the quadrature-axis sub transient reactance.
In [14-16] the authors present models of induction motors in the form of T-type equivalent circuits. The model of an induction motor with a single squirrel cage contains the following resistances: R1, X1 - resistance and reactance of the stator winding, R2, X2 - resistance and reactance of the rotor winding, Rc, X^ - resistance and reactance of
the magnetizing branch, Iln, Uln - harmonics of complex
current and complex voltage of the stator winding, i2n -harmonics of the complex current of the rotor winding, sn - the slip for the nth harmonic.
The rotor of a double-deck induction motor has two squirrel cages: one cage is the starting winding; the other cage is the operating winding. The model of this induction motor contains the following resistances: Rj, Xj - resistance and reactance of the stator winding; X2T - total reactance of the rotor windings dissipation; R1O, X1O - resistance and reactance of the operating stator winding; R1S, X1S -resistance and reactance of the stator starting winding; Rc, X^ - resistance and reactance of the magnetizing branch; I1n, UJ\„ - harmonics of the complex current and complex voltage of the stator winding; i2n - harmonics of complex current of the rotor windings.
Arc devices, whose principle of operation is the electric discharge, are represented in models by piecewise linearized voltage-current characteristics [17] and measured current harmonics [18].
There are several approaches to modeling electronic devices. The most widely used one is the representation of electronic devices by standard spectrum current harmonic sources is [13, 19]. Modeling rectifier converters involves the transfer functions linking the DC and AC sides of the converter [20].
Modeling individual nonlinear devices when analyzing harmonic conditions in electrical networks with a whole
host of nonlinear loads is challenging. In such a situation, one should employ models of equivalent nonlinear loads. To this end, two approaches are proposed in the published research for the adoption. The first approach combines several devices [21-23]. The second approach involves developing equivalent load models based on measurements [24-33]. The models are a set of linear multivariate regression equations that relate measured current harmonics to voltage harmonics and the active power of the modeled load [28, 29].
In [30], a probabilistic model of an equivalent nonlinear load is proposed. Depending on the numbers of current harmonics and the value of the total harmonic distortion of current, the authors suggest considering four types of nonlinear loads. Harmonics of complex current at the common connection point are represented as
l = K £ al ,
i=1
where KE - the share of power of nonlinear loads in
the total power of the equivalent load; ai - the weight coefficient of a nonlinear load of /th type in the total power of nonlinear loads; In), j - the value and phase of the harmonic current of a load of the /th type. All parameters are random variables. This approach requires measuring the current harmonics of each individual load to build the models, which cannot be done because of the large number and complexity of nonlinear loads.
In [31, 32] it is proposed to represent a distribution network with an unknown configuration by the equivalent load, which is modeled by the complex impedance of the distribution network of the «th harmonic (Z_„L) and the «th
harmonic of the complex current source (I L). Both of these values are defined as
Zl = (Unl - Un 2)/(/ 2 - 4),
1'nL 1nl + Un1 ^ <ZnL '
where U ,, U „, i „,, I „, are the measured harmonics of
n2' ' n2'
complex voltage and complex current at the node to which the modeled distribution network is connected in two different modes of the electricity supply network.
In [33], the authors propose decomposing the measured current In into its deterministic ID and random InC components. The expression for the description of InD is determined by fitting a polynomial function to the measured current I n using the method of least squares. The component InC is defined as the difference between the measured current In and the deterministic component InD, i.e., InC = In - InD. The expression for the random component is based on its probability density function. The authors of [33] advise to use the simplest functions, linear and quadratic, to describe InD and InC. The authors of [34] note the disadvantage of the latter method of modeling a nonlinear load. It lies in the fact that the random component of the current InC is rarely captured by one of the known
probability distribution laws. This leads to difficulty in describing the measured current In and requires the use of special ways of description.
The method of modeling nonlinear loads based on measurements of harmonic parameters, despite the noted disadvantages of the above modeling methods, has significant advantages [35, 36]:
a) it is simpler than the method based on reducing to equivalents and combining individual nonlinear devices because according to the results of measurements, the equivalent load model is obtained directly for the particular node of the network to which it is connected;
b) it can capture changes in the equivalent load over time if measurements are made over a long time;
c) it can be applied to any load and at any node of a network of different voltages;
d) it is best suited for estimating the parameters of the equivalent nonlinear load at busbars with mass loading.
Modeling the nonlinear loads connected to high voltage electrical networks as evidenced by measurements of harmonic parameters is the most effective way. High voltage networks are large-distance networks. They cover a significant area. Many nonlinear loads of highpower, as well as lower voltage networks, are connected to their nodes. Each of the loads is an enterprise with its lower voltage electrical network. The main process electrical equipment of the enterprise generates current harmonics into the network. A large amount of auxiliary equipment, including that with nonlinear voltage-current characteristics, receives electric power from the enterprise network. Current harmonics from these networks flow into the high voltage networks. In networks of all voltages, there are switchings of elements in the network and changes in its configuration. As a result of all processes, the harmonic parameters in the networks of all voltages are random, and are probabilistic in nature. Thus, the modeling cannot factor in all the individual nonlinear devices and loads connected to the high voltage electrical networks. Therefore, the most appropriate way of modeling for high voltage networks is based on measurements of harmonic parameters. Measurement is the only way to obtain accurate information. In what follows, we present a methodological approach, algorithms, and a computational program for modeling nonlinear loads based on the measurements of harmonic parameters.
III. Methodological approach, algorithms, and
COMPUTATIONAL PRoGRAM
The studies of the information on harmonic parameters that was obtained from the measurements [37-42] show that the harmonic conditions are random and their parameters are probabilistic in nature. It is determined by the configuration of the network, the composition of its elements, the number of nonlinear loads, the wave and frequency properties of the network, the phase of harmonic source currents, the values of voltage at the terminals of nonlinear electrical equipment, the changes in operating
Network
InN
T"
Measurement point
InL
©
—►
In
Un r
-T.
Nonlinear load
I
I
(!)
i i
i______'_ J-1______J
Fig. 1. Diagram of the network and nonlinear load at the nth harmonic.
conditions and power loads, and many other factors. Because of their unpredictability, nonlinear loads in reallife electrical networks can be modeled only relying on measurements of the harmonic parameters.
Harmonic parameters are calculated by solving a system of equations [16, 43]
U = Zn In, (1)
where Un - a column matrix of the nth harmonic nodal voltages to be determined, Zn - the square matrix of the self-and mutual impedances of the nth harmonic nodes of the network, In - a column matrix of the nth harmonic currents of nonlinear loads connected to the network nodes. The square matrix Zn is calculated as a result of the inversion
of the square matrix of nodal admittance Yn, i.e., Zn = Y-[44], which is formed based on the known parameters of the network elements. Elements of the matrix In are unknown, they have to be determined based on the results of measurements of the harmonic parameters. Each element of the matrix In is a complex number ini = Iani + jIrni, where i is the number of the network node, Iani and Irni are the active and reactive components of the nth harmonic current. The methodological approach and algorithms presented below are meant for their determination.
The harmonic parameters required for modeling the nonlinear load are measured at the node of its connection to the network, marked as a measurement point in the diagram in Fig. 1, where the network and the nonlinear load are represented as current sources. The diagram denotes the following: InN - the vector of the nth current harmonic of the network, InL - the vector of the nth current harmonic of the nonlinear load connected to the node, Un - the vector of the nth harmonic voltage at the node of the nonlinear load connection to the network. The current vector InN is the resultant vector of the nth current harmonic of all nonlinear
20 n
15 -^ 10 5 -i
0
000000000000000000000000000000 50505050505050505050505050505
t, min ^^^^^^^^^^
Fig. 2. Scatter chart of the 3rd current harmonic.
vector InL of the load modeled.
Thus, current In equal to the vector sum of currents InN and InL runs through the connection node, i.e.,
In = Ln + InL. (2)
The following parameters should be measured at the connection node of a nonlinear load: the nth harmonic factor of voltage (KU(n)) and current (KI(n)), the effective values of voltage and current harmonics: Un, In and the values of phase angles of voltage and current harmonics: 9Un, 9In. The arrays of all measured harmonic parameters are series of random variables.
To solve the system of equations (1) it is necessary to determine with a probability of 95% the values of active Ian and reactive Irn currents. The parameter values with a probability of 95% are used in [5] to estimate the deviation of voltage from a sine wave. Values of Ian and Irn are determined based on the measured effective magnitudes of current In and the phase angle 9UI(n) between the harmonics of voltages and currents, according to the expressions
Ian = In cos ^UKny (3)
Irn = In sin ^(n» (4)
where the phase angle is determined using the
measured angles 9Un and %n as
9w(n) = 9Un - (5)
Calculation of the current harmonic values with a given probability requires the knowledge of their cumulative distribution functions that can be obtained with the aid of probability density functions identified after appropriate processing of measurement results. Thus, the nonlinear load model is a set of current harmonic cumulative distribution functions that allow calculating with a given probability the harmonic values of active and reactive currents.
A. Methodological approach to modeling the nonlinear loads based on the measurements of harmonic parameters.
The methodological approach consists of seven steps listed below.
1. Check the measured series of random values of harmonic parameters for the presence of outliers. Outliers are elements that differ significantly in value from the other elements of the series. They make the series of measured parameters non-stationary, do not allow constructing a real histogram and correctly determining the distribution law. Their presence can be determined visually by a scatter chart.
b
prior to switching on the CB after switchin on the CB / > V V g ......... phase A
-I \ / M TTx" phase B
i *< ...... V ' •V i-/ ^___ _ phase C
--->I f
t, min
m <N
a)
<N
<N
<N
40 -, 35 -30 -25 -20 -15 -10 -5 0
prior to switching on the CB
after switching on the CB
phase A phase B
- - phase C
'•♦•■••»■.......................
^ — -» *
in 1> cn <n
<N
t, min b)
m <N
<N
<N
0\ <N
Fig. 3. Changes in a) KU(5) and b) I5 when the capacitor bank is on.
Fig. 2 shows a scatter chart of the 3rd currentharmonic. The diagram clearly shows the four outliers. During processing, outliers can be replaced by neighboring elements or by their average value, or the methods presented in [45, 46] can be employed to replace them.
2. Analyze the harmonic parameters measured to determine the time intervals, which are characterized by parameters that differ significantly in value.
Analysis of measured data revealed that the values of harmonic parameters differ significantly from each other at certain measurement time intervals. These differences result from the changes in the network layout, the parameters of its elements, powers, the nature of loads, and others arising in the operation of electrical networks. In such cases, the measurement time should be divided into intervals, in which operating conditions characterized by certain values of measured parameters take place. Models of nonlinear loads should be developed for each of these intervals. Fig. 3 shows voltage and current harmonics at the network node with capacitor bank (CB) on. The charts demonstrate that after turning on the capacitor bank, there are sharp increases in the 5th harmonic factor of voltages (K^) in the three phases and in the value of the 5th harmonic of currents (I5) in B and C phases. The parameter changes shown in the Figure are due to the harmonic resonant circuits that occur between the capacitor bank capacitance and the network node inductance [47]. Thus,
the entire measutement time interval should be represented by two intervals characterized by operating parameters markedly different in their value: the network parameters before switching on the capacitor bank and the network parameters after switching it on.
3. Establish correlation relationships between current harmonics and corresponding voltage harmonics.
Correlation relationships are analyzed to determine the mathematical methods to be used when modeling current harmonics by measured random variables. The analysis is performed for each of the time intervals established at the second step of the methodological approach. If it is assumed that the random values of current and voltage harmonics are governed by the normal distribution law, then correlation coefficients (r) are calculated to find out the strength of the relationship between them, which are then compared with the coefficients on the Cheddock scale [48]. When the distribution of random variables differs from the normal distribution, methods based on the application of order statistics or the replacement of measured values by their ranks are employed, using the Kenuy test [45, 49], the quadrant test [45, 50], the sign correlation Nelson test [45], and other tests. Table I shows the correlation coefficients rUnIn between current harmonics and voltage harmonics at the nodes of the connection of the aluminum smelter unit (ASU), pulp and paper mill (PPM), and railway substation (RWS) to the electricity supply network.
0
Table 1. Correlation coefficients rUnln.
n 3 5 7 9 11 13 17 19 23 25
runln ASU -0.04 0.02 -0.05 -0.05 -0.43 0.15 0.45 0.20 -0.06 -0.03
PPM 0.61 0.47 -0.06 0.62 0.74 0.56 0.18 0.40 0.59 0.53
RWS 0.29 -0.04 0.37 -0.03 0.12 0.05 0.06 0.10 0.75 0.41
The comparison of rUnIn and the coefficients on the Chaddock scale shows that in most cases, except for the coefficients highlighted in bold in Table I, there is either no correlation between the harmonics of currents and voltages or it is insignificant, which attests to the independence or weak dependence of voltage harmonic values on current harmonics, i.e., the probabilistic nature of both parameters, which makes it possible to apply corresponding mathematical methods when modeling current harmonics. If in some cases, there are correlation relationships between current harmonics and voltage harmonics, it is necessary to identify the existing relationships between them through the regression analysis.
4. Analyze the nth harmonic factors of voltage to select harmonic numbers for modeling.
Modeling can be performed not for all harmonics from the range of the 2nd to the 40th recommended in [5]. The numbers of harmonics whose KU(n) exceed the standard values are chosen for modeling [5]. For this purpose, it is necessary to compare the measured values of KU(n) with the standard values for 95% and 100% of the measurement time. Fig. 4 shows as an example the charts of the measured KU(n) at the nodes of the connection to the electricity supply network of the aluminum smelter unit, pulp and paper mill,
U(3,5)S, KU(7,11)
railway substation, and standard values - Kt
, K,
, K,
S> JVM;13)S> U(17)S> U(9,19,23,25)S
K,
, for 95% of the measurement
time. The Figure shows that at the node of connection of the aluminum smelter unit, the standard values of KU(11)
s, KU(13)S, KU(23)S, KU(25)S are exceeded, at the node of
connection of the pulp and paper mill, the same applies to
KU(5)s, KU(7)s, KU(9)s, KU(11)s, KU(13)s, KU(25)s, and at the node
of the connection of the railway substation, this holds for
KU(3)s, KU(5)s, KU(9)s, KU(11)s, KU(13)s, KU(17)s, KU(19)s, KU(23)s,
KU(25)s. Thus, at the node of connection of the aluminum smelter unit, one needs to model the 11th, 13th, 23rd, and 25th current harmonic; at the node of connection of the pulp and paper mill: the 5th, 7th, 9th, 11th, 13th, and 25th current harmonic; at the node of connection of the railway substation: the 3rd, 5th, 9th, 11th, 13th, 17th, 19th, 23rd, 25th current harmonic.
5. Analyze active and reactive current harmonic directions through the node of the connection to the nonlinear load network in order to select options for modeling.
The values of current harmonics for n > 1 are much smaller than the values of the first current harmonic. Their direction of flow can change for the first harmonic current for many minor reasons. At each time interval defined at the second step of the methodological approach, it is necessary to analyze the directions of the flow of harmonics of active and reactive currents through the node of the connection of the nonlinear load, i.e., relative to the measurement point in Fig. 1. The analysis is performed using the phase angle between the harmonics of voltages and currents 9UI(n). The active and reactive currents are directed from the network to the load if the angle 9UI(n) is within the range of 0 to n/2, i.e., in the 1st quadrant of the complex plane [51]. The active current is directed from the load to the network, and the reactive current is directed from the network to the load if the angle 9UI(n) is within the n/2 to n range, i.e., in the 2nd quadrant. Both currents are directed from the load into the network if the angle 9UI(n) is within the n to 3n/2 range, i.e., in the 3rd quadrant. The active current runs from the network to the load and the reactive current runs in the
2,5
£ 1,5
» ^
1
0,5 0
ASU PPM RWS
K
U(3, 5)S
K
U(7, 11)5
K
KU
U(13)S
U(H)S
KU
U(9,19,23,25)S
Harmonic number
Fig. 4. Measured values of KUin) at the nodes of connection of three enterprises to the network and the standard values for ten harmonics.
2
^ ox
§ ®
3 5
e
Table 2. Distribution of Harmonic Directions for Active and Reactive Currents of the Railway Substation
by the Quadrants of the Complex Plane.
eS° &
3 5 7 9 11 13 17 19 23 25
1 10.3 81.7 0.9 26.2 21.7 19.5 19.7 17.5 60.3 45.7
2 14.7 5.1 19.5 22.4 32.1 26.9 23.3 21.9 26.8 23.0
3 37.5 0.3 68.5 21.9 29.4 33.7 31.7 34.8 4.5 11.0
4 37.5 12.9 11.1 29.5 16.8 19.9 25.3 25.8 8.4 20.3
n
opposite direction if the 9UI(„) is within the 3n/2 to 2n range, i.e., in the 4th quadrant. Thus, the time of measurements of the harmonic parameters at each of the time intervals determined at the second step of the methodological approach can be divided into four more time sub-intervals, which correspond to the locations of the angle 9UI(n) in the 1st, 2nd, 3rd, and 4th quadrants of the complex plane. Sub-ranges (variants) prevailing over others in terms of the number of measurements, i.e., containing the greatest number of measurements out of their total number are selected for modeling. Table II presents the results of the analysis of the directions of the active and reactive currents of the railway substation using the angle 9UI(n). Each cell of
the table corresponds to a quadrant of the complex plane. The cell of the table shows the number of measured 9UI(n) belonging to the corresponding quadrant as a percentage of the total number of measurements 9UI(n), which is 1 440. To model harmonics of active and reactive currents, one should use the options highlighted in bold in the Table II.
6. Develop algorithms for modeling harmonics of active and reactive currents.
Modeling harmonics of active and reactive currents aims to identify their cumulative distribution functions, and then determine the values of Ian and Irn with a probability of 95% or other probability and to solve the system of equations (1). Based on the analysis of information from
The first algorithm
r
The second algorithm
Fig. 5. Flowchart of single current harmonic modeling algorithms.
C End ) Fig. 6. Flowchart of the computational program algorithm.
numerous measurements of current harmonics, we have arrived at the conclusion that two different approaches should be used to identify their cumulative distribution functions. These approaches are presented in Fig. 5 as two algorithms.
The notations used in the Figure are F(x) - cumulative distribution function, f (x,0) - probability density function, x95 - current harmonic value with a 95% probability. The first algorithm, including Blocks 2, 3, 4, 5, 6, and 7, is designed to identify the cumulative distribution functions of random variables of current harmonics corresponding to the known distribution laws. The second algorithm, including Blocks 8, 9, 10, 11, and 12, is designed to identify the cumulative distribution functions of random variables of current harmonics corresponding to mixtures of known distribution laws. Based on the algorithms, we have developed the computational program presented at step B.
7. Check the correctness of calculations with a given accuracy of the values of active and reactive currents.
After calculating the values of Ian and Irn it is necessary to check the correctness of the obtained values. To do this, one has to calculate the magnitude of current In from the values of Ian and Irn, and then compare the obtained value of In with the measured values of current. The calculated In value must not exceed the measured current values. If in some cases the calculated In value is out of range, the reasons for it should be found out and, if necessary, the modeling should be performed again. The reasons for the calculated value of In to exceed the measured current values can be the rounding of the values of Ian, Irn, In upwards, insufficiently accurate identification of the cumulative distribution function to describe harmonics of active and reactive currents, even though the identified cumulative distribution function corresponds to the measured harmonic values of active and reactive currents according to the goodness-of-fit tests. In the latter case, it is necessary to choose another distribution law that most accurately describes the distribution of harmonics of active and reactive currents.
B. Computational program for modeling harmonics of currents of nonlinear loads by measured parameters.
The flowchart of the computational program algorithm is shown in Fig. 6. The program identifies the cumulative distribution functions of harmonics of active and reactive currents and calculates their values with a given probability.
Block 1. Input the series of random values of active or reactive current harmonic X: x1, x2, ..., x, ..., xm, where m is the number of elements of the X series.
Block 2. Determine the number of histogram intervals S, their boundary values xs, the number of random variables of the series X that fall into the sth interval ms, and construct the current histogram. Based on numerous operations of processing the harmonic parameters measurement results and taking into account suggestions available in the
published research [52, 53], this study assumes that the number S is 19 if m > 1000; S is 15 if 500 < m < 1000; S is 11 if m < 500.
Blocks 3-5. Compare the current histogram with the curves of twelve known probability density functions used in the computational program to make hypotheses about the law governing the distribution of current (normal, Weibull, exponential, gamma, Rayleigh, lognormal, beta, minimum value, maximum value, logistic, Maxwell and Cauchy distributions [45, 46, 54, 55]).
Block 6. Calculate the parameters describing the assumed current distribution law of the series of random variables X, according to the expressions presented in [45, 46, 54, 55].
Block 7. Select the test to check the series of random variables X for the presence of outliers depending on the assumed distribution law: 1) if the assumed distribution law is normal, the Irwin test [45, 56] is used to check the series X for outliers; 2) if the assumed distribution law differs from normal, a modification of the Darling test [45] is used to check it.
Blocks 8-9. Check the largest x(m) and smallest x(1) values of a series of random variables X (x^ < X(2) < ... < x(i) < ... < x(m)) arranged in ascending order for outliers with the Irwin test. The statistics t, t* are calculated using the expressions
I 1 m 2
given in Fig. 6, where a = J-(xi - x) . Then they
are compared with the critical value t95 given in the reference handbooks.
Blocks 10-11. Check the largest x{ni) and smallest x^ values of a series of random variablesX (x^ < x{1) < ... < x(i) < •• < x(m)) arranged in ascending order for outliers with a modification of the Darling test.
The statistics L, L* are calculated using the expressions given in Fig. 6, where, F(xim)), F(x(i=l)), F(x(V)) - values of the cumulative distribution function for corresponding x. The statistics are then compared with the a-quantile value of the x2-distribution %2a [2 (m -1)] given in the reference handbooks.
Block 12. Replace outliers with neighboring elements
of the series X as x(m) = x(m - x(1) = x(2).
Blocks 13-14. Check the fit of the assumed current distribution law to the histogram using Pearson's goodness-of-fit test %2ex <xCR [45]. Pearson's goodness-of-fit test statistic is calculated using the expression
JEX
= ^ (ms - mps )2
(6)
.=i mps
where ps is the theoretical probability of the random variables of the series X falling into the sth interval, which is calculated as the difference between the values of the cumulative distribution function of the series X at the end and at the beginning of the sth interval
ps = Fx+1) - F X), (7)
where F(x) is the cumulative distribution function of the series X, calculated as
0,18
0,12
*" 0,06
0
1st component
y^T^Bv ^^^^^ 2nd component
11111ll1111..
<-
4a,
66
I, A
OO o
o <N
2, 2, 2,
4a2
Fig. 7. Current histogram with two peaks.
F ( x) = J f ( x) dx,
mps is the theoretical frequency of occurrence of random variables of series X in the sth interval. Then the calculated
value %2ex is compared with the critical value %2CR given in
[57].
Blocks 15-21. Identify the current probability density function using the methods of decomposing the mixtures of distributions.
First, based on the shape of the current histogram, a hypothesis is made about the number of components in the mixture and the types of their probability density functions. The twelve known probability density functions used in the computational program are proposed as probability density functions for the components of the mixture. Next, a current probability density function is formed as a weight sum of k components of the mixture [55, 58-62]
f (x, ©) = £ Cj Fj (x, 0j), (9)
j=i
where f (x,0) is the desired current probability density function; k > 1 is a natural number; ij is the known probability density function of the jth component of the mixture; 0 = (cj,..., ck, 0j,..., 0k) is the vector of parameters of the mixture components to be determined; cj is the weight of thejth component, Cj > 0, j = 1, ..., k; c1 + .+ ck = 1; 0j is the parameter vector of the jth component of the mixture.
The vector 0 is defined using a series of random variables X as follows. The initial approximations of the vector 0 parameters are set based on the histogram analysis. For example, Fig. 7 shows a current histogram with two peaks. We put forward a hypothesis that the histogram has two components, and their probability density functions are captured by the normal distribution law, which is characterized by two parameters: the mean
and standard deviation (c).
The vector 0j, in this case, consists of two parameters, i.e., 01 =(^,0^ and 02 =(^2,c2), and the vector 0 consists of six ones, i.e., 0 = (c1, c2, o1, c2). The initial values ^, o1, c2 are determined from the histogram as follows. To determine c1, c2, it is assumed that the total area of the
histogram under the curves of the mixture components is
(8)
1. Values of c1, c2 are set by values equal to the fractions of the areas under the corresponding curves of the mixture components. In the case of the normal distribution, the position of the peaks of the histogram relative to the >>-axis allows determining the mean for each of the components. Half the width of the shape under the component probability density curve makes it possible to determine the standard deviation since under the normal distribution, the interval of "four standard deviations" covers 99.99% of the values of the series of random variables X [63].
If the histogram components have distribution laws different from the normal law, the initial approximations of the vector 0 parameters are calculated using censored samples obtained from the series of random variables X [46, 64]. Then 0 is refined by solving the optimization problem
0 = argmin V
(ms - mps )2
s=1
mPs
subject to the following constraints:
-1 = 0,
j=1
- C < 0,
cj - 1 < 0,
- 0j <
(10)
(11)
(12)
(13)
(14)
The constraints correspond to the following: (11) - the sum of all weights must equal 1; (12) - all weights are positive values; (13) - the range that must encompass the values of the weights; (14) - all parameters of the vector 0j are positive values.
Expression (6) is related to coefficients Cj and parameters of the vector 0j by expressions (7)-(9), which show that jEx is a nonlinear function of 0 = (c1, ..., ck, 01, ..., 0k), hence, the optimization problem is nonlinear as well. In the computational program developed, the optimization problem is solved by the method of the generalized reduced gradient. As a result of obtaining the solution, the refined parameters of the vector 0 are calculated to describe the law of the current distribution.
After refining the parameters of the vector 0 and substituting them into (9), we check whether the obtained
x
A B C E F a H I J K L « ■ o P a H S T J v w I X I Y I AA AB
1 Ia5 After combining intervals Frequency in p.ii,
2 1,27075 Characteristics x, ms Interval Interval P. ms mys On, - mp, Y Interval ms mp, (m, - uij)s Interval ph Pt
3 0,79657 ¡1! = 1176 1.9E-16 1 2.3E-29 number from to «V. from to "P. from to
4 1,77933 ■^ntin ~ 1.9E-16 0,20554 41 0,03369 1 1.9E-16 0,20554 0,03389 42 39,65539 0,115401554 2E-16 0,2055 42 39,8 554 0,115401554 2E-1I 0,2055 0,03571 0,03392
5 0,63595 xmx 3,90531 0,41106 77 0,11364 2 0,20554 0,41106 0,07975 77 93,78416 3,003769558 0,2055 0,4111 77 93,7842 3,003789556 0,2055 0,4111 0,06546 0,07982
6 1,34170 s 19 0,61663 130 0,22196 3 0,41108 0,61663 0,10632 130 127,3622 0,05379915 0,4111 0,8166 130 127,362 0,05379915 0,4111 0,6166 0,11054 0,10841
7 1,72906 AK= 0,20554 0,62217 159 0,3443 4 0,81663 0,62217 0,12235 159 143,8605 1,566814614 0,6166 0,8222 159 143,68 1,586814614 0,6166 0,8222 0,1352 0,12245
8 1,39656 1,02771 150 0,46626 5 0,82217 1,02771 0,12396 150 145,8004 0,120963383 0,8222 1,0277 150 145,8 0,120963363 0,6222 1,0277 0,12755 0,12409
9 0,31092 1,23325 144 0,56443 6 1,02771 1,23325 0,11614 144 136,5634 0,402723852 1,0277 1,2333 144 136,563 0,402723652 1,0277 1,2333 0,12245 0.11624
10 2,07614 1,4366 126 0,66655 7 1,23325 1,4386 0,10213 126 120,1002 0,289817922 1,2333 1,4388 126 120,1 0,289817922 1,2333 1,4388 0,10714 0,10221
11 1,32632 Parameters FM 1,64434 №0 0,77161 8 1,4388 1,64434 0,08505 100 100,0222 4.92303E-06 1,4388 1,6443 100 100,022 4.92303E-06 1,4386 1,6443 0,08503 0,08513
12 1,91531 a = 1,60664 1,64966 67 0,63909 9 1,64434 1,64986 0,06748 67 79,357 1,924158695 1,6443 1,8499 67 79,357 1,924158695 1,6443 1,6499 0,05697 0,06754
13 1,57715 ß- 1,32525 2,05542 64 0,6903 10 1,84968 2,05542 0,05121 54 60.22395 0,643225796 1,8499 2,0554 54 60,224 0,643225796 1,8499 2,0554 0,04592 0,05126
14 0,22616 2,26097 33 0,92756 11 2,05542 2,26097 0,03728 33 43,84669 2,683227804 2,0554 2,261 33 43,8467 2,683227804 2,0554 2.261 0,02806 0,03732
15 1,55313 2,46651 30 0,95366 12 2,26097 2,46651 0,0261 30 30,69524 0,015746841 2,261 2,4665 30 30,6952 0,015746841 2,261 2,4665 0,02551 0,02612
16 0,92551 Note1. Formulas in cells D3:D5 use a 2,67205 23 0,97126 13 2,46651 2,67205 0,0176 23 20,69888 0,255617706 2,4665 2.6721 23 20,6969 0,255617706 2,4665 2,6721 0,01956 0,01762
17 1,74214 2,67759 13 0,96273 14 2,67205 2,67759 0,01145 13 13,46456 0,016028588 2,6721 2,8776 13 13,4646 0,016026588 2,6721 2,6776 0.01105 0,01146
18 0,72695 In cells A2:A(m +1]. In the general case, the numberof random variablesmcanbe arbitrary. In the example, m = 1176, so a series of random variables is placed in cells A2:A1177. 3,06314 11 0,96993 15 2,87759 3,08314 0,00719 11 6,459147 0,763169518 2,8776 3,0631 11 8,45915 0,763189518 2,8776 3,0831 0,00935 0,0072
19 0,1044 3,26666 7 0,9943 16 3,08314 3,28866 0,00437 7 5,137931 0,674644076 3,0631 3,2667 7 5,13793 0,674844076 3,0831 3,2887 0,00595 0,00437
20 0,7976 3,49422 4 0,99666 17 3,26666 3,49422 0,00257 4 3,019623 0,318297522 3,2687 3,4942 4 3,01962 0,318297522 3,2687 3,4942 0,0034 0,00257
21 0,0065 3,69976 4 0,99632 18 3,49422 3,69976 0,00146 4 1,718487 3,029000848 3,4942 3,6998 3,4942 3,6998 0,0034 0,00146
22 0,27334 3,90531 2 0,99913 19 3,69976 3,90531 0,00061 2 0,947667 1,166559317 3,6998 3,9053 6 2,66615 4,16874672 3,6996 3,9053 0,0017 0,00081
23 1.53044 Sum 1176 1174,978 17,067 Sum 1176 1174,98 17.039 | Sum 1,00 1,000
24 1,04246
25 1,17709 Critical value of statistics with a probability of 0,95 and ■ Histogram —
2fl 1,07446 the number of degrees of freedom/
27 1,87392
f=S-y 1= 16 fcR 26,29623
28 0,613553 ■0.14 -
23 1,14563 17,039 0.12 -
38 0,71506 iill ^ 0.1 -¿0.08 -\ 0.06 - j
31 1,56259 / \
32 1,6647 /Vote: fj;- theoretical probability of the random variables of the series falling into the ith interval of the histogram; mp,-theoretical \
33 0,71342 \
34 1,99524 0.04 - J r S
35 2,35061 frequency of occurrence of random variables of series in the ith 0.02 -1 >
15 1,59762 interval of the histogram; jjJ.- relative theoretical frequencies of 0 -*- 1 1 li 11- -
37 0,09157 i- n ■i in a CO
38 1,39653 randomvariables falling into the histogram intervals;^' - relative O O T- r- I- ^ n < n
39 0,30546 experimental and theoretical frequencies of random variables falling I, A
40 1,66002 into the histogram intervals.
41 1.11755
Fig. 8. Sheet 6 of the program - Blocks 13-14.
A | B | C | O ' E F | G | H i J | K [M | N | 0 | P | Q | R
1 Ia5 I f alculation of the current value/95 inthe case of a known probability distribution la 19 Calculation of the current value /95 in the case of using the
2 1.27075 1.738619
3 4 5 0.79657 1.77933 0.881156 0.949028 1.026702
0.83595 Parameters F(jc) Calculation /95 Constituent n a
6 7 8 9 10 1.34178 a - 1.806843 2.432319 0.976837 1.303678 1.14 1.60128 1.663351 1 1.19638 0.21062 0.683631
1.72906 ß = 1.325254 F(x) = 0.95 2 1.80796 0.19322 0.316369
1.39656 0.31092 2.07814 3
Calculation/^ 4
5
11 12 1.32632 1.98531 1.078572 0.919175
Constituent Calculation /95
13 14 1.57715 0.22616 0.827141 0.868184 Fl« = 0.68359 'k = 2.001777
F2« = 0.26641 FM = 0.95
15 16 1.55313 0.92551 1.438836 1.205923 Fj« =
F4« = Calculation/^
17 1.74214 1.512808 F5« =
18 0.72695 1.567007
Fig. 9. Sheet 8 of the program - Block 22.
probability density function corresponds to the histogram form using Pearson's goodness-of-fit test [45, 52, 65]. If this test fails, we return to Block 15 to propose other number and/or types of components of the mixture of distributions. If the condition is satisfied, the algorithm proceeds to Block 22 to determine the current value with a 95% probability by solving equation F(x) = 0.95.
The computational program developed based on the algorithm shown in Fig. 6 is implemented in MS Excel and Visual Basic programming environment for Windows applications. The program presents probability density function curves of twelve distributions most commonly used in various analyses referred to above. A dedicated program was written in the Visual Basic programming environment for Windows applications to calculate parameters of twelve distributions. This program consists of twelve subprograms, each intended to compute parameters of one particular distribution.
The computational program consists of eight MS Excel sheets. The first sheet implements Blocks 1-2 of the algorithm, the second - Blocks 3-5, the third - Block 6. Blocks 7-11 of the algorithm are implemented on Sheet 4, Block 12 - on Sheet 5, and Blocks 13-14 - on Sheet 6. Sheet 6 is shown in Fig. 8. Blocks 15-21 of the algorithm are implemented on Sheet 7 of the program, and Block 22 is implemented on Sheet 8. The last sheet of the program is presented in Fig. 9.
Iv Application of the developed methodological approach, algorithms, and computational program to modeling current harmonics of railway traction transformers.
Relying on the developed methodological approach, algorithms, and computational program, we have studied and modeled the 3rd and 5th harmonics of active and reactive currents at the nodes of the connection of traction transformers to the electricity supply network for the four substations of the East Siberian Railway: Mysovaya, Tataurovo, Zaigrayevo, and Novo-Ilyinsk. Electrical power is supplied to the traction network through these transformers to cover high-power nonlinear loads, i.e., electric locomotives. The calculations and analysis have indicated that current and voltage harmonics can be considered as independent random variables, and, therefore, the corresponding mathematical methods can be used. The distribution laws of most harmonics of active and reactive currents do not obey the known distribution laws but have more complex forms, i.e., mixtures of known laws, such as normal, Weibull, exponential, and gamma distributions. The distribution laws of the 3rd and 5th harmonics of active and reactive currents are shown in Table III.
The notations employed in the Table are as follows: W - Weibull distribution, G - gamma distribution, E - exponential distribution, W + W - a mixture of two
Table 3. Harmonic Distribution Laws for Active and Reactive Currents.
Currents
Substation
Mysovaya
Tataurovo
Zaigrayevo
Novo-Ilyinsk
W + W
W + G
E
N + G
Irn W + G G W + G N + W
Ian W + W W + W + W W + E W + W
N + W
N + W
W
N + W
n
I,
an
I
Table 4. Harmonic Values of Active and Reactive Currents Calculated with a 95% Probability, A.
n Substation
Mysovaya Tataurovo Zaigrayevo Novo-Ilyinsk
3 Ian 3.7 - 71.0 3.8 - 9.2
Irn - 2.8 - 73.5 2.5 - 7.2
5 Ian 2.4 - 15.6 2.2 1.7
Irn - 2.5 - 55.2 2.9 3.1
Table 5. Comparison of Ku(3), KU(5) calculated before and after the "Harmonics" software upgrade.
Substation n Ku(n)-P, % Ku(n)-A, % AKu(„), %
Mysovaya 3 1.12 1.10 1.79
5 1.16 1.06 8.62
Posolskaya 3 1.25 1.28 2.40
5 1.19 1.07 10.08
Kamensk 3 1.27 1.31 3.15
5 1.22 1.09 10.66
Selenginsky pulp and paper mill 3 1.29 1.32 2.32
5 1.23 1.10 10.57
Tataurovo 3 1.37 1.44 5.11
5 1.29 1.14 11.63
Severnaya 3 1.36 1.41 3.68
5 1.38 1.26 8.70
Rayonnaya 3 1.37 1.42 3.65
5 1.39 1.29 7.19
Motornaya 3 1.38 1.43 3.62
5 1.41 1.32 6.38
Zaigrayevo 3 1.49 1.56 4.70
5 1.76 1.73 1.70
Novo-Ilyinsk 3 1.66 1.58 4.82
5 2.07 2.02 2.42
Kizha 3 1.69 1.74 2.96
5 2.27 2.24 1.32
3 1.70 1.73 1.76
Petrovsk-Zabaikal'sky ---
J 5 2.56 2.54 0.78
Weibull distributions, W + G - a mixture of Weibull distribution and gamma distribution, N + W - a mixture of normal distribution and Weibull distribution, W + W + W - a mixture of three Weibull distributions, N + G - a mixture of normal distribution and gamma distribution, W + E - a mixture of Weibull and exponential distribution.
Table IV presents the values of the 3rd and 5th harmonics of active and reactive currents calculated with a 95% probability at the substation connection points with regard to the direction of flow. In the Table, the minus sign "-" means that the active or reactive currents are directed from the load to the network. Unsigned (positive) currents are directed from the network to the load. The correctness of the calculations of the harmonic values of active and reactive currents has been verified by checking as described at Step B.
The developed method of modeling nonlinear loads is implemented in the "Harmonics" software package [66] developed at the Melentiev ESI, SB RAS, for the purpose of improvement of the original package.
Table V shows the results of the comparison of KU(3), KU(5) for a segment of the network supplying twelve substations of the East Siberian Railway as calculated before (Ku{nyp) and after (Km„>A) upgrading the "Harmonics" software.
In the previous version of the software, the harmonic parameters were modeled considering the results of Ku(n) measurements at the four substations listed above. We used the harmonic values of active and reactive currents given in Table V to calculate Ku(n)_p, KU(n)_A. The
last column of Table V shows the differences AKu(n) of indices Ku(n)-P from indices Ku(n)_A, calculated as AKU(n) = 100 |Ku(n)-A - Ku(n)-P| / Ku(n)-P. The table shows that the differences in the Ku(3), Ku(5) indices are insignificant.
V. Conclusion
The methodological approach, algorithms, and computational program developed for modeling nonlinear loads to analyze non-sinusoidal conditions in high voltage electrical networks based on measurements of harmonic parameters, as evidenced by long-term modeling experience, can significantly reduce the time of modeling the nonlinear loads and harmonic parameters and improve the accuracy of harmonic values determined for active and reactive currents of nonlinear loads. With the help of the tools developed, we have performed the analysis and modeling of the 3rd and 5th current harmonics at the nodes of the connection of transformers of four substations of the East Siberian Railway to the electrical network that supplies electric power to highpower nonlinear loads, i.e., electric locomotives. The modeling results obtained indicate that the values of current harmonics are distributed according to different laws, including normal, Weibull, exponential, and gamma distributions, and the laws described by mixtures of distributions. The calculation results demonstrate that the approach developed to model nonlinear loads is effective and appropriate. The obtained results will be instrumental in developing the theory of non-sinusoidal operating conditions of electric power systems.
Acknowledgment
The research was carried out under State Assignment
Project (No. FWEU-2021-0001) of the Fundamental
Research Program of Russian Federation 2021-2030.
References
[1] Civil Code of the Russian Federation No. 14-FZ of October 21, 1994. (in Russian)
[2] "On Protection of Consumers' Rights," Federal Law of the Russian Federation No. 2300-1 of February 7, 1992. (in Russian)
[3] "On the Electric Power Industry," Federal Law of the Russian Federation No. 35-RF of March 26, 2003. (in Russian)
[4] L.I. Kovernikova, VV Sudnova, R.G. Shamonov et al., Electric power quality: current state, outstanding issues, and proposals for their solution. Novosibirsk: Nauka, 2017, 219 p. (in Russian)
[5] Electric power. Electromagnetic compatibility of technical means. Quality regulation of electric power in general-purpose electric power supply systems. GOST 32144-2013, Moscow: Standardinform, 2014. (in Russian)
[6] J. D. Greene, C.A. Gross, "Nonlinear Modeling of transformers," IEEE Transactions on Industry Applications, vol. 24, no. 3, pp. 434-438, 1988.
[7] E.F. Fuchs, A.S.M. Masoum, Power Quality in Power Systems and Electrical Machines. Academic Press, 2015.
[8] M.A.S. Masoum, E.F. Fuchs, D. Roesler, "Large signal nonlinear model of anisotropic transformers for nonsinusoidal operation, part II: magnetizing and core loss currents," IEEE Transactions on Power Delivery, vol. 6, no. 4, pp. 1509-1516, 1991.
[9] M.A.S. Masoum, E.F. Fuchs, "Transformer magnetizing current in harmonic power flow," IEEE Transactions on Power Delivery, vol. 9, no. 1, pp. 10-20, 1991.
[10] W. Dugui, X. Zheng, "Harmonic Model of Power Transformer," in Proc. International Conference on Power System Technology, POWERCON'98, August 1998, vol. 2, pp.1045-1049.
[11] T.J. Densem, P.S. Bodger, J. Arrillaga, "Three Phase Transmission System Modelling for Harmonic Penetration Studies," IEEE Transactions on Power Apparatus and Systems, vol. PAS-103, no. 2, pp. 310317, 1984.
[12] A.C. Williamson, "The Effects of System Harmonics upon Machines," in Proc. International Conference on Harmonics in Power Systems, UMIST, England, 1981.
[13] IEEE Power Engineering Society, Tutorial on harmonics modeling and simulation. 1998, 81 p.
[14] S.S. Waters, R.D. Willoughby, "Modeling Induction Motors for System Studies," IEEE Transactions on Industry Applications, vol. Ia-19, no. 5, pp. 875-878, 1983.
[15] F. Forcoles, J. Pedra, M. Salichs, L. Sainz, "Analysis of the induction machine parameter identification," IEEE
Transactions on Energy Conversion, vol. 17, no. 2, pp. 183-190, 2002.
[16] I.V Zhezhelenko, Higher harmonics in electric power supply systems of industrial enterprises, 4th edition, revised and expanded. Moscow: Energoatomizdat, 2000, 331 p. (in Russian)
[17] T. Zheng, E.B. Makram, A.A. Girgis, "Effect of Different Arc Furnace Models on Voltage Distortion,"
in Proc. IEEE International Conference Harmonics Quality of Power (ICHQP), Athens (Greece), Oct. 1998, pp. 1079-1085.
[18] G.W. Chang, Y.J. Liu, C.I. Chen, "Modeling Voltage-Current Characteristics of an Electric Arc Furnace Based on Actual Recorded Data: A Comparison of Classic and Advanced Models," IEEE Power and Energy Society General Meeting - Conversion and Delivery of Electrical Energy in the 21st Century, Pittsburgh (USA), pp. 1-6, 20-24 July 2008.
[19] Task Force on Harmonics Modeling and Simulation, "Modeling and simulation of the propagation of harmonics in electric power networks: Part I. Concepts, models, and simulation techniques," IEEE Transactions on Power Delivery, vol. 11, no. 1, pp. 452-465, 1996.
[20] Task Force on Harmonics Modeling and Simulation, "Characteristics and Modeling of Harmonic Sources - Power Electronic Devices," IEEE Transactions on Power Delivery, vol. 16, no. 4, pp. 791-800, 2001.
[21] A.J. Collin, J.L. Acosta, I.H. Gil, S.Z. Djokic, "Component-based Aggregate Load Models for Combined Power Flow and Harmonic Analysis," in Proc. 7th Mediterranean Conference on Power Generation, Transmission Distribution and Energy Conversion, Agia Napa (Cyprus), 2010.
[22] A.J. Collin, J.L. Acosta, I.H. Gil, S.Z. Djokic, "An 11kV Steady State Residential Aggregate Load Model. Part 1: Aggregation Methodology," in Proc. IEEE PES PowerTech, Trondheim (Norway), June 2011.
[23] A.J. Collin, G. Tsagarakis, A.E. Kiprakis, S. McLaughlin, "Simulating the time-varying harmonics of the residential load sector," in Proc. 16th IEEE International Conference on Harmonics and Quality of Power (ICHQP), Bucharest (Romania), 25-28 May 2014.
[24] A. Miranian, K. Rouzbehi, "Nonlinear Power System Load Identification Using Local Model Networks,"
IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2872-2881, 2013.
[25] H. Renmu, M. Jin, D.J. Hill, "Composite Load Modeling via Measurement Approach," IEEE Transactions on Power Systems, vol. 21, no. 2, pp. 663-672, May 2006.
[26] B.H. Kim, H. Kim, B. Lee B, "Parameter Estimation for the Composite Load Model," Journal of International Council on Electrical Engineering, vol. 2, no. 2, pp. 215-218, 2012.
[27] D. Karlsson, J.D. Hill, "Modelling and identification of nonlinear dynamic loads in power," IEEE Transactions on Power System, vol. 9, no. 1, pp. 157-166, February 1994.
[28] J. Balcells, M. Lamich, E. Griful, M. Corbalan, "Influence of Data Resolution in nonlinear Loads Model for Harmonics Prediction," in Proc. IECON 2016 - 42nd Annual Conference of the IEEE Industrial Electronics Society, Florence (Italy), 23-26 Oct. 2016.
[29] M. Lamich, J. Balcells, M. Corbalan, E. Griful, "Nonlinear Loads Model for Harmonics Flow Prediction, Using Multivariate Regression," IEEE Transactions on Industrial Electronics, vol. 64, no. 6, pp. 4820-4827, June 2017.
[30] M.T. Au, J.V Milanovic, "Development of Stochastic Aggregate Harmonic Load Model Based on Field Measurements," IEEE Transactions on Power Delivery, vol. 22, no. 1, pp. 323-330, January 2007.
[31] S.A. Ali, "A Norton Model of a Distribution Network for Harmonic Evaluation," Energy Science and Technology, vol. 2, no. 1, pp. 11-17, 2011.
[32] E. Thunberg, L. Soder, "A Norton approach to distribution network modeling for harmonic studies,"
IEEE Transactions on Power Delivery, vol. 14, no. 1, pp. 272-277, January 1999.
[33] Probabilistic aspects task force of the harmonics working group subcommittee of the transmission and distribution committee, "Time-varying harmonics: part
I - characterizing measured data," IEEE Transactions on Power Delivery, vol. 13, no. 3, pp. 938-944, July 1998.
[34] Probabilistic aspects task force of the harmonics working group subcommittee of the transmission and distribution committee, "Time-varying harmonics: part
II - Harmonic Summation and Propagation," IEEE Transactions on Power Delivery, vol. 17, no. 1, pp. 279-285, January 2002.
[35] GIGRE Working Group C4.605, "Modelling and Aggregation of Loads in Flexible Power Networks," CIGRE Technical Brochure 566, February 2014.
[36] GIGRE Working Group C4.605, "Modelling and Aggregation of Loads in Flexible Power Networks," Electra, no. 272, pp. 63-69, February 2014.
[37] VC. Luong, L.I. Kovernikova, "Analysis of the harmonic mode at the node of the connection of a traction substation to the supply network," in Proc. VI All-Russian Scientific and Practical Conference "Scientific Initiative of Foreign Students and Postgraduate Students of Russian Universities, " Tomsk (Russian Federation), April 24-26, 2013, vol. 1, pp. 286-291. (in Russian)
[38] L.I. Kovernikova, "Some results of research into harmonics in the high voltage networks with distributed nonlinear loads," Przeglad Elektrotexchniczny, no. 11, pp.239-243, 2013.
[39] L.I. Kovernikova, "Results of the research into the harmonics of loads connected to the nodes of high voltage network," Renewable Energy and Power Quality Journal, vol. 1, no. 12, pp. 553-558, April 2014.
[40] L.I. L.I. Kovernikova, VC. Luong, "Assessment of the steadiness of time series of measured parameters," in
Proc. All-Russian scientific and practical conference with international participation "Improving the efficiency of energy production and use in Siberia", Irkutsk, Russian Federation, April 22-26, 2014, vol. 2, pp. 260-264. (in Russian)
[41] L.I. Kovernikova, "Analysis of probabilistic properties of harmonic currents of loads connected to the high voltage networks," Energy Renewable and Power Quality Journal, vol. 1, no. 13, pp. 575-580, April 2015.
[42] L.I. L.I. Kovernikova, VC. Luong, VK. Ngo, "Special aspects of harmonic mode parameters at the node of the connection of a pulp and paper mill to the supply network," Intellektualnaya elektrotekhnika, no. 3, pp. 15-33, 2018. (in Russian)
[43] J. Arrilaga, D. Bradley, P. Bodger, Power System Harmonics: Translated from the English, Moscow: Energoatomizdat, 1990, 320 p. (in Russian)
[44] L.A. Zhukov, I.P. Stratan, Steady-state operation modes of complex electrical networks and systems: Calculation methods. Moscow: Energiya, 1979, 416 p. (in Russian)
[45] A.I. Kobzar, Applied mathematical statistics. For engineers and academics, 2nd ed., rev. Moscow: FIZMATLIT, 2012, 813 p. (in Russian)
[46] B.Yu. Lemeshko, S.B. Lemeshko, S.N. Postovalov, E.V. Chimitova, Statistical analysis of data, modeling of probabilistic patterns and their study. A computerized approach. Novosibirsk: Novosibirsk State Technical University, 2011, 888 p. (in Russian)
[47] L.I. Kovernikova, "Resonance Modes at Harmonics Frequencies in Electrical Networks," in Proc. E3S Web Conf., 2020, ENERGY-21 - Sustainable Development & Smart Management, 2020, vol. 209, URL: https:// doi.org/10.1051/e3sconf/202020907006.
[48] Shumetov V G., O.A. Kryukova, Methodology and practice of data analysis in management: Methods of univariate and bivariate analysis. Orel: Publishing house of the Orel Branch of the Russian Presidential Academy of National Economy and Public Administration, 2013, 177 p. (in Russian)
[49] M. G. Kenuy, Quick statistical calculations. Simplified methods of evaluation and verification: Handbook. Moscow: Statistika, 1979, 69 p. (in Russian)
[50] Ya. Gayek, Z. Shidak, The theory of rank-order test: Translated from the English. Moscow: Nauka, 1971, 376 p. (in Russian)
[51] A.E. Emanuel, Power definitions and physical mechanism of power flow. John Wiley & Sons, 2010, 264 р.
[52] Applied statistics. Rules for checking the goodness of fit of the theoretical distribution to the experimental one. Part 1. Chi-square type tests, GOST R 50.1.033-2001. Moscow: Standardinform, 2006. (in Russian)
[53] P.V. Novitsky, I.A. Zograf, Estimation of errors in measurement results, 2nd ed., rev. and exp. Leningrad: Energoatomizdat, Leningrad branch, 1991, 304 p. (in Russian)
[54] N. Hastings, B. Peacock, Statistical Distributions: Translated from the English by A. K. Zvonkin. Moscow: Statistika, 1980, 95 p. (in Russian)
[55] R.N. Vadzinsky, Handbook on probability distributions. Saint Petersburg: Nauka, 2001, 295 p. (in Russian)
[56] J.O. Irwin, "On a criterion for the rejection of outlying observations," Biometrika, vol. 17, no. 3-4, 1925.
[57] L.N. Bolshev, N.V. Smirnov, Mathematical statistics tables. Moscow: Nauka. Chief Editorial Board of Mathematics and Physics Publications, 1980, 416 p. (In Russian)
[58] S.A. Ayvazyan, VM. Bukhshtaber, I.S. Enyukov, L.D. Meshalkin; Ed. by S.A. Ayvazyan, Applied statistics: classifications and dimensionality reduction: A handbook. Moscow: Finansy i Statistika, 1989, 607 p. (in Russian)
[59] O.K. Isayenko, VYu. Urbach, "Decomposing mixtures of probability distributions into their components,"
Itogi nauki i tekhn. Ser. Teor. veroyatn. Mat. stat. Teor. kibernet., 13, All-Russian Institute for Scientific and Technical Information, Moscow, 1976, vol. 13, pp. 3758. (in Russian)
[60] VYu. Korolyov, The EM algorithm, its modifications, and their application to the problem of decomposition of mixtures of probability distributions: A theoretical review. Moscow: Publishing House of the Institute of the Russian History, RAS, 2007, 94 p. (in Russian)
[61] K.V Vorontsov, Lectures on statistical (Bayesian) classification algorithms, Computing Centre of the Russian Academy of Sciences, Moscow, 2008. Available at: http://www.ccas.ru/voron/download/ Bayes.pdf (Accessed on September 20, 2017). (in Russian)
[62] E.S. Wentzel, L.A. Ovcharov, Probability theory and its engineering applications, 2nd ed., reprint. Moscow: High School, 2000, 480 p. (in Russian)
[63] L. Zaks, Statistical estimation. Moscow: Statistika, 1976. (in Russian)
[64] G. Kulldorff, Contributions to the theory of estimation from grouped and partially grouped samples. Moscow: Nauka, 1966, 176 p. (in Russian)
[65] VI. Denisov, B.Y. Lemeshko, S.N. Postovalov,
Applied statistics. Rules for checking the goodness of fit of the theoretical distribution to the experimental one. Methodological guidelines. Part I. x2-type tests. Novosibirsk: NSTU Publishing House, 1998, 126 p. (139 p. with Addendum). (in Russian)
[66] L.I. Kovernikova, "The "Harmonics" software package for analysis and normalization of harmonic modes in high-voltage networks," in Proc. of the International Scientific and Practical Conference "Management of electric power energy quality", Moscow, Russian Federation, November 23-25, 2016, printed in 2017, Moscow: OOO "Centr poligraficheskikh uslug "Raduga", pp. 155-164. (in Russian)
Lidiia Kovernikova is a Senior Researcher at the Melentiev Energy Systems Institute of the Russian Academy of Sciences, Irkutsk, Russia. She graduated from Novosibirsk Electrotechnical Institute in 1976. L.I. Kovernikova received her Ph.D. degree in Engineering from the Siberian Energy Institute SB RAS in 1995. Her areas of interest include modeling the non-sinusoidal conditions of electric power systems, analyzing the power quality of electrical networks, and developing technical means to improve the power quality.
Luong Van Chung is a Junior Researcher at the Institute of Ship Design, Hanoi, Vietnam. He graduated from Irkutsk National Research Technical University, Irkutsk, Russia in 2014. Luong Van Chung received his Ph.D. degree in Engineering from the Melentiev Energy Systems Institute SB RAS in 2020. His areas of interest include the development of the power supply system for ships, and power quality.
