An Extreme Learning Machine Based Adaptive VISMA for Stability Enhancement of Renewable Rich Power Systems

: Maintaining power system stability in renewable-rich power systems can be a challenging task. Generally, the renewable-rich power systems suffer from low and no inertia due to the integration of power electronics devices in renewable-based power plants. Power system oscillatory stability can also be affected due to the low and no inertia. To overcome this problem, additional devices that can emulate inertia without adding synchronous machines can be used. These devices are referred to as virtual synchronous machines (VISMA). In this paper, the enhancement of oscillatory stability of a realistic representative power system using VISMA is proposed. A battery energy storage system (BESS) is used as the VISMA by adding an additional controller to emulate the inertia. The VISMA is designed by using Fruit Fly Optimization. Moreover, to handle the uncertainty of renewable-based power plants, the VISMA parameters are designed to be adaptive using the extreme learning machine method. Java Indonesian Power Grid has been used as the test system to investigate the efﬁcacy of the proposed method against the conventional POD method and VISMA tuning using other methods. The simulation results show that the proposed method can enhance the oscillatory stability of the power system under various operating conditions.


Introduction
Renewable-based power plants are becoming more popular over the year due to the need for clean and sustainable energy. Among numerous types of renewable-based power plants, the wind and solar photovoltaic (PV) power plants are exhibiting the leading growth rate in the world [1,2]. Research efforts in [3,4] show the benefits of integrating renewable-based power plants from an economic and environmental point of view. The significant reduction of fuel and environmental costs could be achieved by using renewablebased power plants. Moreover, it also reported that CO 2 and other pollutants are also reduced when renewable-based power plants are integrated into power systems. This condition brings a benefit to both health and climate condition. Despite all the advantages, minimal error by using the DEA method. Research efforts in [16] proposed a method for designing a power system stabilizer based on the Bat algorithm. The application of a cuckoo search algorithm to design a power system stabilizer is also reported in [17]. Moreover, fruit fly optimization (FFO) is reported to give an excellent performance in handling load frequency control [18]. However, a particular operating condition is used to tune the control parameters by using the metaheuristic algorithm. Furthermore, it is essential to adapt operating condition variations in the controller design. Therefore, a control method with adaptability is sought. This paper proposed a method for designing VIC on BESS named virtual synchronous machine (VISMA), which is adaptable to the fluctuation of operating conditions. The extreme learning machine (ELM) method is used to design the adaptive VISMA. Before the VISMA parameters are trained with ELM, the VISMA parameters are tuned using FFO. This process aims to get the best parameters for ELM training. From the comprehensive numerical assessments, it is evident that the proposed method is superior compared to other damping and tuning techniques of VISMA. The rest of the paper is organized as follows: Section 2 presents the dynamic modeling of the conventional power plant, PV, wind, BESS, and VISMA. The Fruit fly optimization, ELM, and the method for designing adaptive VISMA are presented in Section 3. Section 4 represents the results and discussions of the paper. While Section 5 highlights the contribution and the conclusions of this work.

Dynamic Model of Synchronous Power Plant
Transforming the non-linear model of the synchronous generator into a linear model is essential for oscillatory stability studies. The synchronous generator model can be linearized using Park transformation [19]. The details of the park transformation and the synchronous generator non-linear model can be found in [19] and the references therein.
The state-space model of the synchronous generator based on the park transformation for the oscillatory stability can be expressed as Equation (1). (1) The excitation system (e.g., AVR) model is required for power system oscillatory stability analysis. An excitation system is used to regulate the generator output variables, such as voltage, power factor, and also can be used to integrate power system stabilizer (PSS). The variables are set by setting the field flux on the generator. In this research, a fast exciter system is considered. The mathematical representation of the fast exciter can be described by Equation (2).
In Equation (2), K A and T A are the gain constant and time delay constant, respectively. Furthermore, the output of the exciter is limited by the saturation block, as shown in  [20]. The parameters of the fast excitation system are taken from [19,20], and the references therein.

=
( − ) 1 − (2) In Equation (2), KA and TA are the gain constant and time delay constant, respectively. Furthermore, the output of the exciter is limited by the saturation block, as shown in Figure 1 [20]. The parameters of the fast excitation system are taken from [19,20], and the references therein. Figure 1. Schematic diagram of fast exciter [20].
The mechanical torque of the generator depends on the speed droop constant, the time constant of the governor, and the energy source. If there is a change in the generator rotation, the governor will provide feedback to achieve a new speed (new condition). In this research, the governor has been modeled as a first-order differential equation with a gain constant, as shown in Figure 2. The first-order governor model is sufficient enough for the oscillatory stability studies due to the slow time-scale operation of the governor system [19,20]. The initial governor parameters are taken from the [19] and the references therein. In this model, the GSC value is zero and the effect of the combination of turbine system and speed governor produces mechanical power which can be described in Equation (3). In Equation (3), Kg and Tg are gain constant and time constant of the governor [20].

Dynamic Model of WECS
The permanent magnet synchronous generator (PMSG) based wind energy conversion systems (WECS) are considered in this work. The schematic diagram, which includes the PMSG, wind turbine, converter, associated controllers given in Figure 3, is used for the simulation. The power output as a function of wind velocity can be expressed as Equation (4). The mechanical torque of the generator depends on the speed droop constant, the time constant of the governor, and the energy source. If there is a change in the generator rotation, the governor will provide feedback to achieve a new speed (new condition). In this research, the governor has been modeled as a first-order differential equation with a gain constant, as shown in Figure 2. The first-order governor model is sufficient enough for the oscillatory stability studies due to the slow time-scale operation of the governor system [19,20]. The initial governor parameters are taken from the [19] and the references therein.
In Equation (2), KA and TA are the gain constant and time delay constant, respectively. Furthermore, the output of the exciter is limited by the saturation block, as shown in Figure 1 [20]. The parameters of the fast excitation system are taken from [19,20], and the references therein. Figure 1. Schematic diagram of fast exciter [20].
The mechanical torque of the generator depends on the speed droop constant, the time constant of the governor, and the energy source. If there is a change in the generator rotation, the governor will provide feedback to achieve a new speed (new condition). In this research, the governor has been modeled as a first-order differential equation with a gain constant, as shown in Figure 2. The first-order governor model is sufficient enough for the oscillatory stability studies due to the slow time-scale operation of the governor system [19,20]. The initial governor parameters are taken from the [19] and the references therein. In this model, the GSC value is zero and the effect of the combination of turbine system and speed governor produces mechanical power which can be described in Equation (3). In Equation (3), Kg and Tg are gain constant and time constant of the governor [20].

Dynamic Model of WECS
The permanent magnet synchronous generator (PMSG) based wind energy conversion systems (WECS) are considered in this work. The schematic diagram, which includes the PMSG, wind turbine, converter, associated controllers given in Figure 3, is used for the simulation. The power output as a function of wind velocity can be expressed as Equation (4). In this model, the GSC value is zero and the effect of the combination of turbine system and speed governor produces mechanical power which can be described in Equation (3). In Equation (3), K g and T g are gain constant and time constant of the governor [20].

Dynamic Model of WECS
The permanent magnet synchronous generator (PMSG) based wind energy conversion systems (WECS) are considered in this work. The schematic diagram, which includes the PMSG, wind turbine, converter, associated controllers given in Figure 3, is used for the simulation. The power output as a function of wind velocity can be expressed as Equation (4).
In Equation (4), P w represents the power available from the wind, c p is the power coefficient, A presents the rotor swept area of the wind turbine, while, ρ and v represent the air density, and upstream wind speed in the rotor, respectively. The aerodynamic mechanical torque (τ w ) associated with the wind turbine can be obtained as Equation (5) [21].
In Equation (5), ω w presents the rotor speed. For this work, the WECS system without a gearbox has been considered. Therefore, no gearbox model has been developed. Hence, the values of aerodynamic and mechanical torque transmitted to the generator are considered to be the same. Furthermore, the power coefficient of the wind turbine can be expressed by the blade pitch angle and the tip-speed ratio. The power coefficient can be expressed as in Equation (6).
In Equation (4), represents the power available from the wind, is the power coefficient, A presents the rotor swept area of the wind turbine, while, and represent the air density, and upstream wind speed in the rotor, respectively. The aerodynamic mechanical torque ( ) associated with the wind turbine can be obtained as Equation (5) [21].
In Equation (5), presents the rotor speed. For this work, the WECS system without a gearbox has been considered. Therefore, no gearbox model has been developed. Hence, the values of aerodynamic and mechanical torque transmitted to the generator are considered to be the same. Furthermore, the power coefficient of the wind turbine can be expressed by the blade pitch angle and the tip-speed ratio. The power coefficient can be expressed as in Equation (6).
In Equation (6), 1 − 6 are the coefficient, x represents the function related to the wind turbine rotor, and is the constant term. The values of 1 − 6 and x are obtained from Equation (7) [21].
In Equation (7), represents the pitch angle, is the tip ratio. The tip ratio ( ) can be estimated as Equation (8) [21].
The power train model of WECS is required for oscillatory stability studies. The power train model consists of the rotor shaft, generator, blade pitching, and hub with a In Equation (6), c 1 − c 6 are the coefficient, x represents the function related to the wind turbine rotor, and β is the constant term. The values of c 1 − c 6 and x are obtained from Equation (7) [21]. 1 In Equation (7), ϑ represents the pitch angle, µ is the tip ratio. The tip ratio (µ) can be estimated as Equation (8) [21].
The power train model of WECS is required for oscillatory stability studies. The power train model consists of the rotor shaft, generator, blade pitching, and hub with a blade. The power train model used in this paper can be expressed as Equation (9).
In Equation (9), ω g is the mechanical speed of the generator, B m denotes the damping coefficient, aerodynamic torque is represented by τ w_g , equivalent inertia is given by J eq , and τ e is the electromechanical torque. Furthermore, the parameter of the generator side is described by sub-index, g. The comprehensive mathematical model of PMSG can be expressed by Equations (10)-(13) [21].
In Equations (10)- (13), poles are given by p, ω e is the electrical speed, stator resistance is denoted by R s , L d , and L q are the generator inductances, ψ f is the magnetic flux, and L id are L iq the leakage inductances. Moreover, the electromechanical torque as given in Equation (14) is also employed in the PMSG model [22].
The generator side AC/DC converter incorporates the outer and inner controllers. The outer control loop measured the terminal voltage, DC link voltage and compared those measured signals with their reference values. The signal errors are then synthesized using the conventional PI controller to derive reference values of d-axis and q-axis currents. By considering β dgen and β qgen as auxiliary state variables of the outer control loop of the generator side converter, state equations of the controller can be stated as in Equation (15) [23].
The d-axis and q-axis reference currents in the generator side can be expressed as Equations (16) and (17).
Output variables from the outer control loop are then fed into the inner current control loop as reference values and compared with the actual values of generator currents (i dgen , i qgen ). Auxiliary state variables of γ dgen and γ qgen are used to express the state equations of the inner current controllers and expressed as in Equation (18).
A similar method is implemented to the current control loop to determine the modulation indices m * dgen , m * qgen for the generator-side converter. These modulation indices are afterward employed as control variables for the PWM scheme of the converter. Furthermore, the algebraic equations of modulation indices reference signal can be obtained as in Equations (19) and (20) [23]: Similar to the generator-side control, the grid-side inverter control of WECS consists of outer and inner control loops. The active and reactive power reference values are compared to the measured active and reactive output power in the outer loop. State equations of grid side outer control loop can be derived by considering β dgrid and β qgrid as auxiliary state variables as in Equation (21).
The obtained error from the outer control loop is then regulated by PI controllers, yield the reference values for the inner current control loop as follows: Moreover, the inner current control loop generated modulation indices (m * dgrid , m * qgrid ) to provide the switching signal for the grid side inverter. Auxiliary state variables of γ dgrid and γ qgrid are required to provide state equations of the inner current control loops as given in Equation (24).
The algebraic equations of the reference signal for modulation indices of the grid-side inverter are given by Equations (25) and (26) [23].

Dynamic Model of PV Generation
The transmission level PV generator model has been used in this paper. The generic model of the PV, developed by NERC, WECC, and IEEE Task Force, has been used in this study. Based on the NERC's special report, the PV model can be consisting of the grid-side converter and control model of the type-4 wind turbine. Figure 4 shows the dynamic model of a transmission level PV system [24]. The figure shows that the power order has been used to represent the generation of PV array power.
The obtained error from the outer control loop is then regulated by PI controllers, yield the reference values for the inner current control loop as follows: Moreover, the inner current control loop generated modulation indices ( * , * ) to provide the switching signal for the grid side inverter. Auxiliary state variables of and are required to provide state equations of the inner current control loops as given in Equation (24).
The algebraic equations of the reference signal for modulation indices of the gridside inverter are given by Equations (25) and (26)

Dynamic Model of PV Generation
The transmission level PV generator model has been used in this paper. The generic model of the PV, developed by NERC, WECC, and IEEE Task Force, has been used in this study. Based on the NERC's special report, the PV model can be consisting of the gridside converter and control model of the type-4 wind turbine. Figure 4 shows the dynamic model of a transmission level PV system [24]. The figure shows that the power order has been used to represent the generation of PV array power.

Converter Current Limit
Network Interface  The first-order transfer function with unity gain has been considered to represent the converter dynamics. This facilitates the current injection into the network with respect to the electrical controller's active and reactive power command. The active power order to the system depends on the model's power order (PPV), while the reactive power generation depends on the comparative signal generation from the reactive power control loop (e.g., q-axis loop). Three different reactive controls can be implemented in the PV system (e.g., voltage control, power factor, reactive power control) and can be interchangeable based on the flag settings in the controller [24].

Bus Terminal
The controller is comprised of PI reactive power control and converter current limits. The MPPT of PV embraces the logic algorithm to track the maximum power of the PV array. Hence, the MPPT can be neglected due to the no dynamics. The detailed presentation of the PV plant for the electromechanical time-frame analysis including the dynamic parameters can be found in [24].

Dynamic Model of Virtual Synchronous Machine
The battery energy storage system (BESS) including virtual inertia controller is used in this paper to represent the VISMA. The BESS model consists of a second-order model with battery cells. This also includes a first-order model of converter dynamics, and a second-order model of active and reactive power control.
The active and reactive power control signal of a BESS can be obtained as Equations (27) and (28) [25].
In Equations (27) and (28), K BP represents the control loop gain and T BP is the time constant of the active power loop. The K BQ is the reactive power loop gain and T BQ is the time constant of the reactive power controller. The converter dynamic can be obtained as Equations (29) and (30) [25].
In Equations (29) and (30), K R denotes the gain of the converter and T R is the time delay associated with the firing angle, while, K M and I BESS are used to stabilize the BESS during the constant current operation so that the BESS could able to release more power from the battery banks. Moreover, P * BESS and Q * BESS are active and reactive power output of the converter controller. Furthermore, the battery cell dynamics are also required for this study. The dynamics model of the battery cell can be obtained as Equations (31)-(33) [25].
In Equations (31)-(33), R BP and C BP are the self-discharging resistance and capacitance of the battery, while R B1 and C B1 represent the energy and voltage in charging and discharging. Moreover, V BOC and V B1 represent the open-circuit voltage and overvoltage of the battery, respectively. Furthermore, V BT , R BS , and R BT are the terminal voltage, the series resistance of the battery, and the equivalent resistance (e.g., parallel/series connection) of batteries. Figure 5 shows the block diagram of the BESS dynamic model used for the simulation. This is a fifth-order model of battery with converter and the controller. It should be worth noting that the given fifth-order model of battery also consists of POD to compare the proposed VISMA with respect to the auxiliary damping control in the BESS for oscillation damping. To make the BESS operate as VISMA, an additional control scheme has been added. Figure 6 shows the dynamic model of the virtual inertia controller [26]. The input of this controller is a frequency deviation, while the output of this controller is a reference signal of the converter.
ulation. This is a fifth-order model of battery with converter and the controller. It should be worth noting that the given fifth-order model of battery also consists of POD to compare the proposed VISMA with respect to the auxiliary damping control in the BESS for oscillation damping. To make the BESS operate as VISMA, an additional control scheme has been added. Figure 6 shows the dynamic model of the virtual inertia controller [26]. The input of this controller is a frequency deviation, while the output of this controller is a reference signal of the converter.

Regulator Amplifier
Regulator Amplifier Figure 5. Schematic diagram of fifth-order model BESS with POD.

Extreme Learning Machine
Extreme learning machine (ELM) has become popular in recent years in addressing various engineering and other science problems. It was introduced by G. B. Huang in 2006. The ELM structure is constructed by a single hidden layer feed-forward network (SLFN). The ELM method is developed to overcome the weakness in the conventional SLFN, especially in terms of the learning process. The parameters of ELM in its hidden layer are selected randomly and do not have to be tuned. Furthermore, it can produce good generalizations. Thus, it makes the correlation between output and input data is extraordinary fast compared to the other traditional learning algorithms such as SLFN. In addition, the hidden nodes could be set up before the training samples are gained. There are two reasons behind the slow learning process of SLFN [27]. These are: 1. The application of a gradient-based learning algorithm to be conducted for the training process. 2. Parameters on the network are determined via iteration using this learning method.
The manual learning process is used in the gradient-based learning algorithm. The parameters in the conventional gradient-based method are the input weight and hidden ulation. This is a fifth-order model of battery with converter and the controller. It should be worth noting that the given fifth-order model of battery also consists of POD to compare the proposed VISMA with respect to the auxiliary damping control in the BESS for oscillation damping. To make the BESS operate as VISMA, an additional control scheme has been added. Figure 6 shows the dynamic model of the virtual inertia controller [26]. The input of this controller is a frequency deviation, while the output of this controller is a reference signal of the converter.

Regulator Amplifier
Regulator Amplifier

Extreme Learning Machine
Extreme learning machine (ELM) has become popular in recent years in addressing various engineering and other science problems. It was introduced by G. B. Huang in 2006. The ELM structure is constructed by a single hidden layer feed-forward network (SLFN). The ELM method is developed to overcome the weakness in the conventional SLFN, especially in terms of the learning process. The parameters of ELM in its hidden layer are selected randomly and do not have to be tuned. Furthermore, it can produce good generalizations. Thus, it makes the correlation between output and input data is extraordinary fast compared to the other traditional learning algorithms such as SLFN. In addition, the hidden nodes could be set up before the training samples are gained. There are two reasons behind the slow learning process of SLFN [27]. These are: 1. The application of a gradient-based learning algorithm to be conducted for the training process. 2. Parameters on the network are determined via iteration using this learning method.
The manual learning process is used in the gradient-based learning algorithm. The parameters in the conventional gradient-based method are the input weight and hidden

Extreme Learning Machine
Extreme learning machine (ELM) has become popular in recent years in addressing various engineering and other science problems. It was introduced by G. B. Huang in 2006. The ELM structure is constructed by a single hidden layer feed-forward network (SLFN). The ELM method is developed to overcome the weakness in the conventional SLFN, especially in terms of the learning process. The parameters of ELM in its hidden layer are selected randomly and do not have to be tuned. Furthermore, it can produce good generalizations. Thus, it makes the correlation between output and input data is extraordinary fast compared to the other traditional learning algorithms such as SLFN. In addition, the hidden nodes could be set up before the training samples are gained. There are two reasons behind the slow learning process of SLFN [27]. These are: 1.
The application of a gradient-based learning algorithm to be conducted for the training process. 2.
Parameters on the network are determined via iteration using this learning method.
The manual learning process is used in the gradient-based learning algorithm. The parameters in the conventional gradient-based method are the input weight and hidden bias, which are also interconnected from one layer to another. Therefore, it requires a slow learning speed process, and it is often trapped in the local minima.
The formulation of the output function of ELM for generalized SLFN is defined as follows [27]: where ϕ = [ϕ 1 , . . . , ϕ L ] T is the vector of output weight between the hidden layer of L nodes and the output node. The ψ(x) = [ψ i (x), . . . , ψ L (x)] is the output vector of the hidden layer with respect to x. The input data of ELM in d-dimension is mapped by ψ(x) to L-dimensional hidden-layer feature space H. In terms of the learning process, the ELM algorithm is different from the traditional learning algorithm where the goal of the ELM is to reach not only the smallest error value but also the smallest norm of output weights at the same time as described in Equation (35).
where H is the hidden-layer output matrix.
where I, T, and λ r are identity matrix, target matrix, and regularization factor, respectively. The target matrix can defined as T = [t 1 , . . . , t N ] T . By substituting (36) into (34), the output of ELM can be obtained as defined in Equation (37) [28]:

Fruit Fly Optimization
The Fruit Fly Optimization (FFO) was introduced by Pan in 2011 [29]. In terms of smell and vision, fruit flies are superior compared to other similar species. Fruit flies can detect the food from 40 km away [30]. Fruit flies move toward the food according to the smell. Furthermore, they could also reach their food by using their vision capability [31].
FFO has two important steps namely the smell step and vision step. In the smell phase, the fruit flies move toward the food by using their smelling ability. Once they reached near to the food, the vision phase was established [32]. This second step is used to identify their nearest food. The fruit flies will repeat this process until they find the food. The FFO step can be described as follows [33]:

1.
Initialization of parameters (number of particles, number of iterations).

2.
Initialization of the initial position of particles.

3.
Update the location by using smell sensing. The mathematical representation of smell sensing can be described in Equations (38) and (39).

4.
Calculate the distance of the fly from the origin as well as smell concentration judge value (S). The mathematical representation of these procedures can be described using Equations (40) and (41).
Repeat step 3 to step 6 until the algorithm meets the specified criterion.

Tuning VISMA Using FFO
This section describes the procedure of optimizing VISMA parameters (K w and H) using FFO. The objective function of FFO can be mathematically described using a comprehensive damping index (CDI) as given in Equation (47) [34].
The procedure of tuning VISMA parameters based on FFO includes the following steps: Step 1. Linearize the Java-Indonesian grid with all generations and BESS around the chosen operating point.
Step 2. Start the FFO by initializing the number of fruit flies, number of iteration, number of optimized parameters.
Step 3. Generate random position of fruit flies as the initial position.
Step 4. Update the location by using smell sensing.
Step 5. Calculate the distance of the fly from the origin and estimate the smell concentration.
Step 6. Obtain the smell concentrate value.
Step 7. Determine the best fly based on the objective function.
Step 8. Save the best concentration value and the best location of fruit flies.
Step 9. Repeat step 8 until the algorithm meets the criterion (i.e., minimum iteration or objective function).
Step 10. Print the results (K w and H value).

Implementation of ELM and FFO
This section illustrates the process of implementing the proposed method (adaptive VISMA based on Hybrid ELM and FFO). The FFO is used here for the initial parameter tuning. The resulted parameters are then used to obtain the VISMA parameters through the ELM training process. The purpose of this process is to make the VISMA adaptive against the uncertainty of the system operating conditions. The input parameters for training are active power (P) and reactive power (Q). While the output (predicted) parameters for VISMA are K w and H, respectively.
The procedure of adaptive VISMA includes the following steps: Step 1. Java-Indonesian system linearized model used in the FFO should be considered.
Step 2. Initiate to optimize the VISMA parameters.
Step 3. Facilitate the training data from the FFO process.
Step 5. Compute the accurateness of the predicted parameter using mean absolute error as Equation (48) where y i is the predicted data obtained by the proposed learning algorithm, while y i is actual data, and N is the number of data utilized in the training phase. Step 6. If the MSE value is minimum (meet the threshold settings), print out the results.
Otherwise, go to step 4 until the minimum MSE is obtained.

Results and Discussions
The Java-Indonesian electric grid (i.e., 500 kV) has been used in this paper to demonstrate the efficacy of the proposed oscillation damping method. The parameters of the synchronous generators are taken from the realistic model of the Java-Indonesian grid (see Appendix A for system data). The peak load of the realistic system has been used for the simulation. Renewable power generations are considered into the Java grid to simulate the future Indonesian system-2030 scenario (60% from the conventional generators and 40% from renewable-based generations). A large-scale PV plant of 300 MW size is considered to replace a synchronous generator in area 1 (at bus 4). Moreover, an aggregated wind farm is also integrated into area 1. Furthermore, a large-scale BESS of 100 MW size is added to bus 9 in area 2, as shown in Figure 7. The PV and wind power generations are assumed to be operated at the maximum operating point. Table 1 gives the damping performance of a 500 kV Java-Indonesian power network with PV, wind, and BESS. It is also noticeable that there are two weak modes impending from inter-area mode and local mode 2. Thus, this paper's investigation focuses on enhancing the damping performance of these two weak modes in the system.
where ̃ is the predicted data obtained by the proposed learning algorithm, while is actual data, and N is the number of data utilized in the training phase.
If the MSE value is minimum (meet the threshold settings), print out the results. Otherwise, go to step 4 until the minimum MSE is obtained.

Results and Discussions
The Java-Indonesian electric grid (i.e., 500 kV) has been used in this paper to demonstrate the efficacy of the proposed oscillation damping method. The parameters of the synchronous generators are taken from the realistic model of the Java-Indonesian grid (see Appendix A for system data). The peak load of the realistic system has been used for the simulation. Renewable power generations are considered into the Java grid to simulate the future Indonesian system-2030 scenario (60% from the conventional generators and 40% from renewable-based generations). A large-scale PV plant of 300 MW size is considered to replace a synchronous generator in area 1 (at bus 4). Moreover, an aggregated wind farm is also integrated into area 1. Furthermore, a large-scale BESS of 100 MW size is added to bus 9 in area 2, as shown in Figure 7. The PV and wind power generations are assumed to be operated at the maximum operating point. Table 1 gives the damping performance of a 500 kV Java-Indonesian power network with PV, wind, and BESS. It is also noticeable that there are two weak modes impending from inter-area mode and local mode 2. Thus, this paper's investigation focuses on enhancing the damping performance of these two weak modes in the system.

Training Phase
This section describes the training phase of the proposed method. The data used in this paper for the training process are given in Table 2. The active power and reactive power is given in Table 2 are utilized as the input of ELM. Furthermore, the K w and H (virtual inertia parameters) are the output of ELM. The values of active and reactive power are obtained from various operating conditions. While the virtual inertial parameters are the optimized parameters. Figures 8 and 9 illustrate the comparison between actual K w and H parameters with K w and H based on the ELM. It is noticeable that ELM can predict the value of K w and H optimally as indicated by the pattern given in the Figures Moreover, Tables 3 and 4 show the detailed features of Figures 8 and 9. It is found that ELM produces a minimum error for the actual parameter and predicted parameter. Moreover, it is found that the maximum error for K w parameter is 0.0088 and the minimum error is 0.0036. In addition, it is also noticeable that the mean square error (MSE) value for the K w parameter is 0.005817. Furthermore, the maximum and minimum errors for the H parameter are 0.0060 and 0.0047, respectively. While the MSE value for H is 0.005117.     Table 4. Detailed features of the results are given in Figure 9.

Testing Phase
The efficacy of the proposed method is tested, compared, and discussed in this section. The aggregated active power demand of 14,022.4 MW and 8208.16 MVar reactive power demand are considered for the testing of the proposed method. Table 5 shows the damping comparison of the system with BESS (no POD), the system with VISMA, and the system with POD at BESS (optimized using Bat Algorithm). In this paper, only inter-area and local area 2 are used to test the efficacy of the controller as both are identified as weak modes. It is found that a system with VISMA can enhance the damping performance of weak modes significantly. It is noticeable that the inter-area mode damping enhanced from 2.24% to 19.74%, while the local mode 2 dampings increased from 2.11% to 7.42%. The results given in Table 5 show that the system experienced better damping performance for the proposed control compared to the system with BESS (no POD) and the system with POD at BESS. To validate the results given in Table 5, non-linear time-domain simulations are carried out. The time-domain simulations are carried out by applying a small perturbation of 0.01 pu load change. Figures 10 and 11 show the rotor speed responses of G3 and G4. It is found that the system with only BESS exhibits higher oscillation compared to the system with VISMA and POD at BESS. It is also observed that the oscillatory condition of the system with POD at BESS method is almost similar to the system with VISMA. Hence, this non-linear time-domain simulation validated the damping performance results given in Table 5.   Tables 6 and 7 show the detailed features of rotor speed responses of G3 and G4. It is found that systems with only BESS experience longer oscillatory conditions compared to the other scenarios. Moreover, it is noticeable that systems with actual VISMA have smaller overshoot compared to the POD at BESS. Furthermore, Table 8 shows the damping performance of the system under different operating conditions. It is found that under different operating conditions, the damping ratios investigated modes are maintained above the damping standard.

Index
System with BESS System with VISMA System with POD at BESS   Tables 6 and 7 show the detailed features of rotor speed responses of G3 and G4. It is found that systems with only BESS experience longer oscillatory conditions compared to the other scenarios. Moreover, it is noticeable that systems with actual VISMA have smaller overshoot compared to the POD at BESS. Furthermore, Table 8 shows the damping performance of the system under different operating conditions. It is found that under different operating conditions, the damping ratios investigated modes are maintained above the damping standard.

Index
System with BESS System with VISMA System with POD at BESS Figure 11. Rotor speed responses of G4. Tables 6 and 7 show the detailed features of rotor speed responses of G3 and G4. It is found that systems with only BESS experience longer oscillatory conditions compared to the other scenarios. Moreover, it is noticeable that systems with actual VISMA have smaller overshoot compared to the POD at BESS. Furthermore, Table 8 shows the damping performance of the system under different operating conditions. It is found that under different operating conditions, the damping ratios investigated modes are maintained above the damping standard.

P (MW) Q (MVar)
Inter-Area (%) Local Mode 2 (%) Figure 12 shows the damping ratios for different scenarios (i.e., Scenario 1: System without BESS; Scenario 2: System with BESS; Scenario 3: System with VISMA-FFO; Scenario 4: System with VISMA-proposed). The results show that the proposed method could provide better damping to the critical mode. It is worth noting that the performance of VISMA by FFO is lower compared to the proposed method.   Figure 12 shows the damping ratios for different scenarios (i.e., Scenario 1: System without BESS; Scenario 2: System with BESS; Scenario 3: System with VISMA-FFO; Scenario 4: System with VISMA-proposed). The results show that the proposed method could provide better damping to the critical mode. It is worth noting that the performance of VISMA by FFO is lower compared to the proposed method.  Figures 13 and 14 show the non-linear simulation results under three-phase fault at area 2 for the given scenarios in Figure 12. The results in Figures 13 and 14 show that the proposed method could settle the rotor oscillation due to severe fault faster. It is worth noting that the performance of VISMA by FFO is lower compared to the proposed method.  Figures 13 and 14 show the non-linear simulation results under three-phase fault at area 2 for the given scenarios in Figure 12. The results in Figures 13 and 14 show that the proposed method could settle the rotor oscillation due to severe fault faster. It is worth noting that the performance of VISMA by FFO is lower compared to the proposed method.  The comparison between the proposed method (VISMA based ELM) with other tuning methods (VISMA based ANN and VISMA based PSO) is given next. Damping ratio comparison and non-linear time-domain simulation show the proposed method's performance compared with ANN and PSO methods. Table 9 shows the damping comparison between VISMA based ELM, VISMA based ANN, and VISMA based PSO. It is noticeable that VISMA based ELM gives a better damping performance than the other method.  16 show the non-linear time-domain simulation results comparison between the proposed method, VISMA based ANN and VISMA based PSO. It is found that the system gives the best response with the proposed method. It is found that the system with the proposed method has the lowest overshoot and the fastest settling time.  The comparison between the proposed method (VISMA based ELM) with other tuning methods (VISMA based ANN and VISMA based PSO) is given next. Damping ratio comparison and non-linear time-domain simulation show the proposed method's performance compared with ANN and PSO methods. Table 9 shows the damping comparison between VISMA based ELM, VISMA based ANN, and VISMA based PSO. It is noticeable that VISMA based ELM gives a better damping performance than the other method.  16 show the non-linear time-domain simulation results comparison between the proposed method, VISMA based ANN and VISMA based PSO. It is found that the system gives the best response with the proposed method. It is found that the system with the proposed method has the lowest overshoot and the fastest settling time. The comparison between the proposed method (VISMA based ELM) with other tuning methods (VISMA based ANN and VISMA based PSO) is given next. Damping ratio comparison and non-linear time-domain simulation show the proposed method's performance compared with ANN and PSO methods. Table 9 shows the damping comparison between VISMA based ELM, VISMA based ANN, and VISMA based PSO. It is noticeable that VISMA based ELM gives a better damping performance than the other method.  Figures 15 and 16 show the non-linear time-domain simulation results comparison between the proposed method, VISMA based ANN and VISMA based PSO. It is found that the system gives the best response with the proposed method. It is found that the system with the proposed method has the lowest overshoot and the fastest settling time. Hence, the system would have a better damping margin when the system is added with VISMA based on ELM. Hence, the system would have a better damping margin when the system is added with VISMA based on ELM.

Robustness Analysis
Two indicators are used to assess the robustness of the system with the proposed method. The structured singular value (SSV) is used as the first indicator. The SSV can be used to examine the robustness of the system against uncertainty. The paper states that the system consists of large-scale PV and wind power plants with intermittency in output. Therefore, it is essential to evaluate the robust stability against uncertainty. The generalized technique for robustness assessment using the linear fractional transformation (LFT) is depicted in Figure 17, where M indicates the complex transfer matrix. The complex transfer matrix can be described as described in Equation (49).  Hence, the system would have a better damping margin when the system is added with VISMA based on ELM.

Robustness Analysis
Two indicators are used to assess the robustness of the system with the proposed method. The structured singular value (SSV) is used as the first indicator. The SSV can be used to examine the robustness of the system against uncertainty. The paper states that the system consists of large-scale PV and wind power plants with intermittency in output. Therefore, it is essential to evaluate the robust stability against uncertainty. The generalized technique for robustness assessment using the linear fractional transformation (LFT) is depicted in Figure 17, where M indicates the complex transfer matrix. The complex transfer matrix can be described as described in Equation (49).

Robustness Analysis
Two indicators are used to assess the robustness of the system with the proposed method. The structured singular value (SSV) is used as the first indicator. The SSV can be used to examine the robustness of the system against uncertainty. The paper states that the system consists of large-scale PV and wind power plants with intermittency in output. Therefore, it is essential to evaluate the robust stability against uncertainty. The generalized technique for robustness assessment using the linear fractional transformation (LFT) is depicted in Figure 17, where M indicates the complex transfer matrix. The complex transfer matrix can be described as described in Equation (49).
In Equation (49), ∆ ∈ C (q 1 +p 1 ) is another complex matrix. Moreover, w and z are the vectors of perturbation output and input. The upper LFT with respect to ∆ can be described as in Equation (50).
In Equation (50), the unperturbed system and the nominal transfer matrix knowledge due to the perturbation is labeled as M 22  : Figure 17. Structure of robustness analysis using SSV [35]. following steps can be used for the SSV analysis in the power system: Step 1. Non-linear system should be linearized around the selected operating points.
Define the uncertainty (this could be the number of operating conditions, penetration levels of renewable generations, and others).
Obtain the matrix A with different uncertainty by a certain interval.
Get the s d (that will also vary within the interval of different uncertainty. Hence the s d can be used to capture the variation of uncertainty).
Estimate the  by using Equation (51).
Step 6. Obtain   The details about the SSV calculation can be found in [35] and the references therein. The second indicator is the transient kinetic energy (TKE) assessment. The TKE assessment estimates the amount of mechanical energy used by the system in the transient condition. The smaller the value, the less the mechanical energy being used. The TKE can be estimated as in Equation (53).
In Equation (53), the angular momentum of the rotor at synchronous speed and speed deviation of i th generator in a system with n machines are described as Ji and Δωi. Table 10 shows the SSV values comparison between the system with BESS only, the system with In Equation (51), the number of full blocks, repeated scalar blocks, and identity matrix is denoted by ∆ F , δ s and I rs , δ s ∈ R, ∆ s ∈ C m s ×m s . Moreover, the smallest structured uncertainty ∆ is defined as SSV. Furthermore, the SSV or µ can be expressed as in Equation (52). In Equation (52), µ(M) determines the robustness of the system. According to LFT, the system can be specified as a robustly stable system if and only if maxµ(M) < 1. The following steps can be used for the SSV analysis in the power system: Step 1. Non-linear system should be linearized around the selected operating points.
Step 2. Define the uncertainty (this could be the number of operating conditions, penetration levels of renewable generations, and others).
Step 3. Obtain the matrix A with different uncertainty by a certain interval.
Step 4. Get the δ s (that will also vary within the interval of different uncertainty. Hence the δ s can be used to capture the variation of uncertainty).
The details about the SSV calculation can be found in [35] and the references therein. The second indicator is the transient kinetic energy (TKE) assessment. The TKE assessment estimates the amount of mechanical energy used by the system in the transient condition. The smaller the value, the less the mechanical energy being used. The TKE can be estimated as in Equation (53).
In Equation (53), the angular momentum of the rotor at synchronous speed and speed deviation of ith generator in a system with n machines are described as J i and ∆ω i . Table 10 shows the SSV values comparison between the system with BESS only, the system with the proposed controller method, and the POD at BESS. It is found that the system with BESS has SSV above 1, while the system with the proposed controller method exhibits an SSV value well below 1. Moreover, the system with POD at BESS and VISMA by FFO shows an SSV value lower than 1. However, the SSV values of the system with POD at BESS and VISMA by FFO are higher than the system with the proposed method. Hence, it can be stated that the proposed controller can make the system robust against uncertainty. Moreover, Table 11 shows the TKE value comparison for G3 and G4. It is noticeable that the TKE value for both G3 and G4 are smaller when the proposed controller method is added to the system.

Conclusions
This paper proposed an optimal method to design an adaptive VISMA based on ELM to enhance the small-signal stability of a renewable-rich realistic representative power system. From the simulation results, it is found that the ELM can estimate the optimal parameters for VISMA. In addition, by using ELM, VISMA can automatically adjust the parameters when there is a sudden change in the operating condition. Therefore, it exhibits better performance than the VISMA by FFO. It is worth noticing that the proposed method is also superior compared to other tuning methods for VISMA (e.g., VISMA based PSO and VISMA based ANN). Moreover, the VISMA design using ELM can handle the uncertainty in the power system as indicated by the value of the SSV (the SSV value is below 1). Furthermore, it is also found that by using VISMA based ELM, the kinetic energy produced by the system is smaller compared to the other method for both generators 3 and 4 (0.0000000000304 for G3 and 0.00000000011309 for G4).   Table A2. Java-Indonesian grid generator data [34].