Continuous state space model of the ITER pulsed power electrical network for stability analysis

Fusion Engineering and Design 139 (2019) 62–73

Contents lists available at ScienceDirect

Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes

Continuous state space model of the ITER pulsed power electrical network for stability analysis

T

Claudio Finottia, Elena Gaioa, , Ivone Benfattob, Inho Songb, Jun Taob ⁎

a b

Consorzio RFX (CNR, ENEA, INFN, Università di Padova, Acciaierie Venete SpA), Corso Stati Uniti 4, 35127 Padova, Italy ITER Organization, Route de Vinon-sur-Verdon, CS 90 046, 13067, St. Paul Lez Durance Cedex, France

ARTICLE INFO

ABSTRACT

Keywords: ITER Thyristor Controlled Reactors (TCR) ac/dc converters Power systems State space model Eigenvalue analysis Stability analysis Small signal analysis Selective Modal Analysis (SMA)

The ITER large power supplies of the superconducting magnets and heating and current drive systems, fed by the Pulsed Power Electrical Network, can absorb from the grid active and reactive power up to 500 MW and 950 Mvar respectively. In this paper, a new analytical approach based on the state space formulation is proposed to investigate the dynamic stability of such complex system, including the main power components and their control, with the aim of identifying possible instability phenomena due to interactions among the ac/dc converters, reactive power compensation systems and related controllers. This model describes the dynamics of the main components of the ITER pulsed power supply systems and it is based on small signal approach to tackle the non linearity of the ac/dc conversion and reactive power compensation systems; the discrete phenomena have been approximated by continuous transfer functions, so the linear control theory can be applied for the stability analysis, making this method very effective and fast. Moreover a modular approach is adopted; thus this analytical model can be also easily adapted to other similar electrical power systems.

1. Introduction ITER experiment, under construction in Cadarache (France) [1], aims at demonstrating the feasibility of nuclear fusion as viable energy source. Many large and complex plants are necessary for the ITER operation, and one of the most powerful is the pulsed power supply system, unique for complexity and size, which supplies the superconducting coils and the heating and current drive systems. It is composed of many ac/dc converters with a total installed power of about 2 GVA, which can absorb during an experimental pulse a total active and reactive power from the grid up to 500 MW and 950 Mvar respectively [2]. A reactive power (Q) compensation system rated for 750 Mvar and harmonic filters are provided to reduce the reactive power and current harmonics injected into the grid in order to comply with the requirements of the French Transmission System Operator (TSO) [3]. The ratings guarantee the compensation in steady state; however, in transient conditions instability phenomena could occur, due to interactions between the ac/dc conversion, reactive power compensation systems and their controllers. This issue was first studied via a numerical code (using SimPowerSystem™ of Matlab Simulink toolbox) capable of reproducing the instantaneous current and voltage profiles of the whole ITER power supply system [3]. However this approach, besides requiring very long computational times, can show the instability ⁎

evolution only once it is known, since sensitivity tools to identify critical conditions cannot be applied. For this reason, an analytical approach based on the state space formulation, similar to the studies carried out for the High Voltage Direct Current (HVDC) systems [4–6], was selected to analyze these possible instability phenomena. As a first step, a model to investigate power and voltage stability and harmonic analysis in case of quasi static or steady state conditions has been worked out [7]. Then, dynamic analytical models describing all the ITER PS main components, included their control systems were developed and finally joined in a complete model of the whole ITER PS system. To this analytical model, the linear control theory tools, for example the eigenvalue and Selective Modal Analyses (SMA) [8–10], can be applied. This approach is very effective for the identification of instability causes and very helpful to find solutions; moreover, the computation times are very short, thus allowing fast analyses of different ITER scenarios, in terms of current and voltage profiles versus time to be supplied by the converters. The detailed descriptions of the dynamic analytical models of the ac/dc converters and of the reactive power compensation system are given in [11] and [12], while this paper gives an overview of the overall model and shows an example of application for the stability analysis of the ITER PS system. Section 2 introduces the ITER pulsed power supplies. In Section 3,

Corresponding author. E-mail address: [email protected] (E. Gaio).

https://doi.org/10.1016/j.fusengdes.2018.12.058 Received 31 October 2017; Received in revised form 17 December 2018; Accepted 18 December 2018 0920-3796/ © 2018 Elsevier B.V. All rights reserved.

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

the main assumptions for the setup of the state space model are discussed. A modular approach has been followed; for the ac/dc converters and the TCR models, already described in [11] and [12] respectively, the main aspects only are recalled in this paper. The models of the other parts (as harmonic filters, thyristor firing synchronization and electrical network) of the ITER pulsed power supplies are described in detail and the integration of the single modules in the whole model is explained. The analysis is carried out at small signal level, therefore NewtonRaphson power flow method is used to calculate the equilibrium point. In the last part of the paper, a stability analysis of the combined operation of the ac/dc converters, reactive power compensation systems and related controllers has been carried out using the state space model and the results have been validated by numerical simulations of a detailed model of the whole power supply system built using SimPowerSystem™ Matlab toolbox [13], taken as reference. Finally, possible unstable operation conditions of the ITER power supply system have been identified by the state space model and solutions have been proposed to avoid them.

Table 1 Nominal dc output voltage and current of the ac/dc converter base unit.

No Load Voltage [V] Nominal Current Idc_nom [kA]

PF

CS

TF

CC

VS

± 1420

± 1350

± 900/450

± 130/ ± 440

± 1400

± 55

± 45

68

± 10

± 22.5

Table 2 ac/dc converters for each busbar with the indication of the number of units connected in series. Busbar 1 66 kV PF3 (3 units) CS3U (2 units) CS2L (2 units) CS1 (4 units) TF (1 unit)

2. The ITER pulsed power supply system The ITER pulsed power electrical network is shown in Fig. 1. The ITER experiment is supplied by the High Voltage grid at 400 kV with a short circuit power of about 10–12 GVA. The ac/dc conversion systems of the superconducting magnets and of the heating and current drive systems are equally distributed among three independent 66 kV and 22 kV distribution systems supplied by three identical three windings step-down transformers (nominal power 300/250/150 MVA and voltage 400/66/22 kV, % impedance at 300 MVA: xps% = 12%, xpt% = 27%, xst%=15%, where the subscripts indicate the windings: primary, secondary and tertiary). For each 66kV busbar, a 250 MVA reactive power compensation system based on Thyristor Controlled Reactor (TCR) technology and harmonic tuned filters are installed [2,3]. In the first operational phase, ITER will operate with a limited set of ac/dc converters and then an upgrade is foreseen to reach the full performance. In this paper, the configuration with the final set the power supplies has been assumed for the validation of the state space model and its first application for stability analysis. Below a short description of the main components is given; it is underlined that the last ones are still under development and can be modified or improved in the future; therefore their description has to be considered like a possible solution for the control system design. As another preliminary remark, it is specified that for the scope of the work the ac/dc converters of the superconducting coils only have been considered in the state space model, since they require the majority of the power from the grid, thus causing the major impact.

Busbar 2

Busbar 3

22 kV

66 kV

22 kV

66 kV

22 kV

CCU1 (1 unit) CCS1 (1 unit) CCL1 (1 unit)

PF2 (3 units) PF4 (3 units) VS1 (6 units)

CCU2 (1 unit) CCS2 (1 unit) CCL2 (1 unit)

PF5 (3 units) CS2U (2 units) CS3L (2 units) PF (2 units) PF6 (2 units)

CCU3 (1 unit) CCS3 (1 unit) CCL3 (1 unit)

is divided in several groups: Poloidal Field (PF) coils, Central Solenoid (CS) coils, Toroidal Field (TF) coils, Correction (CC) coils (Upper, Side, Lower) and Vertical Stabilization (VS) coils. Table 1 shows the ratings of the ac/dc converter base unit for each coil group and Table 2 indicates the number of basic units to be connected in series to supply each coil. The base units of the PF, CS, and CC are composed of two back to back thyristor bridges (called here subunits) connected in parallel in twelve pulse configuration (Fig. 2). For TF coil, unidirectional current only is requested, so anti-parallel thyristor bridges are not necessary; for VS ac/dc converter, the subunits in parallel are in 6 pulses configuration. The ac/dc converter base unit can operate in different modes depending on the absolute value of the coil current ICoil with respect to the nominal current Idc_nom:

• 12 pulse operation: when the coil current is high (for |I • •

2.1. The ac/dc converters of the superconducting magnets The main features of the ac/dc conversion system are shortly recalled below; further details can be found in [11]. The magnetic system

Coil| > 0.4 · Idc_nom); it is equally distributed between the two subunits connected in parallel Circulating current mode: when the coil current is around zero (for |ICoil| ≤ 0.2 · Idc_nom); an additional circulating current is imposed between the thyristor bridges of the subunits connected in antiparallel to control continuously the coil voltage even during the sign change of the coil current 6 pulse operation: in this phase, only one subunit supplies the coil and this operation mode allows the transition between the 12 pulse and circulating current mode without thyristor misfiring (for 0.2 · Idc_nom < |ICoil| ≤ 0.4 · Idc_nom).

2.2. The reactive power compensation system The maximum reactive power demand of ITER ac/dc converters can reach 950 Mvar, while the maximum reactive power consumption allowed by the TSO is 200 Mvar. For this reason, a reactive power compensation system rated for 250 Mvar will be installed for each 66 kV busbar to compensate for the reactive power demand and to reduce the current harmonics injected by the ac/dc converters into the grid [3]. The reactive power compensation system is composed of TCR and Harmonic Filters connected in parallel (Fig. 3). The branch inductance Lr of the TCR is 112 mH, while the parameters of the filters are given in the Table 3. Moreover for the 3rd and 23rd harmonic filters, a

Fig. 1. Schematic representation of the ITER electrical network. 63

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

Fig. 2. ac/dc converter base unit (left), operation modes versus current (right).

resistance RHP (1 kΩ) is connected in parallel to the RL elements for the filtering of high order harmonics. The scheme of the reactive power compensation system of one busbar is shown in Fig. 3.

Table 3 Main parameters of the Harmonic Filters.

3. The state space model of the whole ITER pulsed power supply Some simplified assumptions are necessary to represent the ITER power supply system by state space formulation, similarly to [5,6]:

• a non-linear algebraic-differential equation system based on the • • •

equations related to the fundamental components for ac systems and mean ones for the dc systems have been used the ac variables have been represented in the dq frame [14], so the ac and dc systems have the same frequency frame the small signal analysis around an equilibrium point has been applied for nonlinear systems the discrete phenomena have been approximated with equivalent continuous transfer functions.

harmonic order h

Capacity [μF]

Inductance [mH]

Resistance [Ω]

3 5 7 11 13 23

16.7144 50.1433 33.4289 33.4289 26.1983 17.4656

68.7219 8.2466 6.3112 2.5558 2.3349 1.1189

0.2699 0.037 0.0496 0.0401 0.0489 0.0703

x = A x + B inp y = C x + D inp

(1)

where Δx is the state vector, Δinp are the inputs (dq components of the busbar voltages or input control signals), and Δy are the outputs (dq components of the currents absorbed by the subsystems or output control signals). The matrixes A,B,C,D are derived from the calculation of the Jacobian matrixes around an equilibrium point of the nonlinear algebraic-differential equations describing the dynamic behavior of the subsystem. In the following, the state space models of the different subsystems are described; for the state space models of ac/dc converters and reactive power compensation system, an introductory section is dedicated to their control system, thus completing the description given in [11,12]

With this last assumption, the continuous state space model can describe phenomena with characteristic frequency lower than half switching frequency of the main power electronic devices (TCR and ac/ dc converters) according to Nyquist–Shannon sampling theorem. Therefore considering that the switching frequency of the thyristors is 300 Hz, the state space model can describe phenomena with maximum frequency lower than 150 Hz. The equilibrium point necessary for the small signal analysis is defined by ITER scenario in terms of voltage and current in the coils. A modular approach has been adopted, so the state space model of the ac/dc converters, of the TCR, of the Harmonic Filters, of the synchronization block to the MV busbars voltage and of the calculation of the relevant rms value have been separately developed, and then integrated in the whole model. All the subsystems state space models can be represented with the small signals approach by the standard notation:

3.1. The state space model of the ac/dc converters 3.1.1. Control system of the ac/dc converters Each unit is equipped with a feedforward voltage controller to quickly apply to the coil the voltage requested by the Plasma Control System. The voltage drop is calculated by multiplying the measurement of the dc output current of each bridge and the equivalent impedance of the step down transformer and of the interphase reactor, and the output of the controller is normalized with respect to the measurement of the phase-to-phase busbar voltage Urms. The derivative of the voltage

Fig. 3. The reactive power compensation system. 64

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

Fig. 4. Block scheme of the ac/dc converter base unit in circulating current mode and its control system.

reference Uconv_unit_ref is limited to avoid fast variations of the thyristor firing angles, which could cause current unbalance between the two subunits in parallel operation. Inner current regulator is implemented to regulate the current sharing between the two parallel subunits in 12 pulse operation or to control the circulating current between the antiparallel thyristor bridges in circulating current mode. The bandwidth of the closed loop of this controller has been set to 24 Hz. Here, the control system in circulating current mode is shown as an example in Fig. 4 (the control systems of the unit in 6 and 12 pulses configuration are described in detail in [11]); the normalized reference ref_n1 for subunit 1 is:

ref _n1 =

exceeded and 0.05 is the margin necessary for the small adjustments carried out by the feedback regulators. Idc_max and Idc_min are functions of the coil current ICoil and they are:

• In 12 pulses operation: I = I = I • In 6 pulses operation: I = I = I • In circulating current mode: I = I + I dc_max

dc_max

2 ktra

3.1.2. Development of the state space model of the ac/dc converter The main aspects of state space model of the ac/dc converter described in [11] are recalled here. The dynamic response of the ac/dc converter can be approximated with a fixed time delay equal to half a switching time, that in a continuous state space model can be represented by a low pass filter with a time constant from 0.833 ms to 1.667 ms depending on the operation mode of the unit, i.e. 12 or 6 pulses or circulating current modes. The model is based on the wellknown equations of the thyristor ac/dc converter. The magnet system is not included in the model; the load is represented by a constant current generator which value is equal to the current circulating in the coil, at the given equilibrium point of the scenario. Since the small signal analysis has been adopted, the derivative of the voltage reference is very low, so the derivative limiter of the voltage reference does not have any effect and it is not included in the model. The equations describing the three operation modes are slightly different due to the dissimilar power and control configurations, therefore three state space models have been developed. Nevertheless, the relevant inputs and outputs for the development of the ITER whole model are the same, making easier their construction. The input variables Δinpconv are the voltage reference ΔUconv_unit_ref, the dq components of the busbar voltage Δed_bus and Δeq_bus, the PLL signal ΔθPLL and the rms grid voltage measurement ΔUrms, while the output variables are ac/dc converter currents Δiconv_unit_d and Δiconv_unit_q. The state space model at the small signals around an equilibrium point can be written as:

Urms (2)

and for the subunit 4, ref_n4 results:

ref _n 4 =

Uconv _ unit _ ref + (Ri + 2 Rtra ) Idc _ m4 + 6 f Ltra Idc _ m4 + out _PI 3

2 ktra

Urms (3)

where Idc_m1 and Idc_m4 are the dc output current measurements of the thyristor bridges, f is the fundamental frequency of the grid, Rtra and Ltra are the equivalent winding resistance and leakage inductance of the step down transformer respectively (secondary side), ktra is the turn ratio of the step down transformer, Ri is the resistance of the interphase reactors and out_PI is the output of the inner current regulator. The normalized references are converted in thyristor firing angles by an arccos function. Finally, in ITER the voltage control strategy of the ac/dc converters composed of more than one unit in series is based on the sequential control, i.e. one unit is used for the voltage regulation, while the other units provide the maximum or minimum available voltages, thus reducing the reactive power demand. The equations for the calculation of the maximum Udc_Max and minimum Udc_Min voltages are:

Udc _ Max =

3

2 Urms [cos ( ktra

Min )

0.05]

(6 f Ltra + Ri + 2 Rtra)

x conv _ unit = A conv _ unit

(4)

Idc _ max (ICoil )

PLL

Udc _ Min =

3

2 ktra

Urms [cos (

Min )

0.05]

Urms

(6 f Ltra + Ri + 2 Rtra)

Idc _ min (ICoil ) + 2 6 f Ltra Idc _ max (ICoil )

Coil/2

Coil

dc_max Coil circ_ref, Idc_min = Icirc_ref where Icirc_ref = 0.2 · Idc_nom. Moreover in this case the maximum and minimum voltage must have the same absolute value, so it yields: Udc_max_ccm = min(|Udc_max|, | Udc_min |) and Udc_min_ccm = -min (|Udc_max|, | Udc_min |)

Uconv _ unit _ ref + (Ri + 2 Rtra) Idc _ m1 + 6 f Ltra Idc _ m1 + out _PI 3

dc_min

dc_min

x conv _ unit + Bconv _ unit _ edq + Bconv _ unit _ ref

i conv _ unit _ d = Cconv _ unit x conv _ unit iconv _ unit _ q

(5)

ebus _ d + Bconv _ unit _ PLL ebus _ q

Uconv _ unit _ ref

(6)

where the states and the calculation of the Jacobian matrixes are defined in [11]. It is highlighted that some inputs and outputs given in [11] are not reported here, because they are internal variables without interaction

where αMin (15°) and γMin (25°) are the minimum firing (in rectifier mode) and margin (in inverter mode) angles respectively to be not 65

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

Fig. 5. The reactive power compensation system (bold line) and its control system (dashed line).

with the other systems; in particular with the selected input and output variables, the matrix Dconv_unit is a zero matrix. As described in the previous section, more units in series can be used in sequential control; due to the resulting different operating points, the state space model (6) for each unit type (regulation, maximum or minimum voltage) has to be calculated individually. Therefore the state space model of the coil ac/dc converter is the union of the state space models of the unit working in regulation (with subscript “Reg”), at the maximum voltage (with subscript “Max”), and at the minimum voltage (with subscript “Min”). The related matrixes become:

ACONV =

BCONV _ PLL =

BCONV _ ref =

CCONV =

kRe g =

=

0 kRe g Bconv _ unit _ Re g _ ref Bconv _ unit _ Re g _ PLL Bconv _ unit _ Max _ PLL + 0 kMax Bconv _ unit _ Max _ ref Bconv _ unit _ Min _ PLL 0 kMin Bconv _ unit _ Min _ ref

(9)

3

2 ktra 3

[cos (

2 ktra

Min )

[cos (

nconv _ unit _ Max kMax

0.05]

Min )

0.05]

nconv _ unit _ Min kMin

(12)

If there are no units working either at maximum or minimum voltage, it is sufficient not to include in the construction of the matrixes. The state space model of the coil ac/dc converter results:

x CONV = ACONV

x CONV + BCONV _ edq + BCONV _ ref

ebus _ d + BCONV _ PLL ebus _ q

UCONV _ ref

iCONV _ d = CCONV x CONV iCONV _ q

(7)

(8)

Cconv _ unit _ Re g nconv _ unit _ Max

kMin =

Udc _ Min Urms

=

Urms

Bconv _ unit _ Re g _ edq Bconv _ unit _ Max _ edq Bconv _ unit _ Min_ edq

Bconv _ unit _ Re g _ ref 0 0

Udc _ Max Urms

PLL

Aconv _ unit _ Re g 0 0 0 Aconv _ unit _ Max 0 0 0 Aconv _ unit _ Min

BCONV _ edq =

kMax =

(13)

3.2. The state space model of the reactive compensation system 3.2.1. Control of the reactive power compensation system The state space model of the TCR is described in detail in [12]; here, the main aspects only are recalled. The description is done with reference to the block scheme of Fig. 5, showing the reactive power compensation and relevant control system. The TCR control strategy aims at satisfying in real time the equation:

QRPC = Qconv

(10)

(14)

where QRPC is the reactive power provided by the reactive power compensation system and Qconv is the total reactive power demand of the ac/dc converter connected at the same busbar. The reactive power measurement of the ac/dc converters is the TCR control input. It is obtained applying the dq transformation [14] to the sum of the three phase currents of the ac/dc converters (iconv_bus_d_m and iconv_bus_q_m) and phase to ground voltages of the related busbars (ebus_d_m and ebus_q_m); then, the dq components are filtered by a low pass filter (obtaining iconv_bus_d_f, iconv_bus_q_m, ebus_d_m and ebus_q_m), and

nconv _ unit _ Min Cconv _ unit _ Min

Cconv _ unit _ Max (11)

where the nconv_unit_Max and nconv_unit_Min are the number of units operating at the maximum and minimum voltage respectively, and the constant kReg, kMax and kMin are derived from (4) and (5): 66

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

the reactive power Qconv_f of the ac/dc converter is calculated by the equation (see also block scheme of Fig. 5):

3 (ebus _ q _ f iconv _ bus _ d _ f 2

Qconv _ f =

ebus _ d _ f iconv _ bus _ q_ f )

xTCR _ L Qconv _ f

PLL

(15)

Urms iTCR _ d = CTCR iTCR _ q

The control signal refTCR is a normalized susceptance, which can vary from 0 (the TCR are not conducting) to 1(the TCR is conducting at its maximum conduction angle at full power), and it is calculated by:

refTCR = 2

f

2 Urms

TCR)

=

2

2

TCR

+ sin (2

TCR )

diLhr _ d dt diLhr _ q dt

deChr _ d dt

(17)

deChr _ q dt

where the angle αTCR is related to the zero cross of the respective phase to phase voltage on the TCR and it can vary from π/2 rad (maximum consumption of reactive power with refTCR = 1) to π rad (consumption of reactive power equal to 0 with refTCR = 0).

xTCR _ L + BL _ edq

ebus _ d + BL _ PLL ebus _ q

PLL

=

= =

ebus _ q

eChr _ d + eChr _ q

ihr _ d +

(ebus _ q iconv _ bus _ d

ebus _ d iconv _ bus _ q )

Lhr

Lhr iLhr _ d

Rhr iLhr _ d Rhr iLhr _ q

(21)

ihr _ q

Chr eChr _ q Chr

Chr eChr _ d

(22)

Chr

ihr _ q = iLhr _ q +

3 2

Lhr iLhr _ q

Lhr

ihr _ d = iLhr _ d +

(23)

ebus _ d eChr _ d RHP ebus _ q eChr _ q

(24)

RHP

From the previous equations the matrixes of the state space model of the harmonic filters can be derived:

where the input variables are the normalized reference ΔrefTCR, the dq busbar voltages Δebus_d and Δebus_q and the PLL signal ΔθPLL. The output variables are the TCR current variations ΔiTCR_d and ΔiTCR_q. Also in this case DL is a zero matrix. To complete the model of the TCR, the control system has to be considered. As shown in Section 3.3, the low pass filters in the dq measurements of Fig. 5 can be moved at the end of the calculation of the reactive power and it can be demonstrated that the synchronization angle θPLL has no impact in the calculation of the reactive power measurement Qconv_f, thus the state equation can be written:

dt

=

ebus _ d

• for high pass filter (for index hr equal to 3 and 23):

+ BL _ ref

(18)

(

=

ihr _ d = iLhr _ d ihr _ q = iLhr _ q

refTCR

Qconv _ f +

xTCR _ L Qconv _ f

• For the tuned filters (for index hr from 5 to 13):

iTCR _ d = CL xTCR _ L iTCR _ q

dQconv _ f

ibus _ conv _ d ibus _ conv _ q

where Chr and Lhr are the capacitance and the inductance of the filter at the hr harmonic order, and eChr_d and eChr_q are the dq voltage components on the capacitor, while iLhr_d and iLhr_q are the dq current components in the inductor and ω is the fundamental angular frequency 314 rad/s. The total dq currents ihr_d and ihr_q adsorbed by each harmonic filter are (see Fig. 5):

3.2.2. State space model of the TCR and control Similarly to ac/dc converters, the time delay of firing angle assumed as the average time between two switching instants is represented in the state space model by a low pass filter with time constant equal to 1.667 ms. Moreover a set of band-pass and band-stop filters describes the TCR switching operation accounting for the three-phases arrangement. The resulting state space model is [12]:

xTCR _ L = AL

+ BTCR _ idq

ebus _ d + BTCR _ PLL ebus _ q

3.2.3. State space model of the harmonic filters The Harmonic Filters equations can be derived by the state equations of the capacitor and inductance in dq frame, which are

(16)

where BC is the equivalent susceptance of the harmonic filters at the fundamental frequency. A look up table converts the control signal refTCR into the thyristor firing angle αTCR by reversing the function h (αTCR) [15]:

h(

+ BTCR _ edq

(20)

Qconv _ f

Lr BC 3

xTCR _ L Qconv _ f

= ATCR

xHR _ FIL = AHR _ FIL xHR _ FIL + BHR _ FIL _ edq

ebus _ d ebus _ q

iHR _ FIL _ d = CHR _ FIL xHR _ FIL + DHR _ FIL _ edq iHR _ FIL _ q

ebus _ d ebus _ q

(25)

where ΔxHR_FIL are the states (voltages and currents of the capacitors and inductors respectively), and ΔiHR_FIL_d and ΔiHR_FIL_q are the sum of the dq current components of all harmonic filters Δihr _d and Δihr _q. 3.3. The state space model of the thyristor synchronization with the busbar voltage and rms voltage calculation

)

TQconv

The thyristor firing of the ac/dc converters and of the TCR have to be synchronized with the instantaneous voltages of the respective supply busbars. As it is shown later, this part of the control system is very important in the stability analysis. In this study, a Phase Locked Loop (PLL) method [16] is used for the thyristor synchronization, where the three input phase to ground voltages ebus_abc are transformed into dq components [14], then the angle error θerr is calculated by the arctan function and it is corrected by the PI controller. The Voltage Controlled Oscillator (VCO) integrates the angular frequency providing the reference angle locked to the fundamental component of the

(19) where iconv_bus_d and iconv_bus_q are the dq components of current adsorbed by the ac/dc converter connected to the busbar, and TQconv is the time constant of the low pass filter. Calculating the new Jacobian matrixes of (16) and (19) respect with the inputs Urms, ebus_d, ebus_q, iconv_bus_d, iconv_bus_q and with the new state Qconv_f, and with few algebraic manipulations, the whole TCR state space model including the control system can be written in the form:

67

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

Fig. 6. Block scheme of the Phase Locked Loop and of the rms voltage calculation in dq frame.

measured busbar voltages. The block scheme of the PLL in dq frame is shown in Fig. 6. The rms value of the phase to phase busbar voltages is given by:

Urms =

3 2 2 (ePLL _ d _ f + ePLL _q_f ) 2

d

dt

cos ( PLL) sin ( sin ( PLL) cos (

PLL )

kp _ PLL = 2

ebus _ d ebus _ q

Ti _ PLL =

(27)

where ebus_d and ebus_q are the dq components of the busbar voltage and ePLL_d and ePLL_q are the results of the rotation. Defining the amplitude ebus and phase angle θV of the busbar voltage:

ebus = V

(

ebus _ q ebus _ d

)

(28)

Ubus =

the angle error θerr results: err

= arctan = arctan

= arctan

ebus cos ( V ) sin ( PLL) + ebus sin ( V ) cos ( PLL ) ebus cos ( v ) cos ( PLL) + ebus sin ( V ) sin ( PLL )

sin ( cos (

ebus _ q ebus _ d

PLL

+ w1 _ PLL

2

fPLL 20 fPLL

(33)

3 2 2 (ePLL _ d + ePLL _ q) 2

(34)

and substituting ePLL_d and ePLL_q with (27), it can easily obtained:

ePLL _ q ePLL _ d

Ubus = (29)

V

PLL )

V

PLL )

=

V

PLL

= arctan

ebus _ q ebus _ d

PLL

dUrms = dt

(30)

(

(

ebus _ q

)

PLL

)

(35)

Urms + Ubus TUrms

(36)

where TUrms is the time constant of the low pass filter. For the construction of the PLL state space model, the input signals are the dq components of the busbar voltage ebus_d and ebus_q (lumped in the vector inpPLL), the synchronization angle θPLL and the rms voltage Urms are the outputs and the states in (32), (36), while the remaining state w1_PLL of the PI integrator is given in (31) (the three states are

Finally the set of the equations can be completed considering the PI regulator and the VCO:

kp _ PLL arctan e _ kp _ PLL err bus d dw1 _ PLL = = dt Ti _ PLL Ti _ PLL

3 2 2 (ebus _ d + ebus _ q) 2

It can be noted that the synchronization angle θPLL does not affect the rms voltage calculation. Finally, the equation of the low pass filter results:

and simplifying by Werner equations, it yields: err

+ w1_ PLL = kp _ PLL arctan

Concerning the rms calculation of the phase to phase voltage Urms, it can be observed that two low pass filters are present in the block scheme, but after the linearization foreseen for the small signal approach this structure can be simplified with only one low pass filter downstream of the square root. So the calculation of the rms voltage Ubus is equal to:

2 2 ebus _ d + ebus _q

= arctan

err

where w1_PLL is the state related to the PI integrator, and kp_PLL and Ti_PLL are the proportional gain and the time constant of the integrator of the PI regulator respectively. It is highlighted that the response of the PLL to a variation of the phase angle of the busbar voltage depends principally on the proportional gain kp_PLL of the PI controller, while the response to the variation of the fundamental frequency depends on the time constant Ti_PLL. Once selected the bandwidth (fPLL) in Hertz of the phase angle synchronization and assuming a factor 20 for the decoupling of the frequency estimation in order to avoid interactions between the two control functions, the parameters of the PI regulator result:

(26)

PLL )

= kp _ PLL

(32)

where ePLL_d_f and ePLL_q_f are the dq voltages at the output of the PLL and the subsequent low pass filter (Fig. 6). In the dq frame, the fundamental frequency corresponds to 0; the same states for the ω0 angular frequency in dq frame, in the PLL model of Fig. 6. This means that the synchronization angle θPLL, which estimates the phase angle of the measured busbar voltage, is without the sawtooth component. The abc→dq transformation can be represented by a rotation matrix of the dq frame and it yields:

ePLL _ d ePLL _ q =

PLL

(31) 68

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

lumped in the vector xPLL=[w1_PLL,θPLL, Urms]). The previous equations can be lumped in a differential system:

where ΔxPS is the vector of the state variables which lumps together all the state variables of the subsystems installed in the different busbars, Δed_source and Δeq_source are the dq components of the supply voltage (Fig. 7), ΔinpPS is the vector of the voltage references of the ac/dc converters installed for each busbar and ΔiPS_d and ΔiPS_q are the dq components of the currents adsorbed by the whole ITER power supply system. The number of the states ΔxPS of the ITER power supply state space model depends on the operation point and it is around 500 on average.

(37)

xPLL = f (xPLL , inpPLL )

and applying the small signal analysis around the equilibrium point p= (xPLL0, inpPLL0), the Jacobian matrixes of the state space model result:

APLL = CPLL

f xPLL

p

BPLL =

f inpPLL

p

0 1 0 0 0 = DPLL = 0 0 1 0 0

(38)

3.5. Calculation of the equilibrium point for a given scenario

3.4. The setting up of the whole state space model

The calculation of the matrixes A, B, C and D for each subsystem is performed around an equilibrium point, and the necessary information for their calculation are the amplitude and phase of the busbar voltages, the power flow at each node and the dc output voltage and current of

The state space model related to one busbar at 66 kV is given by the union of the matrixes of the single subsystems: xPLL xHR _ FIL xTCR = x CONV 1

APLL

0 0 0 AHR _ FIL 0 0 BTCR _ PLL CPLL 0 ATCR BTCR _ idq CCONV 1 0 0 BCONV 1 _ PLL CPLL ACONV 1

x CONVn

BCONVn _ PLL CPLL

0

0

0

0 0

BTCR _ idq CCONVn 0 0 0 ACONVn

xPLL xHR _ FIL xTCR + x CONV 1

BPLL _ edq BHR _ FIL _ edq BTCR _ edq BCONV 1 _ edq

xCONVn

BCONVn _ edq

0 0 ebus _ d 0 + B ebus _ q CONV 1 0

0 0 0 0 0

UCONV 1 _ ref UCONVn _ ref

0 BCONVn

xPLL

ibus _ d = [ 0 CHR _ FIL CTCR CCONV 1 ibus _ q

CCONVn ]

xHR _ FIL xTCR + [DHR _ FIL _ edq] x CONV 1

ebus _ d ebus _ q

(39)

x CONVn

each ac/dc converter. The dc output voltages and currents of the ac/dc converters can be derived easily from voltage and current waveforms in the superconducting coils of the ITER scenario, while the busbar voltages and the power flows at each node are unknown. Therefore, the calculation of equilibrium point has been carried out using the Newton-Raphson power flow method [17]. A short description of the method is given below. The impedance network is shown in Fig. 7, composed of the impedance of the grid and the equivalent impedances of the three stepdown transformers. The total number of the nodes of the network is 11. The active Pinj(k) and reactive Qinj(k) power injected at each node k are equal to:

where Δx… are the state variables, Δebus_d and Δebus_q are dq components of the busbar voltage, ΔUCONV…_ref are the input voltage references of the ac/dc converters, and Δibus_d and Δibus_ q are the dq components of the total current absorbed by all subsystems connected in parallel. The subscripts of the matrixes mean:

• “PLL”: the Phase Locked Loop for the synchronization of the grid • • •

(only one PLL is assumed for the synchronization of all the subsystems connected to the same bus) “HR_FIL”: the Harmonic Filter subsystem, “TCR”: the TCR subsystem (when the TCR thyristors are switched off because of high Q demand of the ac/dc converter, this subsystem is not included in the model) “CONV”: the ac/dc converters connected to the busbar, from 1 to n

xPS = APS

xPS + BPS _ edq_ source

iPS _ d = CPS xPS + DPS _ edq _ source iPS _ q

N = 10 j=0

Pinj(k ) =

It can be noted that there are some cross matrixes between the different subsystems (see not diagonal elements of the state matrix in (39)), because some outputs are the inputs of other subsystems. The state space model of the bus at 22 kV is the same of the 66 kV bus, but without the part related to the reactive power compensation system (i.e. the Harmonic Filters and the TCRs). The state space model of the whole ITER power supply system is obtained by manipulating the state space matrixes of each busbar considering the impedance in series of the step down transformers and, finally, the impedance of the grid; this manipulation introduces a coupling between parallel buses, so the stability of each bus is strictly dependent on the others. The final structure of the analytical model results:

|Vk ||Vj| (Gkj cos(

+ Bkj sin( Qinj(k ) =

N = 10 j =0

k

k

j)

for k = 1, 2, ,

j ))

|Vk ||Vj |(Gkj sin(

k

j)

Bkj cos(

k

j ))

(41) where V and θ are the amplitudes and the phases of the voltages at each node, Gkj and Bkj are the conductances and the susceptances of the admittance matrix. It can be noted that the k index starts from 1, because the source voltage is known, with amplitude and phase equal to 400 kV and 0 electrical degrees respectively. The active Pload and reactive Qload load powers for the nodes 2,3,4 is equal to zero, while for the other nodes is equal to the power sum of the ac/dc converters and of the reactive power compensation system (if any). The quantities Pload and Qload are functions of the voltage and current on the coils of the ITER plasma scenario (known) and of the node voltage amplitude (unknown). The Jacobian matrix for the Newton-Raphson solution has been analytically calculated:

ed _ source + BPS _ inp inpPS eq _ source ed _ source + DPS _ inp inpPSs eq_ soure (40)

69

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

J=

(Pload (k ) + Pinj (k ) )

(Pload (k ) + Pinj (k ) )

Vk

k

(Qload (k ) + Qinj (k ) )

(Qload (k ) + Qinj (k ) )

Vk

k

scenario (see Table 4). Applying the eigenvalue analysis [8–10] to the state matrix related to this operating point, an unstable mode is observed and its eigenvalue is equal to 2.73 ± 377i. For the identification of the instability cause, the participation factor pki of the states to the unstable mode has been calculated by the equation:

(42)

The Newton-Raphson procedure is an iterative process and the injected power and load power for each node are calculated by a first guess of the amplitude Vk(0) and phase angle θ k(0) of the node voltages. The error in the estimation of the power for each node k is calculated:

Pk = Pload (k ) + Pinj (k ) Qk = Qload (k ) + Qinj (k )

pki =

(43)

(44)

where

Vk k

=

J

Pk Qk

1

(45)

Then, the new voltage guess is used for a new calculation of the injected and load power values. The iterative process is stopped when the error in the estimation of the power for each node k is lower than a certain threshold. Finally the amplitudes and phase angles of the node voltages are obtained. The maximum error in the estimation of the power and voltages has been set at about 1% of the nominal values. 4. Validation of the state space model with the Power System Transient Model A first example of application of the State Space Model (StSpM) of the whole ITER power supply system has been a stability analysis in dependence of the PLL bandwidth. In this example, the time constants of the low pass filter in the rms calculation TUrms in (36) and in the reactive power measurement TQconv in (19) have been set to:

TUrms = TQconv =

2

1 fPLL

(47)

where li is the left eigenvector of the state matrix of the unstable mode and the index k indicates the element of li related to kth state of the state space model [10]. It can be noted in Fig. 8 that there are three peaks in the participation factors, which correspond to the state θPLL of the PLLs of the three busbar at 66 kV. Reducing the bandwidth of the PLL this mode becomes stable. The StSpM has been validated by a comparison of the results of the stability analysis above described, with the results obtained by a numerical simulation of the same power system in the same operating conditions, carried out with the Power System Transient Model (PSTM). The PSTM is described in detail in [2] (some modifications have been introduced in order to be aligned with the up-to-date design inputs and control systems described in this paper) and it is implemented in SimPowerSystemTM blockset of Matlab Simulink [12], which reproduces the instantaneous currents and voltages waveforms and represents in detail all the nonlinear and discrete components (thyristor switching, etc…). The simulation with the PSTM model has been carried out to verify the findings from the StSpM; first, the PLL bandwidth for a stable operation (fPLL = 15 Hz) has been assumed to reach the equilibrium condition; then the fPLL bandwith has been increased to 40 Hz at 0.65 s and the system has become unstable, as predicted by StSpM. At the same time, pulsed variations (100 ms long) have been applied to the voltage references of the ac/dc converters both in the StSpM (run with the State Space toolbox of Matlab Simulink [13]) and in the PSTM to compare the waveforms obtained with the two different models. Fig. 9 shows the effect of the perturbation of the equilibrium point in the dq components of the current adsorbed by the whole ITER power supply system for both PSTM and StSpM. Then, the same simulation has been repeated maintaining the stable PLL bandwidth; the results are shown in Fig. 10. It can be noted in both figures that there is a good agreement

The new guess for the node voltages is equal to

Vk (1) = Vk (0) + Vk k (1) = k (0) + k

( Re {lki })2 Re {li} (Re {li})T

(46)

where, as a first tentative, the PLL bandwidth (fPLL) for all 66 kV and 22 kV busbars has been selected equal to 40 Hz. An operating point in terms of coil voltages and currents has been taken from the ITER

Fig. 7. Impedance representation of the ITER electrical network.

70

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

between the PSTM and the State Space Model (StSpM) waveforms. There are only small differences in the amplitudes and frequency oscillations of the current due to the simplified assumptions in the StSpM. The accuracy of the StSpM has been evaluated comparing the growth rate (Fig. 11) and the frequency (Fig. 12) of the unstable mode for different value of the PLL bandwidth for both models. Running the State Space Model in different operation points, it yields that the PLL bandwidth threshold between stable and unstable operation is identified by the StSpM with a maximum error lower than 10% and the error in the estimation of the mode frequency is lower than 5%. Hence, the StSpM has been proved to be useful to quickly find the critical operating conditions of the ITER power supply system, since the study of the eigenvalues provides some sensitivity on the stability of the different operation conditions, Moreover, applying the linear control theory, for example the Selective Modal Analysis, the State Space Model easily allows identifying solutions to solve the instability causes, which can be better investigated in a second moment and in detail by the PSTM, making the analysis faster and more efficient.

Fig. 8. Participation factors pki of the states to the unstable mode.

the ac/dc converters, even if it is not straightforward. No unstable eigenvalues is detected with the PLL bandwidth lower than 15 Hz; the best compromise between the PLL dynamic response and the stability of the whole system is with a PLL bandwidth of 10 Hz. This example of application of the state space model developed has shown how it allows quickly identifying the critical control configuration (PLL bandwidth) and potential unstable conditions during an operating scenario of the ITER power supply system.

5. ITER scenario analysis by state space model Another example of application of the StSpM is the verification of critical phases of the ITER scenario. It lasts 1 hour, but for sake of simplicity in this example, a short scenario of 11 s composed of points sampled every 10 ms has been produced with the most salient phases relevant for ITER power supply system operation. The total power consumptions (active Pconv, reactive Qconv, and apparent Sconv) of the ac/dc converters in this short scenario are shown in Fig. 13. An eigenvalue analysis of the state matrix has been carried out for each operation point of this scenario for different values of the PLL bandwidth. The stability analysis highlights three modes (one for each 66 kV busbar) which can be stable or not depending on the operation point and on the bandwidth of the PLL, and with a frequency in the range 60–65 Hz (forward sequence) in dq frame. In the alternating frame, the frequency of these modes is 110–115 Hz, which is exactly in the range of the resonance introduced by the harmonic filters. Therefore the explanation of this instability phenomenon could be that the ac/dc converters (the TCR contribution is very low) can inject into the grid non integer current harmonics during transient conditions, which trigger the resonance. As a consequence, an oscillation occurs in the busbar voltage, and finally it appears in the PLL output reentering in the ac/dc converters. If the PLL bandwidth is high enough, this oscillation is not filtered, giving the instability. The maximum real part of the potential most unstable eigenvalue (i.e. the highest real part among the three identified modes) is shown in Fig. 14 and, comparing with Fig. 13, it can be noted that the real part of the instability is correlated to the total apparent power consumption of

6. Conclusions This paper proposes a new analytical approach for the stability analysis of the ITER pulsed power supply system. This method based on state space formulation describes the dynamic behavior of the main components and is capable to identify possible instability phenomena due to interactions among the ac/dc converters, reactive power compensation systems and related controllers. An example of instability has been found by the eigenvalue analysis of the state matrix and its cause has been easily identified by the calculation of the participation factors, and a corrective action has been proposed. The analysis results obtained by the developed State Space Model (StSpM) of the whole ITER pulsed power supply system have been confirmed by comparison with a numerical simulation of the same power system in the same operating conditions, made with another model (PTSM), capable to reproduce the instantaneous circuit waveforms. The comparison has given very good results thus validating the StSpM results and proving its capability to quickly identify critical operating conditions and possible solutions. This sensitivity cannot be provided by programs like PTSM, which are also very much heavier and time consuming and should be used for verification of the results only.

Table 4 Voltages and currents in the superconducting coils in one operating points of a ITER scenario. Busbar 1

Busbar 2

Busbar 3

66 kV

22 kV

66 kV

22 kV

66 kV

22 kV

IPF3 = 33kA VPF3=-25V ICS3U = 15kA VCS3U = 99V ICS2L = 16kA VCS2L = 26V ICS1 = 36kA V CS1 = 4V ITF = 68kA VTF = 29V

ICCU1 = 9.9kA VCCU1 = 24V ICCS1 = 9.9kA VCCS1 = 24V ICCL1 = 9.9kA VCCL1 = 24V

IPF2 = 22kA VPF3 = 54V IPF4 = 26 kA VPF4=-3V IVS1 = 7.6kA VVS1=-0.7kV

ICCU2 = 9.9kA VCCU2 = 24V ICCS2 = 9.9kA VCCS2 = 24V ICCL2 = 9.9kA VCCL2 = 24V

IPF5 = 37kA VPF5 = 21V ICS2U = 10kA VCS2U = 84V ICS3L = 9.8 VCS3L=-10V IPF1 = 34kA VPF1 = 26V IPF6 = 38kA VPF6=-17V

ICCU3 = 9.9kA VCCU3 = 24V ICCS3 = 9.9kA VCCS3 = 24V ICCL3 = 9.9kA VCCL3 = 24V

71

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al.

Fig. 12. Frequency of the unstable mode for different value of the PLL bandwidth for both models. Fig. 9. Variation of the dq currents adsorbed by the whole ITER power supply after the perturbation of the equilibrium condition for an unstable case (fPLL = 40 Hz) for both models.

Fig. 13. The total power consumptions (active Pconv, reactive Qconv, apparent Sconv) of the ac/dc converters during the short scenario.

Fig. 10. Variation of the dq currents adsorbed by the whole ITER power supply after the perturbation of the equilibrium condition for a stable case (fPLL = 15 Hz) for both models.

Fig. 14. Real part of the unstable eigenvalue.

Acknowledgment and disclaimer The views and opinions expressed herein do not necessarily reflect those of the European Commission and of the ITER Organization. Fig. 11. Growth rate of the unstable mode for different value of the PLL bandwidth for both models.

References [1] ITER experiment web site, www.iter.org. [2] J. Tao, I. Benfatto, J.K. Goff, A. Mankani, F. Milani, I. Song, H. Tan, J. Thomsen, ITER Coil Power Supply and Distribution System, IEEE/NPSS 24th Symposium on Fusion Engineering (2011), https://doi.org/10.1109/SOFE.2011.6052201. [3] A.D. Mankani, I. Benfatto, J. Tao, J.K. Goff, J. Hourtoule, J. Gascon, D. CardosoRodrigues, B. Gadeau, The ITER reactive power compensation and harmonic filtering (RPC & HF) system: stability & performance, IEEE/NPSS 24th Symposium on Fusion Engineering (2011), https://doi.org/10.1109/SOFE.2011.6052341. [4] S. Chiniforoosh, J. Jatskevich, A. Yazdani, V. Sood, V. Dinavahi, J.A. Martinez, A. Ramirez, Definitions and applications of dynamic average models for analysis of power systems, IEEE Transactions on Power Delivery 25 (4) (2010) 2655–2669, https://doi.org/10.1109/TPWRD.2010.2043859.

The state space models can be easily updated with further details coming from the final design and manufacturing of the ITER pulsed power supply system and can be used as an aid to verify different operating scenarios in terms of stable operation of the converters and reactive power compensation systems. Moreover, they can be integrated with new features as an aid for the design phase of new fusion power plants like DEMO.

72

Fusion Engineering and Design 139 (2019) 62–73

C. Finotti et al. [5] D. Jovcic, N. Pahalawaththa, M. Zavahir, Analytical modelling of HVDC-HVAC systems, IEEE Trans. Power Deliv. 14 (2) (1999) 506–511, https://doi.org/10. 1109/61.754095. [6] D. Jovcic, N. Pahalawaththa, M. Zavahir, Small signal analysis of HVDC-HVAC interactions, IEEE Trans. Power Deliv. 14 (2) (1999) 525–530, https://doi.org/10. 1109/61.754098. [7] C. Finotti, E. Gaio, I. Benfatto, A.D. Mankani, J. Tao, Analytical model for stability analysis in high power ac/dc converters applied to the ITER case, IEEE 15th Int. Conf. on Harmonics and Quality of Power (2012) 209–215, https://doi.org/10. 1109/ICHQP.2012.6381207. [8] L. Rouco, Eigenvalue-based methods for analysis and control of power system oscillations, IEEE Colloquium on Power System Dynamics Stabilisation (1998), https://doi.org/10.1049/ic:19980031. [9] F.C. Schweppe, Selective modal analysis with applications to electric power systems, PART I: heuristic introduction, IEEE Transactions on Power Apparatus and Systems (9) (1982) 3117–3125, https://doi.org/10.1109/TPAS.1982.317524 PAS101. [10] W.A. Hashlamoun, M.A. Hassouneh, E.H. Abed, New Results on Modal Participation Factors: Revealing a Previously Unknown Dichotomy, IEEE Transactions on Automatic Control 54 (7) (2009) 1439–1449, https://doi.org/10.1109/TAC.2009.

2019796. [11] C. Finotti, E. Gaio, I. Benfatto, I. Song, J. Tao, Continuous state-space model in dq frame of the thyristor AC/DC converters for stability analysis of ITER pulsed power electrical system, IEEE Trans. Plasma Sci. 44 (11) (2016), https://doi.org/10.1109/ TPS.2016.2608947. [12] C. Finotti, Elena Gaio, Continuous model in dq frame of Thyristor Controlled Reactors for stability analysis of high power electrical systems, Int. J. Electr. Power Energy Syst. 63 (2014) 836–845, https://doi.org/10.1016/j.ijepes.2014.06.045. [13] MATLAB Software by MathWorks: www.mathworks.com. [14] P. Kundur, N.J. Balu, M.G. Lauby, Power System Stability and Control, McGrawHill, (1994). [15] J. Dixon, L. Moran, J. Rodriguez, R. Domke, Reactive power compensation technologies: state-of-the-art review, Proc. IEEE 93 (12) (2005) 2144–2164, https://doi. org/10.1109/JPROC.2005.859937. [16] F. Blaabjerg, R. Teodorescu, M. Liserre, A.V. Timbus, Overview of control and grid synchronization for distributed power generation systems, IEEE Trans. Ind. Electron. 53 (5) (2006) 1398–1409, https://doi.org/10.1109/TIE.2006.881997. [17] B. Stott, Review of load-flow calculation methods, Proc. IEEE 62 (7) (1974) 916–929, https://doi.org/10.1109/PROC.1974.9544.

73