Improving Small-Signal Stability of Power Systems With Significant Converter-Interfaced Generation

This work presents a hierarchical control strategy to improve the stability of electrical networks with significant converter-interfaced generation (CIG). Due to the lack of inertia of CIG systems, these networks can undergo a high rate of change of frequency, compromising the frequency stability. In a first level control, a local controller based on the virtual synchronous generator (VSG) concept is used to emulate inertia and provide short-term frequency regulation. However, the inclusion of significant VSG units can have a negative impact on the damping of inter-area oscillations. Therefore, in a second level control, a centralized controller is proposed to damp these low-frequency electromechanical oscillations affected by VSGs. Several practical issues such as the identification of a system model for the control design, the compensation of communication delays, and the discrete-time implementation of the controller are particularly analyzed. The introduced supplementary controls allow increasing the penetration of renewable energy sources without jeopardizing the frequency and small-signal stability. Eigenvalue analysis and nonlinear hybrid simulations combining DIgSILENT and Python are performed to validate the proposed control strategy.


I. INTRODUCTION
F REQUENCY measurements indicate that the continuous growth of nonconventional renewable energies is reducing the system inertia, consequently degrading the frequency stability [1]. To deal with this situation, system operators are demanding converter-interfaced generation (CIG) systems able to emulate inertia and support the short-term frequency regulation [2]. In the last years, several techniques have been proposed to provide such supplementary services. They can be broadly classified into two approaches. In the first approach, the converter is synchronized with the grid using a phase-locked loop and exchanging an additional power proportional to the variation and time derivative of the frequency [3]; in the second approach, the converter is controlled to emulate the dynamic behavior of a synchronous generator (e.g., virtual synchronous generators (VSGs) [4], synchronverters [5], power-synchronization control [6], and synchronous power control [7]- [9]). In weak systems, the second approach has shown to be very effective to improve the frequency control in terms of both the rate of change of frequency (RoCoF) and the frequency nadir [10]. The emulation of a synchronous generator entails some challenges, similar to the installation of a new conventional generator in a multi-machine power system [11]. A VSG introduces an electromechanical mode [9], and the virtual rotating mass affects existing inter-area modes. The creation and interconnection of areas based on VSGs modify the small-signal characteristics of the system and change the oscillation frequencies in which power system stabilizers were originally tuned. As a consequence, the higher the virtual inertia units are dispatched, the higher the risk of adverse effects on the damping of power system electromechanical oscillations. Additional control loops based on local measurements have been proposed in [12]- [14] to improve the stability and damping properties of VSGs. However, the damping term of these schemes mainly impacts on the local electromechanical mode of the VSG, and it is not fully effective to mitigate inter-area oscillations.
Different decentralized control approaches using local measurements have also been considered in [15]- [19] to improve the damping of power system oscillations. Controls based on local measurements can damp certain modes, but they lack sufficient observability to properly damp inter-area oscillations [20]- [22]. On the other hand, the increasing number of phasor measurement units (PMUs) in power systems and the advances in wide-area measurement systems have enabled to implement control strategies using remote measurements in nearly real time [22]- [24]. Remote measurements transmit information related to the overall network dynamics and have a much higher observability of inter-area modes. In addition, controls based on remote measurements achieve the damping effect with lower control actions [20] and are less prone to destabilize local modes than controls based on local measurements [25]. Although current communication networks face several challenges such as data packet dropouts, random communication delays, and false data injection attacks, different solutions are being proposed to deal with these issues (e.g., see [26]- [28]), making remote communications more reliable and failure resilient. The study of effective strategies to damp inter-area oscillations (e.g., those based on global information) will enable the large-scale integration of renewable energy sources and a better use of their limited energy storage for frequency support services.
In this paper, a hierarchical control strategy is proposed to improve the frequency stability and the electromechanical oscillation damping of multi-machine power systems with large numbers of CIG plants. The control strategy manages the output power of these plants using both local and remote control signals. In a first level control, a local controller provides inertial response with a well-damped local mode, whereas in a second level control, a centralized controller takes advantage of the high inter-area mode observability of wide-area measurements. Two major challenges arise in the practical implementation of centralized controllers, namely, to obtain a system model for the control design and to compensate the time delay of communication channels [29]. This work approaches the first issue by using an identification method to estimate a multi-input multi-output (MIMO) model of the system from time-domain data. These data can be obtained either from field measurements or from software packages that transmission system operators employ to analyze large power systems (e.g., PSS/E and DIgSILENT PowerFactory) [30]- [32]. On the other hand, the second issue is solved by designing a discrete-time controller that considers both the transmission delay and the discrete nature of the sampling and control update processes.

II. FREQUENCY STABILITY AND SMALL-SIGNAL ANALYSIS OF POWER SYSTEMS WITH CIG
Power systems are currently facing the challenge of the continuous displacement of conventional generation by CIG. A case study based on the Argentinian power system is used to illustrate this situation (see Fig. 1). More details about this case study will be given in Section VI-A. Figure 2 shows the response of the system frequency for a loss of infeed of 800 MW (corresponding to 4% of the total generated power). Three cases are considered: the first one is a current scenario with conventional generation and negligible CIG (CONV case); the second one is a future scenario with a penetration of 25% of CIG working as constant power sources (PQ-CIG case); and the third one is the previous case, but with CIG working as VSGs [8] (VSG-CIG case). When the CONV and PQ-CIG cases are compared, a frequency stability degradation in terms of both the frequency nadir (from 49.8 Hz to 49.71 Hz) and the maximum deviation of the RoCoF (from −0.11 Hz/s to −0.16 Hz/s) is observed. The VSG-CIG case has a better frequency response than the PQ-CIG case, but a poorly-damped oscillation arises in its RoCoF curve (see blue line in the bottom subplot of Fig. 2). To analyze this oscillation, the eigenvalues of the three cases are shown in Fig. 3. In the VSG-CIG case (blue diamonds), a poorly-damped inter-area mode appears around 0.5 Hz. Consequently, the VSG is able to damp its local electromechanical mode, but it is not able to properly damp the inter-area modes appearing after a significant VSG-based CIG is included in the power generation share. To address this problem, we propose to add a remote damping signal to virtual power plants to effectively improve the damping ratio of the affected low-frequency electromechanical modes.

III. VIRTUAL POWER PLANTS
The CIG units in a power plant (or area) are controlled using the VSG structure shown in Fig. 4. The active power loop is based  on the classical swing equation of a synchronous machine, and the reactive power loop regulates the desired reactive power. The virtual inertia and damping coefficients are given by H and D, respectively, whereas the current reference to be tracked by the inner current control loop is calculated using the virtual impedance R + jX (see [7] and [8] for further details).
In this work, a virtual power plant is an aggregation area of distributed generation and energy storage systems. An example of such a heterogeneous power plant can be observed in Fig. 5. The total active power delivered to the rest of the system is defined as p poi = p 0 + Δp, where p 0 is the steady-state power and Δp is an additional transient power. This additional power will work either as a probing signal (injected by the identification method) or as a damping signal (generated by the centralized controller). These two modes of operation are selected by the switch SW (see Fig. 5); they will be described in the following sections.
A plant-level control is considered for aggregation purposes. This is achieved by using a power plant controller (PPC) that controls a cluster of elements as a single power unit. The task of the PPC is to coordinate the additional power of both the distributed generation and the energy storage system (ESS) to satisfy that the power Δp = Δp i + Δp ess follows its reference value Δp .

IV. SYSTEM IDENTIFICATION
The proposed centralized control is designed using an identified model of the system. With this aim, a procedure based on the subspace state-space system identification (N4SID) method is described below [31]. In this section, local (first-level) controls are activated; in this way, the centralized (second-level) control is designed considering the dynamics introduced by local controls [33].
The identified model is calculated from time-domain inputoutput data obtained by injecting probing signals to the power system and recording the simulated response. The system is modeled considering the most detailed models available in DIgSILENT to represent power converters, synchronous generators, automatic voltage regulators, power system stabilizers, and turbine-governor systems. Damping controllers based on identified models obtained from simulation software were also studied in [32] and [34]- [36].
The identification procedure can be divided into three stages. In a first stage, the control actions and measurements of the centralized control are selected as the inputs and outputs of the identified model. As previously mentioned, a supplementary signal added to the active power delivered by CIG plants is used to act on the power system. In the case study, the CIGs labeled NOA, COM, CRU, and BAH (see Fig. 1) are chosen because they have the highest rated power, resulting in the following input vector To obtain an acceptable control level, the total power managed by the selected CIG has to be enough to impact on the power system dynamics. On the other hand, voltage angle differences between buses from different areas are considered as measurements due to their high observability of inter-area modes [37]. The output vector is chosen by selecting bus voltage angles from the main areas. Subscripts in (2) indicate the selected buses and determine where the PMUs must be placed. In a second stage, the input-output data required for the system identification are obtained by injecting probing signals from the input vector (1) and recording the response of the output vector (2). To focus the probing signal energy on the frequency range of interest (i.e., inter-area oscillations in the range of 0.2 to 1.0 Hz), the following band-limited signal is used to excite the system between frequencies f L and f H [30] where the constant K id defines the maximum amplitude of the signal, and T p is the time when it is applied. The considered input signals and the measured outputs are shown in Fig. 6. Finally, the probing signals and the recorded measurements sampled at 25 Hz (i.e., a sampling time of h = 40 ms) are introduced to the N4SID algorithm that computes the state-space matrices A d , B d , C d , and D d of the identified model [see (4) and (5) in the next section]. Figure 7 shows a block diagram to illustrate the above methodology.
The N4SID problem can be defined as follows. Given the measurement sequence of input and output vectors {u 0 u 1 u 2 · · · u M −1 } and {y 0 y 1 y 2 · · · y M −1 }, find the matrices A d , B d , C d , and D d so that the identified model has a good fit to these measurements [31]. This procedure involves matrix operations, such as computing block-Hankel and projection matrices, performing singular value decomposition and solving a linear equation system. A step-by-step description of the N4SID algorithm is summarized in [31], and the full development can be found in [38]. In this study, the MATLAB System Identification Toolbox is used to perform the N4SID. Figure 8 compares the eigenvalues of the case-study continuous-time system (a 252nd-order model) and the resulting identified discrete-time system (a 20th-order model). The  identification could capture the low-frequency electromechanical modes from the recorded time-domain data (see the three critical inter-area modes between 0.5 and 0.75 Hz). Because it is more common to see the eigenvalues plotted in the s-plane, and for comparison purposes, the eigenvalues of the discrete system are transformed from the z-plane to the s-plane using the transformation λ z = e hλ s [39].

V. CENTRALIZED CONTROL STRATEGY
The steps of the centralized control design are described in this section.

A. Identified Model
The identified discrete-time model can be written in a statespace form as follows where x ∈ n x ×1 , u ∈ n u ×1 , and y ∈ n y ×1 are the state, input, and output vectors of the identified model, respectively. The matrices A d , B d , C d , and D d are obtained from the N4SID method as explained in the previous section. The input vector has the control actions (i.e., the additional active powers delivered by CIG plants), and the output vector has the remote measurements (i.e., the bus voltage angle differences).

B. Time Delay of Communication Channels
The time delay between the moment a measurement is obtained by the PMU and the moment it is available at the central control station can be modeled as a discrete transfer function z −M (i.e., an Mth-order delay model q k+M = y k ). For a jth measurement, this delay can be written in a state-space form as follows where and h j ∈ M ×1 is the state vector. The signal y kj represents the jth measurement at the PMU location, and the signal q kj represents the jth measurement at the central control location.
Considering that the centralized controller has n y measurements, the complete MIMO delay model can be obtained by combining the single-input single-output delay model of each communication channel, resulting in where with h = [h 1 · · · h n y ] T ∈ Mn y ×1 , y = [y 1 · · · y n y ] T ∈ n y ×1 , and q = [q 1 · · · q n y ] T ∈ n y ×1 . This approach allows considering different time delays for each communication channel.
In a similar way, a model is obtained for the time delay between the moment the control action is computed by the centralized controller and the moment the control action is actually applied by CIG plants. Considering an Nth-order model for the time delay of each control signal, it results where g ∈ Nn u ×1 is the state vector and the matrices A i , B i , C i , and D i are obtained as in the previous MIMO delay model (9)- (14). The vectors v ∈ n u ×1 and u ∈ n u ×1 represent the control signals at the central control location and at the CIG location, respectively. The measurement sampling and the control signal update are performed at a rate of 25 Hz, and a time delay of 200 ms is considered for both measurement and control communication channels; therefore, time-delay models with N = M = 5 are implemented. This approach in the discrete-time domain is able to represent communication time delays more accurately than other approaches in the continuous-time domain based on Padé approximations [40]- [42].

C. Signal-Conditioning and Filtering Stage
Measurements are filtered before computing the control law. The filter consists of a high-pass (washout) filter to guarantee a zero control signal in steady state and a low-pass filter to avoid interactions with high-frequency modes. The pre-filter is added to each input of the centralized controller and limits the control signals within the frequency range of interest. The lower and upper cutoff frequencies of the filter are 0.016 Hz and 3 Hz, respectively. The filter model in a state-space form is given by where s ∈ n f n y ×1 , q ∈ n y ×1 , and w ∈ n y ×1 are the state, input, and output vectors of the filter, respectively, and n f is the filter order.

D. System Extension
An extended model is considered in the control design to compensate both the time delay and the filtering effects (a similar approach can be found in [41]). This model is obtained by combining the identified model (4)-(5), the communication delay models (9)-(10) and (15)- (16), and the filter model (17)-(18), yielding with an extended state vector x e = [ x g h s ] T ∈ n e ×1 , and The above matrices can be simplified because in most of the cases the feedforward matrix D is null (i.e., the models have a strictly proper transfer function). The extended model order is given by n e = n x + Nn u + Mn y + n f n y = 61, where n x = 20 and n f = 2 are the orders of the identified and pre-filter models, n u = 4 and n y = 3 are the numbers of input and output signals of the centralized controller, and the values N = M = 5 are defined by the communication delay. A controller of order n e = 61 can be managed by current CPUs; however, this number can be increased in cases with more input and output signals or higher orders in the identified and prefilter models. For this reason, in the next subsection we reduce the extended model order to a predefined number n r set by the designer.

E. Reduction of the Extended Model Order
Model reduction is often used to simplify the design and to reduce the computational cost of inter-area oscillation damping controllers (e.g., see [33], [34] and [43]). In this study, the MATLAB Control System Toolbox is used to reduce the extended model order to n r < n e via the balanced truncation method, obtaining the following system where x r ∈ n r ×1 andw ∈ n y ×1 are the state and output vectors of the reduced-order model, respectively. The model reduction provides the degree of freedom to choose the final order of the controller (in the case study, it is set to n r = 30), whereas the values n x , n y , n u , n f , N , and M can be chosen to meet the system characteristics and the desired performance. In case n r = n e , no reduction is performed, and the matrices used in the following control design steps (i.e., A r , B r , C r , and D r ) are directly chosen as the extended matrices (A e , B e , C e , and D e ).
The reduction technique guarantees (within the selected frequency range) that, for a given input vector v, the outputw approximates the output w of the extended model (20). However, the reduced-order model (25) cannot be used to design a state-feedback control law because the states x r do not have physical meaning and they cannot be measured. To overcome this problem, a state observer is designed, as described in the following subsection. Note that, in case no reduction is performed, the state observer is still required because the states x of the identified model have the same problem.

F. State Estimation
A Luenberger observer is proposed to estimate the states of the reduced-order model from the known signal ŵ The 'hat' indicates an estimated value, and G r ∈ n r ×n y is the observer gain matrix. The observer consists of a prediction term based on the model (25) and a correction term proportional to the error w −ŵ (see [39]).

G. Control Law
The control signal v k is obtained from a state-feedback control law as follows v k = −K rxrk (29) where the estimated statesx rk are fed back and K r ∈ n u ×n r is the control gain matrix. An overview of the different blocks, matrices, and signals considered in the centralized control is shown in Fig. 9.

H. Gain Selection
The following procedure is performed to choose both the control gain K r and the observer gain G r so that the closed-loop system is stable with a desired dynamic response. For control tuning purposes, it is assumed that the input-output behavior of the extended model has been perfectly captured by the reducedorder model; this means thatw k = w k .
Substituting the control law (29) into (25) and defining the estimation error e k = x rk −x rk , it is obtained Considering (25) and (27), after some mathematical arrangements, the dynamics of the estimation error can be written as Then, the closed-loop system is obtained by combining (30) and (31) By the separation theorem, the dynamics of the system (32) is given by the eigenvalues of the matrices (A r − B r K r ) and (A r − G r C r ). Linear techniques can be used to design K r and G r (as two separate problems) so that these matrices have eigenvalues in a desired region. In our work, the optimal quadratic technique (LQR) is used (see [39] for further details).

I. Stability Verification
The gain selection is achieved by assuming that the outputs of the reduced and extended models are equal (i.e.,w k = w k ). Because the reduced model is an approximation of the extended model, some mismatches are expected. Therefore, the closedloop stability is calculated without the assumption made in the control tuning. With this aim, the dynamics of both the extended states x e and the controller statesx r are combined. Substituting the control law (29) into (19), and using (20) and (28) in the observer (27), the complete closed-loop system is obtained as follows By computing the eigenvalues of the system (33), we can verify the final location of the closed-loop eigenvalues; if the reducedorder model is properly designed, they must be near the ones set in the tuning procedure. This condition will be checked for the case study in the following section.

A. Case Study and Co-Simulation Structure
A case study based on the Argentinian power system is considered in this work (see single-line diagram in Fig. 10 and parameters in [44]). This system has hydro, wind, and solar resources located far from the main load centers that, together with a weakly meshed network, make inter-area oscillations a critical issue. We analyze scenarios with different virtual inertia and damping coefficients in VSGs. In all cases, similar levels of improvement in the inter-area oscillation damping are observed after including the centralized controller. Due to space limitations, we show the most representative results using VSGs with coefficients H = 7 s and D = 35 pu/pu. A new case called VSG-CIG WAC is considered in this section. This case is based on the VSG-CIG case (previously defined in Section II) but including the centralized controller.
A co-simulation able to combine continuous-and discretetime components is conducted to verify the performance of the proposed control strategy. The power system is modeled in DIgSILENT PowerFactory, and the centralized control loop and time delays are implemented in Python. The control signals of the centralized controller are updated in the DIgSILENT simulation every 40 ms (i.e., a data transfer rate of 25 Hz); during this time interval, the DIgSILENT model is run in an open-loop fashion until the next data package arrives from the Python script. Figure 11 shows the structure of this hybrid simulation.

B. Frequency and Small-Signal Stability
To analyze the impact of the centralized controller on the frequency stability, the introductory test of Fig. 2   considering the VSG-CIG WAC case (see Fig. 12). The variation and time derivative of the frequency are similar in the VSG-CIG and VSG-CIG WAC cases. However, the poorly-damped oscillation of the RoCoF observed in the VSG-CIG case is now more damped in the VSG-CIG WAC case (compare blue and purple lines in the bottom subplot of Fig. 12). To assess this damping improvement, the eigenvalues of the VSG-CIG and VSG-CIG WAC cases are shown in Fig. 13. The low-frequency electromechanical modes associated with inter-area oscillations  are significantly damped by the centralized control action (compare blue diamonds with purple squares).

is repeated
In Fig. 13, besides the eigenvalues of the closed-loop system (33) (i.e., the VSG-CIG WAC case), we show the eigenvalues of the closed-loop system (32) (i.e., assuming that the reduced model perfectly captures the behavior of the extended model). It is confirmed that the eigenvalues set in the tuning procedure (light-blue dots) practically match the final closed-loop eigenvalues (purple squares), particularly in the frequency range of the inter-area oscillations from 0.2 to 1.0 Hz. In addition, a comparison of the time-domain responses of the system (32) and the detailed system implemented in DIgSILENT is shown in Fig. 14. The trajectories are very similar in terms of both oscillation frequency and damping ratio, which also agrees with the small-signal characteristics compared in Fig. 13. Figures 13  and 14 show that the system (32) and the detailed system have a similar dynamic behavior justifying the tuning procedure performed in Section V-H.

C. Transient Stability
The transient stability and control performance are evaluated applying a 100-ms short circuit at the GBA bus. Tests considering different levels of both CIG penetration and system load are shown in Figs. 15 and 16, respectively, where speed deviations of synchronous generators for the cases with and without the centralized controller are plotted. In agreement with the previous small-signal analysis, an oscillation damping improvement is observed when comparing the VSG-CIG WAC case with the VSG-CIG case for the different operating conditions and dispatches.
The additional powers of CIG plants demanded by the centralized controller are shown in Fig. 17 (for the sake of clarity, only the case with CIG penetration level of 25% and medium system load is plotted). The maximum power variation is limited to 10% of the rated power of the plants; the discrete behavior of the transmitted control signals is seen in the zoom window of Fig. 17. A total time delay of 440 ms is observed between the short circuit is applied at t = 1 s and the damping power signals start rising. This delay is the sum of the measurement channel delay (200 ms), the sampling time of the centralized controller (40 ms), and the control channel delay (200 ms).
A performance comparison of the centralized controller with another power oscillation damping (POD) controller is introduced in Fig. 18. In this test, the previous fault is applied and a CIG penetration level of 25% and medium system load are used. The considered POD controller is based on the phasecompensation method, similar to standard power system stabilizers (PSSs), commonly found in the literature (e.g., see [40] and [43]). Both control approaches are able to improve the damping of the inter-area oscillations; however, a higher damping ratio is achieved with the proposed control strategy.

D. Failures in Communication Channels
Several tests have been conducted to assess the impact of communication link failures on the centralized controller Fig. 15. Comparison of the system response with the centralized controller (VSG-CIG WAC case, purple lines) and without the centralized controller (VSG-CIG case, blue lines). Scenarios with CIG penetration levels of 15%, 25%, and 35% are shown in each subplot, considering a medium system load. Fig. 16. Comparison of the system response with the centralized controller (VSG-CIG WAC case, purple lines) and without the centralized controller (VSG-CIG case, blue lines). Scenarios with low, medium, and high system loads are shown in each subplot, considering a CIG penetration level of 25%.  control station and a CIG plant, and (c) the complete loss of communication between a PMU and the central control station. The oscillation damping capability of the centralized controller is reduced by the communication link failure but, in all cases, the resulting performance is better than the scenario without the centralized controller (compare generator speed deviations of cases with and without the centralized controller in the top subplots of Fig. 19).

VII. CONCLUSION
Virtual synchronous generators are emerging as a promising solution to overcome the frequency stability degradation in power systems dominated by renewable energy sources and dc interconnections. A potential drawback of the installation of VSG-based plants in multi-machine power systems was analyzed in this work. It was found that these plants, when installed at large scale, could have an adverse effect on the damping of inter-area oscillations. To solve this problem, a centralized control approach was proposed using wide-area measurements. These measurements have higher observability of inter-area oscillations than local measurements and give more effective damping signals to the VSG controller. Several practical issues such as model identification, communication delays, and discrete-time implementation were discussed, and the proposed strategy was verified in different operating conditions and communication failure scenarios. The implementation of these supplementary controls provides a reliable integration of large numbers of CIG plants, and both small-signal stability and frequency stability are enhanced.