Curve-fitting-based method for modeling voltage–current characteristic of an ac electric arc furnace

Electric Power Systems Research 80 (2010) 572–581

Contents lists available at ScienceDirect

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

Curve-fitting-based method for modeling voltage–current characteristic of an ac electric arc furnace Yu-Jen Liu, Gary W. Chang ∗ , Rong-Chin Hong Department of Electrical Engineering, National Chung Cheng University, 168 University Rd., Min-Hsiung, Chia-Yi 621, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 23 April 2009 Received in revised form 22 September 2009 Accepted 24 October 2009 Available online 24 November 2009 Keywords: Cubic spline interpolation Electric arc furnace Harmonics Voltage–current characteristic

a b s t r a c t Field measurements of voltage and current are the most effective way for characterizing the electric response of an ac electric arc furnace that describes its nonlinear behavior. Sufficient measured information can be adopted to determine the background harmonic level in a power system, to characterize specific sources of harmonics, and to develop an appropriate nonlinear voltage–current characteristic. In this paper, a curve-fitting-based method called cubic spline interpolation is proposed to model the voltage–current characteristic of an ac electric arc furnace in the steady state. Meanwhile, the actual measured data are collected for modeling use. Two classic methods, harmonic current injections and equivalent harmonic voltage sources, for modeling the electric arc furnace load are reviewed and used to evaluate the performance of proposed model. Results obtained from the measured data and computer simulations of the three electric arc furnace models are then compared according to the voltage and current waveforms, as well as the voltage–current characteristic. It is shown that the proposed model is more accurate than the two classic approaches for harmonic assessment of electric arc furnaces and can be used for modeling similar types of nonlinear loads in the harmonic penetration study. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Electric arc furnaces (EAFs) have been widely used in steelmaking industries. The EAF employs high temperatures produced by the low-voltage and high-current electric arc existing between the electrodes and melting material to smelt scrap iron raw materials. Nowadays, EAFs are designed for very large power input ratings. Due to the nature of the electric arc and the meltdown processes, these devices cause significant waveform distortions such as harmonics, interharmonics and flickers in the supply network. Therefore, electric utilities and their customers always pay much attention to mitigate power quality (PQ) problems associated with EAF loads. For an effective mitigation of PQ problems occurred, obtaining the time response of the EAF becomes of great importance to investigate the impact of the nonlinear and time-varying load to the connected power system. Meanwhile, an accurate EAF model is necessary for assessment of the usefulness of different solutions for PQ problems associated with EAFs. Usually, the steel-making from EAFs is a batch type of operations, which generally needs to add raw materials two to three times during the complete process which takes about 2–3 h. Fig. 1 shows a measured power variation of the EAF during a typical pro-

∗ Corresponding author. Tel.: +886 5 2428302; fax: +886 5 2720862. E-mail address: [email protected] (G.W. Chang). 0378-7796/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2009.10.015

cess. There are three distinct periods included in the steel-making process: striking, melting, and refining [1]. In the striking the threephase electrodes of the furnace are lowered and contact with the scrap steel, which lead to the electric arc build-up. The melting is then started. Finally, the whole process ends in stably refining. However, considering the complexity of the EAF operation and the randomness associated with each operation period, it is very difficult to develop an accurate deterministic model for describing the dynamic behavior of an EAF during the striking and melting periods. Figs. 2 and 3 show example waveforms of the measured 20 cycle arc voltages and currents during striking and melting, respectively, where a nonstationary stochastic phenomenon can be observed. Traditionally, most research efforts of EAF models focus on steady-state modeling used to describe the operation of EAFs during refining. Factors that affect EAF operations usually include the melting or refining materials, the electrode position, the electrode arm control scheme, and the supply system voltage and impedance. Modeling the EAF load depends on several parameters (e.g. arc voltage, arc current, and arc length) which are determined by the positions of electrodes [2]. Thus, it is of importance to know the electric responses of EAFs to develop an appropriate EAF model for PQ study. In reality, some physical responses of EAFs are not easy to derive from the theoretical study and assumptions are often made to simplify the modeling task. The most effective and simplest manner for obtaining the electric responses of EAF loads is the field

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

573

Fig. 1. Measured power variation of the ac EAF.

Fig. 3. (a) Arc voltage and (b) arc current during melting stage.

voltage is defined as a nonlinear function of the arc length. Classic models based on nonlinear differential equations are adopted to model the EAF in [10,11], where Mayr’s and Cassie’s equations are solved to obtain the nonlinear conductance of the EAF. A harmonicdomain solution method is used in [12,13], where the EAF model is developed from the energy balance equation and is a nonlinear differential equation of the arc radius and the arc current. However, some physical parameters are not easy to obtain and to be used for the aforementioned EAF models. In practice, the field measurement of EAF voltages and currents is the most effective way to obtain the required parameters for modeling the EAF. In this paper, a nonlinear conductance EAF model based on a curve-fitting technique is proposed. To show the accuracy of the proposed model, two classic EAF models based on actual measured data are also introduced for comparison. The following briefly reviews the two classic methods: harmonic current injection and harmonic voltage source models. Fig. 2. (a) Arc voltage and (b) arc current during striking stage.

2.1. Harmonic current injection model (HCIM) measurements of EAF voltages and currents. Based on the measured EAF voltages and currents, this paper proposes a dynamic EAF conductance model which is developed by the cubic spline interpolation method (CSIM). To evaluate the performance of the proposed model, two classic EAF models: harmonic current injection model (HCIM) and harmonic voltage source model (HVSM), are also reviewed and compared. An actual power system of a steel plant for field measurements and simulations is then illustrated to test the usefulness of the three EAF models. 2. Review of classic models for voltage–current characteristic of the ac electric arc furnace Voltage or current distortions associated with EAFs result from nonlinear voltage–current (i.e. vi) characteristics of such loads. Many steady-state EAF models have been developed for PQ studies. The simplest one is to model the EAF by an inductor in series with a resistor [3]. A nonlinear resistance model based on using the piecewise linearization method to describe the furnace with the assumed vi characteristic has been proposed in [4,5]. In [6–9], the EAF is modeled as a time-varying voltage source and this arc

The HCIM is a commonly used method for modeling nonlinear loads in the harmonic study. In this model, harmonic components of the measured arc current of the EAF load are obtained by performing the fast Fourier transform (FFT) and then are injected into the power network through the secondary side of the furnace transformer. The harmonic current source for describing EAF load is represented by its Fourier series as: iarc (t) =

H  h=1

ah sin(hωt) +

H 

bh cos(hωt)

(1)

h=1

where H is the highest order of harmonics under considerations and the Fourier coefficients may change randomly during each fundamental cycle and the coefficients are selected as a function of the measured arc current. Fig. 4 shows the equivalent circuit of the HCIM used in the study. Once the currents of the EAF load are measured, FFT is performed to analyze the arc current and then inject all considered harmonic components into the EAF transformer secondary side. The EAF harmonic voltages thus can be determined.

574

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

Fig. 4. Harmonic current injection model of an EAF.

Fig. 5. Harmonic voltage source model of an EAF.

2.2. Harmonic voltage source model (HVSM) In the HVSM, as shown in Fig. 5, the EAF is modeled as a steadystate harmonic voltage source behind series impedance, Zcable , which is composed of the EAF transformer secondary cable to the electrode. The measured arc voltage, vmeasured (t), is taken at the EAF transformer secondary side and the voltage drop across the cable, vcable (t), is calculated by the arc current and cable impedance. The simulated arc voltage, varc (t), is derived from the difference between vmeasured (t) and vcable (t). The mathematic expressions of these voltages are given in (2)–(4).

 H

vmeasured (t) =

Vh sin(hω0 t + h ),

(2)

h=1

vcable (t) = Lc

di(t) + Rc i(t), dt

varc (t) = vmeasured (t) − vcable (t).

(3) (4)

After (4) is calculated, the FFT is applied to obtain harmonic voltages used to represent the EAF. The EAF harmonic currents are then can be determined by the voltage source model. 3. Proposed model for voltage–current characteristic of the electric arc furnace Though both of the HCIM and HVSM provide simple manners to model the EAF load, however, the drawback is that the EAF model is less accurate when the harmonic orders under considerations are insufficient. In this paper, the cubic spline interpolation model (CSIM) based on a curve-fitting technique is proposed to model the EAF as a function of nonlinear conductance [14]. The development of the proposed EAF model is considered as a problem of function approximations. A set of mathematical polynomial is obtained for modeling that allows one understanding the derivation of voltage–current characteristic of the EAF load. The CSIM is based on building up the cubic polynomials of the nonlinear arc conductance to describe the behavior of the EAF. A cubic spline is a spline constructed of piecewise third-order polynomials which pass through a set of data points without knowing slopes. The following presents the procedure of modeling EAF by using the CSIM.

Fig. 6. (a) Measured conductance of the EAF with its interpolated curve for one fundamental cycle and (b) the first five points of (a).

Consider a measured waveform shown in Fig. 6 of one fundamental cycle with 64 sampling points (i.e. a sampling frequency of 3840 Hz) for the EAF conductance. Give the measured waveform defined on [a,b] and a set of data points a = x0 < x1 < · · · < xn = b, where n is the number of sampling points minus 1 (here n = 63). It is seen that there is an interval between two adjacent measured points, however, the values of the EAF conductance (i.e. the slope) within each interval is unknown. Thus, an unknown conductance function, Gi (x), is defined for each interval. The CSIM is then used to derive from a three-order polynomial, as shown in (5), which can be used to express the unknown conductance function, Gi (x). Gi (x) = ai + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3

(5)

In (5), i = 0, 1, 2, . . . , n − 1 and the coefficients ai , bi , ci , and di are unknown. The next task is thus to determine these coefficients so that the value of Gi (x) can be exactly estimated. For this purpose, a cubic spline interpolated Gi (x) is a function that must satisfies the following requirements: • Each spline must pass through the given data points yi . Gi (xi ) = yi

and

Gi (xi+1 ) = yi+1 ,

i = 0, 1, 2, . . . , n − 1

(6)

• Interior data points between each spline must be continuous. Gi+1 (xi+1 ) = Gi (xi+1 ),

i = 0, 1, 2, ..., n − 1

(7)

• The first and the second derivatives of the splines must be continuous across the interior data points. Therefore, the spline forms a smooth function.  Gi+1 (xi+1 ) = Gi (xi+1 ),

i = 0, 1, 2, ..., n − 2

(8)

 Gi+1 (xi+1 )

i = 0, 1, 2, . . . , n − 2

(9)

=

Gi (xi+1 ),

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

Furthermore, an additional boundary condition also must be considered in the CSIM. Depending on different types of boundary conditions imposed at x0 = a and xn = b, two widely used boundary conditions, natural and clamped boundaries, are defined as follows: Gi (x0 ) = Gi (xn ) = 0 Gi (x0 ) = D0

(natural)

and Gi (xn ) = Dn

(10) (clamped)

(11)

where D0 and Dn are the values of the first derivatives of the unknown functions. The natural boundary condition generally gives less accurate results than the clamped condition near the endpoints of the interval [x0 , xn ], unless the function Gi nearly satisfies Gi (x0 ) = Gi (xn ) = 0 [14]. Though the clamped boundary leads to more accurate approximation, however, it is necessary to have either the values of the derivatives at the endpoints satisfying the boundary condition. An alternative to the natural boundary condition that does not require knowledge of the derivative of the function is the nota-knot condition, which is also adopted in the study. Such condition requires that G (x) must be continuous at x1 and at xn−1 . By applying aforementioned requirements and boundary conditions in the definition to the cubic polynomials of (5), the calculation of the coefficients of the cubic polynomials representing each interval of Fig. 6 is simplified by a suitable choice of these algebraic representations of the equations. Let yi be the value of each measured data point for i = 0, 1, 2, . . . , n. A set of mathematical expressions for using the CSIM with the not-a-knot condition to estimate the function Gi (x)of (5) is given by (i) Set ai = yi for i = 0,1,2„. . .,n. (ii) Solve: (h0 + 2h1 )c1 + (h1 − h0 )c2 =





3 h1 (a2 − a1 ) − (a1 − a0 ) , h1 + h0 h0

for i = 1

(12)

and hi−1 ci−1 + 2ci (hi−1 + h1 ) + hi ci+1 =

3 3 (a − ai ) − (a − ai−1 ), hi i+1 hi−1 i

for i = 2, . . . , n − 2 (13)

and (hn−2 − hn−1 )cn−2 + (2hn−2 + hn−1 )cn−1 =

3 hn−1 + hn−2

h

n−2

hn−1

where hi = xi+1 − xi , (iii) Set c0 = c1 +

cn = cn−1 +

hn−1 (cn−1 − cn−2 ). hn−2

bi =

1 h (a − ai ) − i (ci+1 + 2ci ), 3 hi i+1

for i = 0, 1, 2, . . . , n − 1. (17)

(v) Set di =

1 (c − ci ), 3hi i+1

for i = 0, 1, 2, . . . , n − 1.

(18)

Eqs. (12)–(14) of step (ii) form a linear system with n − 1 equations for the n − 1 unknown coefficients c1 , c2 ,. . ., and cn−1 . In a matrix form of Ac = q, the following matrices are obtained:

⎡

(h0 + 2h1 ) h2 ⎢ A=⎣ . . . 0

⎡

⎡

⎤

c1 c2 .. .

⎢ c=⎢ ⎣

⎥ ⎥ ⎦

cn−1

⎢ ⎢ ⎢ ⎢ q=⎢ ⎢ ⎢ ⎢ ⎣

(h1 − h0 ) 2(h1 + h2 ) . . . 0

··· ··· ··· ···

(hn−2

0 0 . . . − hn−1 )

⎤

0 0 ⎥ , . ⎦ . . (2hn−2 + hn−1 ) n×n

(19)

,

n×1

(20)





h1 3 (a2 − a1 ) − (a1 − a0 ) h1 + h0 h0





⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

3 3 (ai+1 − ai ) − (a − ai−1 ) hi hi−1 i .. .  3 hn−2 (an − an−1 ) − (an−1 − an−2 ) hn−1 + hn−2 hn−1

.

(21)

n×1

The coefficients ai can be found in step (i) and ci are given in (15), (16) and (20). Then the coefficients bi and di can be derived from (17) and (18), respectively. Therefore, a completive spline function for describing the nonlinear conductance of Fig. 6(a) for the EAF can be expressed by

G(t) =

⎧ G0 (t) ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩

G1 (t) .. . Gn−1 (t)

t0 ≤ t ≤ t1 t1 ≤ t ≤ t2 , .. . tn−1 ≤ t ≤ tn

(22)

(14)

4. Case study

(15)

By comparing the performances of the three described methods for modeling the EAF, Fig. 7 shows an one-line diagram of the actual power system of a steel plant used for illustration. In the system, the plant is fed directly from the utility power network at 161 kV, a delta-wye main transformer is connected to transform 161–11.4 kV and a 50-ton three-phase ac EAF connected to an 11.4 kV/460 V transformer rated at 33 MVA is investigated. In addition, to observe

for i = 1, 2, . . . , n − 1.

h0 (c1 − c2 ), h1

(iv) Set

where Gi (t) is a piecewise cubic polynomial and i = 1,2,. . .,n−1.



(an − an−1 ) − (an−1 − an−2 ) ,

575

(16)

Fig. 7. One-line diagram of the power system under study.

576

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

Fig. 8. MATLAB/SIMULINK-based model of the system in Fig. 7.

Table 1 Measured current and voltage harmonic components used for HCIM and HVSM at the later stage during refining period. HCIM

Fig. 9. Equivalent HCIM for EAF load in Fig. 8.

the mitigation of voltage/current distortions caused by the EAF, devices such as the static VAR compensator (SVC) are also modeled. A PQ meter is located at the EAF transformer secondary side for measurements. Also, the power system under study modeled by MATLAB/SIMULINK is shown in Fig. 8, where the functional block of the EAF load can be replaced by HCIM, HVSM, and the proposed model, respectively.

HVSM

Harmonics

Current (%)

Harmonics

Voltage (%)

2 3 4 5 6 7 8 9 10 11 12 13 14 15

2.5 6.5 0.8 8 1 3.5 0.5 2 0.3 1.75 0.6 0.7 0.5 0.6

2 3 4 5 6 7 8 9 10 11 12 13 14 15

0.5 18 0.5 2.5 1 2 0.5 1.5 0.5 1.2 0.6 1.6 0.4 0.8

It is known that the operation of EAF is always a nonstationary stochastic behavior. Even though the EAF is operated during the refining period, it is still not easy to find a long-lasting steady-state response. Therefore, in this study three different vi characteristics at the earlier, middle, and later stages of the refining period will be tested to demonstrate the usefulness of the proposed model

Fig. 10. Equivalent HVSM for EAF load in Fig. 8.

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

by measuring the instantaneous voltage waveforms, current waveforms, and vi characteristics at phase A with a fixed sampling rate of 3840 Hz for one fundamental cycle. It must be noted that the measured voltage and current data both at the earlier and middle stages during the refining period still present dynamic behavior. Thus the HCIM and HVSM cannot be adopted to model such vi characteristic since these two methods are only applicable for modeling nearly steady-state behavior. On the contrary, the proposed method is cable of modeling different vi characteristic at different stages during the refining period. 4.1. Simulations for HCIM and HVSM The equivalent EAF models developed by HCIM and HVSM for the EAF load block of Fig. 8 are shown in Figs. 9 and 10, respec-

577

tively. Based on the measured currents and voltages obtained at the later stage during the refining period, modeling EAF by the HCIM and HVSM methods which consider harmonics up to the 15th order are sufficient for simulations. Table 1 lists the contents of measured harmonic voltage and current used for representing the two classic models. Figs. 11 and 12 display the modeling results of voltage waveforms, current waveforms, and vi characteristics under simulations for the HCIM and HVSM. It is seen that the classic methods cannot precisely track the dynamic characteristics of the actual measured EAF current and voltage due to the lack of complete measured voltages and currents to represent the vi characteristic and ignoring some high-order harmonics in both models. Therefore, either the HCIM or HVSM cannot produce an accurate function approximation for the EAF vi characteristic. Both classic models lead to larger modeling errors.

Fig. 11. HCIM modeling results obtained at the later stage during refining period: (a) voltage waveform, (b) voltage spectra, (c) current waveform, (d) current spectra, and (e) vi characteristic.

578

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

Fig. 12. HVSM modeling results obtained at the later stage during refining period: (a) voltage waveform, (b) voltage spectra, (c) current waveform, (d) current spectra, and (e) vi characteristic.

4.2. Simulations for the proposed model In the proposed model, as shown in Fig. 13, the CSIM is adopted to fit the measured EAF conductance and then an approximated EAF conductance can be produced in the block of conductance model and is connected to the system under study. Figs. 14 and 15 show

test results of the EAF operated at the earlier and middle stages during the refining period, respectively. It is found that, as the EAF enters the refining period, the EAF vi characteristics are still highly nonlinear. Voltages and currents per cycle in these stages are not identical and the vi characteristic presents as a multi-valued function. Therefore, it is difficult to directly model EAF based on its vi

Fig. 13. Proposed model for EAF load in Fig. 8.

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

Fig. 14. CSIM modeling results obtained at the earlier stage during refining period: (a) voltage waveform, (b) current waveform, and (c) vi characteristic.

Fig. 15. CSIM modeling results obtained at the middle stage during refining period: (a) voltage waveform, (b) current waveform, and (c) vi characteristic.

579

580

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

Fig. 16. CSIM modeling results obtained at the later stage during refining period: (a) voltage waveform, (b) voltage spectra, (c) current waveform, (d) current spectra, and (e) vi characteristic.

characteristic since the instantaneous electric response cannot be easily determined by the multi-valued relationship. The proposed method thus models the EAF load by its dynamic conductance. Once the EAF voltage and current are obtained from the conductance model, the vi characteristic can be determined. To observe the results obtained at the later stage during the refining period, the sampling frequency of the PQ meter used for field measurements is set to be 3840 Hz. It means that a one fundamental cycle with 64 sampling points of the measured EAF conductance is used for modeling input. The harmonic components to be observed are up to the 32nd order with the given sampling frequency. However, some measured high-frequency harmonic components of the EAF arc currents and voltages are relatively small in comparison to their fundamental components and they are not shown in the harmonic spectra. Therefore, the harmonic spectra of currents and voltages to be observed for the EAF are also up to the 15th order. The accuracy of the results derived from the proposed model depends only on a numerical tolerance set

for the model. By observing Fig. 16, it is found that the proposed model presents a relatively accurate fitting than those of the classic ones. 4.3. Solution accuracy Table 2 indicates the average relative error of each waveform in Figs. 11, 12, and 14–16, where the actual measurement is adopted as the reference. It is noted that both the EAF voltage and current determined by the HCIM yield the largest errors among the three compared methods. Despite the use of HVSM can effectively reduce errors, results are still less accurate than those obtained by the proposed model. Literature surveys show that the harmonic current injection and equivalent harmonic voltage source methods are still commonly used in harmonic penetration studies. Thus, it is concluded that equivalent harmonic voltage source method is suggested for EAF modeling when only considering the HCIM and HVSM. Overall, by observing Table 2, it is seen that a better agree-

Y.-J. Liu et al. / Electric Power Systems Research 80 (2010) 572–581

581

Table 2 Average relative errors between HCIM, HVSM and CSIM compared to the measured data. Average error (%)

Models

Stage

Electric quantity

Harmonic current injection model (HCIM)

Earlier

EAF voltage EAF current

7.46 11.66

Middle

EAF voltage EAF current

6.19 8.84

Later

EAF voltage EAF current

8.26 8.48

ment has been reached by using the proposed model. To model the EAF by the proposed method at the earlier and middle stages during the refining period, the errors are larger than those obtained in the later stage since the number of used sampling points for the EAF model is not sufficient to fully describe the EAF dynamics in the prior two operating stages. 5. Conclusions In this paper, the proposed and two classic EAF modeling methods based on actual measured data are compared. Field measurements taken from the power system of a steel plant is used to show the accuracy of the compared models through computer simulations. Results have shown that the HCIM and HVSM are less accurate when modeling the vi characteristic of the EAF load. On the other hand, the proposed model for the EAF yields more accurate and satisfactory results. The proposed technique is also suitable for modeling other nonlinear devices in a simple and sufficient manner and can be easily incorporated with a harmonic power flow algorithm for harmonic studies. However, if the used sampling points obtained by actual measurements for modeling are not sufficient, the accuracy of test results may be highly affected. In addition, to model EAF characteristics at different operating stages requires corresponding field measurements or experiments. Future work will focus on the development of more accurate models of the EAF with advanced methods such as neural networks to include different operating stages.

Harmonic voltage source model (HVSM)

6.86 7.45

Proposed model (CSIM)

4.23 5.12

References [1] G.T. Heydt, Electric Power Quality, 1st ed., Stars in a Circle Publications, 1991. [2] T. Zheng, E.B. Makram, A.A. Girgis, Effect of different arc furnace models on voltage distortion, in: Proceedings of the IEEE Eighth International Conference on Harmonics Quality of Power (ICHQP), vol. 2, 1998, pp. 1079–1085. [3] GIGRE Working Group, Harmonics, characteristic parameters, methods of study, estimates of existing values in the network, Electra 36 (5) (1981) 35–54. [4] S. Varadan, E.B. Makram, A.A. Girgis, A new time domain voltage source model for an arc furnace using EMTP, IEEE Trans. Power Deliv. 11 (3) (1996) 1685–1691. [5] M.A.P. Alonso, M.P. Donsion, An improved time domain arc furnace model for harmonic analysis, IEEE Trans. Power Deliv. 19 (1) (2004) 367–373. [6] G.C. Montanari, M. Loggini, A. Cavallini, L. Pitti, D. Zaninelli, Arc-furnace model for the study of flicker compensation in electrical networks, IEEE Trans. Power Deliv. 9 (4) (1994) 2026–2036. [7] L. Tang, S. Kolluri, M.F. McGranaghan, Voltage flicker prediction for two simultaneously operated AC arc furnaces, IEEE Trans. Power Deliv. 12 (2) (1997) 985–992. [8] R. Horton, T.A. Haskew, R.F. Burch IV., A time-domain ac electric arc furnace model for flicker planning studies, IEEE Trans. Power Deliv. 24 (3) (2009). [9] V. Terzija, V. Stanojevic, Power quality indicators estimation using robust Newton-type algorithm, IET Proc. Gen. Transm. Distrib. 151 (4) (2004) 477–485. [10] K.J. Tseng, Y. Wang, D.M. Vilathgamuwa, An experimentally verified hybrid Cassie–Mayr electric arc model for power electronics simulations, IEEE Trans. Power Electron. 12 (3) (1997) 429–436. [11] H. Mokhtari, M. Hejri, A new three phase time-domain model for electric arc furnace using MATLAB, in: Proceedings of the IEEE Asia-Pacific Transmission and Distribution Conference and Exhibition, vol. 3, 2002, pp. 2078–2083. [12] E. Acha, A. Semlyen, N. Rajakovic, A harmonic domain computational package for nonlinear problems and its application to electric arcs, IEEE Trans. Power Deliv. 5 (3) (1990) 1390–1397. [13] E. Acha, M. Madrigal, Power Systems Harmonics: Computer Modeling and Analysis, John Wiley & Sons, England, 2001. [14] E.M. Gisela, U. Frank, Numerical Algorithms with Fortran, Springer, 1996.