Modeling of loads dependent on harmonic voltages

Electric Power Systems Research 152 (2017) 367–376

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

Modeling of loads dependent on harmonic voltages Marcelo Brunoro a,b,∗ , Lucas Frizera Encarnac¸ão b , Jussara Farias Fardin b a b

IFES — Federal Institute of Espírito Santo, Av. Vitória, 1729-Jucutuquara, Vitória, ES 29040-780, Brazil UFES — Federal University of Espírito Santo, Av. Fernando Ferrari, 514-Goiabeiras, Vitória, ES, 29075-910, Brazil

a r t i c l e

i n f o

Article history: Received 24 March 2017 Received in revised form 26 July 2017 Accepted 31 July 2017 Keywords: Load modeling Harmonics Parameter estimation Power system harmonics Power system analysis computing

a b s t r a c t Electrical system models are important to allow the accomplishment of several studies aiming at reducing losses, and improving power quality, among others. In the current context, carrying out harmonic analysis becomes necessary due to the increasing insertion of nonlinear loads in the system. Since the networks voltage and current profiles are strongly affected by the load behavior, their modeling is essential to perform such analyses. Different from the other harmonic models found in the literature, this paper presents a new proposal of load modeling, combining the ZIP model and a cross admittance matrix. This combination brings together the benefits of load characterization with the traditional ZIP model, which provides some physical knowledge about the load, as well as the frequency crossing given by an admittance matrix. In this way, it is possible to accurately determine the harmonic power injection in the load bus, by knowing its harmonic voltage. The present proposal still considers a limitation for the ZIP coefficients, in order to identify the power ratio in terms of constant impedance, constant current and constant power. In addition to discussing the method for determining the load model parameters, using exhaustive search and multiple linear regression, a real case study with data obtained with a power quality meter exemplifies the model application in an electronic load. The results show that the proposed harmonic model was able to represent the load with high accuracy, and the parameters found provide information about the type of the modeled load. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Several factors have made the electric power system planning and control increasingly challenging to increase its robustness and compliance with normative criteria. Among these factors, we can mention: unbalance between generated and demanded power; distributed generation intensification with renewable energy resources using static power converters; nonlinear loads quantity increase with switched devices and saturated magnetic devices. Static converters are the largest nonlinear loads used in the industry for different purposes, such as adjustable speed drives and uninterrupted power supplies. Therefore, in order to maintain the electric power system with reduced energy losses and with the desired voltage and current profile, providing acceptable levels of quality, safety, availability, reliability, and economic viability, various studies and analyses about the system are needed [1–4].

∗ Corresponding author at: IFES — Federal Institute of Espírito Santo, Av. Vitória, 1729-Jucutuquara, Vitória, ES, 29040-780, Brazil. E-mail addresses: [email protected], [email protected] (M. Brunoro), [email protected] (L.F. Encarnac¸ão), [email protected] (J.F. Fardin). http://dx.doi.org/10.1016/j.epsr.2017.07.030 0378-7796/© 2017 Elsevier B.V. All rights reserved.

In this context, considering that the loads have a great influence on the busbar voltage profile and circulating network current, the load representation has a significant impact on the power system analysis and control functions [2,5]. However, the load models are the least known models among the various components of the system [6]. The different and more complex characteristics of modern loads in relation to older loads require detailed studies of their actual behavior in power flow studies and in harmonic analyzes. Many problems of power system analysis consider constant value loads ignoring its dependence on voltage, current and frequency. Thus, in the absence of appropriate models for the representation of actual loads, power system analyses may produce erroneous results [6,7]. All of these factors make it challenging to model the load for applications in the power system. Due to the increase in the use of nonlinear loads and electronic devices in the power system, the harmonic currents injection in the distribution networks has increased significantly, causing distortion in the voltage drops along the network. Consequently, this change in the sinusoidal voltage profile causes distortions even in the current of linear loads connected to the bus [1,4,6,8]. Currently, in addition to the nonlinear loads, the distribution systems are submitted to harmonic currents and voltages due to the presence of power inverters to interface solar and wind power plants [9]. The

368

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

existence of harmonic voltages and currents in the system causes serious problems such as transformers and conductors heating, electronic equipment malfunction and power demand increase. Consequently, harmonic monitoring has become an important task to ensure the power quality [1,6,8]. Therefore, it is fundamental to perform electric system harmonic analysis, adopting harmonic models of loads, aiming to increase the system efficiency through the application of techniques for the reduction of harmonic distortion [6,8,10,11]. Static models are widely used for the representation of loads in various analyses of the power system. The characterization of static loads may be done by the ZIP model, which represents the active and reactive powers of the load by components of constant impedance, constant current and constant power, considering voltage variations [6], and thus providing some physical significance on the load [12]. This model is widely used to represent modern residential and industrial loads in several studies on the power system as in stability analyses [5], power flow studies and power system planning [6], voltage and reactive power control [13], Conservation Voltage Reduction (CVR) [14], evaluation of the Available Transfer Capability (ATC) [15], islanding detection for inverter-based Distributed Generation (DG) [5] and fault location [16]. In Ref. [17], a time variant ZIP model was applied in the allocation of DGs in a mesh distribution system. A quasi real-time ZIP model was proposed in Ref. [18] in a CVR study. A history of load modeling including the ZIP model can be seen in Ref. [14]. On the other hand, models that consider harmonic voltage and current components are used to estimate the nonlinear loads impact on the system. Among the harmonic load modeling we mention: fixed current source model, where each harmonic current component of the load is represented by a constant current source [10,19]; Norton model, in which each load harmonic current component is represented by the Norton equivalent circuit, containing an admittance in parallel with a current source [20–22]; coupled Norton model, whose equation considers a constant current for each harmonic order of load current and an admittance matrix responsible for the harmonic coupling between components of several harmonic orders of voltage and current of the load [19,23]; and harmonically coupled admittance matrix model, which contains an array of admittances that relates voltages and currents of several harmonic orders [19]. This paper focuses on the proposition and parameters determination of a novel model for harmonic loads considering voltage variations. Thus, after determining the load model parameters, following the methodology proposed in this article, it is possible to accurately determine the active and reactive power harmonic injection in the load bus as its voltage function. This proposal associates elements widely used in load modeling, which is the ZIP model and a cross-admittance matrix. Unlike the other harmonic models presented previously, this combination associates the benefits of the characterization of loads by the model ZIP, as well as the frequency crossing given by the admittance matrix, necessary in nonlinear loads. The ZIP model provides some physical knowledge about the load, besides flexibility in the representation of different types of loads, considering voltage variations. Thus, the proposed model may represent the load in a widely diffused form in studies on the system through the components of constant impedance, constant current and constant power. The cross-admittance matrix provides the cross-interaction between several harmonic voltage components in the composition of the power and current of the modeled load, providing the representation of nonlinear loads. Some characteristics of the coupled Norton model are also incorporated since in this proposal there is a constant current component, but with a constant power factor, in addition to the cross coupling of frequencies given by the admittance matrix. This model can subsidize several harmonic analyses

in systems with the presence of numerous nonlinear loads. Such analyses can be used to detect loads that contribute significantly to the system voltages distortion through the injection of harmonic currents, as well as to the application of measures to mitigate unwanted harmonic components. The parameters of the ZIP portion of the fundamental component of the model are determined by exhaustive search, while the other parameters are obtained by multiple linear regression. All parameters are found from measured voltage and current data of the load. Such data should consider different operating points, covering the whole range of values in which the system is to be modeled. This study also proposes the use of the coefficients of the ZIP model in a predetermined range, and it can help in the interpretation of the type of load modeled in terms of the proportions of constant impedance, constant current and constant power. A case study, using actual data, is presented in this paper showing the results of this approach for a three-phase full-bridge diode rectifier. The proposed model was able to represent the tested nonlinear load with high accuracy, although it contains a lower number of parameters than the coupled Norton model for analyses that consider more than three harmonic orders. 2. Load modeling Load models must be compatible with the intended analysis. Thus, several models are proposed in the literature. Static models are widely used in load steady state studies, whereas dynamic models are typically used in load transient studies. However, in order to consider the load harmonic components, it is necessary to use harmonic modeling. This section shows the load models considered in the proposed modeling. 2.1. Polynomial model or ZIP model Static models express the characteristics of the load, at any instant of time, as a function of the load bus voltage, where the active and reactive power components are considered separately. This model is only indicated for situations in which the voltage and frequency variations are small or slow, where the steady state is reached quickly. The static model can also be applied when the interest is focused only on steady state analyses. Most composite loads have this behavior and static modeling can be employed [15,24]. The static model is usually represented in the literature as an exponential and polynomial model [25–27]. It is possible to represent loads with constant power, constant current or constant impedance characteristics, with exponents equal to 0, 1 or 2, respectively. The ZIP denomination can be used due to the model composition by constant impedance, constant current and constant power components. This representation is described in (1) and (2). P = P0 [pZ (

|V | 2 |V | ) + pI ( ) + pP ] |V0 | |V0 |

(1)

|V | 2 |V | ) + qI ( ) + qP ] |V0 | |V0 |

(2)

Q = Q0 [qZ (

where subscript 0 represents the initial operation condition of the respective variable and the coefficients pZ , pI and pP , whose sum is equal to 1, define the proportion of each component in the load active power characterization. The same reasoning can be applied to the reactive power coefficients [6,18,25]. In Ref. [3], the substation load modeling of the electric power distribution system is performed using the exponential and ZIP model with dynamic estimation of its parameters by the weighed least squares method in the recursive form. In Ref. [13], however,

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

an important optimization problem is addressed in smart grids, proposing an algorithm for voltage and reactive power control in distribution systems in which the load is represented by the ZIP model, whose parameters were selected based on urban load data. The values of the ZIP model coefficients for several loads can be seen in Refs. [14,18].

369

Optimization with least squares method are applied. The most significant characteristics of these methods are shown in Ref. [24]. Practical circumstances that cause disturbances are addressed in Ref. [30], such as capacitor switching, transformer tap change, load variation and current injection. It is important that the collected data from the system contain information about the characteristics to be modeled.

2.2. Harmonic modeling 3. Proposed model With the increasing presence of nonlinear loads connected to the power system, modeling that considers harmonic components becomes essential to estimate the impact of current loads or new clients in the system, as well as to investigate the effectiveness of harmonic mitigation techniques. In this case, it is indispensable to characterize the harmonic behavior of loads in studies on losses and on the voltage and current harmonic components [28,29]. When a linear load undergoes a harmonic voltage, it will only absorb the harmonic current of the same order. In contrast, nonlinear loads can absorb currents with different harmonic orders in relation to the applied voltage [28]. Thus, models that consider the effect of several harmonic voltage components on a given current harmonic component are more accurate depending on the load type to be represented [19,23]. The constant current source model suggests that they are independent of the load voltage, consequently disconsidering the harmonic interaction between different voltage and current frequencies, causing errors in the modeling of certain loads [19]. In the Norton model, each load harmonic component is represented as a constant current source in parallel with an impedance. Thereby, the load harmonic currents are related only to their same order harmonic voltages [30,31]. While the classical Norton model neglects the coupling between different harmonic frequencies between the voltage and current of the load, the coupled Norton model promotes such coupling using an admittance matrix. Disregarding this interaction between frequencies can lead to significant errors in the modeling of load behavior [23,31,32]. In the harmonically coupled admittance matrix model, it is assumed that each load current harmonic component is dependent on the various harmonic voltage components considered in the analysis, with the relation between the current and the voltages given by admittances. The load current components up to the N harmonic order are described in (3).

⎡

I1

⎤

⎡

Y11

Y12

Y13

···

Y1N

⎤ ⎡

V1

⎤

⎢ I ⎥ ⎢ Y1 Y2 Y3 · · · YN ⎥ ⎢ V ⎥ ⎢ 2⎥ ⎢ 2 2 2 2 ⎥ ⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ I ⎥ ⎢ Y1 Y2 Y3 · · · YN ⎥ ⎢ V ⎥ ⎢ 3⎥=⎢ 3 3 3 3 ⎥·⎢ 3 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢ . .. .. .. ⎥ ⎢ .. ⎥ .. ⎣ .. ⎦ ⎣ .. . . . . ⎦ ⎣ . ⎦ IN

YN1

YN2

YN3

···

YNN

(3)

VN

where Ih and Vh are the load current and voltage harmonic phasor components of order h, respectively; and Yhv is the admittance submitted to the harmonic voltage component of order v, used to determine the load current harmonic component of order h [28].

A widely used model called the polynomial or ZIP model, given by (1) and (2), represents satisfactorily loads for several steady state and transient studies. However, in the latter case, the ZIP model must be combined with models that represent the dynamic part of the load [33,34]. According to Ref. [30], the ZIP model is not able to consider the harmonic effects in the load representation. On the other hand, some approaches to load modeling consider the use of an admittance matrix, as in the Norton model where each load current harmonic component is dependent on the harmonic voltage component of the same order. Meanwhile, in the harmonically coupled admittance matrix model, given by (3), there is harmonic voltage dependence of several orders in each harmonic current component. Considering that the approach with the ZIP model does not accurately represent a load subject to harmonic voltage components, it is proposed, alternatively, that the active and reactive load powers be characterized by the combination of the ZIP model and Crossed Frequency Admittance Matrix model, providing the effect of the harmonic coupling. So, for each harmonic component of order h, the per phase load apparent power can be obtained as described in (4). In this proposal, the ZIP portion of the model will represent the load physical behavior, and its dependence. Concerning the harmonic voltage components, it is represented by the portion containing the cross admittance matrix. It is important to mention that the ZIP portion of the proposed model assumes greater relevance in the load characterization, and the portion containing the cross admittance matrix complements the load representation taking several harmonic voltage components into account. SLh = Sh + Sh

where the real and imaginary part of the apparent power Sh are obtained from (1) and (2), respectively, which were rewritten in (5) and (6); and Sh is described in (7) up to the N order. Ph = P0h [ pZ V h 2 + pI V h + pP ]

(5)

Qh = Q0h [ qZ V h 2 + qI V h + qP ]

(6)

There are different types of load models and methods for determining their parameters. In the measurement-based approach, indicated for situations where the model structure is known, techniques such as evolutionary strategies, artificial intelligence and systems identification can be used. Genetic Algorithms and Particle Swarm Optimization are also used to estimate the models parameters. Combinations of Genetic Algorithms and Particle Swarm

 N

S h =

|Vv |2 Yhv , withv = / h

(7)

v=1

where V h is given by (8); the relation between the ZIP coefficients are given by (9) and (10); and the admittance Yhv decomposed into real and imaginary part is given by (11). |Vh | |V0h |

(8)

pZ + pI + pP = 1

(9)

qZ + qI + qP = 1

(10)

Yhv = Ghv + jBhv

(11)

Vh = 2.3. Determining model parameters

(4)

where Ghv is a conductance and Bhv , a susceptance subject to the

voltage harmonic component of order v, used to determine S of order h.

370

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

In this proposal, the ZIP model coefficients pZ , pI , pP , qZ , qI and qP will be obtained for h equal to 1, which is, from the fundamental component of the current and voltage load data. These coefficients will also be used to determine the apparent power of the further harmonic orders. For each harmonic component h the values of P0h and Q0h will be estimated, as well the Yhv values. By choosing the measure that will be used as initial value, consequently the apparent power initial value can also be established. However, we decided to estimate this power initial value to increase the model accuracy, due to the presence of noise in the voltage and current measurements used to determine such power. The values of V0h must be established as one of the load voltage measures Vh used in the model parameter estimation process. In this case, the first measure of the estimation process data set was adopted, as suggested by Refs. [15,26]. Observing (5) and (6) it is noted that the various harmonic components magnitudes of the initial voltage V0h cannot be equal to zero. Thus, harmonic orders whose voltage magnitudes are zero or close to this value must be excluded from the model. For the case where there is symmetry in the load voltages and currents, where the even order harmonic components are equal to zero, the load model can represent its characteristics only with the components of odd order. In situations where the multiples of third order harmonic components are very close to zero, these components must also be excluded from the model. When applicable, the mentioned harmonic components elimination does not significantly compromise the proposed model response and it is advantageous because it results in the reduction of the number of parameters and consequently the number of points necessary to the identification process. 4. Proposed model parameters determination In order to identify the proposed model parameters, the least squares method was used. In this method, the measured values must cover the entire range of values in which the system is to be modeled. In this case, it is necessary for the system to receive disturbances that can be natural or artificial, such as load changes and transformer tap changes [30]. In a distribution system, the voltage variations cause changes in the load currents, which modify the conductors voltage drop, resulting again in the load voltage change. Measurement requirements must be observed so that no data problems occur in the identification of model parameters. Among them, we can mention: set measurements that represent different operating conditions, involving the situations that are intended to model and ensuring that the system of equations is linearly independent; Linear independence between the harmonic voltages components to exist only one parameters set; Number of measurements greater than the number of parameters to be determined; and harmonic voltage components amplitude sufficiently high to permit measurement and low enough to avoid problems due to nonlinearities. Energy quality meters can be used to determine the measurements for the model estimation and validation process. However, the meter must be able to provide the magnitude and angle values of the load harmonic voltage and current components to be modeled up to the highest harmonic order considered in the model. In this way, it is possible to calculate the load apparent power for each harmonic order h, as described in (12). SLh = PLh + jQLh = Vh · Ih ∗

(12)

where PLh and QLh are the active and reactive load powers, respectively; and Vh and Ih are the load voltage and current phasors, respectively, for the harmonic order h. Aiming at the application of the least squares method to identify the model parameters, the system of equations needs to be overde-

termined for obtaining a representative response in the presence of noises in the measurements [35]. Therefore, the number of samples from the parameter estimation process must be greater than the number of parameters estimated. Thus, the quantity M of measurements must be greater than the quantity of harmonic components considered in the model. The determination of the minimum quantity of samples considered in the model will be exemplified for two situations in which there are components to be excluded due to the reduced magnitude of their voltages. The first situation is for the case where the even-order components must be excluded, and the quantity of M measures that should be collected for the system of equation construction must meet the M > Modd relation, with Modd defined in (13). The second situation is for the case where the even and multiples of third order components are disregarded from the model, and the quantity of M measures to be collected is given by M > Mpos , with Mpos defined in (14). Modd =

N+1 2

(13)

Mpos =

2N + 3 − (−1)N 6

(14)

where N is the highest harmonic order considered in the model where M must always be an integer greater than 3, where the values of Modd and Mpos are not always integers. For example, for the load representation up to the 10th harmonic order (N = 10), considering only the odd-order components, the quantity of samples required is given by M > Modd . Thus, (13) has Modd = 5.5, so that M > 5.5, that is, the quantity of samples required for the estimation process must be an integer value greater than or equal to 6. Typically, the load representation is improved by increasing the amount of harmonics considered in the model, but this implies an increase in the number of measurements necessary for the parameter identification process. However, a high number of measures helps increase the immunity of the model parameters identifying process in the presence of measurement noises. On the other hand, if the quantity of measures is limited, (13) or (14) can help determine the largest harmonic order that the load model can represent. It should be emphasized that for the validation process of the parameters model, other measures are needed besides those which have already been used in the identification process. The process of identifying the parameters was divided into two stages. In the first stage, the parameters of the model for the fundamental component are estimated, including the ZIP model coefficients. Therefore, in this stage the following parameters will be determined: P01 , pZ , pI , pP , Q01 , qZ , qI , qP and Y12 , Y13 , · · ·, Y1N being N the highest harmonic order considered in the model. Thus, the interaction of the various voltage harmonic components in the fundamental component is provided by the admittances Y1v , with v assuming values of 2 to N. In the second stage, the model parameters for higher order components than the fundamental are estimated considering the ZIP coefficients determined in the first stage. Here the following parameters are obtained for each harmonic order h with values from 2 to N: P0h , Q0h and Yh1 , · · ·, Yhv , · · ·, YhN where v is different from h. These steps are detailed below and Fig. 1 shows the summary sequence for obtaining the parameters, highlighting the steps to determine the parameters of the fundamental component and the other harmonic components. 4.1. Step for determining parameters of the fundamental component The process of obtaining the fundamental component parameters of the load apparent power was divided into two parts. The

Higher order parameters

Set maximum and minimum limits for P01, Q01 and ZIP coefficients Set step size for P01, Q01 and ZIP coefficients Determine P01, Q01 and ZIP coefficients by exhaustive search

and reactive powers occur for |V1 | = |V01 |, which is, for V 1 = 1. Then, once the initial voltage value V0h is defined as the first load voltage sample of the estimation process, the maximum and minimum limits of the initial powers can be obtained for this condition. In a simplified way, in order to consider the entire admissible range of the active power initial value, the linear regression was performed

Determine vector Y1v by multiple linear regression, with v = 2, 3, ..., N Determine by multiple linear regression P0h, Q0h and the vector Yhv with v = 1, 2, ..., N, for each harmonic order h = 2, 3, ..., N, where v ≠ h.

between the active power measurements and the respective V 1 values. Then, the maximum limit for the active power initial value

Fig. 1. Summary sequence for parameter determination.

first part obtains the parameters of the Sh term and in the second part, the parameters of Sh , both of (4) with h equal to 1. In (4) for the Sh portion to be preponderant in the proposed model and its ZIP coefficients to represent the proportion of each load component as impedance, current and constant power, the term Sh is initially represented as an error that encompasses the measurement noises and the modeling error. Thus, from (4)–(6) and the relationships described in (9) and (10), the measure m of the real and imagim = V m · I m ∗ from nary part of the load apparent power, given by SL1 1 1 (12), can be represented as described in (15) and (16), respectively. Therefore, initially, the Sh term parameters with h equal to 1 is determined, where P01 , pZ , pI and pP are the active power parameters and Q01 , qZ , qI and qP are the reactive power parameters. m

m

was determined by making V 1 = 1 on the superior parallel curve of the regression line containing the farthest active power measured point. Similarly, the minimum limit of the initial value of the active power was determined using the inferior parallel curve to the regression line that contains the farthest active power measured point. The same process was adopted for reactive power. Fig. 3 illustrates the achievement of these limits in a case study. The regression equation can be seen in (17) and (18). ˆ PL1 = X · ˇ P1 + εP1

(17)

ˆ QL1 = X · ˇ Q1 + εQ1

(18)

where

⎡

(15)

m

m m QL1 = Q01 [qZ (V 1 2 − 1) + qI (V 1 − 1) + 1] + eQ 1

⎤

PL1 = ⎢ ⎢

(16)

M PL1

m

m and em are the errors for where V 1 = |V1m |/|V01 | from (8); and eP1 Q1 the m-th measure. Considering that the ZIP model coefficients indicate the load proportion with constant impedance, constant current and constant power, their values are expected to be limited between 0 and 1. Nevertheless, this limit is not considered in Refs. [3,25,34], where the ZIP coefficients come to assume negative values and, in some cases, values greater than 1. With multiple linear regression, it is possible to determine these coefficients, as well as the initial powers, but this method does not allow the parameter values to be restricted in the expected range. Then the exhaustive search method was used with the objective of minimizing the quadratic m and em of (15) and (16), respectively. In this sum of the errors eP1 Q1 way, the exhaustive search can find the values of the ZIP coefficients always limited to the range between 0 and 1, providing some physical information about the load from the portion of constant impedance, constant current and constant power given by the coefficients. As the determination of the parameters of the model is not performed in real time in which the measurements are performed, the use of a slower method does not compromise the result. In order to make the exhaustive search faster, it is essential that the search space be as narrow as possible. Thus, in the i-th iteration of the ZIP coefficients search algorithm, its limits are defined as follows: 0 ≤ piZ ≤ 1; 0 ≤ piI ≤ (1 − piZ ); 0 ≤ qiZ ≤ 1 and 0 ≤ qiI ≤ (1 − qiZ ). The pI maximum limit considers the fact that pZ + pI ≤ 1 since the

1 PL1

⎢ P2 ⎥ ⎢ L1 ⎥ ⎥, ⎥ ⎣ ... ⎦

m

m m PL1 = P01 [pZ (V 1 2 − 1) + pI (V 1 − 1) + 1] + eP1

371

relation described in (9) must be followed. The same is true for qI for the relation described in (10). It is known that the initial value of the fundamental portion of the active power P01 and of the reactive power Q01 is within the range of measured values. This fact can be used to define the initial powers limits of the exhaustive search. However, simply establishing P01 and Q01 between the maximum and minimum values of the active and reactive power measurements, respectively, would lead to the use of still high limits. Thus, by observing (5) and (6), for h equal to 1, it can be seen that the initial values of the active

Exhaustive search

Fundamental component parameters

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

⎡

1 QL1

⎡

1

⎤

⎡ 1 ⎤ εP1 ⎢ ⎥ ⎢

2 ⎥ ⎢ ⎥ ⎢1 V ⎥ ˆ0 ˇ ⎢ ε2P1 ⎥ P1 ⎢ 1 ⎥ ˆ ⎢ ⎥ X=⎢ , ˇ = e ε = ⎥ P1 ⎢ . ⎥ ⎢ . . ⎥ P1 ˆ1 ˇ . ⎣ P1 ⎢ .. .. ⎥ . ⎦ ⎣ ⎦ M 1

V1

1

V1

εP1

M

⎤





⎡ ε1 ⎤ Q1

⎢ Q2 ⎥ ⎢ ε2 ⎥ ˆ0 ˇ ⎢ L1 ⎥ ⎢ Q1 ⎥ ⎥ , ˇˆ Q1 = Q 1 e εQ1 = ⎢ ⎥ Q L1 = ⎢ ⎢ . ⎥ ⎢ . ⎥ ˆ1 ˇ ⎣ .. ⎦ ⎣ .. ⎦ Q1 M QL1

εM Q1

ˆ and ˇ ˆ are the vectors with the estimates for the regreswhere ˇ P1 Q1 sion lines coefficients of the active and reactive power, respectively; and X is full column rank. The ˆ symbol indicates that the variable is an estimate. The adjusted regression models are described in (19) and (20). 

ˆ ˆ L1 = X · ˇ P P1

(19)



ˆ ˆ L1 = X · ˇ Q Q1

(20)

where the estimates for the regression lines coefficients are described in (21) and (22). −1 T ˆ ˇ · XT · P L1 P1 = [X · X]

(21)

−1 T ˆ ˇ · XT · Q L1 Q1 = [X · X]

(22)

372

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

Thus the minimum and maximum limits of the active power are determined by (23) and (24), respectively, and the limits for the reactive power by (25) and (26). min P01

=

ˆ0 ˇ P1

ˆ1 +ˇ P1



ˆ L1 ) + min(PL1 − P



ˆ 1 = ␸1 T · ␸1

−1

· ␸1 T · S1

(37)

(23)



max ˆ0 + ˇ ˆ 1 + max(PL1 − P ˆ L1 ) P01 =ˇ P1 P1

way all parameters of the model, for the fundamental component can be estimated.

(24)

4.2. Step for determining parameters of high order harmonic components



min ˆ0 + ˇ ˆ 1 + min(QL1 − Q ˆ L1 ) Q01 =ˇ Q1 Q1

(25)



max ˆ0 + ˇ ˆ 1 + max(QL1 − Q ˆ L1 ) Q01 =ˇ Q1 Q1

(26)

For the exhaustive search algorithm it is necessary to establish the limits for the search space and the step size used in the iterative process to determine the parameters. The step size of each parameter is given by the parameter limits and by the amount of parts into which its search space should be divided. At the end of the iterative process, the parameters will be defined as the set of values that minimizes the sum of the squared errors between the measured and estimated values. From (9), (10) and the parameters already found, the coefficients pP and qP are given by (27) and (28), respectively.

With the parameters of the fundamental component, it is possible to determine the parameters of the other harmonic orders. From (4) the measure m of the load active power and reactive power, for the harmonic component of order h, greater than the fundamental, are presented (38) and (39), respectively. m m PLh = ϕPh · ˆ Ph + εm Ph

(38)

m m QLh = ϕQh · ˆ Qh + εm Qh

(39)

m and ϕm where the vectors of measures ϕPh Qh are described in (40) and (41), respectively, considering the ZIP model coefficients already determined for the fundamental component; and the model parameter vectors ˆ Ph and ˆ Qh are given by (42) and (43), respectively.

pP = 1 − pZ − pI

(27)

qP = 1 − qZ − qI

(28)

m ϕPh = [ pZ V 2 + pI V + pP h h

Then, the estimates of the active power and reactive power of the Sh term, for a given fundamental voltage value, are described by (29) and (30), respectively.

m ϕQh = [ qZ V 2 + qI V + qP h h

Pˆ 1m

m

= P01 [pZ V 1

2

m

+ pI V 1 + pP ]

m

(29)

m

Qˆ 1m = Q01 [qZ V 1 2 + qI V 1 + qP ]

(30)

where the estimated apparent power is given by (31). Sˆ 1m = Pˆ 1m + jQˆ 1m

(31)

In the second part the parameters of Sh of (4) are obtained considering the estimates of the apparent power with h equal to 1. Thus, from (4) and (31), the measure m of the Sh term for the fundamental component can be written as in (32). m ˆ S1m = ϕ1 ·  1 + εm 1

(32)

m −S ˆ m ; ϕm is the vector of measures described in where S1m = SL1 1 1 ˆ (33);  1 is the vector of estimated parameters of apparent power

represents the modeling, measurement described in (34); and εm 1 and noise error. m ϕ1 = [ |V2m |2

ˆ 1 = [ Y12

Y13

|V3m |2 ···

|VNm |2 ]

··· Y1N ]

T

(33) (34)

Once the vector of the measures is defined for the fundamental component, it is possible to write the observation matrix ␸1 . For M measures, from (32), it is possible to write the overdetermined system according to (35).

⎡

S11

⎢ S2 ⎢ 1 ⎢ ⎢ . ⎣ .. S1M

⎤

⎡

1 ϕ1

⎤

⎡

ε11

⎥ ⎢ ϕ2 ⎥ ⎢ ε2 ⎥ ⎢ 1⎥ ⎢ ⎥ = ⎢ ⎥ · ˆ 1 + ⎢ 1 ⎥ ⎢ . ⎥ ⎢ . ⎦ ⎣ .. ⎦ ⎣ .. M ϕ1

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(35)

m

m

m

|V1m |2 · · · |Vvm |2 · · · |VNm |2 ]

(40)

|V1m |2 · · ·|Vvm |2 · · ·|VNm |2 ]

(41)

ˆ Ph = [ P0h

Gh1

···

Ghv

···

GhN ]

T

ˆ Qh = [ Q0h

Bh1

···

Bhv

···

BhN ]

T

(42) (43)

/ h and 1 < h ≤ N. where v = Once the vector of the measurements is defined, it is possible to generically represent the observation matrices for the load active power and reactive power. For M measurements, from (38) and (39), it is possible to write the overdetermined system according to (44) and (45).

⎡ P1 ⎤ Lh

⎡

1 ϕPh

⎤

⎡ ε1 ⎤ Ph

⎢ P 2 ⎥ ⎢ ϕ2 ⎥ ⎢ ε2 ⎥ ⎢ Lh ⎥ ⎢ Ph ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ · ˆ Ph + ⎢ Ph ⎥ ⎢. ⎥ ⎢. ⎥ ⎢. ⎥ ⎣ .. ⎦ ⎣ .. ⎦ ⎣ .. ⎦ M PLh

⎡ Q1 ⎤ Lh

M ϕPh

⎡ ϕ1 ⎤ Qh

εM Ph

⎡ ε1 ⎤ Qh

⎢ Q 2 ⎥ ⎢ ϕ2 ⎥ ⎢ ε2 ⎥ ⎢ Lh ⎥ ⎢ Qh ⎥ ⎢ Qh ⎥ ⎢ ⎥=⎢ ⎥ · ˆ Qh + ⎢ ⎥ ⎢. ⎥ ⎢. ⎥ ⎢. ⎥ ⎣ .. ⎦ ⎣ .. ⎦ ⎣ .. ⎦ M QLh

M ϕQh

(44)

(45)

εM Qh

The matrix representation of (44) and (45) can be seen in (46) and (47), respectively. P Lh = ϕPh · ˆ Ph + E Ph

(46)

Q Lh = ϕQh · ˆ Qh + E Qh

(47)

The estimates of the parameter vectors can be obtained by the multiple linear regression method as described in (48) and (49). In this way, all model parameters for the component of order h greater than the fundamental can be obtained.

εM 1

The matrix representation of (35) is written in (36). S1 = ␸1 · ˆ 1 + ␧1

m

(36)

The vector with the estimated parameters ˆ 1 can be obtained by multiple linear regression method as described in (37) [36]. In this

−1 T T ˆ Ph = [ϕPh · ϕPh ] · ϕPh · P Lh

(48)

−1 T T ˆ Qh = [ϕQh · ϕQh ] · ϕQh · Q Lh

(49)

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

ZS

373

Table 1 ZIP model coefficients for the case study.

I

Load

Fig. 2. System block diagram for the case study.

Fig. 3. Regressor line and initial reactive power limits.

5. Case study and results In order to apply the proposed model, a circuit was built to carry out the tests and thus obtain the measured values necessary for the estimation and validation processes of the model parameters. A distribution system was represented by a Thévenin equivalent circuit, composed of a voltage source VS and an system equivalent impedance ZS , as shown in Fig. 2 [19], where V and I are the load voltage and current, respectively. Elements such as transformers and nonlinear loads connected to the system distort the voltage applied to the load to be modeled. The circuit of Fig. 2 seeks to represent such influence through the VS source. In addition, the equivalent impedance ZS considers the harmonic voltage drop in the network caused by the load harmonic current. In this case study the ZS impedance is composed of an inductance in each phase, with nominal value equal to 0.5 mH. To represent the presence of nonlinear loads in the system, the voltage VS contains several harmonic components, and the nominal effective value of its fundamental component is equal to 100 V with frequency of 60 Hz for the three phases. The three-phase power supply used was the Pacific Power Source model 360AMX. The implemented nonlinear load is a three-phase full-bridge diode rectifier, supplying a RL load, which is, constant impedance, whose nominal resistance and inductance values are equal to 45  and 1 mH, respectively. In this case study, the load representation was performed considering up to the 25th harmonic order, which is, N equal to 25. With the power quality meter Voltech model PM3000A, the magnitude and angle values of the voltage and current harmonic components were obtained up to the 25th order in one input phase of the three-

pZ

pI

pP

qZ

qI

qP

0.9915

0.0085

0.0000

1.0000

0.0000

0.0000

phase rectifier. The circuit was considered balanced and the model to be obtained will be applied to each system phase. The disturbances required to identify the load model parameters were obtained by random amplitude variations of the Thévenin equivalent circuit voltage source VS . Such source voltage variations occur in both the fundamental component and the higher order harmonic components, each within a range that assumes different percentage values for each harmonic order. The variations of the fundamental voltage component were approximately 3%. These disturbances aim at the circuit operation at different points with linear independence between the collected samples and their harmonic components, allowing the model parameters identification. Random variations in the source voltage amplitude were realized in the fundamental component and in the components of order 5, 7, 11, 13, 17, 19, 23 and 25. In this way, the load voltage V will be influenced by the source voltage VS and by the voltage drop on the ZS impedance due to the load current I which contains several harmonic components. By previously evaluating the values of the load measured harmonic voltages it was found that the even-order and multiples of third order components were not relevant and, therefore, could be excluded from the modeling. Thus, for the value of N, defined previously equal to 25 from (14), it was verified that the quantity of measures M must be greater than 9. Considering the presence of noises in the measurements, we adopted the quantity of measures M equal to 18 for the overdetermined system construction of the model parameters estimation process. In the tests performed, 26 load voltage and current measurements were collected. Approximately 70% of the collected data, 18 samples, were used in the parameter estimation process and the other data, 8 samples, were used in the model validation process. 5.1. Model parameters According to the procedure already presented, in the first part of obtaining the parameters for the fundamental component, it is necessary to establish the minimum and maximum limits for the active power and reactive power to be used in the exhaustive search. Exemplifying the way of obtaining these limits, Fig. 3 shows the min and Q max as well as the regressor line reactive power limits Q01 01 ˆ L1 used to determine these limits from the collected samples. Q Fig. 3 also shows the points of the fundamental component of the load reactive power, obtained from the data used in the parameter estimation process. The step size used in the parameters search pZ , pI , qZ and qI was equal to 5 × 10−4 and for P01 and Q01 their respective search spaces were divided into 1000 parts, which is, the step size for P01 max − P min )/1000 and for Q max was equal to (P01 01 was equal to (Q01 − 01 min Q01 )/1000. At the end of the procedure, we obtained the ZIP model coefficients as well as the initial values of the active power and reactive power for the fundamental component and higher order components, which are described in Tables 1 and 2, respectively. The powers of the other harmonic orders are considered equal to zero. In addition, it was determined from the cross-admittance matrix, whose magnitudes are shown in Fig. 4, where it is possible to verify that the order v harmonic voltages contribute in a different way to the composition of the load power of order h, as described by (7).

374

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

Table 2 Initial active power and reactive power for the case study.

6

h

P0h (W)

Q0h (VAR)

1 5 7 11 13 17 19 23 25

404.622 1.490 −0.312 −0.207 −0.045 −0.049 −0.009 0.006 0.008

25.753 0.252 −0.232 −0.176 −0.067 −0.111 −0.034 −0.039 −0.037

4 2 0 -2 -4 -6

Measured Estimated 0

5

10

15

20

Fig. 5. Current waveforms for the case study.

Current Magnitude (A)

4 Measured Estimat ed

3 2 1 0

1

3

5

7

9 11 13 15 17 19 21 23 25 Harmonic Compo nent Order

200 100 0 -100

Fig. 4. Admittance magnitudes obtained for the case study.

-200

5.2. Results The present case study aims to apply the new proposed model, as well as the methodology to obtain its parameters, to represent the electronic load constructed: in this case, a three-phase, full-bridge diode rectifier, feeding an RL load. According to the methodology presented previously, the values of the following model parameters were obtained for the studied load: initial powers (P0h and Q0h ) and ZIP coefficients (pZ , pI , pP , qZ , qI and qP ), as well as the elements of the admittance matrix (Yhv ).Thus, considering the parameters already determined, the load representation can be performed using Eq. (4)–(8). From (4) it is possible to determine the estimate for the harmonic components of the complex apparent power of the modeled load SLh from the harmonic components of the voltage applied on this load Vh for each harmonic order h. With the estimated values of SLh and the voltages Vh , it is also possible to determine the estimates of the harmonic current components of the load by Ih = (SLh /Vh )* for each harmonic order h. Thus, the load current as a function of time t, considering only N harmonic components, is given by (50). i(t) =

√

N 

2

|Ih | sin (hωt + h )

(50)

h=1

where |Ih | and h are the rms value and angle of the component Ih of the load current, respectively; and ω = 2 f. For this case study, f is equal to 60 Hz and N is equal to 25. One way of verifying the proximity between an estimated and measured value of the load current is by comparing the graphs of the waveforms of these currents. Thus, the waveforms of approximately one period of a measured and estimated load current sample, reconstituted from the harmonic components up to the

Measured Estimated 1

3

5

7

9

11 13 15 17 19 21 23 25

Fig. 6. Measured and estimated current: (a) magnitude; (b) phase.

25th order as described by (50), can be seen in Fig. 5. As the worst case, the sample used refers to the current that provided the largest Total Harmonic Distortion (THD) difference between the measured values of the validation data and their estimated values [4]. The measured and estimated magnitudes and angles of the same load current sample are shown in Fig. 6, where graph (a) allows comparison between current magnitudes and graph (b) allows comparison between current angles. It can be seen in Figs. 5 and 6 that the estimated values are very close to the measured values. One indication of this proximity is that the greatest THD difference found between the measured and estimated currents of the validation samples was only 2.3 × 10−3 . Another way of verifying the proximity between the measured and estimated values of the validation data is by determining the largest relative rms current error between the measured and estimated values, which in this case was only 0.39%. The errors between the active and reactive power measured and estimated can be visualized in Fig. 7. The load active power and reactive power are determined by the measured currents, and voltages of the validation data and the respective estimated active power and reactive power are obtained by the model.

N Once defining S0 = P0 + jQ0 = P + jQ0h it is possible to h=1 0h determine |S0 | from the data in Table 2. Thus, considering the errors of the active power and reactive power of the validation data, shown in Fig. 7, it is possible to verify that the largest magnitudes of these errors in relation to |S0 | are equal to 0.21% and 0.08%, respectively. This shows the small contribution of the active power

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

0.5

0

-0.5

-1

1

2

3

4

5

6

7

8

Sample 0.4 0.2

375

constant current and constant power. It was possible to verify that for the tested rectifier some parameters can be suppressed, due to their small contribution to the harmonic power of the load. Such an evaluation allows the model to be reduced without significantly compromising the load representation. The results of the case study showed that for the modeled load, several harmonic voltage components contributed to define the value of a given harmonic power and current component of the load. In addition, the various load harmonic power components were affected differently by voltage variations, indicating the influence of the model on the harmonic flow and consequently the harmonic injection in the network. This confirms the applicability of the model for harmonic studies in systems subject to variable harmonic voltages.

0

References

-0.2 -0.4 1

2

3

4

5

6

7

8

Fig. 7. Error between measured and estimated power: (a) active; (b) reactive.

and reactive power errors in the composition of the apparent load power. When the power values P0h and Q0h provided by the model are very low in relation to |S0 |, it means that the voltage component of order h does not contribute significantly to the respective load power in the same order. In this case, it was verified that the absolute values of P0h /|S0 | and Q0h /|S0 | for h equals 19, 23 and 25 are less than 0.01% so that the terms of the ZIP portion of the model of these harmonic orders do not contribute significantly to the composition of the apparent power of the tested load and can be neglected in some analyses of this load. Another important result is that several admittance magnitudes shown in Fig. 4 are nonzero, indicating that only the ZIP model would not satisfactorily represent the tested nonlinear load. For a linear load the values of the cross admittance matrix would be zero. Furthermore, for the tested rectifier, it was possible to determine that its load is of the constant impedance type, based on the estimated ZIP model coefficient values, described in Table 1. This analysis can help characterize the load or set of loads, showing its contribution to the distortion of the voltage and current in the system. 6. Conclusions A proposed load model that combines the traditional ZIP model with a cross-admittance matrix was presented so as to allow the modeling of nonlinear loads that underwent variable voltages with harmonic components. The model is applicable for the modeling of one or more harmonic loads in systems with the presence of other nonlinear loads, allowing us to obtain some physical knowledge about the represented load. The procedure for using the methodology was described, as well as the mathematical development of the model. A case study was presented to show the application of the methodology proposed in this article. The proposed new model was able to model with high accuracy a three-phase full-bridge diode rectifier with RL load, characterizing its behavior even when it undergoes variable harmonic voltages. The proposed methodology allowed the modeling of the tested load even with the presence of measurement noise in the data collected by the power quality meter. For the load modeled in the case study the ZIP coefficient limitation did not affect the model accuracy and allowed to determine the proportion of each portion in relation to constant impedance,

[1] IEEE Recommended Practice and Requirements for Harmonic Control in Electric Power Systems, IEEE Std 519-2014 (Revision of IEEE Std 519-1992), 2014, 1–29. [2] V. Terzija, G. Valverde, C. Deyu, P. Regulski, V. Madani, J. Fitch, S. Skok, M.M. Begovic, A. Phadke, Wide-area monitoring, protection, and control of future electric power networks, Proc. IEEE 99 (2011) 80–93. [3] L.T.M. Mota, A.A. Mota, Load modeling at electric power distribution substations using dynamic load parameters estimation, Int. J. Electr. Power Energy Syst. 26 (2004) 805–811. [4] R. Arghandeh, A. Onen, J. Jung, D. Cheng, R.P. Broadwater, V. Centeno, Phasor-based assessment for harmonic sources in distribution networks, Electr. Power Syst. Res. 116 (2014) 94–105. [5] A. Samui, S.R. Samantaray, An active islanding detection scheme for inverter-based DG with frequency dependent ZIP–exponential static load model, Int. J. Electr. Power Energy Syst. 78 (2016) 41–50. [6] N. Sasidharan, N. Madhu M, J.G. Singh, W. Ongsakul, An approach for an efficient hybrid AC/DC solar powered Homegrid system based on the load characteristics of home appliances, Energy Build. 108 (2015) 23–35. ´ M. Todorovski, Optimal capacitor placement in distorted [7] J. Vuletic, distribution networks with different load models using Penalty Free Genetic Algorithm, Int. J. Electr. Power Energy Syst. 78 (2016) 174–182. [8] I.D. Melo, J.L.R. Pereira, A.M. Variz, P.A.N. Garcia, Harmonic state estimation for distribution networks using phasor measurement units, Electr. Power Syst. Res. 147 (2017) 133–144. [9] L.S. Xavier, A.F. Cupertino, J.T. de Resende, V.F. Mendes, H.A. Pereira, Adaptive current control strategy for harmonic compensation in single-phase solar inverters, Electr. Power Syst. Res. 142 (2017) 84–95. [10] M. Lamich, J. Balcells, M. Corbalán, E. Griful, Nonlinear loads model for harmonics flow prediction, using multivariate regression, IEEE Trans. Ind. Electr. 64 (2017) 4820–4827. [11] D.D. Ferreira, E.A. Nagata, S.C. Ferreira, J.M. de Seixas, C.A. Duque, C.A.G. Marques, J.D.S. Guedes, A.S. Cerqueira, Method based on independent component analysis for harmonic extraction from power system signals, Electr. Power Syst. Res. 119 (2015) 19–24. [12] Measurement-Based Load Modeling, EPRI Technical Report 1014402, EPRI, Palo Alto, CA, 2006. [13] E. Feinberg, J. Hu, E. Yuan, A stochastic search algorithm for voltage and reactive power control with switching costs and ZIP load model, Electr. Power Syst. Res. 133 (2016) 328–337. [14] A. Bokhari, A. Alkan, R. Dogan, M. Diaz-Aguiló, F. d. León, D. Czarkowski, Z. Zabar, L. Birenbaum, A. Noel, R.E. Uosef, Experimental determination of the ZIP coefficients for modern residential, commercial, and industrial loads, IEEE Trans. Power Deliv. 29 (2014) 1372–1381. [15] A. Kumar, J. Kumar, ATC with ZIP load model—a comprehensive evaluation with third generation FACTS in restructured electricity markets, Int. J. Electr. Power Energy Syst. 54 (2014) 546–558. ˜ An impedance relation [16] A. Herrera-Orozco, J. Mora-Flórez, S. Pérez-Londono, index to predict the fault locator performance considering different load models, Electr. Power Syst. Res. 107 (2014) 199–205. [17] V.V.S.N. Murty, A. Kumar, Mesh distribution system analysis in presence of distributed generation with time varying load model, Int. J. Electr. Power Energy Syst. 62 (2014) 836–854. [18] M. Manbachi, H. Farhangi, A. Palizban, S. Arzanpour, Quasi real-time ZIP load modeling for conservation voltage reduction of smart distribution networks using disaggregated AMI data, Sustain. Cities Soc. 19 (2015) 1–10. [19] R. Senra, W.C. Boaventura, E.M.A.M. Mendes, Assessment of the harmonic currents generated by single-phase nonlinear loads, Electr. Power Syst. Res. 147 (2017) 272–279. [20] C.A. Canesin, L.C.O. de Oliveira, J.B. Souza, D.D.O. de Lima, R.P. Buratti, A time-domain harmonic power-flow analysis in electrical energy distribution networks, using Norton models for non-linear loading, IEEE 16th International Conference on Harmonics and Quality of Power (ICHQP) (2014) 778–782.

376

M. Brunoro et al. / Electric Power Systems Research 152 (2017) 367–376

[21] A. Javadi, A. Hamadi, L. Woodward, K. Al-Haddad, Experimental investigation on a hybrid series active power compensator to improve power quality of typical households, IEEE Trans. Ind. Electr. 63 (2016) 4849–4859. [22] F. Karimzadeh, S. Esmaeili, S.H. Hosseinian, A novel method for noninvasive estimation of utility harmonic impedance based on complex independent component analysis, IEEE Trans. Power Deliv. 30 (2015) 1843–1852. [23] A.S. Folting, J.M.A. Myrzik, T. Wiesner, L. Jendernalik, Practical implementation of the coupled norton approach for nonlinear harmonic models, Power Systems Computation Conference (PSCC) (2014) 1–7. [24] P. Regulski, D.S. Vilchis-Rodriguez, S. Djurovic, V. Terzija, Estimation of composite load model parameters using an improved particle swarm optimization method, IEEE Trans. Power Deliv. 30 (2015) 553–560. [25] J.R. Marti, H. Ahmadi, L. Bashualdo, Linear power-flow formulation based on a voltage-dependent load model, IEEE Trans. Power Deliv. 28 (2013) 1682–1690. [26] J.V. Milanovic, On unreliability of exponential load models, Electr. Power Syst. Res. 49 (1999) 1–9. [27] M. Dös¸o˘glu, A.B. Arsoy, Modeling and simulation of static loads for wind power applications, Neural Comput. Appl. 25 (2014) 997–1006. [28] M. Fauri, Harmonic modelling of non-linear load by means of crossed frequency admittance matrix, IEEE Trans. Power Syst. 12 (1997) 1632–1638.

[29] Y. Sun, G. Zhang, W. Xu, J.G. Mayordomo, A harmonically coupled admittance matrix model for AC/DC converters, IEEE Trans. Power Syst. 22 (2007) 1574–1582. [30] R. Bagheri, J. Moghani, G. Gharehpetian, S. Mirtalaei, Determination of aggregated load power consumption, under non-sinusoidal supply using an improved load model, Energy Convers. Manage. 50 (2009) 1563–1569. [31] A.S. Koch, J.M.A. Myrzik, T. Wiesner, L. Jendernalik, Evaluation and validation of Norton approaches for nonlinear harmonic models, PowerTech (POWERTECH) 2013 IEEE Grenoble (2013) 1–6. [32] P.F. Ribeiro, Time- Varying Waveform Distortions in Power Systems, John Wiley & Sons, 2009. [33] P. Ju, C. Qin, F. Wu, H. Xie, Y. Ning, Load modeling for wide area power system, Int. J. Electr. Power Energy Syst. 33 (2011) 909–917. [34] B.-H. Kim, H. Kim, Measurement-based estimation of the composite load model parameters, J. Electr. Eng. Technol. 7 (2012) 845–851. [35] L.A. Aguirre, Introduc¸ão à Identificac¸ão de Sistemas—Técnicas Lineares e Não-Lineares Aplicadas a Sistemas Reais, 3rd ed., UFMG, Belo Horizonte, 2007. [36] G.R. Rey, L.M. Muneta, Power Quality Harmonics Analysis and Real Measurements Data, InTech, 2011.