Optimal Load Restoration in Active Distribution Networks Complying with Starting Transients of Induction Motors
H. Sekhavatmanesh, J. Rodrigues, C. L. Moreira, J.A.P. Lopes, R. Cherkaoui
1 Abstract โ Large horsepower induction motors play a critical role as industrial drives in production facilities. The operational safety of distribution networks during the starting transients of these motor loads is a critical concern for the operators. In this paper, an analytical and convex optimization model is derived representing the starting transients of the induction motor in a semi-static fashion. This model is used to find the optimal energization sequence of different loads (static and motor loads) following an outage in a distribution network. The optimization problem includes the optimal control of the converter-based DGs and autotransformers that are used for the induction motor starting. These models together with the semi-static model of the induction motor are integrated into a relaxed power flow formulation resulting in a Mixed-Integer Second Order Cone Programming (SOCP) problem. This formulation represents the transient operational limits that are imposed by different protection devices both in the motor side and network side. The functionality of the proposed optimization problem is evaluated in the case of a large-scale test study and under different simulation scenarios. The feasibility and accuracy of the optimization results are validated using I) off-line time-domain simulations, and II) a Power Hardware-In-the-Loop experiment.
Index Terms โAutotransformer, Convex optimization, Distribution network, DG converter set points, Induction motor starting, Load energization sequence, Relaxed power flow formulation, Semi-static model.
NOMENCLATURE A. Parameters (unless mentioned, all are in p.u.) ๐ท ๐ Importance factor of the load at bus ๐ . ๐น ๐๐๐กโ Square of the nominal thermal ampacity limit of line ๐๐ . ๐ป ๐ Inertial constant of the motor m (sec.) ๐ ๐๐๐ฅ Number of steps assigned to each starting motor ๐๐ ๐ (๐๐ ๐ ) Active (Reactive) load voltage sensitivity at bus ๐ . ๐ ๐,๐ก0 (๐ ๐,๐ก0 ) Active (Reactive) nominal load at bus ๐ at time ๐ก . ๐ ๐๐๐ฅ,๐๐ท๐บ Ampacity limit of the DG converter at node ๐ . ๐พ๐ ๐ Coefficient of the friction and windage loss in the motor m. ๐ฟ ๐,๐ก0 Load energization status at node i and time t according to the steady-state analysis (1/0). ๐ ๐๐ (๐ฅ ๐๐ ) Resistance (Reactance) of line ๐๐ . ๐ ๐,๐ Slip value of the motor ๐ at step k. โ๐ ๐ Slip interval between two successive steps for the motor ๐ . ๐ค ๐๐ , ๐ค ๐๐ Weighting factors of the objective function terms. B. Variables (unless mentioned, all are in p.u.) ๐น ๐๐,๐,๐ Square of current flow magnitude in line ๐๐, at step ๐ during the starting of the motor ๐ . ๐น๐ ๐,๐๐ท๐บ (๐น๐ ๐,๐๐ท๐บ ) Square of the active/reactive current references of the DG converter at node i for the starting of the motor ๐ . ๐ฟ ๐,๐ก Binary decision variable indicating if the load at node i is energized at time t or not (1/0) ๐ ๐,๐,๐๐ท (๐ ๐,๐,๐๐ท ) Active (Reactive) load power at bus ๐ , at step ๐ during the starting of the motor ๐ . ๐ ๐๐,๐,๐ (๐ ๐๐,๐,๐ ) Active (Reactive) power flowing in line ๐๐ , at step ๐ during the starting of the motor ๐ . ๐ ๐,๐,๐๐๐ข๐ (๐ ๐,๐,๐๐๐ข๐ ) Active (Reactive) power from the substation node ๐ , at step ๐ during the starting of the motor ๐ . ๐ ๐,๐,๐๐ท๐บ (๐ ๐,๐,๐๐ท๐บ ) Active (Reactive) power injection from the DG node ๐ , at step ๐ during the starting of the motor ๐ . ๐ฅ๐ ๐ Integer variable representing the tap position of the autotransformer at the motor ๐ . ๐กฬ ๐,๐ An approximation for the acceleration time of the motor m until step k (sec.). โ๐กฬ ๐,๐ An approximation for the time lengths of the step k in the acceleration period of motor m (sec.). ๐ ๐,๐๐๐๐ Electrical torque of the motor m , at step k . ๐ ๐,๐๐๐๐ Mechanical torque of the motor m , at step k . ๐ ๐,๐,๐ Square of voltage magnitude at bus ๐, at step ๐ during the starting of the motor ๐ . C. Indices ๐, ๐
Index of nodes ๐๐ Index of branches ๐ Index of slip step ๐ Index of nodes hosting motor loads in the off-outage area ๐ก Index of time step D. Sets ๐ Set of nodes ๐ ๐ Set of nodes hosting motor loads ๐ ๐ Set of protected nodes ๐ ๐ Set of nodes hosting static loads ๐ โ Set of nodes in the off-outage area ๐ ๐โ Set of nodes hosting motor loads in the off-outage area ๐ set of all the time samples in the restorative period ๐ Set of lines ๐ ๐ Set of protected lines I. I NTRODUCTION
Large-scale induction motors have been used for more than a century as industrial drives for compressors, pumps, fans, or blowers[1]. The large and reactive current driven by induction motors during the starting period can impose high risks both in the motor side and in the network side [2]. The analysis of the motor acceleration transients has been studied in the literature in the context of many power system applications. For example, the effect of motor reacceleration is studied in power system restoration [3]โ[5], voltage and frequency stability of islanded microgrids [6]โ[10], and protection settings within industrial facilities [1], [11]. In these papers, the behavior of the induction motor is evaluated using only time-domain simulations. These approaches would be very time-consuming if they are used for decision-making problems with huge and complex solution spaces. For such problems, it is
Optimal Load Restoration in Active Distribution Networks Complying with Starting Transients of Induction Motors
H. Sekhavatmanesh,
Member, IEEE , J. Rodrigues, C. L. Moreira, J.A.P. Lopes,
Fellow, IEEE , R. Cherkaoui,
Senior, IEEE Fig. 1. The discretized toque-speed curve of the induction motor needed to formulate the motor starting dynamics in an analytical way such that it can be integrated to an optimization problem. Such analytical approaches are proposed in [12]โ[15]. The authors of [12] developed nonlinear differential algebraic formulations to estimate the voltage dip during the motor starting. This voltage dip is predicted in [13] using neural network for an induction motor with a certain kVA capacity installed on a bus with a certain short circuit capacity. The developed formulations are all nonlinear and non-convex. Therefore, they cannot be integrated into the convex optimization problems. The presented formulation in [14] is incorporated into a maximum restorable load problem and solved iteratively using a heuristic approach. This approach is applicable only in the case of simple problems involving only one decision variable. The authors in [15] presented a quadratic optimization problem for the minimization of the voltage deviation with respect to the nominal value during the motor starting. In that paper, the motor reactive power during the acceleration period is approximated as a simple algebraic function of the terminal voltage magnitude. This simplified model does not account for the dynamics of the motor starting. In consequence, it cannot evaluate correctly the feasibility of the motor starting with respect to the operational safety constraints. In this paper, we study the induction motor starting within the context of the load restoration problem. When a fault occurs in a distribution network, once it is isolated, the area downstream to the fault place remains unsupplied which is called the off-outage area . This area is re-energized by the healthy neighboring feeders using switching operations. This new configuration remains for a so-called restorative period until the faulted element is repaired. During this period, it is aimed to restore the loads in the most optimal way such that the total energy not supplied is minimized. In this regard, the network security constraints must be respected especially in case of starting large motor loads. The problem of finding the optimal load energization sequence given a new network configuration is referred as the load restoration problem. Apart from the load energization sequence, another decision variable that can support the network security constraints during the motor acceleration is the control of converter-interfaced generation units. In [16]โ[19], different control strategies are proposed for the control of the DG converter under grid fault condition. The aim is to support the low-voltage right through (LVRT) capability of the DG. For achieving this aim, these control strategies focus either only on the active power or only s R ls jX m jX lr jX S r R i U A B Fig. 2. The equivalent circuit of the induction motor. on the reactive power injection of the DG. This weakness is addressed in [20] and [21] by considering the resistive characteristics of the distribution line impedances and controlling both active and reactive power set points of the DG converter. However, these control methodologies are aimed to support the voltage only at a single node. Therefore, the network is simplified using its Theveninโs equivalent seen from that node. In this paper, an analytical optimization model is derived for the load restoration problem. The main goal is to derive the optimal energization sequence of different type of loads. The problem is formulated as MISOCP and solved using Gurobi solver. Compared with the state-of-the-art algorithms proposed so far for the restoration problem, the major contributions of this paper are the following: -
The proposed load restoration problem incorporates a convex semi-static model of induction motor loads representing their starting transients. -
A convex model is derived for the converter-interfaced DGs working in constant current mode following voltage sags induced by motor load startings. This model is integrated into the optimal load restoration problem in order to obtain the optimal current set points of these DGs. Using these set points, the DGs support the electrical safety constraints in an optimal fashion during motor starting transients. -
The developed optimization problem includes a convex model of the autotransformer that is used for the starting of the induction motor. Therefore, the optimal tap setting of this autotransformer is derived from the proposed optimization problem. -
The transient operational limits imposed by under-voltage and over-current protection devices are integrated into the developed optimization problem. The aim is to guarantee that the starting transients of motor loads will not trigger the protection devices that exist in the distribution network. II.
MODELLING OF THE INDUCTION MOTOR STARTING
In this section, a semi-static model is proposed for the starting dynamics of the induction motor such that it can be integrated into the power flow formulation in a convex fashion. The aim is to formulate the operational safety constraints in the whole distribution network during the motor starting period. Fig. 1 shows the general trajectory of the electrical torque ( ๐ ๐๐๐ ) generated by an induction motor during its acceleration assuming that the motor voltage is fixed. The load torque ( ๐ ๐๐๐๐ ) and the acceleration torques ( ๐ ๐๐๐ ) are also shown in this figure. The former is defined as the summation of the mechanical torque ( ๐ ๐๐๐ ) on the shaft and the friction and 3 windage loss ( ๐พ๐ ๐ (1 โ ๐ ๐,๐ ) ). The latter is the difference between the electrical and load torques. The electrical torque of the motor load at node ๐ is formulated in (1) using the equivalent circuit of the induction machine shown in Fig. 2, where ๐ ๐ and ๐ ๐๐ represent the resistance and reactance of the stator; ๐ ๐ and ๐ ๐๐ are the resistance and reactance of the rotor; and ๐ ๐ is the magnetization reactance [12]. ๐ ๐,๐๐๐๐ = ๐ ๐ ๐ ๐กโ /๐ ๐,๐ (๐ ๐ /๐ ๐,๐ + ๐ ๐กโ ) + (๐ ๐๐ + ๐ ๐กโ ) (1) where, ๐ ๐,๐ = (๐ ๐ โ ๐ ๐ )/๐ ๐ ๐ ๐กโ = ๐ ๐ ๐ ๐2 ๐ ๐ 2 +(๐ ๐๐ 2 +๐ ๐2 ) ๐ ๐กโ = ๐ ๐ ๐ ๐2 ๐ ๐ 2 +(๐ ๐๐ 2 +๐ ๐2 ) ๐ ๐กโ = ๐ ๐ 2 ๐ ๐ +๐ ๐๐ ๐ ๐ (๐ ๐๐ +๐ ๐ )๐ ๐ 2 +(๐ ๐๐ 2 +๐ ๐2 ) ๐ ๐ and ๐ ๐ are the synchronous speed and the rotor speed of the motor at node ๐ , respectively. ๐ ๐ and ๐ ๐๐ represent the resistance and reactance of the stator. ๐ ๐ and ๐ ๐๐ are the resistance and reactance of the rotor. ๐ ๐ is the magnetization reactance. ๐ ๐กโ , ๐ ๐กโ and ๐ ๐กโ are the Theveninโs voltage square, resistance and reactance seen from the rotor terminals (AB in Fig. 2), respectively. Inserting the expressions of these equivalent parameters in (1) results in a formulation for ๐ ๐,๐๐๐๐ that is a linear function of the square of the voltage terminal ( ๐ ๐ ) and non-linear function of the slip ( ๐ ). As illustrated in Fig. 1, in order to make a convex model, we discretize the slip range from the standstill ( ๐ = 1 ) until the stable point into ๐พ ๐๐๐ฅ fixed steps with equal length of โ๐ . In this regard, the slip value during each step k is assumed fixed and equal to a given value ๐ ๐,๐ , which is considered as a parameter. This assumption is justified according to the discussion provided in [2]. In this paper, the terminology of slip step (or shortly step ) refers to each of these discretized intervals within the transient acceleration period of each energized motor load . As discussed in [2], the step length โ๐ should be small enough depending on the total inertia of the induction motor. The next step is to derive the time duration of a given step ๐ indicated by โ๐ก ๐,๐ . The inverse of this time length is obtained using the dynamic motion equation given in (2) [22]. ๐,๐ = โ12๐ป ๐ . ฮ๐ ๐ (๐ ๐,๐๐๐๐ โ ๐ ๐,๐๐๐๐ โ ๐พ๐ ๐ (1 โ ๐ ๐,๐ )) (2) ๐กฬ ๐,๐ = โ โ๐กฬ ๐ โ ,๐๐๐ โ =1 (3) Since the acceleration torque during each step is fixed, the time derivative of the slip is represented by โ๐โ๐ก . In order to derive โ๐ก ๐ , the piece-wise linear approximation method is used as explained in Appendix A. In this regard, we add (13)-(16) to the set of constraints, replacing ๐ฅ and ๐(๐ฅ) with ๐,๐ and ฮ๐ก ๐,๐ , These slip steps ( ๐ ) should not be confused with time steps ( ๐ก ), which refer to low-resolution intervals within the whole restorative period. respectively. Therefore, โ๐กฬ ๐,๐ is obtained in a linear way as an approximation for ฮ๐ก ๐,๐ . In (3), we obtain an approximation for the acceleration time of motor m until step k ( ๐กฬ ๐,๐ ) by adding the approximated time lengths of all previous steps to step k ( โ๐กฬ ๐,๐ ). The obtained variable ๐กฬ ๐,๐ will be used in section III.C to derive the transient voltage and current limits. During each single step, since the slip and therefore all the parameters of the motor equivalent circuit shown in Fig. 1 are fixed, the electrical state variables can be represented in the phasor domain. The aim of the next section is to obtain the values of these state variables for a given step ๐ and in the whole distribution network using AC power flow equations. In this paper, we neglect the DC term imbedded in the starting current of the induction motor. This assumption is justified because of the low X/R ratio in distribution networks. Therefore, the DC term of the starting current disappears shortly. III. R ESTORATION P ROBLEM F ORMULATION
In this section, the load restoration problem is presented as an example to show how the proposed semi-static model of the motor starting could be integrated into the power flow formulation. Due to the difficulties of integrating the transient constraints related to the motor load starting in the full restoration problem formulation, a two-stage approach is exploited. In the first stage, we solve the restoration problem according to the formulation presented in [23]. This restoration problem contains the model of passive and active elements only in steady-state conditions. Solving this optimization problem provides I) the optimal configuration of the network (line switching variables) and II) the optimal load restoration sequences during the restorative period. This stage of analysis is referred in this paper as the steady-state analysis . In the second stage, which is studied in this paper, we take the restoration solution obtained from the first stage and we modify it concerning the transient constraints of motor load starting. In this stage , we modify only the obtained load restoration sequence, while considering the starting dynamics of induction motors, which were neglected in the first stage. It means that the line switching variables are fixed to the ones obtained from the steady-state analysis. We assume that considering the starting dynamics of motor loads does not affect the optimal network configuration that was obtained during the steady-state analysis. In this second stage, we just change the energization sequence of the loads in the most optimal way such that the transient constraints are respected during the starting of each induction motor in the off-outage area The main objective is to minimize the resulting increase of energy not supplied referring to its value obtained from the steady-state analysis. The main decision variables of the proposed optimization problem are three folds, namely, I) the energization sequence of 4 different loads during the time ( ๐ฟ ๐,๐ก ), II) the optimal tap setting of the autotransformer that is used for the starting of the induction motor ( ๐ฅ๐ ๐ ), and III) the optimal active/reactive current injections by dispatchable DGs ( ๐น๐ ๐,๐๐ท๐บ /๐น๐ ๐,๐๐ท๐บ ). The optimization problem is formulated in the form of a MISOCP. The optimization problem in the first stage (steady-state analysis) includes power flow formulation for each time step t during the restorative period. However, in the second stage (studied in this paper), the optimization problem includes the power flow formulation for each step ๐ in the acceleration period of each motor load ๐ in the off-outage area. In this paper, it is assumed that only one motor load can be started at time and only once the starting transient of any other motor disappeared. As suggested in [1], the motor loads in an industrial plant are categorized into groups mainly based on their functional processes. These groups of motor loads are considered to be restored in successive time steps with certain intervals. Only the motor loads that are in the same group are restored simultaneously. In this regard, a given motor load in the proposed optimization problem can represent a group of motor loads in the LV network. The dynamic parameters of this aggregated motor load are specified according to the strategy given in [8]. According the assumption mention above, to each motor load in the off-outage area a specific set of steps k ( ๐ โ{1,2, . . , ๐ ๐๐๐ฅ } ) is assigned. Therefore, all the electrical state variables are indexed with ๐ and ๐ . The constraints involving these indices should hold for all the steps and all the motors in the off-outage area. Minimize : ๐น ๐๐๐ = ๐ ๐๐ . ๐น ๐๐ + ๐ ๐๐ . ๐น ๐๐ (4) ๐น ๐๐ = โ โ ๐ท ๐ . (๐ฟ ๐,๐ก0 โ ๐ฟ ๐,๐ก ). ๐ ๐,๐ก0๐กโ๐๐โ๐ (5) ๐น ๐๐ = โ โ โ ๐ ๐๐ . ๐น ๐๐,๐,๐ijโ๐๐ ๐๐๐ฅ ๐=1๐โ ๐ ๐โ (6) Subject to: {๐ฟ ๐,๐ก โค ๐ฟ ๐,๐ก0 , ๐ฟ ๐,๐ก โค ๐ฟ ๐,๐ก+1 โถ ๐ โ N โ ๐ฟ ๐,๐ก = 1 โถ ๐ โ N\N โ โ๐ก โ ๐ (7) ๐ณ๐๐๐ ๐ด๐๐ ๐๐๐๐๐ (8)
๐จ๐ช ๐๐๐๐๐ ๐๐๐๐ ๐๐๐๐๐๐๐๐๐๐๐ (9)
๐ซ๐ฎ ๐ด๐๐ ๐๐๐๐๐ (10)
๐จ๐๐๐๐๐๐๐๐๐๐๐๐๐๐ ๐ด๐๐ ๐๐๐๐๐ (11)
๐ป๐๐๐๐๐๐๐๐ ๐ช๐๐๐๐๐๐๐๐๐๐ (12)
The objective function ( ๐น ๐๐๐ ) is formulated in (4) as the weighted sum of the reliability ( ๐น ๐๐ ) and operational ( ๐น ๐๐ ) objective terms. As mentioned earlier, the main objective is to minimize the unsupplied energy of loads due to the shifting in their energization times, which is referred in this paper as the reliability objective term. This energy is calculated in (5) summing the power of the all the loads that are not restored ( ๐ฟ ๐,๐ก = 0 ), whereas they were commanded to be restored according to the steady state analysis ( ๐ฟ ๐,๐ก0 = 1 ). Unlike the reliability term, the operational term has a very small weighting coefficient. This term is expressed in (6) as the total active line power losses in the distribution network. This term is included in the objective function just to satisfy the exactness condition stated in [24] for the relaxed AC power flow formulation. According to this condition, the objective function of a minimization (/maximization) problem should strictly increase (/decrease) with the total line power losses in the network. Therefore, the squared current variables will be bounded at optimal values, ensuring the exactness of the optimal solution. Constraint (7) enforces the load at node ๐ and time ๐ก to remain unrestored ( ๐ฟ ๐,๐ก = 0 ) if it is commended to be so according to the results of the steady-state analysis ( ๐ฟ ๐,๐ก0 = 0 ). Moreover, according to (7), once a load is restored ( ๐ฟ ๐,๐ก = 1 ), it should remain supplied during the rest of the restorative period ( ๐ฟ ๐,๐ก+1 = 1 ). The loads that are not in the off-outage area should remain always supplied. A. Load Modeling
Equations (8 .a) and (8 .b) extract, respectively, ๐ ๐,๐0 and ๐ ๐,๐0 , as the nominal active and reactive powers of the load at node ๐ and at the time instant ๐ก , when the motor ๐ is started (๐ฟ ๐,๐ก โ ๐ฟ ๐,๐กโ1 = 1) . At this time t , if the load at node i is not restored ( ๐ฟ ๐,๐ก = 0 ), then ๐ ๐,๐0 and ๐ ๐,๐0 will be zero. The product of binary variables in (8 .a) and (8 .b) introduce non-linear terms. There terms are linearized according to the reformulation technique proposed in [25]. ๐ ๐,๐0 = โ ๐ ๐,๐ก0 . ๐ฟ ๐,๐ก . (๐ฟ ๐,๐ก โ ๐ฟ ๐,๐กโ1 ) ๐กโ๐ โ๐ โ ๐, โ๐ โ ๐ ๐ (8.a) ๐ ๐,๐0 = โ ๐ ๐,๐ก0 . ๐ฟ ๐,๐ก . (๐ฟ ๐,๐ก โ ๐ฟ ๐,๐กโ1 ) ๐กโ๐ โ๐ โ ๐, โ๐ โ ๐ ๐ (8.b) In the following, we use ๐ ๐,๐0 and ๐ ๐,๐0 to formulate the active and reactive power consumptions of different type of loads. We start with the static loads. The active and reactive power of the static load at node ๐ and at step ๐ , during the starting of motor load ๐ are expressed in (8.c) and (8.d), respectively. In this regard, the exponential model is used. Assuming that ๐ ๐,๐,๐ is close to , the model is linearized using the binomial approximation approach as proposed in [23]. The products of the binary variables ๐ฟ ๐,๐ก (in the formulation of ๐ ๐,๐0 and ๐ ๐,๐0 ) and the positive continuous variable ๐ ๐,๐,๐ introduce non-linear terms in (8.c) and (8.d). In order to preserve the linearity, these terms are re-formulated as given in Appendix B. ๐ ๐,๐,๐๐ท = ๐ ๐,๐0 (1 + (๐ ๐,๐,๐ โ 1)) ๐๐ ๐ /2 โ โ ๐ ๐,๐0 (1 + ๐๐ ๐ ๐,๐,๐ โ 1)) โ๐ โ ๐ ๐ , โ๐ โ ๐ ๐ , โ๐ (8.c) ๐ ๐,๐,๐๐ท = ๐ ๐,๐0 (1 + (๐ ๐,๐,๐ โ 1)) ๐๐ ๐ /2 โ โ ๐ ๐,๐0 (1 + ๐๐ ๐ ๐,๐,๐ โ 1)) โ๐ โ ๐ ๐ , โ๐ โ ๐ ๐ , โ๐ (8.d) Now, we move to the formulation of the motor load powers. The active and reactive power of the motor load at node ๐ and at step ๐ , during the starting of motor load ๐ are expressed in the following. For the motor load that is starting ( ๐ = ๐ ), (8.e) and (8.g) express the active and reactive power consumptions according to the equivalent circuit shown in Fig. 2. The other motor loads ( ๐ โ ๐ ) that are already energized are modeled in 5 (8.f) and (8.h) as PQ constant loads during the starting of the motor load ๐ . ๐ ๐,๐,๐๐ท = { ๐ ๐,๐,๐ ( ๐ ๐,๐๐กโ ๐ ๐,๐๐กโ 2 + ๐ ๐,๐๐กโ 2 ) โถ ๐ = ๐๐ ๐,๐0 โถ ๐ โ ๐ โ๐, ๐ โ ๐ ๐ โ๐ (8.e) (8.f) ๐ ๐,๐,๐๐ท = { ๐ ๐,๐,๐ ( ๐ ๐,๐๐กโ ๐ ๐,๐๐กโ 2 + ๐ ๐,๐๐กโ 2 ) โถ ๐ = ๐๐ ๐,๐0 โถ ๐ โ ๐ โ๐, ๐ โ ๐ ๐ โ๐ (8.g) (8.h) where, ๐ ๐,๐๐กโ and ๐ ๐,๐๐กโ represent the Theveninโs equivalent resistance and reactance seen from the terminals of motor ๐ at a given step ๐ , respectively. The value of these Theveninโs equivalent impedances depend on the slip value at each step k . B. AC power flow formulation
Constraints (9.a)-(9.d) represent the second-order cone relaxation of the branch flow model proposed in [26] for each step ๐ in the acceleration period of each motor load ๐ . The aim is to extract the electrical state variables (i.e. voltage, current, and power flow variables). These variables are needed for the optimal control of the converter-interfaced generation units and for checking the transient constraints as will be discussed in sections III.C and III.E, respectively. ๐ ๐,๐,๐ = ๐ ๐,๐,๐ โ 2(๐ ๐๐ . ๐ ๐๐,๐,๐ + ๐ฅ ๐๐ . ๐ ๐๐,๐,๐ ) + ๐น ๐๐,๐,๐ (๐ ๐๐2 + ๐ฅ ๐๐2 ) (9.a) โ๐๐ โ ๐, โ๐ โ ๐ ๐ , โ๐ ๐ ๐๐,๐,๐ = โ ๐ ๐๐ โ ,๐,๐๐ โ โ ๐ +๐ ๐๐ . ๐น ๐๐,๐,๐ + ๐ ๐,๐,๐๐ท โ ๐ ๐,๐,๐๐๐ข๐ โ ๐ ๐,๐,๐๐ท๐บ (9.b) โ๐๐ โ ๐, โ๐ โ ๐ ๐ , โ๐ ๐ ๐๐,๐,๐ = โ ๐ ๐๐ โ ,๐,๐๐ โ โ ๐ + ๐ฅ ๐๐ . ๐น ๐๐,๐,๐ + ๐ ๐,๐,๐๐ท โ ๐ ๐,๐,๐๐๐ข๐ โ ๐ ๐,๐,๐๐ท๐บ (9.c) โ๐๐ โ ๐, โ๐ โ ๐ ๐ , โ๐ ๐,๐,๐ โค ๐ฃ ๐๐๐ฅ2 โ๐ โ ๐, โ๐ โ ๐ ๐ , โ๐ (9.e) ๐น ๐๐,๐,๐ โฅ ๐ ๐๐,๐,๐2 + ๐ ๐๐,๐,๐2 ๐ ๐,๐,๐ โ๐๐ โ ๐, โ๐ โ ๐ ๐ , โ๐ (9.d) Constraint (9.a) expresses the nodal voltage equation as given in [26]. The last term in the right hand side of (9.a) is usually neglected, since it is much smaller than the other terms. Constraints (9.b) and (9.c), respectively, concern with the active and reactive power balances at the extremities of each line. The first term in the right hand side of (9.b) and (9.c), represent, respectively, the sum of active and reactive power flows in lines that are connected to bus ๐ except the line ๐๐ . According to these two equations, the power flow from bus ๐ to bus ๐ is equal to the sum of power flows in other lines connected to bus ๐ plus the net injection power from bus ๐ and the power loss in line ๐๐ . Constraint (9.e) imposes the maximum voltage limit at the nodes of the network . Constraint (9.d) is the relaxed version of the current flow equation in each line according to [26]. This constraint is implemented in the form of the following second order cone constraint. โ 2๐ ๐๐,๐,๐ ๐๐,๐,๐ ๐น ๐๐,๐,๐ โ ๐ ๐,๐,๐ โ โค ๐น ๐๐,๐,๐ + ๐ ๐,๐,๐ โ๐๐ โ ๐, โ๐ โ ๐ ๐ , โ๐ C. Converter control of generation units
In this paper, we consider only dispatchable generation units including DGs, storage systems, and/or static synchronous compensators that are interfaced with the grid via full-bridge power converters. In this paper, the term of DG is used to refer to all these generation units. The re-energization of the DGs after the fault is accounted for in the steady-state analysis. In the optimization problem studied in this paper, it is assumed that the DGs are already connected to the grid and they work in normal state . It means that the DG converter controls the active and reactive power injections (or active power and voltage) following pre-determined set points. When a motor load starts, the voltage at the DG hosting node drops. This voltage drop is detected at the DG hosting node and launches the proposed voltage support control scheme, referred in this paper as current saturation mode . In this mode, instead of the active and reactive powers, the injection current of the DG is controlled. This control mode is achieved through already existing Fault-Right-Through (FRT) strategies in the converter interface while assuming that their current set points can be modified. According to the following formulation, we compute the optimal values of current set points in case of starting each of the motor loads during the restorative period. ๐น๐ ๐,๐๐ท๐บ + ๐น๐ ๐,๐๐ท๐บ = ๐ ๐๐๐ฅ,๐๐ท๐บ 2 (10.a) ๐ ๐,๐,๐๐ท๐บ 2 โค ๐น๐ ๐,๐DG . ๐ ๐,๐,๐ (10.b) ๐ ๐,๐,๐๐ท๐บ 2 โค ๐น๐ ๐,๐DG . ๐ ๐,๐,๐ (10.c) The magnitude of the current injection by the DG converter at node i is forced in (10.a) to be equal to its maximum current limit ( ๐ ๐๐๐ฅ,๐๐ท๐บ ). Constraints (10.b) and (10.c) derive the square of the active and reactive current components, respectively. These constraints are relaxed versions of the original formulations that are equalities instead of inequalities. The aim is to build a convex model of the DG control in current saturation mode. These relaxations will be exact according to the discussion provided in Appendix A. Constraints (10.b) and (10.c) are implemented in the form of second order cone constraints as expressed in the following: โ 2๐ ๐,๐,๐๐ท๐บ ๐น๐ ๐,๐DG โ ๐ ๐,๐,๐ โ โค ๐น๐ ๐,๐DG + ๐ ๐,๐,๐ โ 2๐ ๐,๐,๐๐ท๐บ ๐น๐ ๐,๐DG โ ๐ ๐,๐,๐ โ โค ๐น๐ ๐,๐DG + ๐ ๐,๐,๐ D. Autotransformer tap setting
In some cases, the applied voltage to the induction motor terminal is reduced during the starting period using an autotransformer. This leads to reduce the starting current, power loss and radiated heat during the acceleration period. However, while reducing the voltage terminal, the starting torque will be reduced as well. It causes to lengthen the acceleration period. Therefore, there is a trade-off in setting the tap position of the autotransformer. In what follows, the autotransformer is modeled and incorporated into the optimization problem. The tap position of the autotransformer will be set according to the obtained optimal solution before starting the motor. 6 b) S iโjโ S iโjโ Z p Z s Z p Z s i r ๏ณ๏ c) .ฮr a) 'i U 'j U 'i U 'j U i U j U 'i U 'j U i U j U Fig. 3. Modelling of the autotransformer. a) schematic, b) standard equivalent circuit, c) linearized model.
Fig. 3 shows the equivalent circuit of an autotransformer, where ๐ ๐ and ๐ ๐ denote the impedances at the primary and secondary sides, respectively. ๐ ๐โฒ๐โฒ is the apparent power flowing from the auxiliary node ๐โฒ to the auxiliary node ๐โ . Based on this equivalent circuit, the constraints related to each autotransformer installed at the terminals of motor ๐ are formulated as follows: ๐ ๐,๐,๐ = ๐ ๐,๐,๐ . (1 + ๐ . ฮ๐ ๐ ) (11.a) ๐ ๐,๐,๐ โ ๐ ๐,๐,๐ . (1 + 2 ๐ . ฮ๐ ๐ ) (11.b) where, is the ratio change between two consecutive taps of the autotransformer. Assuming that the turn ratio of the autotransformer is close to 1, the non-linear term in (11 .a ) is linearized in (11 .b ) using the binomial approximation. The product of the integer variable โ๐ ๐ and the positive continuous variable ๐ ๐,๐,๐ makes a non-linear term. This non-linear term is reformulated according to the linearization strategy given in [23]. For this aim, first the integer variable ฮ๐ ๐ should be expressed as the weighted sum of auxiliary binary variables. Then, the product of each of these auxiliary binary variables and the positive continuous variable ๐ ๐,๐,๐ is linearized using the approach given in Appendix B. Finally, Fig. 3.c shows how the derived model in (11 .b ) together with the primary and secondary impedances of the autotransformer are incorporated into the AC-power flow equations. E. Transient Constraints
In this section, transient constraints are expressed regarding the safe starting of the induction motor in a distribution network. ๐ ๐,๐๐๐๐ โฅ ๐ ๐,๐๐๐๐ + ๐พ๐ ๐ (1 โ ๐ ๐,๐ ) (12.a) ๐ฬ ๐,๐๐๐๐ โค ๐ ๐,๐,๐ โ๐ โ ๐ ๐ (12.b) ๐น ๐๐,๐,๐ โค ๐นฬ ๐๐,๐,๐๐๐๐ฅ โ๐๐ โ ๐ ๐ (12.c) In order to avoid the induction motor to stall, (12.a) enforces the electrical torque to be larger than or equal to the load torque. For a given slip ๐ ๐,๐ , the electrical torque of the motor ๐ is obtained using (1) and the mechanical load torque is determined according to the torque-speed curve of the mechanical load. This curve is assumed to be given for a specific load on the shaft. In order to avoid the tripping of the under voltage relays and the over-current relays, (12.b) and (12.c) are added for every step k of the starting of each motor load m. Fig. 4 shows the typical protection curves reported in [5] and [10]. These curves min k U (0.8) (0.917) (0.88) max, ij k F thij F thij F thij F thij F thij F (sec) k t max, ij k F min k U (a) (b) (sec) k t Fig. 4. Typical protection curve of a) over-current and b) under-voltage relays used in this paper represent the values of the under-voltage ( ๐ ๐,๐๐๐๐ ) and over-current ( ๐น ๐๐,๐,๐๐๐๐ฅ ) limits as functions of the acceleration time ( ๐กฬ ๐,๐ ), which was formulated in (3). In order to preserve the linearity in terms of variable ๐กฬ ๐,๐ , ๐ ๐๐๐๐ and ๐น ๐๐,๐๐๐๐ฅ are approximated by ๐ฬ ๐๐๐๐ and ๐นฬ ๐๐,๐๐๐๐ฅ , respectively, according to the piecewise linear approximation method explained in Appendix A. In this regard, for obtaining ๐ฬ ๐๐๐๐ , we add (13)-(16) to the set of constraints, replacing ๐ฅ and ๐(๐ฅ) with ๐กฬ ๐,๐ and ๐ ๐๐๐๐ , respectively. For deriving ๐นฬ ๐๐,๐๐๐๐ฅ , we augment the set of constraints by (13)-(16), replacing ๐ฅ and ๐(๐ฅ) with ๐กฬ ๐,๐ and ๐น ๐๐,๐๐๐๐ฅ , respectively. IV. N UMERICAL R ESULTS
In this section, the functionality of the proposed optimization model is evaluated using a test distribution network shown in Fig. 5. This network is based on a 11.4kV distribution network in Taiwan. The base power and energy values are set to 1MW and 1 MWh, respectively. The details regarding the nodal and line data are given in [27]. Except the motor loads that are indicated in Fig. 5, the rest of loads are assumed as static loads. The nameplate power ratings of the induction motors at nodes {13}, {20, 29}, and {41} are equal to 805, 435, and 80.5 horsepower, respectively. The parameters of the equivalent circuit (see Fig. 2) of these motor loads are adopted from the real data given in [1]. The static loads are assumed to be of constant-impedance type. Therefore, their load-voltage sensitivity coefficients are set to ๐พ ๐ = 2 and ๐พ ๐ =2 . It is assumed that each node in the network shown in Fig. 5 is equipped with a load breaker. The load variation data along time is according to the practical data reported in [27]. In Fig. 5, the critical loads are identified with โ*โ. The priority factor ( ๐ท ๐ ) of these loads is equal to 10 and for the other loads is equal to 1. It is assumed that a fault occurs on the substation 86-84 shown in Fig. 5. The faulted substation is isolated by opening its nearest breakers at both sides. The restorative period is assumed from 17:00 P.M. to 03:00 A.M [28]. As mentioned in section III, in order to obtain the optimal restoration strategy, first, the steady-state analysis is performed according to [23]. The reconfiguration results are shown in Fig. 5. The areas in the off-outage area that are colored with red, green, blue, and yellow are energized through closing tie-switches T1, T2, T4, and T7. The 7 Fault OLTC2OLTC1 B1 F ee d er A B2 F e e d e r B
15 2218 23 B3 Feeder C
25 29 B4 F ee d er D B5 F ee d er E B6 F ee d er F
85 8747 B7 F ee d er G B8 B9 F ee d er H F ee d er I B10 F e e d e r J B11
T3 T4T5T8 T9 T11 F e e d e r K T10
Sectionalizing switchCircuit BreakerTie-switchCritical load T6
630 42
Motor Load M M MM
T1 T7T2 M
20 505864 6783
Fig. 5. The test distribution network under the post-fault configuration. optimal load restoration sequence obtained from the steady-state analysis is depicted on a time axes in Fig. 6.a. Then, the proposed optimization problem in this paper is solved in order to shift the energization instants of certain loads to account for the starting dynamics of the motor loads. In this regard, three sets of simulation scenarios are studied. The algorithm is implemented on a PC with an Intel(R) Xeon(R) CPU and 6 GB RAM; and solved in Matlab/Yalmip environment, using Gurobi solver. Branch-and-Bound method is used to handle the developed mixed-integer optimization problem. The slip step size ( ฮ๐ ) is assumed to 0.05 p.u. A. Simulation scenario I: linear load toque
The first set of simulation scenarios are studied assuming that the mechanical loads of all the motor loads have linear torque-speed characteristics. The magnitude of the nominal mechanical torques (at the synchronous speed) are equal to 0.4, 0.3, and 0.05 p.u. for the loads at nodes {13}, {20, 29}, and {41}. Scenario I.a is defined assuming the critical and sensitive loads are set according to the default conditions shown in Fig. 5. The optimization problem is solved and the resulting optimal decisions are shown in Fig. 6.b. This figure shows the loads whose energization times should be shifted with reference to the results obtained from the steady-state analysis (see Fig. 6.a). The time resolution of this study is chosen to be 1 hour. As it can be seen in Fig. 6, the motor loads are energized in sequential steps with sufficient time intervals such that the starting transients of different motors do not overlap [1]. This shifting of the load energization times causes 29.7 p.u. additional energy not supplied (while considering the priority factors ๐ท ๐ ). This simulation scenario includes a large-scale off-outage area and 4 unsupplied motor loads resulting in 134 binary and 24282 continuous variables. However, the solution for such a large study case is found just in 12.28 second. In order to see the effects of the critical loads on the optimal results, scenario I.b is studied. In this scenario, the loads at nodes {24, 30, 33, 40} are considered as additional critical loads. The optimal results obtained for the load restoration sequence is depicted in Fig. 6.c. Compared to the results of scenario I.a (see Fig. 6.b), it can be seen that the energization times of the new critical loads are not postponed in scenario I.b. The reliability objective value is obtained equal to 29.95 p.u. The computation time is 11.64 second. a) steady state analysis results {41}
17 319 23 224 1 {12,14,30}{24,33} {21}{22}{40}{13}{20}{29} b) Scenario I.a
17 319 23 224 1 {12,14,30}{24,33} {21}{22}{40}{13} {29} {41} {41}
17 22 318 23 224 1 {12,14,18} {21}{22}{13} {20} {29}{27} c) Scenario I.b {41}
17 319 23 224 1 {12,14,30,18}{33} {21}{22}{40}{13}{20}{29} d) Scenario I.c e) Scenario II.a {41}
17 319 23 24 {12,14,30}{33} {21}{22}{40} {13}{20}{29} f) Scenario II.b
22 22
Load restoration sequence:
17 22 3 Time(hour)
A B C
18 19 D E F
24 1 G H A={2,7,8,10,11,12,14,18,23,30}B={4,5,27}C={3,6,24,33}D={9,40} E={13,22} F={21}G={20} H={29,41}Restoration time shifting of the static load at node i from time t1 to time t2 t1 t2{i}{i} Restoration of the motor load at node t at time t t {20} Fig. 6. The optimal load energization sequences.
In the next step, scenario 1.c is defined such that node 18 is added to the protected nodes. The rest of the simulation conditions are the same as in scenario I.a. The results in Fig. 6.d shows that the restoration of the load at node 18 is shifted to the time when all the motor loads are already energized. The reason is that the under voltage limit at node 18 (during the motor starting transients) cannot be respected without shifting the energization of many of other loads. This decision leads to increase the reliability objective value from 29.07 p.u. in scenario I.a to 29.87 p.u. The computation time in scenario I.c is 11.69 second. B. Scenario II: fixed load torque
In scenario II.a, it is assumed that the mechanical load on the shaft of the induction motor at node 20 has a fixed torque-speed characteristic equal to 0.06 p.u. The rest of the simulation conditions are the same as in scenario I.a. According to the results shown in Fig. 6.e, for the safe starting of the motor load at node 20, the same amount of loads should be shifted with respect to the ones in scenario I.a, although the mechanical load power in this scenario is so much less than the one in scenario I.a. The reason is that under a fixed-torque mechanical load, the induction motor can accelerate only if the starting torque (electrical torque at the standstill, s=1) is larger than the mechanical torque. In order to generate this starting toque, the voltage at the motor terminal should be large enough. This is obtained by shifting the energization of loads in the off-outage area, which results in a reliability objective value equal to 29.07 p.u.. The computation time is 16.98 second. In Scenario II.b, the test case of the scenario II.a is studied while assuming an autotransformer at the terminals of the motor load at node 20. It is located in series between node 20 and the induction motor terminals during its acceleration period. Once the motor reaches to 80% of its nominal speed, this autotransformer is taken out from the circuit and the motor is directly connected to the grid [29]. This autotransformer enables 8
Fault B5 F ee d er E B6 Feeder F T8 T9
31 42
M DG
T10 OLTC1
Fig. 7. Part of the test distribution network under the post-fault configuration in simulation scenario III. Table I. The parameters of the induction motor at INESC TEC.
Size (W) ๐ (ฮฉ) ๐ (ฮฉ) ๐ (ฮฉ) ๐ (ฮฉ) ๐ ๐ (ฮฉ) ๐ป(sec) ยฑ20% voltage regulation range in 5 steps ( = 10%, ๐ = 4 ). According to the results shown in Fig. 6.f, with the optimal setting of the autotransformer, there is no need any more to shift the energization time of loads at nodes {20, 24}. The optimal setting of the autotransformer is obtained with the tap position equal to -1. It means that the autotransformer reduces the voltage at motor terminals, which in turn reduces the starting current magnitude of the induction motor. Therefore, the magnitudes of the line voltage drops are reduced. In this regard, this optimal tap setting improves the quality of the restoration solution by increasing the margins of the nodal voltage magnitudes with respect to the transient under-voltage limits. The optimal value of the reliability objective is 25.175 p.u. and the computation time is 14.57 second. It will be validated in the next section that all the transient operational limits are respected with the obtained restoration solutions. Among all the constraints, the voltage magnitude at the accelerating motor node has a very narrow margin with respect to the transient minimum voltage limit. C. Scenario III: DG control in current saturation mode
Scenario III is defined, where a fault occurs on line 84-30. Fig. 7 shows the part of the network that is affected by the post-fault configuration. As it can be seen, the only motor load in this part of the network is at node 41. The parameters of this induction motor are given in Table I. The mechanical load torque on the shaft of the induction motor changes linearly with speed and equals to 2.2 N.m. at the synchronous speed. There is also a DG at node 42, with the ampacity limit equal to 3.65 A. These parameters are according to the parameters of the physical induction motor and DG in the test setup at INESC TEC. The proposed optimization problem is solved in case of this simulation scenario. The optimal load restoration sequence is to energize the static loads at nodes {32,35,36,37} together with the motor load at node {41} at the beginning of the restorative period (at 17:00 P.M). Once the transients of the motor starting disappear (at 17:15 P.M.), we will restore the loads at nodes {31,33,39,40}. The control of the DG converter at node 42 enters to the saturation mode during the motor acceleration period. The optimal active and reactive components of the converter current references (square of ๐น๐ ๐,๐๐ท๐บ and ๐น๐ ๐,๐๐ท๐บ used in (10)) are obtained as 2.106 and 2.970 ampere, respectively. Fig. 8. Simulation results for scenario III. a) slip, b) electrical and mechanical torques, c) motor starting current, d) voltage at node 41. V. F EASIBILITY V ALIDATION R ESULTS
In this section, it is aimed to validate the solution feasibility of the proposed semi-static optimization model using an off-line simulation and a physical test experiment. It will be illustrated that the dynamics of the induction motor starting are represented into the optimization problem with a sufficient degree of accuracy. For this aim, the results obtained from scenario III and given in section IV.C are applied on the distribution network shown in Fig. 7. A. Time-domain simulation results
In this section, an off-line model of the distribution network shown in Fig. 7 is built in Matlab/Simulink. The motor loads and static loads are represented by the Simulink model of the induction motor, and the impedance loads, respectively. The DG is modeled with a controllable dynamic load. First, we apply the obtained optimization results of scenario III on the off-line simulation model and then we start the motor. Fig. 8 shows the simulation results. In this figure, the electrical state profiles obtained from the optimization problem are compared with the ones obtained from the time-domain simulation. As it can be seen, the proposed semi-static optimization model does not represent the initial overshoot transients of the motor inrush current. As mentioned in section II, in deriving the semi-static model, we neglect the
Fig. 9. Simulation results for scenario II.b. a) slip, b) electrical and mechanical torques, c) motor starting current, d) voltage at node 20. Triphase PM15Power Amplifier PMSG Induction MotorResistive Load WorkstationOPAL-RT OP5600Real-Time SimulatorDC/AC InverterDC Power Supply
Fig. 10. The experimental PHIL test setup in the smart grid laboratory at INESC TEC.
Motor Load
Controller BatteryCable
Real-Time SimulatorCurrent SensorPower AmplifierWorkstation Motor VoltageMotor Current Power Hardware
Fig. 11. Block diagram of the PHIL test setup.
DC term in the motor starting current. Disregarding these very fast transients, Fig. 8 shows that the proposed semi-static model represents accurately the behavior of electrical state variables during the motor acceleration at the motor and network sides. In a similar fashion, the optimal results obtained from scenario II.b are tested using a time-domain simulation in Matlab/Simulink. The results at t=01:00, when the motor load at node 20 is started (see Fig. 6.f), are provided in Fig. 9. These results validate the feasibility of the optimal setting that was found in scenario II for the tap position of the auto-transformer. B. Experimental test
As the next step of the validation study, a Power Hardware In the Loop (PHIL) experiment is performed. The test setup is implemented as shown in Fig. 10 in the smart grid laboratory at INESC TEC [30]. The block diagram of the laboratory test setup is depicted in Fig. 11. This PHIL test setup consists of three main parts. I) The first part includes the Matlab/Simulink network model of the whole network shown in Fig. 7 except the nodes 41 and 42. This model is compiled, and then executed by a Real-Time Simulator, namely OPAL-RT OP5600. II) The second part includes the Power Amplifier which magnifies the motor voltage signal received from the Real-Time Simulator. In this test setup TriPhase PM15 is used as the power amplifier with the nominal voltage equals to 400V. It can tolerate up to 30A of peak current per phase. In order to respect this ampacity limit during the motor starting transients, the base power and voltage of the test network are reduced to 320 VA and 150 V, respectively. III) The third part of the PHIL test setup is the physical hardware, including the motor load at node 41 and the DG at node 42 that are connected through a co-axial cable. In order to emulate a mechanical load with a linear torque-speed characteristic, the induction motor is loaded with a permanent magnet synchronous generator connected to a resistive load. The DG is emulated using an AC-DC inverter that is supplied by a DC power supply and controlled in current saturation mode. We start the induction motor and measure the voltages at the motor terminals during the motor acceleration period. As shown in Fig. 12, this waveform is always above the transient under-voltage limit. It verifies that the under-voltage relay installed at the motor terminals will not trip in case of starting the motor load. In Fig. 13, the voltage magnitude measured at the motor
Fig. 12. The voltage magnitude measured at the motor terminals during its acceleration period.
Fig. 13. The voltage at the motor terminal during its acceleration period obtained from the experiment and from the optimization model. terminal during the experiment is compared with the voltage profile obtained from the proposed semi-static optimization model. This figure shows that except at the first instants after the motor energization, the results of the semi-static model are sufficiently accurate with respect to the experimental results. At the first instants following the motor starting, the voltage at the motor terminal experiences a voltage dip that cannot be represented using the semi-static model. This voltage dip is caused by the DC term in the motor inrush current, which is neglected in deriving the semi-static model. As it can be seen in Fig. 13, this transient voltage dip disappears very fast and does not affect the feasibility of the solution with respect to the transient under-voltage limit. Fig. 14 shows the voltage, current, and power measured at the terminals of the DG during the motor acceleration period. The active and reactive components of the DG current are almost controlled to the reference values, which are obtained in section IV.C as 2.106 and 2.97 Ampere, respectively. Fig. 14 validates that the DG is working in current constant mode during the period that the DG experiences the voltage dip at its terminals due to the inrush starting current of the motor. VI. C ONCLUSION
This paper presents a semi-static optimization model for the starting dynamics of induction motors. This model is integrated into the power flow formulation in a convex fashion. The 10
Fig. 14. The measured voltage, current and power of the DG during the motor acceleration period in the PHIL test experiment. resulting optimization problem represents the transient states of an active distribution network during the motor acceleration transients. The current saturation mode of the DG converter control is modeled in this formulation in a convex way. In this regard, the optimal references for the active and reactive components of the injection current are obtained. In addition, the proposed optimization formulation includes a convex model for the optimal tap setting of the autotransformer that is used for the starting of the induction motors. The resulting optimization formulation can be applied in any decision-making problem that is facing the dynamics imposed by the starting of induction motors. This model is used in this paper to solve the optimal load restoration problem in a distribution network. In this problem, the optimal energization sequence of different loads (static and motor loads) are obtained such that the total energy of loads that cannot be supplied will be minimized. The functionality of the proposed optimization model in solving large-scale load restoration problems is evaluated in the case of different simulation scenarios. The accuracy of the proposed semi-static optimization model in providing feasible solutions was verified using off-line time-domain simulations and also using a Power-Hardware-In-the-Loop test experiment. VII. R EFERENCES [1] G. S. Grewal, S. Pocsai, and M. M. Hakim, โTransient motor reacceleration study in an integrated petrochemical facility,โ
IEEE Trans. Ind. Appl. , vol. 35, no. 4, pp. 968โ977, Jul. 1999, doi: 10.1109/28.777207. [2] H. Sekhavatmanesh, R. Cherkaoui, J. Rodrigues, C. L. Moreira, and J. A. P. Lopes, โA convex model for induction motor starting transients imbedded in an OPF-based optimization problem,โ to be presented at the PSCC2020, Porto, Portugal, p. 7. [3] I. Beil, A. Allen, A. Tokombayev, and M. Hack, โConsiderations when using utility-scale battery storage to black start a gas turbine generator,โ in , 2017, pp. 1โ5, doi: 10.1109/PESGM.2017.8274529. [4] Z. Zhao and B. Ooi, โFeasibility of fast restoration of power systems by micro-grids,โ
Transm. Distrib. IET Gener. , vol. 12, no. 1, pp. 126โ132, 2018, doi: 10.1049/iet-gtd.2017.0323. [5] H. Qu and Y. Liu, โMaximum restorable load for substation during power system restoration,โ in , 2009, pp. 1โ4, doi: 10.1109/SUPERGEN.2009.5348350. [6] C. L. Moreira, F. O. Resende, and J. A. P. Lopes, โUsing Low Voltage MicroGrids for Service Restoration,โ
IEEE Trans. Power Syst. , vol. 22, no. 1, pp. 395โ403, Feb. 2007, doi: 10.1109/TPWRS.2006.888989. [7] D. Wu, H. Wu, and H. Dongt, โInfluence of Induction Motor Starting on Microgrid,โ in , 2018, pp. 376โ381, doi: 10.1109/APPEEC.2018.8566305. [8] F. K. Tuffner, K. P. Schneider, J. Hansen, and M. A. Elizondo, โModeling Load Dynamics to Support Resiliency-Based Operations in Low-Inertia Microgrids,โ
IEEE Trans. Smart Grid , vol. 10, no. 3, pp. 2726โ2737, May 2019, doi: 10.1109/TSG.2018.2809452. [9] H. Zheng and C. L. DeMarco, โA New Dynamic Performance Model of Motor Stalling and FIDVR for Smart Grid Monitoring/Planning,โ
IEEE Trans. Smart Grid , vol. 7, no. 4, pp. 1989โ1996, Jul. 2016, doi: 10.1109/TSG.2016.2548082. [10] J. C. Gomez, M. M. Morcos, C. A. Reineri, and G. N. Campetelli, โBehavior of induction motor due to voltage sags and short interruptions,โ
IEEE Trans. Power Deliv. , vol. 17, no. 2, pp. 434โ440, Apr. 2002, doi: 10.1109/61.997914. [11] J. Bredthauer and N. Struck, โStarting of large medium voltage motors: design, protection, and safety aspects,โ
IEEE Trans. Ind. Appl. , vol. 31, no. 5, pp. 1167โ1176, Sep. 1995, doi: 10.1109/28.464534. [12] S. Hasan, K. M. Muttaqi, R. Bhattarai, and S. Kamalasadan, โA Coordinated Control Approach for Mitigation of Motor Starting Voltage Dip in Distribution Feeders,โ in , 2018, pp. 1โ6, doi: 10.1109/IAS.2018.8544554. [13] F. G. Limbong, โThe use of neural network (NN) to predict voltage drop during starting of medium voltage induction motor,โ in , 2016, pp. 156โ160, doi: 10.1109/ICITACEE.2016.7892428. [14] H. Qu and Y. Liu, โMaximizing restorable load amount for specific substation during system restoration,โ
Int. J. Electr. Power Energy Syst. , vol. 43, no. 1, pp. 1213โ1220, Dec. 2012, doi: 10.1016/j.ijepes.2012.05.049. [15] M. Falahi, K. L. Butler-Purry, and M. Ehsani, โInduction Motor Starting in Islanded Microgrids,โ
IEEE Trans. Smart Grid , vol. 4, no. 3, pp. 1323โ1331, Sep. 2013, doi: 10.1109/TSG.2013.2271261. 11 [16] P. Rodriguez, A. V. Timbus, R. Teodorescu, M. Liserre, and F. Blaabjerg, โFlexible Active Power Control of Distributed Power Generation Systems During Grid Faults,โ
IEEE Trans. Ind. Electron. , vol. 54, no. 5, pp. 2583โ2592, Oct. 2007, doi: 10.1109/TIE.2007.899914. [17] L. Hadjidemetriou, P. Demetriou, and E. Kyriakides, โInvestigation of different Fault Ride Through strategies for renewable energy sources,โ in , 2015, pp. 1โ6, doi: 10.1109/PTC.2015.7232594. [18] Y. Yang, F. Blaabjerg, and H. Wang, โLow-Voltage Ride-Through of Single-Phase Transformerless Photovoltaic Inverters,โ
IEEE Trans. Ind. Appl. , vol. 50, no. 3, pp. 1942โ1952, May 2014, doi: 10.1109/TIA.2013.2282966. [19] J. Miret, A. Camacho, M. Castilla, L. G. de Vicuรฑa, and J. Matas, โControl Scheme With Voltage Support Capability for Distributed Generation Inverters Under Voltage Sags,โ
IEEE Trans. Power Electron. , vol. 28, no. 11, pp. 5252โ5262, Nov. 2013, doi: 10.1109/TPEL.2013.2246190. [20] A. Camacho, M. Castilla, J. Miret, L. G. de Vicuรฑa, and R. Guzman, โPositive and Negative Sequence Control Strategies to Maximize the Voltage Support in ResistiveโInductive Grids During Grid Faults,โ
IEEE Trans. Power Electron. , vol. 33, no. 6, pp. 5362โ5373, Jun. 2018, doi: 10.1109/TPEL.2017.2732452. [21] A. Camacho, M. Castilla, J. Miret, L. G. de Vicuรฑa, and G. L. M. Andrรฉs, โControl Strategy for Distribution Generation Inverters to Maximize the Voltage Support in the Lowest Phase During Voltage Sags,โ
IEEE Trans. Ind. Electron. , vol. 65, no. 3, pp. 2346โ2355, Mar. 2018, doi: 10.1109/TIE.2017.2736486. [22] M. Tooley,
Plant and Process Engineering 360 . Elsevier, 2009. [23] H. Sekhavatmanesh and R. Cherkaoui, โAnalytical Approach for Active Distribution Network Restoration Including Optimal Voltage Regulation,โ
IEEE Trans. Power Syst. , pp. 1โ1, 2018, doi: 10.1109/TPWRS.2018.2889241. [24] L. Gan, N. Li, U. Topcu, and S. H. Low, โExact Convex Relaxation of Optimal Power Flow in Radial Networks,โ
IEEE Trans. Autom. Control , vol. 60, no. 1, pp. 72โ87, Jan. 2015, doi: 10.1109/TAC.2014.2332712. [25] H. Sekhavatmanesh and R. Cherkaoui, โOptimal Infrastructure Planning of Active Distribution Networks Complying with Service Restoration Requirements,โ
IEEE Trans. Smart Grid , vol. PP, no. 99, pp. 1โ1, 2017, doi: 10.1109/TSG.2017.2716192. [26] J. A. Taylor and F. S. Hover, โConvex Models of Distribution System Reconfiguration,โ
IEEE Trans. Power Syst. , vol. 27, no. 3, pp. 1407โ1413, Aug. 2012, doi: 10.1109/TPWRS.2012.2184307. [27] C.-T. Su, C.-F. Chang, and J.-P. Chiou, โDistribution network reconfiguration for loss reduction by ant colony search algorithm,โ
Electr. Power Syst. Res. , vol. 75, no. 2, pp. 190โ199, Aug. 2005, doi: 10.1016/j.epsr.2005.03.002. [28] S. ฤurฤiฤ, C. S. รzveren, L. Crowe, and P. K. L. Lo, โElectric power distribution network restoration: a survey of papers and a review of the restoration problem,โ
Electr. Power Syst. Res. , vol. 35, no. 2, pp. 73โ86, Nov. 1995, doi: 10.1016/0378-7796(95)00991-4. [29] R. F. McElveen and M. K. Toney, โStarting high-inertia loads,โ
IEEE Trans. Ind. Appl. , vol. 37, no. 1, pp. 137โ144, Jan. 2001, doi: 10.1109/28.903136. [30] C. Gouveia et al. , โExperimental validation of smart distribution grids: Development of a microgrid and electric mobility laboratory,โ
Int. J. Electr. Power Energy Syst. , vol. 78, pp. 765โ775, Jun. 2016, doi: 10.1016/j.ijepes.2015.12.005. [31] H. P. Williams,
Model Building in Mathematical Programming . John Wiley & Sons, 2013. [32] E. M. L. Beale and J. Tomlin, โSpecial facilities in a general mathematical programming system for nonconvex problems using ordered sets of variables,โ
Oper. Res. , vol. 69, pp. 447โ454, Jan. 1969. VIII. A PPENDICES A. Piecewise Linear Approximation
Consider f(x) as a continuous function with the domain of [ ๐ฅ , ๐ฅ ๐ ]. The concept of the piece-wise linear approximation introduced in [31] is shown in Fig. 15, where ๐ฬ(๐ฅ) denotes an approximation function for f(x) . We divide the domain of function ๐ by n break-points ๐ฅ , ๐ฅ , โฆ , ๐ฅ ๐ . Between each two successive breaking points ๐ฅ ๐ and ๐ฅ ๐+1 , the function f is approximated with a straight line connecting the points ๐ฅ ๐ and ๐ฅ ๐+1 . In order to formulate this approximation method, we introduce n non-negative auxiliary variables ๐ ๐ under (13) - (16) . Assume a given point x between two breaking point ๐ฅ ๐ and ๐ฅ ๐+1 . The value of x is expressed in (13) in terms of the weighted sum of breaking points ๐ฅ ๐ and ๐ฅ ๐+1 , with ๐ ๐ and ๐ ๐+1 as the weighting coefficients. Using the same auxiliary variables, function ๐ฬ(๐ฅ) is formulated as in (14) . The sum of all auxiliary variables ๐ ๐ should be one (15) . Therefore, ๐ฬ(๐ฅ) consists of straight lines between successive breakpoints. For the given point ๐ฅ โ [๐ฅ ๐ , ๐ฅ ๐+1 ] , all the auxiliary variables, except ๐ ๐ and ๐ ๐+1 , will be zero. Because as given in (16) , variables ๐ are forced to respect the Special Ordered Set-2 (SOS2) constraint [32]. According to this constraint, out of all ๐ ๐ variables, at most two successive variables can be non-zero. Most of the commercial solvers have the feature to account for these types of constraints during the Branch-and-Bound search algorithm. ๐ฅ = โ ๐ ๐ ๐ฅ ๐๐๐=1 )13( ๐ฬ(๐ฅ) = โ ๐ ๐ ๐(๐ฅ ๐ ) ๐๐=1 )14( โ ๐ ๐๐๐=1 = 1 )15( ๐ ๐ โฅ 0 โถ ๐๐๐2 )16( We can increase the accuracy of this approximation using more number of breaking points ( n ). However, it leads to more number of auxiliary variables, which could increase the computation burden of the optimization problem. 12 (x) f x x x x n x (x ) f ( ) f x (x ) f (x ) n f (x) f Fig. 15. Piece-wise linear approximation of an arbitrary continuous function ๐(๐ฅ) B. Elimination of product of variables
In this section, a method is provided for the linearization of the constraints which incorporates a product of two variables. The product of two variables ๐ฅ and ๐ฅ can be replaced by one new variable y, subject to the constraints given in (17) and (18) . The proof of this linearization is provided in [31]. It is assumed that ๐ฅ is a binary variable and ๐ฅ is a positive continuous variable, for which โค ๐ข holds. (17) ๐ฅ โ ๐ข(1 โ ๐ฅ ) โค ๐ฆ โค ๐ฅ (18) A. Exactness of the relaxation used in the modeling of the DG control in the current saturation mode.
According to [24], it can be proved that the relaxation made in (10.b) and (10.c) will be exact under the following conditions: I)
The maximum voltage limit at the DG hosting node is not binding during the acceleration period of the motor load. II)
The line ampacity limits are not binding during the acceleration period of the induction motor. III)
The objective function strictly decreases with the square of active and reactive power injections of the DG ( ๐ ๐,๐,๐๐ท๐บ 2 and ๐ ๐,๐,๐๐ท๐บ 2 ). Condition I holds, since the voltage at the DG node is already dropped due to the starting of the induction motor. Therefore, there is a large margin for the voltage at the DG node with respect to its maximum limit. Condition II is usually ensured during the planning phase of the distribution networks. For example, in industrial networks, the line ampacity limits are sufficiently larger than the starting current of induction motors. Under these two conditions, the optimal value of the objective function ( ๐น ๐๐ in (4)) will be bounded only by the minimum voltage constraint. In this regard, during the starting of an induction motor, the minimum voltage limit at the hosting node of the starting motor is the bottleneck constraint. Therefore, in order to restore more loads (to decrease the value of objective function ๐น ๐๐ ), the voltage at the motor node should increase. From the other hand, with the starting of the motor load (increasing of power absorption at the motor node), the voltage at the DG terminal decreases. It means that the sensitivity of the voltage at the DG node ( ๐ ๐ท๐บ ) is positive with respect to the power injection at the motor node ( ๐ ๐ : negative of power absorption). Due to the symmetry of the grid admittance matrix, we can intuitively infer that the sensitivity of the voltage at the motor node ( ๐ ๐ ) with respect to the power injection at the DG node ( ๐ ๐ท๐บ ) is also positive. The mathematical proof of this expression is beyond the scope of this paper. Therefore, if the DG power injection increases, the voltage at the motor load increases and therefore the total unrestored energy of loads ( ๐น ๐๐๐๐