WAMS-Supported Power Mismatch Optimization for Secure Intentional Islanding

: It is expected that a coordinated operation of several system integrity protection schemes will become a necessity in the future. This research represents an innovative strategy for coordinating under-frequency load shedding and intentional controlled islanding schemes for improving electric power system stability and resilience. In the great majority of real-world cases, both approaches follow conventional tactics, i.e., disconnecting a ﬁxed number of feeders at predeﬁned frequency thresholds and isolating a predeﬁned area of a power system regardless of the actual conditions. Under the newly arisen network conditions in which weather-dependent distributed energy sources introduce a signiﬁcant level of intermittency, conventional approaches need to be upgraded in order to retain a high level of power system operation security. In this paper, a mixed-integer linear programming approach is used to adjust the island size, including/excluding additional substations according to the available amount of generation in the region. The ﬁne-tuning of the power rebalancing is achieved by potentially blocking selected load shedding stages. This minimizes the power imbalance of the newly formed islands, which helps to reduce the number of partial or even total blackouts and also accelerates the power system’s restoration process. The optimization approach was tested on a generic IEEE 39-bus network and shows promising results along with the capability of coping with real-world applications using wide-area monitoring systems as a source of real-time measurements. The results also indicated the importance of appropriate load modelling since both voltage and frequency dependence are recognized to have a signiﬁcant effect on intentional controlled islanding.


Introduction
By definition, an electric power system (EPS)'s operation is deemed stable when all state variables (such as frequency, rate of change of frequency, voltage magnitudes, etc.) are within the permissible bounds and steady. When these conditions are not met, corrective actions need to be taken to avoid cascade tripping that may result in an unsupervised EPS separation or even a widespread blackout. These actions are referred to as system integrity protection schemes (SIPSs) whose task is to eliminate the risk of severe socio-economic damage caused by a blackout. The most commonly applied SIPSs are under-frequency and under-voltage load shedding (UFLS and UVLS, respectively) schemes [1] that often suffice for EPS stabilization [2]. However, this is not always the case [3,4].
To set up an additional protection barrier, a so-called intentional controlled islanding (ICI) scheme has been presented and researched recently [5][6][7][8][9]. ICI is a corrective measure scheduled for activation only when the remaining set of SIPSs is fully exhausted and not capable of stabilizing the system on its own. Its aim is to split the EPS into multiple stable and self-sustainable islands. To this end, it is of the utmost importance to ensure a sufficient level of active power balance within the island.
In this respect, it has been shown in the existing literature that coordination between ICI and other SIPSs, such as UFLS, needs to be researched [10]. Namely, it turned out that an UFLS scheme operating by means of frequency and rate-of-change-of-frequency (RoCoF) relay tripping is not always able to avoid frequency instability.
Constructional limitations in steam and gas turbines dictate avoidance from operating below 47.5 Hz (in 50 Hz systems) and 57.0 Hz (in 60 Hz systems) [11]. The reason is that the life expectancy of the turbines below these thresholds starts to significantly decrease and they are exposed to potential physical damage due to the excitation of the resonant frequencies of the turbine blades. To avoid such undesired consequences, power plants using steam and gas technology are equipped with under-frequency protection. However, the malfunction or poor parameterization of under-frequency protection can also be the cause of a blackout. A good example to support this claim is the 2003 Italy blackout [4].
The most commonly used criterion by an ICI scheme for the determination of island boundaries is the minimization of a power exchange with the rest of the system. This approach is based on the real-time monitoring of the transmission line power flows [7,12]. In fact, this is very similar to minimizing the active power imbalance within the island ( [8,13]) while other approaches focus on synchronous generator angle instability [12]. In several cases, the coherence criterion was used to determine groups of synchronous generators that do not experience significant electromechanical oscillations against each other [14][15][16] (i.e., finding coherent groups of generators and include them in the same island). The two most evident drawbacks of the latter approach are (i) the dependence of calculation results on the pre-incident equilibrium point and (ii) not considering topology changes. As an improvement of this approach, spectral clustering is often used in combination with the advances of the slow coherency approach [9]. Nevertheless, some of the drawbacks remain, such as ignoring the possibility of transient stability occurrence ( [5,17]). Hence, researchers in [6] built multi-stage constrained spectral clustering, taking into account the simultaneous effect of frequency as well as the active and reactive powers. In [18], authors assume inter-area oscillations as a contingency. Therefore, it can be concluded that the existing literature reveals versatile ICI approaches.
The questions that need to be answered by an ICI concept are where and when the islanding should take place. In this paper, the focus is given to the question of where the most appropriate boundaries of the island are located. On the other hand, islanding will take place once the frequency reaches a predefined threshold and is, therefore, not considered as part of the optimization procedure. For this purpose, we selected a frequency threshold of 57.6 Hz (an equivalent for a 50 Hz system would be 47.9 Hz [19]).
With regard to the cause of critical EPS operation that requires ICI activation, authors in most of the existing publications assume a short-circuit fault whose consequences challenge the angle stability of synchronous generators. This type of contingency was recognized as an initiating event in 44% of the recorded blackouts [20]. This is why searching for coherent groups of generators makes sense. However, a power plant outage could also represent a plausible event with a low probability on the one hand, yet a high impact on the other. To the best of the authors' knowledge, only a few research papers considered a power plant outage as an initial event, which is probably due to the fact that it was detected as an initiating event in only 11% of the cases [20]. In the course of EPS deterioration, generator tripping was indeed one of the consequences, but the operational issues began with a short-circuit fault [10,21].
Only a few researchers have discussed the coordination of UFLS and ICI [10]. Authors in [10] emphasize the need for further research on this topic and claim that this is an open question for further research. The potential of including other SIPS in ICI is also emphasized because voltage instability often occurs prior to transient stability and thus needs to be researched urgently ( [10,22]). Furthermore, by foreseeing possible dynamic stability issues during the restoration process, the existing ICI could be further improved [23,24]. For this purpose, we used mixed-integer linear programming (MILP), which is widely used in various energy management system applications. Its simplicity and efficiency make it applicable in areas such as ICI [10,21], post-disaster utility reaction [25], peak operation of gas-fired generating units [26], using renewable energy sources (RES) for frequency Energies 2021, 14, 2790 3 of 13 stability in an islanded network following a major fault [27] or even minimizing CO 2 emissions from the steelwork industry [28].
The novelty of this paper is the establishment of efficient coordination between UFLS and ICI schemes with the aim of forming stable islands. In addition, further actions are not required after island separation (such as load shedding or generator tripping), which differs from the existing literature according to the authors' best knowledge. The methodology is independent of EPS configuration and adaptable to different operating conditions. This paper is organized as follows: In Section 2, the problem to be solved is described alongside the network used for the case study. The methodology selected and adapted for an ICI based on MILP is discussed in Section 3. In Section 4, the results obtained from the proposed strategy applied on an IEEE 39-bus test system are given, alongside the discussion about the results and our future work. The conclusions are finally drawn in Section 5.

Problem Description
The current trends in EPS development indicate an increasing level of penetration of RES both in the distribution network in terms of distributed generation (DG) as well as the transmission system mostly in terms of wind parks. Many of these belong to the Converter-Interfaced Generation (CIG) type of units not contributing inherently to system inertia. Consequently, this trend results in the decrease and time volatility of EPS inertia masses since conventional power plants based on synchronous machine technology (coal and nuclear) are scheduled for closure, and RES availability significantly depends on weather conditions. In this way, spinning reserves are difficult to predict as are frequency excursions after power imbalance events [29], both in terms of frequency nadir as well as RoCoF. This means that the existing solutions for assuring EPS frequency stability need to be updated in order to adjust to volatile conditions in real time.
We consider that the coordination of multiple SIPS will be crucial in the future for performing a set of ultimate remedial actions involving load curtailment (UFLS) and system splitting (ICI). In the upcoming subsections, all relevant information about such a coordination between SIPS is provided, including the modelling of the electrical network and its elements, UFLS and a wide-area monitoring system (WAMS) as a provider of real-time system measurements.

Network
In this paper, the authors decided to employ one of the most commonly used test networks for such purposes, an IEEE 39-bus test system also known as the New England test system. It consists of 10 synchronous generators, 34 transmission lines (modelled with π sections), 12 transformers and 19 equivalent loads. The synchronous machines are modelled with a sixth-order representation along with respective governors and automatic voltage controllers. As most of the machines in this specific network are part of thermal power plants (five nuclear and three coal), the IEEEG1 governor model is used for all of them. For the sole hydro power plant (denoted as G10 in Figure 1), the IEEEG3 governor model is assumed. As for the excitation controller, all machines are equipped with an IEEET1 excitation model. Since power oscillation damping is not the focus of our work, power system stabilizers are not activated. The synchronous machine G1 represents the remaining part of the network, to which New England is interconnected. More information about the test system parameters can be accessed in [30].

UFLS Protection
Under-frequency relays, performing the task of UFLS protection, are assumed to be located in all substations that include an equivalent load (see Figure 1). The chosen UFLS setting follows the conventional approach with six stages represented by frequency thresholds and the amount of curtailed load in each of them. We made the decision to deliberately not consider any intentional time delays since it has been proven many times to worsen the EPS frequency response and the selectivity of UFLS [10,21]. Namely, intentional time delays may cause lower frequency nadir due to a slower relay trip signal and the overlapping of (otherwise consecutive) stages. The settings, summarized in Table 1, comply with the latest ENTSO-E regulations [31,32] and have been proven many times as robust. In publications [10,21], the authors regarded a 50% de-loading as the maximum UFLS reach, yet we decided to lower this value to 45% according to the recommendations for the continental European interconnection [31,32]. Since UFLS is to safeguard the overall system integrity, it is up to the transmission system operators (TSOs) to request its general requirements, whereas a tight collaboration with distribution system operators (DSOs) is needed to set up a list of specific medium voltage feeders that are to be included in individual UFLS stages. Usually, a selection of feeders is made based on their yearly or few-month average power consumption. However, on a shorter time scale, power fluctuations might significantly deviate from average values and, therefore, UFLS activation might not always yield the expected results. This is why we consider that real-time power measurements have the potential to improve this setback. This is definitely a requirement for a successful and stable transition into island operation.
One of the most important aspects in EPS modelling for UFLS and islanding studies is to carefully select an appropriate load model. In general, power consumption depends Based on the technical reference [30], the rated current of an individual 345-kV transmission line is 1000 A. In terms of active power, this means approximately 600 MW of transmission capacity. G1 located at bus 39 is connected to two neighboring buses (i.e., bus 1 and bus 9) by means of double circuit lines, which corresponds to a doubled transmission capacity of 1200 MW.

UFLS Protection
Under-frequency relays, performing the task of UFLS protection, are assumed to be located in all substations that include an equivalent load (see Figure 1). The chosen UFLS setting follows the conventional approach with six stages represented by frequency thresholds and the amount of curtailed load in each of them. We made the decision to deliberately not consider any intentional time delays since it has been proven many times to worsen the EPS frequency response and the selectivity of UFLS [10,21]. Namely, intentional time delays may cause lower frequency nadir due to a slower relay trip signal and the overlapping of (otherwise consecutive) stages. The settings, summarized in Table 1, comply with the latest ENTSO-E regulations [31,32] and have been proven many times as robust. In publications [10,21], the authors regarded a 50% de-loading as the maximum UFLS reach, yet we decided to lower this value to 45% according to the recommendations for the continental European interconnection [31,32]. Since UFLS is to safeguard the overall system integrity, it is up to the transmission system operators (TSOs) to request its general requirements, whereas a tight collaboration with distribution system operators (DSOs) is needed to set up a list of specific medium voltage feeders that are to be included in individual UFLS stages. Usually, a selection of feeders is made based on their yearly or few-month average power consumption. However, on a shorter time scale, power fluctuations might significantly deviate from average values and, therefore, UFLS activation might not always yield the expected results. This is why we consider that real-time power measurements have the potential to improve this setback. This is definitely a requirement for a successful and stable transition into island operation. One of the most important aspects in EPS modelling for UFLS and islanding studies is to carefully select an appropriate load model. In general, power consumption depends on both the voltage and the frequency. With regard to voltage, the most typical load model representatives are constant impedance (Z), constant current (I), constant power (P) or the exponential model. The latter enables the inclusion of frequency dependence as well. In this way, load modelling plays an important role in ameliorating the consequences of a power imbalance event [33]. It is up to the selected load parameters to determine to what extent this effect is present. One has to approach this task with care since an inappropriate setting could have a direct impact on whether a certain UFLS stage (or even ICI) is activated or not.
There are a lot of differences between the characteristics of a residential, industrial or any other type of load. Additionally, these characteristics are subject to time variation, e.g., seasonal, daily and/or hourly. The most common and, at the same time, the most basic equation representing load voltage and frequency dependence is as follows: in which P(t), Q(t), U(t) and ∆f(t) are actual active power, reactive power, voltage magnitude and frequency deviation at the load bus in relation to time, respectively. P 0 , Q 0 , U 0 and f 0 are their corresponding values prior to the imbalance or, in other words, their steadystate values. The four parameters k pu , k qu , k pf and k qf determine how the load will react to variations in frequency and voltage. In this paper, multiple sets of these four coefficients will be used in order to verify EPS frequency response to a wide range of possible structures of consumers represented by an equivalent load. The selected sets of coefficients will include a constant power model, a constant impedance model with neglected frequency dependence, and a combination of voltage and frequency dependence based on the coefficients taken from [34]. For the residential load type, mean values for k pf and k pf are used. In this paper, feeders controlled by UFLS frequency relays were considered an asset to ICI. Namely, the island power mismatch can be eliminated much more efficiently and accurately by having a possibility to include or exclude certain feeders from the last (or last two in certain cases) ULFS stages.

Phasor Measurement Units and Wide-Area Monitoring System
Phasor measurement units (PMUs) are handy devices to measure and report voltage and current synchrophasor measurements, frequency and RoCoF with a report rate in the range between 10 and 50/60 Hz, depending on the EPS nominal frequency [35][36][37]. A WAMS system integrates data streams provided by an arbitrary number of PMU devices and enables system operators to have almost real-time insight into the EPS operating conditions (during an ENTSO-E split in January 2021, the report rate was determined as 20 Hz [2]). For the purpose of this paper, PMUs need to be located in (i) all buses that might be subject to islanding and (ii) generator buses. One has to be aware that the number of PMUs could be optimized to provide sufficient system observability, yet this was not the focus of our research at this stage.

Methodology Intentional Controlled Islanding
An EPS can be represented as an undirected graph model G = (V, E), where V denotes the system buses {v 1 , v 2 , . . . , v n }, while E denotes the edges e i,j ∈ E, (i, j = 1, . . . , n). The latter is associated with the weight w i,j , representing the active power flow on the edge e i,j . A cut-set (2) removes a set of edges with the aim of splitting G into K islands (in this case, two islands V h and V k ) as described in: In order to illustrate how the graph theory can be applied to the EPS, an illustrative (fictional) example is shown in Figure 2. The exemplary EPS consists of two generating buses (nodes denoted as 1 and 8) and two load buses (denoted as 5 and 7). Each power line or transformer is represented as a graph edge between two nodes. In order to compute the minimum cut, power lines 3-6 and 4-6 need to be disconnected, which is represented by a red dashed curve. With that action, we impose a certain power mismatch in each of the two formed islands, i.e., a power surplus of ∆P = 2 pu in the first and a power deficit of ∆P = −2 pu in the second. Power rebalancing would, therefore, require tripping some generation in the first island and load shedding in the second. In order to illustrate how the graph theory can be applied to the EPS, an illustrative (fictional) example is shown in Figure 2. The exemplary EPS consists of two generating buses (nodes denoted as 1 and 8) and two load buses (denoted as 5 and 7). Each power line or transformer is represented as a graph edge between two nodes. In order to compute the minimum cut, power lines 3-6 and 4-6 need to be disconnected, which is represented by a red dashed curve. With that action, we impose a certain power mismatch in each of the two formed islands, i.e., a power surplus of ∆P = 2 pu in the first and a power deficit of ∆P = −2 pu in the second. Power rebalancing would, therefore, require tripping some generation in the first island and load shedding in the second. The optimization approach used in this paper is based on linear programming, more precisely mixed-integer linear programming (MILP). For the sake of simplicity, a DC power flow is used within the optimization for considering both the transmission line's and the transformer's active power limits. With that, we accept that losses are neglected but, at the same time, we gain some significant calculation speed. Our assumption was that losses do not have an immense impact on the overall power balancing problem, especially having in mind the uncertainties of actual load characteristics. In addition, with loads introducing quite significant uncertainties, we believe ignoring losses does not deteriorate the situation to an unacceptable level, especially when they are estimated to approximately 1-2% for an IEEE 39-bus test system.
In seeking the best possible power balance, we decided to assume that a governor reaction to frequency deviations is too slow to have a noticeable effect [38]. ICI is expected to activate in case UFLS is not able to rebalance the system by itself. Such conditions are without a doubt accompanied by extreme RoCoF values, which means that the frequency can deteriorate to a very low value in a very short period of time.
We designed the islanding concept to run in two subsequent steps. In the first step, the generated and consumed powers within the island are roughly matched. Seeking a minimum power flow disruption (PFD) is used for the splitting criteria as follows: To perform the criteria (3) in the undirected graph G = (V, E) with edge weights , , the partitioning cost could be described by: The optimization approach used in this paper is based on linear programming, more precisely mixed-integer linear programming (MILP). For the sake of simplicity, a DC power flow is used within the optimization for considering both the transmission line's and the transformer's active power limits. With that, we accept that losses are neglected but, at the same time, we gain some significant calculation speed. Our assumption was that losses do not have an immense impact on the overall power balancing problem, especially having in mind the uncertainties of actual load characteristics. In addition, with loads introducing quite significant uncertainties, we believe ignoring losses does not deteriorate the situation to an unacceptable level, especially when they are estimated to approximately 1-2% for an IEEE 39-bus test system.
In seeking the best possible power balance, we decided to assume that a governor reaction to frequency deviations is too slow to have a noticeable effect [38]. ICI is expected to activate in case UFLS is not able to rebalance the system by itself. Such conditions are without a doubt accompanied by extreme RoCoF values, which means that the frequency can deteriorate to a very low value in a very short period of time.
We designed the islanding concept to run in two subsequent steps. In the first step, the generated and consumed powers within the island are roughly matched. Seeking a minimum power flow disruption (PFD) is used for the splitting criteria as follows: To perform the criteria (3) in the undirected graph G = (V, E) with edge weights w i,j , the partitioning cost could be described by: where z i,j is a binary variable determining whether an edge is included in an island or not. By (5), we wanted to avoid single-generator islands in case it generates the lowest possible power by means of a technological minimum. In this step, the algorithm determines to which island the individual buses will be connected. This is the so-called "island-formation" step which, of course, results in one island having a slight surplus of power and the other a deficit of power. In the second step, the remaining power mismatch will be eliminated by ordering either the blocking or the forced activation of certain feeders in the UFLS stages. In this way, we aim to maximize the supplied EPS load in both islands while fulfilling the frequency stability related constraints: where P d,i represents the consumption in bus i (index i is selected among load buses in the selected island, either V h or V k ), P g,i represents the generation at bus i (index i is selected among the generator buses in the selected island, either V h or V k ) and P start d,i represents the consumption at bus i before the main contingency occurs. Factors k min and k max depend on the sign of the power imbalance: in case of a power surplus k min = 0.55 and k max = 1 and in case of a power deficit k min = 0 and k max = 0.55. The value of 0.55 originates from having 55% of the system load in operation after UFLS is fully exhausted (see Table 1). We assumed that each feeder represents a 5% share of the total EPS load and is equipped with a frequency relay that constantly monitors the signal for activation/deactivation. This makes it possible to further decrease the loading in 5% steps.
As an additional constraint, the transmission capacity of the transformers and transmission lines is taken into consideration: where δ i and δ j denote voltage phase angles at buses i and j, and P min i,j and P max i,j denote the minimum and maximum allowable power transfer over the element between buses i and j with a reactance X i,j . The DC power flow is used to check whether these limits are violated in the upcoming island formation. If the DC power flow indicates that the expected islands will not be sustainable, a selected cut-set is denoted as inappropriate.
According to (4)-(9), the calculation begins at the moment a contingency is detected by a WAMS system (see Figure 3). Therefore, the frequency relay settings and cut-sets limits are determined before the frequency reaches critical values that require ICI activation. For extreme RoCoF values, such as the one that appeared in Australia in 2016 (6 Hz/s) [3], one could assume several other robust solutions, e.g., either running the algorithm periodically when the EPS is still in the steady-state operating conditions, or considering RoCoF as an additional optimization parameter that would provide an indication on when to initiate the islanding. the stability of the islands. If there are none, then commands are sent to the load frequency relays carrying new settings. The entire process is performed prior to the first UFLS stage activation. Afterwards, the Python script waits for UFLS or ICI activation.
The approach from Figure 3 was tested using the PowerFactory DIgSILENT 2021 SP 1 software for RMS (root mean square) dynamic simulations in conjunction with Python. Namely, PowerFactory offers a feature of running an arbitrary Python code in each time step of the simulation.

Case Study
In order to invoke extreme frequency fluctuations, the tripping of multiple synchronous generators is the most obvious option (as exemplified in the IEEE 118 test system from [10,21] where two generators were tripped to reach and fully exhaust UFLS and eventually activate ICI). However, in this paper we decided to trip G1 instead, which represents an equivalent of the remaining network. It should be noted that generation tripping does not affect the network topology.
Furthermore, multiple scenarios were observed during the operation of all machines due to high system loading. In those conditions, we examined the effect of having different load models, ranging from constant P, constant Z, and down to the exponential model with/without frequency dependence.  Figure 3. Overall algorithm structure of the proposed approach.
In each time step of the simulation, the measurements from PowerFactory, such as bus frequency and load power, are obtained from the Python script, which is used for making a decision whether any actions are required. When a contingency is detected (i.e., the output power of one of the generators drops to zero), (4) and (5) are triggered to find a set of the most suitable lines to trip once the ICI activation frequency threshold is met. Based on this, two islands are expected. The idea is to maximize the overall consumption in both of the islands combined (6), regardless of whether there is a power deficit or a power surplus in individual islands. The condition that has to be considered in the process is that at the moment of islanding, the power generation in both of the islands must be greater than the consumption (7). The taken actions involve the redistribution of UFLS stages by manipulating each load in order to achieve (8). After that, the DC power flow calculation (9) is performed to verify the occurrence of bottlenecks that might endanger the stability of the islands. If there are none, then commands are sent to the load frequency relays carrying new settings. The entire process is performed prior to the first UFLS stage activation. Afterwards, the Python script waits for UFLS or ICI activation.
The approach from Figure 3 was tested using the PowerFactory DIgSILENT 2021 SP 1 software for RMS (root mean square) dynamic simulations in conjunction with Python. Namely, PowerFactory offers a feature of running an arbitrary Python code in each time step of the simulation.

Case Study
In order to invoke extreme frequency fluctuations, the tripping of multiple synchronous generators is the most obvious option (as exemplified in the IEEE 118 test system from [10,21] where two generators were tripped to reach and fully exhaust UFLS and eventually activate ICI). However, in this paper we decided to trip G1 instead, which represents an equivalent of the remaining network. It should be noted that generation tripping does not affect the network topology.
Furthermore, multiple scenarios were observed during the operation of all machines due to high system loading. In those conditions, we examined the effect of having different Energies 2021, 14, 2790 9 of 13 load models, ranging from constant P, constant Z, and down to the exponential model with/without frequency dependence.
In Figure 4, a frequency response of the observed test system is depicted for applying three different load models. As can be seen, the frequency drops below the selected ICI activation threshold in all three cases (see the red horizontal line), whereas the underfrequency protection of generators is activated only when the constant Z load model was used (with or without frequency dependence; see the black line). This is why in the continuation, only the constant Z model (with and without frequency dependence) is considered as a worst-case scenario. In Figure 4, a frequency response of the observed test system is depicted for applying three different load models. As can be seen, the frequency drops below the selected ICI activation threshold in all three cases (see the red horizontal line), whereas the underfrequency protection of generators is activated only when the constant Z load model was used (with or without frequency dependence; see the black line). This is why in the continuation, only the constant Z model (with and without frequency dependence) is considered as a worst-case scenario. At simulation time t = 1 s, generator G1 is tripped which causes an outage of 2900 MW. Prior to this contingency, the total generation equals 5646 MW. The frequency begins to decay, UFLS is activated and, finally, the system is split in half. The algorithm selected the cut-set on lines 15-16 and 16-17 as the most appropriate (denoted with the red dashed line in Figure 1). Portions of the network that are to operate either in island 1 or island 2 are also clearly marked in Figure 1. Without the coordination between ICI and UFLS, island 1 would initiate island operation with an active power deficit of 2254 MW, whereas island 2 would initiate island operation with a surplus of active power in the amount of 927 MW. Clearly, the mere range or these imbalance values is enough to conclude that both islands have little chance for successful stabilization, which was confirmed by our simulations (depicted with dashed curves in Figure 5 with the legend entry w/o). By applying the presented methodology, the frequency in both islands stabilized successfully (see solid lines in Figure 5 with the legend entry w). The orange and the red curves depict simulation results when the load is modelled with a constant impedance and, therefore, no frequency dependence. The blue and the green curves, on the other hand, additionally consider the load as frequency dependent (legend entry 1 refers to island number 1, legend entry 2 refers to island number 2).
In island 1, apart from UFLS as determined by Table 1, the algorithm dictated an additional disconnection of loads at buses 18, 25, 26, 27, 28, 29 and 31. At the same time, the load at bus 15 was decreased to 30% of its initial value. In island 2, the algorithm suggested to keep loads at buses 20, 21, 22, 23 and 24 in operation. Their tripping in accordance with UFLS was therefore blocked. In addition, the load at bus 16 was scheduled to 85% of its initial value. The analysis clearly showed that our approach has a beneficial impact on EPS operation and increases the probability for avoiding blackouts in critical operating conditions.  At simulation time t = 1 s, generator G1 is tripped which causes an outage of 2900 MW. Prior to this contingency, the total generation equals 5646 MW. The frequency begins to decay, UFLS is activated and, finally, the system is split in half. The algorithm selected the cut-set on lines 15-16 and 16-17 as the most appropriate (denoted with the red dashed line in Figure 1). Portions of the network that are to operate either in island 1 or island 2 are also clearly marked in Figure 1. Without the coordination between ICI and UFLS, island 1 would initiate island operation with an active power deficit of 2254 MW, whereas island 2 would initiate island operation with a surplus of active power in the amount of 927 MW. Clearly, the mere range or these imbalance values is enough to conclude that both islands have little chance for successful stabilization, which was confirmed by our simulations (depicted with dashed curves in Figure 5 with the legend entry w/o). By applying the presented methodology, the frequency in both islands stabilized successfully (see solid lines in Figure 5 with the legend entry w). The orange and the red curves depict simulation results when the load is modelled with a constant impedance and, therefore, no frequency dependence. The blue and the green curves, on the other hand, additionally consider the load as frequency dependent (legend entry 1 refers to island number 1, legend entry 2 refers to island number 2). The impact of considering the frequency dependence of the load model can be observed by comparing the green and the red dashed curve (island 2), as well as the yellow with the blue curve (island 1), in Figure 5. One can observe slight differences in the overall frequency response. In the case without UFLS-ICI coordination, the frequency overshoot was reduced by 0.4 Hz (from 63.6 Hz to 63.2 Hz). With regard to the post-mortem value of the frequency in both islands, no significant differences can be observed.

Conclusions
The presented methodology for the coordination between UFLS and ICI by applying mixed-integer linear programming shows good results in increasing the probability for avoiding blackouts in critical operating conditions. The feasibility of the approach was tested on an IEEE 39-bus test system model, which is most commonly used for this type of research. The combination of minimum power flow disruption and minimum power imbalance is used for determining the most appropriate cut-set. The latter is achieved with the use of flexible UFLS stages.
The concept can be implemented in those EPSs that have WAMS installed with the support from a sufficient number of PMU devices. With regard to EPS modelling, we emphasized the need for putting proper focus on the modelling of the load whenever frequency-related events are simulated. Simulations showed that the selected load model has a huge impact on the EPS frequency response and can, therefore, play a crucial role in how real-life conditions can be represented by the dynamic model.
Results presented in this paper were obtained by means of RMS dynamic simulations. More accurate and realistic EMT dynamic simulations are expected in the near future by using a real-time digital simulator and a hardware-in-the-loop test set-up. Such analysis will consider several realistic limitations that were ignored in this paper, e.g., relay measuring accuracy, relay measuring, tripping delay, etc.    In island 1, apart from UFLS as determined by Table 1, the algorithm dictated an additional disconnection of loads at buses 18, 25, 26, 27, 28, 29 and 31. At the same time, the load at bus 15 was decreased to 30% of its initial value. In island 2, the algorithm suggested to keep loads at buses 20, 21, 22, 23 and 24 in operation. Their tripping in accordance with UFLS was therefore blocked. In addition, the load at bus 16 was scheduled to 85% of its initial value. The analysis clearly showed that our approach has a beneficial impact on EPS operation and increases the probability for avoiding blackouts in critical operating conditions. The impact of considering the frequency dependence of the load model can be observed by comparing the green and the red dashed curve (island 2), as well as the yellow with the blue curve (island 1), in Figure 5. One can observe slight differences in the overall frequency response. In the case without UFLS-ICI coordination, the frequency overshoot was reduced by 0.4 Hz (from 63.6 Hz to 63.2 Hz). With regard to the post-mortem value of the frequency in both islands, no significant differences can be observed.

Conclusions
The presented methodology for the coordination between UFLS and ICI by applying mixed-integer linear programming shows good results in increasing the probability for avoiding blackouts in critical operating conditions. The feasibility of the approach was tested on an IEEE 39-bus test system model, which is most commonly used for this type of research. The combination of minimum power flow disruption and minimum power imbalance is used for determining the most appropriate cut-set. The latter is achieved with the use of flexible UFLS stages.
The concept can be implemented in those EPSs that have WAMS installed with the support from a sufficient number of PMU devices. With regard to EPS modelling, we emphasized the need for putting proper focus on the modelling of the load whenever frequency-related events are simulated. Simulations showed that the selected load model has a huge impact on the EPS frequency response and can, therefore, play a crucial role in how real-life conditions can be represented by the dynamic model.
Results presented in this paper were obtained by means of RMS dynamic simulations. More accurate and realistic EMT dynamic simulations are expected in the near future by using a real-time digital simulator and a hardware-in-the-loop test set-up. Such analysis will consider several realistic limitations that were ignored in this paper, e.g., relay measuring accuracy, relay measuring, tripping delay, etc. Funding: This work was funded by the Slovenian Research Agency through the Electric Power Systems No. P2-0356 research program, the funding mechanism for young researchers and the resource management for low latency reliable communications in smart grids-LoLag, J2-9232 project.

Conflicts of Interest:
The authors declare no conflict of interest. Edge from node i to node j w i,j

Nomenclature
Weight edge from node i to node j z i,j Binary variable showing edge status between node i and node j P i,j Active power flow from node i to node j X i,j Reactance between nodes i and j δ i Voltage angle of ith node P d,i Load active power at node i after UFLS activation P start d,i Load active power at node i in steady state P g,i Generator active power at node i after UFLS activation k pu Exponential coefficient for voltage dependence of active power k qu Exponential coefficient for voltage dependence of reactive power k pf Frequency dependency coefficient for active power k qf Frequency dependency coefficient for reactive power {•}0 Subscript for variables in steady state {•}(t) Time-dependent variables k min Minimum value for a load in pu in island transition k max Maximum value for a load in pu in island transition