Estimation of Rwanda ’ s Power System Inertia as Input for Long-term Dynamic Frequency Response Regulation Planning

Changes in the quota between conventional generation and renewable energy sources constitute major chal- lenges that the modern power systems have to encounter. Conventional power plants are replaced by renewable generation (e.g. wind turbines, photovoltaics) that contribute to the reduction of power system inertia. This may introduce frequency stability issues because frequency is affected by the amount of system inertia, along with the response of controllable frequency reserves and the amount of power imbalance. Therefore, the estimation and analysis of power system inertia and the frequency response assessment is essential to ensure power system stability and security. This paper proposes a novel method to estimates the inertia constant for three different periods in future, namely, 2025, 2035 and 2050 based on the produced future energy scenarios (FES) for Rwandan ’ s power system. In addition, the paper evaluates the frequency response dynamics for each scenario. Results show that the highest progression in renewable energy resources penetration resulted to a larger reduction in the system inertia constant (from 7.2 in control area 1 to 3.83s in control area 3) and the largest frequency drop was observed during the high progression scenario in the year 2050 where the PV and imported power penetration was expected to reach more than 30% of the total installed capacity.


Introduction
As of December 2020, the total installed capacity in Rwanda is around 226.7MW with an import option of 5.5MW. This total capacity of Rwanda includes 45% hydro, 27% diesel, 15% methane gas, 7% peat and 6% solar [1], [2]. The projection shows that the total installed capacity will surpass 1,450MW by 2050. It will also reach 1,383MW and 1, 295MW during medium and high progression scenarios with the current CO 2 reduction policy in place, respectively [3]. These policies include the replacement of diesel-powered power plants with PV and imports [3]. This will be achieved by upgrading the major inter-country power flow in the region and by exporting and/or importing power from Ethiopia through Sudan via several transmission lines (through Uganda) in operation by 2022 [4]. This region will also be connected to the South Africa Power Pool via the Tanzania-Zambia line. With the regional integration and renewable resource penetration into Rwandan's power systems, stability and security issues will also have to be addressed in the future.
Rwanda Utility Regulatory Authority in grid code defines frequency ranges that the customers could expect their supply to stay within [5] . According to this, the sufficient response need to be maintained to ensure that under normal conditions frequency is maintained within 50Hz ±0.5% without interconnection to any power pool, and the control area is considered to be under normal frequency conditions when: (a) the immediate demand can be met with the available scheduled resources, including any expensive contingency resources; (b) the Area Control Error (ACE) deficit does not exceed the available reserves for longer than 15 minutes; (c) the frequency is not less than 49.8Hz for longer than 15 minutes; (d) applicable power pool control performance criteria are not violated; and (e) the frequency is between 49.5 and 50.5Hz. The control area is said to be under abnormal conditions otherwise. Following any system disturbance (i.e., the loss of the largest single unit on the interconnected power supply), the frequency band is extended to operate between 49.0Hz and 51.0Hz (±2%). Under system stress the frequency on the power system could experience variations within the limits of 50Hz ±2.5% (i.e., between 48.75Hz and 51.25Hz). If multiple contingencies occur simultaneously in the system, it enters in a system blackout and its operating condition changes to 'extreme' scenario, where the frequency may drop below 47.5Hz or reach above 51.5Hz (i.e., − 5%/+3%).
Researchers have been trying to find out how the future of Rwanda's power system may look like in the coming years. In [6], the author tracks the possible available and untapped renewable energy resources and outlines credible pathways for Rwanda's energy future in the next 30 years and beyond by considering how much energy the country will need and where this energy could come from by reviewing different energy generation sharing in effective ways to meet the demand. Their results show that fuel-based power plants will be decommissioned, and thus the CO 2 emissions will be reduced. However, about the way forward with power system frequency stability is usually not considered in the literature. In [3], three different scenarios were produced: basic, medium, and high progression, respectively, which reviewed Rwanda's electricity generation capacity, and assessed how the annual output from this generation could change over the next 30 years. This study examines how various technologies could develop between now and 2050, considering the grid integration of renewable generation and electricity interconnectors. It also discusses the development of the mix of generation technologies for each scenario. This paper investigates power system inertia and frequency behaviour under different policies for the period of 30 years from today. Important contextual factors that had not been prominent in the literature for Rwanda's power system have been reviewed and the findings are hoped to create an insightful direction on how these factors can open researchers' involvement including the study and proposition of new frequency control methods. Researchers have proposed different load frequency control methods. In [7], the Frequency Control of Future Power Systems is discussed. The authors considered model representation of a population of the water heater devices for the demand side frequency response. The highlighted model representation of a population of battery energy storage system (BESS)-based distributed energy resources such as smart electric vehicles (EVs) charging, large-scale BESSs, and residential and non-residential BESSs. The simplified Great Britain power system and the 14-machine South-East Australian power system were used to demonstrate the effectiveness of the new methods in controlling power system frequency following a disturbance.
In [8], a novel Adaptive Deadbeat-Based Control for Load Frequency Control of Low Inertia System in Interconnected Zones North and South of Scotland is proposed. An adaptive deadbeat (ADB) controller was developed to investigate its capability in providing a fast frequency response to an electrical power system. This controller was developed to meet the requirements of the National Grid System Operability Framework (SOF), which requires frequency to be accelerated in line with a fast rate of change of frequency (RoCoF) when a high rate of nonsynchronous machines is presented. The controller's parameters were optimized using particle swarm optimization (PSO) to ensure a robust operation and to maintain the proper operation of the power system. The design of the ADB controller was then integrated with the multiarea model of the north and south zones of Scotland. This model was developed to conform to the future energy requirements scenario stated by National Grid whereby regional control can be provided in both the north and south of Scotland. In comparison with the standard PI and Fuzzy-PI controllers used in the four highlighted case studies, it was shown that the ADB controller was able to significantly reduce the RoCoF and deviation of frequency when a sudden loss of generation occurred in a low inertia zone.
Also, the Load Aggregation over a Time of Day to provide Frequency Response in the Great Britain Power System was presented in [9]. The paper discussed the dynamic frequency control to evaluate the capacity that can be gathered from the aggregation of domestic heat pumps and fridges for frequency response. It also estimated the potential of frequency response at a particular time during winter and summer days, and the relationship between both loads (domestic heat pumps and fridges) to provide Firm Frequency Response service. A case study on the simplified Great Britain power system model was developed. Based on this case study, three scenarios of load combination were simulated according to the availability of the load and considering cost savings. Results showed that the aggregation of heat pumps and fridges offered large power capacity and, therefore, an instantaneous frequency response service was achievable. Techniques for estimating system inertia, on the other hand, have been proposed. Modelling-based approaches and standard measure methods are the two groups into which the methodologies have been subdivided. Also, the estimation of inertia at both the individual and system levels is investigated further. Following that, we'll go over the merits and limitations of various approaches for estimating inertia in greater depth. Finally, this paper proposes a potentially useful approach for predicting the predicted inertia constant of the electric power system in the future. This work is intended to contribute in power system in general and in frequency response control in particular. Several difficulties remain unanswered despite the fact that tremendous progress has been made; some of these are discussed in greater depth below. According to [10], the extended Kalman filter determines the system inertia constant by starting with the simple generator model as a starting point and working backwards. When there are large disturbances, there are numerous approaches for determining the inertia constant, including the least-squares [11,12] optimization methodology [13] and several Kalman filter methods [14][15][16]. As a follow-up to their previous work, a detection method with different controllers is used in [17][18][19] to estimate the inertia constant; however, they are only applicable to the fractional derivative generation model. While model-based methods have been used widely to estimate inertia, the quality of the estimation is contingent on how well the model has been constructed and tested. They may also be unsuccessful in estimating the virtual inertia of the system, particularly if the swing equations do not adequately replicate the system behavior and control techniques have not yet been specified. Take note that the computation of virtual inertia for non-synchronous devices is performed using the frequency divider formula [22] in [20,21], which requires the admittance matrix and interior reactance of the system to be determined using Thevenin equivalent as the only inputs. But all power transmission network parameters and transformer data are constantly prone to variations as a result of unknown variables. This causes the network admittance matrix to be less exact than it should be. By analyzing the frequency of the power system's transmissions, the European Network of Transmission System Operators for Electricity develops a method that can be used in the Nordic power system to estimate the amount of kinetic energy present during disturbances, which is described in [23]. It has turned out to be a considerably more challenging task than initially anticipated, owing to the fact that the frequency might behave quite differently in different portions of the system, depending on the operational circumstances and location of the fault, making it difficult to predict where a fault would occur. For estimating the overall kinetic energy in the Nordic power system, as well as its distribution between different price zones, it was discovered that online measurements were an effective technique of gathering data. The information gathered can be utilized to track the evolution of kinetic energy through time, which can then be used for future studies of the data. Several factors affect the accuracy of the inertia constants in this evaluation. For example, generator improvements, particularly for older units, may affect the accuracy of the inertia constants in this assessment. Because of concerns such as the validity of the inertia constants owing to the installation and alteration of power generators, particularly for older units, and other factors, it is difficult to estimate the accuracy of the inertia constants. A proposal to measure the inertia of the GB Power System using Synchro phasor Measurements is presented in [ [24]] , and the method is detailed in greater depth. According to the proposed technique, a suitable event for analysis is first identified for a specific section of the global positioning system network, and then the observed transients are filtered in order to produce an accurate estimate of inertia for a specific region of the network.Despite the various methodologies that have been offered, it is not possible to estimate the system inertia constant over a long period of time since they do not take into account the integration of renewable resources. This paper offers a convenient approach for estimating the inertia constant for future power systems, which is mostly based on the future energy scenarios that have been developed.

Research elaboration
The inertia constant for 2025, 2035 and 2050 is estimated for the Rwandan's power system during basic, medium, and high progression scenarios. The frequency response is also investigated following system disturbances in the network. This studied network consists of 42 interconnected generators, thermal, hydro, and solar power plants, that are controlled through three coherent control areas: North (Area 1), South (Area 2), and Western zones (Area 3), respectively. This research is expected to contribute in power system planning for long term in terms of frequency response dynamic. It is necessary to know the amount of kinetic energy present in the power system in order to investigate many frequency-related issues. As a result, knowing the amount of kinetic energy in the system is important in order to operate the system as safely and efficiently as possible. The goal of the This research is expected to contribute in increasing understanding of power system behavior as it relates to system inertia in the future by proposing a method to determine the current level of kinetic energy in the power system, as well as the levels of kinetic energy that can be expected in the immediate future.

System Inertia Estimation for current power system
Electricity According to [25], the frequency response of a power system can be estimated by examining the electromechanical dynamics of the rotor of a generator (i.e., how the electrical and mechanical part of the rotor interact with each other) [25]. The dynamics of the motion of the rotor of a single generator (G i ) can be described by the 'Swing Equation of Synchronous Generators' represented in Eqn (1) .
Where J i is the combined moment of inertia of the prime mover (turbine and shaft) and generator in kg.m 2 , ω mi is the mechanical rotating speed of the rotor in rad/s, and T mi and T ei are the mechanical and electromagnetic torque in Nm, respectively. Equation (2) is ob- The quantity Jω m can be denoted by the symbol 'M' and represents the inertia constant. It is related to kinetic energy of rotating masses, to be denoted by W k , as shown in Eqn (3).
Substituting (3) into (2) yields the following expression in Eqn (4) : Since ω m does not change by large amount before stability is lost, M is evaluated at the synchronous speed and is assumed to remain as a constant value [26], [27]. According to different studies in the literature, power system analysis is generally expressed as in per unit quantities. The swing equation can therefore be expressed in per unit by dividing it by the base power S B , which results in Eqn (5) as follows: The per unit inertia constant (H) can now be expressed as in Eqn (6): Substituting H in (6) into (5) yields the following expression in Eqn (7): The power system network considering North, South and Western zones with all generators is shown in Fig. 1. These generators respond to sudden changes in demand simultaneously and therefore, it is important to study the electromechanical dynamics of the whole power system where the total inertia constant of the system needs to be derived and determined.
In general, the total inertia constant of a multimachine power system is calculated after choosing a common system base [26], [27].
If H i is the inertia constant of machine 'i' expressed on the common MVA base S B , and H Gi is the inertia constant of machine 'i' expressed on the machine rated MVA (S Gi ), then the expression presented in Eqn (8) can be obtained as in the following.
It is suggested in [23], [28],that S B refers to the total of the rated apparent power of an individual generator, as shown in Eqn (9) as follows: Equations (8) and (9) are hence extended to 'N' number of generators swinging coherently to determine their equivalent inertia constant (H) as shown in Eqn (10).
For some generators, according to Cummins Generator Technologies in [29], [30], [31], the inertia constant of an individual generator (H Gi ) is calculated using the parameter 'WR 2 ' (moment of inertia) in kg.m 2 as shown in Eqn (11). The value of this parameter is taken from the manufacturer's published data sheet and is used with a constant 'c' that considers the running speed and a factor to deal with the kg.m 2 units of 'WR 2 '.
The constant 'c' has a base value of 49.3 based on a 2-pole alternator running at 3000rpm. By using the collected data for each power plant, the current inertia constant for Rwanda's power system is calculated. As the system under consideration is a three-area interconnected system with hydro, solar, and thermal power plants. Table I presents the generation data where the average inertia for each generation type in each control area is determined and calculated using the expressions derived in Eqn (1) to Eqn (11).

Load Frequency Control Model
After determining the total system inertia in Table I, a frequency response model for the simplified system in Fig. 5 with all generators of the same zone and same type aggregated as a single generator unit is studied. The test model in MATLAB Simulink represents a typical Rwandan power system as a three-area, hydro, thermal, and solar system connected by tie-lines as shown in Fig. 2. In addition to data presented in Table I, further calculations and assumptions were made to obtain other relevant network base parameters. These parameters are presented in Tables II to VIII in the following sub-sections.

Governor Speed Regulation, system Stiffness and Load Damping Determination
In multi-generator system, if all generators are assumed to operate at the same synchronous speed, the equivalent generator model representation in [27] can be used. The equivalent load-damping constant (D eq ) is calculated in Eqn (12).
Where n is the number of generators and D i is the individual load response to the frequency deviation for each generator. The composite power/frequency characteristic of a power system depends on the combined effect of the droops of all generator speed governors and the frequency characteristics of all loads in the system. By using this composite power/frequency characteristic, the load damping constant and the governor speed regulation can be determined [26], [27], [32]. For a system with 'n' generators, the load-damping constant (D) is expressed as a percentage change in load divided by percentage change in frequency, the steady state frequency deviation (Δf SS ) following a load change (ΔP L ) is given by Eqn (13).
With R i expressed in per unit represents the individual governor speed regulation and therefore, the composite frequency response characteristic of the system can be obtained in Eqn (14).
Where the equivalent governor speed regulation is determined as in Eqn (15): The composite frequency response characteristic β is normally expressed in MW/Hz and is referred to as the stiffness of the system. The composite power/frequency characteristic of a power system is also illustrated in Fig. 3. An increase of system load by ΔP L (at normal frequency) results in a total generation increase of ΔP G due to governor action and a total system load reduction of ΔP D due to its frequency  sensitive characteristic [25], [27], [33]. The load damping factor of the power system under study is determined based on the data showing the load and frequency variations in different time intervals collected from different control areas. The load and frequency changes collected values are presented in Table II. A sample of five different events in each substation of areas was collected and the average of each change is calculated and shown in Table II.  Table II also presents the calculated load damping factor in each area. The calculated average load damping factors at North, South and West zones are approximated as 1.06, 1.37, and 1.11, respectively.
The system under study is modelled with supplementary control   . 3. Composite governor and load characteristic modes. It is in this regards that both the stiffness of the system and the governor speed control constant are calculated. The governor speed regulation constant (R), and the system stiffness (β) are calculated using the relationships represented in Fig. 3. The calculated values for R and β for generation and installed capacity are shown in Table III. As discussed, the network is a supplementary controlled interconnected system with three different control areas. To form the basis for supplementary control, the Load Frequency Control (LFC) system in each area of network should control the interchange power and local frequency with other control areas [27], [34,35], [36], [37]. Therefore, the described dynamic LFC system model shown in Fig. 1 is analysed by considering the power signals and the reactance values of tie-lines, as shown in Fig. 4.
In normal operation the tie-line power which flows from control area 1 to control area 2 (P tie, 12 ) is obtained by Eqn (16).
Where X 12 is the tie-line reactance between control areas 1 and 2. δ 1 and δ 2 are the angles of end voltages V 1 and V 2, respectively. For small deviations in the angles Δδ 1 and Δδ 2 , the power in the tie-line changes and this change is calculated, as shown in Eqn (17).
The term |V1‖V2 | X12 cos(δ 1 − δ 2 ) is defined as the synchronising coefficient (T 12 ) of the tie-line and hence its power deviation takes the following form in Eqn (18): The frequency deviation (Δf) is related to the reference angle (Δδ), shown by Eqn (19) as follows: Equation (20) can be obtained by expressing the tie-line power deviations in term of Df, as follows: The Laplace transformation of Eqn (20) yields the following expression in Eqn (21): In general, the incremental tie-line power that flows from control area 'i' to 'j' is given in Eqn (22), shown as follows: In addition to the MATLAB model, the Rwandan power system was also modelled and simulated in IPSA+ Power tool by conducting a loadflow to obtain the voltage and angle in each control area. Table IV presents the voltage and angle with the calculated synchronising coefficient in each tie-line and area, respectively.
As the data to calculate the turbine and governor time constants are not available, the turbine and governor time constants are selected based on the standards presented in the IEEE Power & Energy Society technical report (PES-TR1) [38], [39]. The selected parameters and the results of calculations for each control area are presented in Table V.
The base power P Base is set to the value of 250 MVA, while the system load change ΔP L , which occurs at time t=0s, is set as 0.2p.u. in each control area. The definition of symbols employed in Table IV is shown in  Table VI.

Inertia constant estimation for future power system
The novel method to estimate the inertia e for future power system, is  mainly based on the produced future energy scenarios. By using this method, the hydro, and thermal inertia constant can be estimated according to the available capacity for a given year. The total inertia for a control area j in the year i, (H i j ) is estimated using Equation (23). Based on quantitative analysis, this paper proposes an approach for estimating the total system inertia by combining historical and planned installed capacity data with projections for the integration of renewable energy resources. An actual power system's historical energy and renewable energy integration percentage data. The paper develops and evaluates a mathematical model for long-term estimation of system inertia and then its effect on frequency responses. By estimating and reducing the total error in inertia over a certain period of time, the model parameters are found. For error minimization, the sequential quadratic programming (SQP) approach is used in conjunction with a constrained optimization technique [40,41]. The inertia estimation model, is represented using the modified functional-link net (MFLN) [42]. In the expression Eqn (22a), a base inertia constant (H 0 ) is calculated which is again used to find the inertia constant in each power system control area for a given year.
With The base inertia constant presented in Eqn (22a) is then used to calculate the inertia constant of a given control area as in Eqn (23).
S h j , hydro power installed capacity for the control area in MVA; S th j , Thermal power installed capacity for the control area in MVA; S nh j , Total non-inertia installed capacity for the control area in MVA

Frequency Response for Current Power System
Simulation results for frequency response are shown in Fig. 10. With the increase in load of 0.2pu in each control area, the frequency fluctuates and the power outputs of ΔP m1 , ΔP m2 , and ΔP m3 increases providing their primary frequency response at 49.88Hz. This response was enough to supply the load and the frequency stop dropping after around t = 8s (see Fig. 5) for the whole system to be synchronised.
f 1, f 2, and f 3 , represents frequency values for control areas 1, 2 and 3, respectively. The frequency in control area 3 (f 3 ) is synchronised to 49.88Hz after around t = 4s of the disturbance while the same occurred after around t = 6s and t = 8s in control areas 1 and 2, respectively. It is seen that the synchronisation in control areas 2 and 3 occurred faster than in area 1. The slow response for area 1 is due to the contribution and connection of PV generation. This new generation synchronizes the system at the new operating frequency of 49.88Hz (see Fig. 6). According to control loops shown in [27], the change in mechanical power in each area due to governor action is given by ΔP m1 = − ΔF 1 /R 1 and ΔP m2 = − ΔF 2 /R 2 . Thus, the Northern zone control area increases the generation by 19.6MW while the generation increases by 14.3MW and 15.625MW in control areas 2 and 3, respectively.
The total change in generation is 49.525MW (19.6 + 14.3 + 15.25), which is 0.475MW less than 50MW load change. This happens since the frequency is not synchronised with the nominal frequency of 50Hz during the simulation.

Frequency Response Study for Future Power System
The Rwandan three-area system model was developed in MATLAB Simulink environment. The three areas are assumed to be operating in parallel at the nominal frequency of 50Hz. The simulation considered the existing power system following the disturbance applied at t = 2s (increase in load of 50MW). Next subsection presents the future inertia constant estimated for the Rwandan power system under basic, medium, and high progression scenarios, respectively. To conduct simulation cases, the estimated variations of PV and imports during basic, medium, and high progression scenarios have been taken into considerations, and all losses including connection, wiring, and other losses for both PV and interconnectors are neglected throughout the simulation.

Basic progression scenario
The simulation investigates the behaviour of the network model when PV generators and interconnection are considered. Since solar PV and imported power do not provide inertia for frequency regulation, then the overall system inertia will be affected. This reduction in the inertia in each area is computed by using the same expression shown in Eqn (10). The change in the equivalent inertia for basic progression scenario is presented in Table VII.
With a variation in load of 50MW at t = 2s, the frequency deviation response in each area drops from the nominal value of 50Hz as shown in Fig. 7 (control area 1), Fig. 8 (control area 2), and Fig. 9 (control area 3), respectively. This drop in frequency is due to the increase in the supplied load and the corresponding increase in PV generation and power import.   Figure 7 shows the frequency response in the control area 1 with PV and import penetration for case 1, case 2, and case 3 representing the years of 2025, 2035, and 2050, respectively. The frequency drop was highest at 30% PV and import penetration because at this point the value of the equivalent inertia had reduced to 4.29s. As the safe operating frequency is 49.8Hz, at 30% PV and import penetration in the control area 1, the system reaches the operating frequency that is below the safe operating limit of 49.8Hz. In addition, for cases 1 and 2, the frequency response dropped from nominal value; however, these values were recorded to be within safe operating limits of the power system. Therefore, to avoid the adverse effect of severe frequency deviation, a suitable control strategy must be applied to assist in load-generation balance to maintain frequency at safe operation.
As observed in control area 1 for the projected inertia constants in years 2025, 2035, and 2050, the highest frequency drop occurs at 28% penetration of PV and import. However, as the penetration rate is lower compared to the control area 1, the frequency drop is also low (approximately 0.03Hz difference). The frequency response characteristic in control area 2 with PV and import penetration for three cases is shown in Fig. 8.
It is observed that the frequency in control area 2 for cases 2 and 3 returns to its steady state after around t = 3s after the disturbance while the same occurred after around t = 6s in control area 2 for case 1. It is seen that the synchronisation in this control area in 2035 and 2050 occurred faster than in 2025. However, both cases 2 and 3 show a higher frequency drop with a lower operating frequency due to the reduction in the inertia constant. The frequency deviations for the third area are illustrated in Fig. 9, which clearly shows the maximum overshoot of around 0.48Hz in the year 2025, 0.54Hz for the year 2035, and 0.7Hz for the year 2050, i.e., 0.96, 1.08 and 1.4% in frequency drop and a settling time of 5 to 6 seconds for a step change in load after t = 2s.
Case 3 (with 16% PV and import penetration) shows the highest value of maximum frequency deviation, indicating that the reduction in inertia causes the frequency to fall below its operating value (49.64Hz).

Medium progression scenario
The progressions in PV and imports penetration and the resulted reduction in the inertia constant are presented in Table VIII, and the frequency response characteristics for each year and each area are shown in Fig. 10.
In this scenario, the increasing penetration of renewable technologies and thus a reduced inertia constant were observed. For example, as it is estimated, the PV solar installed capacity and imports were projected from 25 to 37% from the year 2025 to 2050 in the Northern region of the country and this reduced the inertia constant from 5 to 3.6s (i.e.,

28% reduction).
It is seen from the results in Fig. 10 that as the integration of renewable energy resources is increased in the system, the frequency deviation also increases slightly. It has been observed that for higher PV and import penetration, the peak undershoots and settling time increases while the operating frequency reduces to lower values compared to the normal operating frequency. At the same time, the system experiences larger frequency deviation oscillations (see bottom results of Fig. 10) for the year 2050 due to higher penetration rate of renewables, whereas the system frequency oscillations are minimum for the year 2025 due to lower integration of PV.

High progression scenario
In this scenario, simulation results in Fig. 11 rather show a higher peak below the operating frequency for each control area compared to other progression scenarios. In addition, in the year 2050, the system experiences larger frequency fluctuations and moves towards instability state as the system restoration time is more than 20s. This is because the highest progression in renewable energy resources penetration resulted to a larger reduction in the calculated inertia constant (from 7.2 in control area 1 to 3.83s in control area 3). By comparing the results in each year, the trajectories of deviation and settling time are almost the same in years 2025 and 2035, whereas overshoots in frequency deviations observed at each time instant are almost the same in control areas 1 and 3 due to the connection and operation of thermal power plants (see Fig. 11).

Conclusion
This study primarily focused on the integration of low carbon technologies and reviewed their impact on Rwanda's power system inertia by analysing the behaviour of the system with dynamic frequency response in MATLAB Simulink. The work conducted in this paper aims and hopes to contribute to frequency regulation planning for future studies. Three distinct scenarios (basic, medium, and high progression) with three cases were carried out to determine the ability of the existing power system to cope with PV and import penetration for the years of 2025, 2035 and 2050, respectively.
The results show that the existing conventional generators could not absorb variations in electricity outputs from PV and import penetration, which as a result may be causing the frequency to deviate beyond its safe operating limits. The calculated inertia constant values may reduce from 7.2s to around 3.83s during high progression scenario which results in frequency deviation to be around 46.5Hz (which is roughly 7% below its normal operating frequency of 50Hz), and the system experiences larger fluctuations and oscillations in the year 2050 due to higher integration of renewable resources and hence, there is also a risk of system instability. With the increased integration of PV and import in Rwanda's power systems, further simulation results showed a reduction in the equivalent inertia constant in each control area. Furthermore, with a variation in 50MW load at the time t = 2s, the frequency deviation response in each area and scenario dropped below the nominal operating value of 50Hz.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.