Optimal Allocation of Renewable Energy Hybrid Distributed Generations for Small-Signal Stability Enhancement

This paper solves the allocation planning problem of integrating large scale renewable energy hybrid distributed generations and capacitor banks into the distribution systems. Extraordinarily, the integration of renewable energy hybrid distributed generations such as solar photovoltaic, wind, and biomass takes into consideration the impact assessment of variable generations from PV and wind on the distribution networks’ long term dynamic voltage and small-signal stabilities. Unlike other renewable distributed generations, the variability of power from solar PV and wind generations causes small-signal instabilities if they are sub-optimally allocated in the distribution network. Hence, the variables related to small-signal stability are included and constrained in the model, unlike what is obtainable in the current works on the planning of optimal allocation of renewable distributed generations. Thus, the model is motivated to maximize the penetration of renewable powers by minimizing the net present value of total cost, which includes investment, maintenance, energy, and emission costs. Consequently, the optimization problem is formulated as a stochastic mixed integer linear program, which ensures limited convergence to optimality. Numerical results of the proposed model demonstrate a significant reduction in electricity and emission costs, enhancement of system dynamic voltage and small-signal stabilities, as well as improvement in welfare costs and environmental goodness.


Introduction
Distributed generations (DGs) are modular power generating systems that are placed in the distribution systems proximate to consumption centers in order to satisfy immediate power needs, defer investment on transmission and distribution upgrade and expansion, reduce production and welfare costs, reduce losses, diversify energy resources, improve system reliability, power quality, and network stabilities, and many other power system benefits accrue from it [1,2]. The technologies adopted in DG integrations may be renewable based (solar photovoltaic (PV), wind, biomass, fuel cells, etc.), or non-renewable based (internal combustion engines, etc.), or a hybrid of both. As renewable energy resources (RES) hybrid, DG provides a sustainable option due to its infinitude quantities, complementarity, environmental goodness, technological advancement, and economic profitability. However, inappropriate type, suboptimal sizing, and improper location of renewable energy hybrid distributed generations (REHDGs) such as PV and wind cause small-signal oscillations in the distribution systems (DSs). This is due to the variability of output power generated and injected from these REHDG, which totally depends on intermittent solar radiation, temperature, and wind speed. Several research works have agreed that the variability of power generated from intermittent renewable resources relative to load or vice versa results in power system oscillations [3][4][5][6][7][8]. That is, power imbalance between the total generations from the REHDGs, other plants and transmission feeder supply, and power demand aggregate at a time results in small-signal instabilities in the distribution system [9,10]. At large scale integration levels, the small-signal instabilities due to the variability of REHDGs power can have significant effects on the distribution networks (DNs) [5,7,10]. All the aforementioned concerns make the formulation and optimization of REHDGs' optimal allocation problems a task to solve with simple mathematical models. To obtain a realistic model, it is very important to represent the network as a dynamic model, use multiple periods for the planning horizon, and include all the pivotal constraints. Therefore, the planning and design of the optimal allocation of renewable energy hybrid DGs (REHDGs) requires a serious consideration of the type of network configurations, types of DG technologies, their intermittency modeling, and the number, capacities, and locations of the units. This is to achieve minimum total costs while the dynamic small-signal stability requirements of the network are simultaneously met and appropriately assessed.
There have been several diverse objectives continuously employed by many researchers on the optimal allocation of DGs in DSs. Some of these objectives are minimization of system losses and enhancement of voltage profile, minimization of line loss, minimization of investment and operation costs, minimization of total penalties for system loss compensations, minimization of renewable DG penetration, maximization of DG capacity, maximization of system reliability, and so on. Minimization of the NPV of the total cost as an objective function is very common in the planning of optimal allocation of REHDGs. Analytical based methods were proposed in [11][12][13] for the planning and operation of the optimal locations and sizes of DGs in DSs. Numerical based methods whose algorithms find numerical solutions for different DG allocation problems were employed in [14][15][16][17][18][19][20]. Linear programming (LP) [14,15], mixed integer non-linear programming (MINLP) [16,17], quadratic programming (QP) [18,19], and optimal power flow (OPF) [20] are some of the common numerical methods applied in solving DGs' allocation problem. Intelligent search (IS) based methods are differently employed to solve the optimal sizing and placement of DGs problems. IS methods utilize artificial intelligence (AI) algorithms like the genetic algorithm (GA) [21,22], particle swarm optimization (PSO) [23,24], simulated annealing (SA) [25][26][27], harmony search (HS) [28,29], big bang crunch (BBC) [30,31], the fireworks algorithm (FA) [32,33], and the water drop algorithm (WDA) [34,35].
However, some current works use mixed integer linear programming (MILP) [36][37][38][39] for the mathematical formulation of the allocation model due to its ability to detail the physics and mechanics of the problems and find global optimal solutions with less computational requirement. In Munoz Delgado et al., MILP was used in the expansion planning of DS to minimize the NPV of total cost of investment, maintenance, production, power losses, and unserved power [36]. A multi-stage and stochastic mathematical model, formulated as a mixed integer linear programming (MILP) problem, was employed in [37] to find the optimal time, sizes, and placement of renewable DGs, compensators, and energy storage systems with a view toward minimizing the NPV of the total cost. In [38], a chance constrained stochastic MILP model was formulated for determining optimal decisions in DGs' investments with operational uncertainties modeling. The model was further optimized with a variant of the evolutionary method, the vertical sequencing protocol, to minimize the total cost of investment and operation. The authors in [39] proposed MILP to optimize DG capacity hosting of a radial distribution network through the reconfiguration of existing tie and smart switches with the objective to maximize total DG capacity deployed into the network at a minimized cost. All the works discussed here evaluated static voltage stability, and few modeled the uncertainties of renewable resources, but were unable to assess the effects of the variability of renewable powers injected on the networks' dynamic voltage and small-signal stabilities. Similarly, the global optimality of their solutions was not reported.
Despite several research studies in the areas of optimal sizing and placement of distributed generation problems, most of these works considered only the optimal location and sizing of a single DG unit, and those that were multiple DGs were mostly conventional sources. Most of the works did not evaluate the network stability, and those that did only evaluated static voltage stability, but not even a dynamic one. They did not evaluate long term dynamic small-signal stability, which is the prerequisite for a power system to be in operation in the first instance. Further improvements introduced in this study include the modeling of uncertainties in REHDGs' allocation expansion planning model to account for the renewable resources' intermittency, the usage of time varying load demands in a dynamic distribution network, and the use of dynamic model of the planning horizon to achieve optimal long term planning where capacity is added as and when needed. The planning of the optimal sizing, timing, and placement of renewable DG units to attain minimum cost and enhance the small-signal stability level during integration in the distribution network is the strength of this paper. Hence, its contributions are as follows: • To the authors' knowledge, no literature has ever presented an assessment of long time dynamic small-signal stability in the planning optimization of REHDGs' allocation in distribution networks.

•
Unlike previous studies, this work incorporates the variables necessary for the network stability and includes pivotal constraints that are necessary in REHDGs' optimization for enhancing the long time dynamic small-signal stability of distribution networks.

•
In this study, dynamic planning is employed for the planning horizon, as opposed to the static model usually applied in most research works. The planning horizon is split into various time periods, which are comprised of a specific number of years and sub-years.

•
Unlike most existing research works, this work models the uncertainty of intermittent RES and implements time varying loads in a dynamic model of the distribution network.
In this paper, a mixed integer linear programming algorithm is proposed as done in [37], to find the optimal time, sizes, and locations of multiple REHDGs in the distribution network systems (DNS). The typical syntax of a mixed integer linear programming formulation is presented as follows: Subject to: where F m (u, v) is the objective function, u and vare the vectors of continuous and discrete variables, respectively, and Expressions (2) and (3) are the inequality and equality constraints, respectively. Consequently, an allocation optimization model that determines the optimal sizes, time, and locations of renewable power units (solar PV, wind, and sugar-cane biomass) and capacitor banks (CB) and constrains variables that are responsible for network stability is developed. The optimal integration of these technologies is to essentially meet the objectives of maximizing the RES power generated and injected into the DNS and enhancing the small-signal stability (SSS) of the distribution systems. The resulting optimization model is formulated as a mixed integer linear mathematical program. In addition, the principle of fast decoupled power flow is applied in the linearization of traditional non-linear AC network model.
The remainder of the paper is arranged as follows: The mathematical modeling of renewable resources and load is presented in Section 2. The output power of PV modules and wind turbines that can be harnessed during the planning time is calculated in Section 3. Section 4 formulates the optimization model for solving the RES allocation planning problem. Section 5 theorizes about the power network model, renewable generators dynamic models, and the eigenvalue analysis approach necessary for evaluating the SSS of the proposed model. The discussion and analysis of the results from the case studies used for the validation of the proposed model are presented in Section 6. Finally, the main conclusions of this paper are drawn in Section 7.

Modeling of Renewable Energy Resources and Load
This section presents the mathematical formulation of models for renewable energy resources and power demands (load). The annual data for solar irradiance and ambient temperature and wind speed used for the estimation of solar PV and wind generators are the measurements collected in the KwaZulu-Natal region of South Africa through the South Africa National Energy Development Institute (SANEDI) and Wind Atlas of South Africa (WASA), respectively. These resources are evenly distributed over a range of distance in this region. The site under study was equally chosen because of the large scale cultivation of sugar-cane, whose waste is an economic quantity for large scale generation of electricity [40].

Processing of Climatic Data
Five years of climatic data of these renewable resources collected from the studied region were processed using the beta and Weibull probability distribution function for estimating hourly solar irradiance and wind speed, respectively. The respective years were grouped into seasons: summer, autumn, winter, and spring. A month was taken to be thirty (30) days, of which each seasonal segment had 450 solar irradiance and wind speed data points, i.e., 5 years × 3 months a season × 30 days a month. From the seasonal data, the shaping (a and k) and scaling (b and c) parameters of the beta and Weibull distribution functions were respectively generated for each seasonal segment using MATLAB function "dfittool". This was done to redistribute the raw irradiance and wind speed data to follow the beta and Weibull distributions, respectively. Thereafter, one year was selected as the study period and split into seasons, of which each season was described by a day in that season segment. The scaling and shaping parameters mentioned were then applied to generate the frequency distribution of irradiance and wind speed data of a typical day for each season, respectively. Thus, the beta and Weibull distribution functions of each hour were generated using (5) and (6), respectively. For clarity, the typical day of a season was the day whose 24 hourly values gave the lowest standard deviation as corresponded to the 24 hourly average values of such a season. The typical day was further divided into 24 h segments, which referred to particular hourly periods of the whole season. That is, a whole year had 96 h periods for 24 h a day per season.

Modeling of Solar Irradiance
The solar irradiance data usually had two modes of distribution functions for the same hour of a typical day of each season. The beta PDF parameters obtained from the seasonal segments were used to obtain a typical day's hourly behavior of solar irradiance data as in the equation below: where: f B (I) = the beta distribution function of I I = solar irradiance a = the scaling parameter of the beta distribution function b = the shaping parameter of the beta distribution function Γ = the Gamma function

Modeling of Wind Speed
The wind speed of the site under study followed a Weibull probability distribution, and its behavior was modeled with the Weibull PDF f W (v) as in (6). This model was acceptable based on the good results obtained from the Anderson and Darling (AD) goodness of fit test.
where: k = the shaping parameter of the Weibull distribution function λ = the scaling parameter of the Weibull distribution function v = wind speed in m/s

Modeling of Sugar-Cane Biomass
The annual electricity output E anl delivered by the biomass power generation system with the rated electricity generator depended on the capacity utilization factor, CUF, as given below: Then, the hourly electricity output of sugar-cane waste/shaft was modeled with this equation: where: E bio,h = the output energy of the biomass/electricity generator per hour P bio,h = the rated power of the biomass/electricity generator per hour η d f g = the efficiency of the biomass/electricity generator

Modeling of Load Demand
The load profiles used in this study were computed by down-scaling the actual annual load profile of South Africa [41] to the IEEE-14 bus and -118 bus systems and implemented as case studies. They presented the hourly load demand level for all the planning years.

Calculation of PV Generator Output Power
The PV module output power corresponding to each state depended on the solar irradiance and temperature of the site, as well as the module characteristics. After generating beta PDF for a specific 24 h period of a typical day, the output power of a PV module during various states was estimated from the PV module characteristics' power performance curve for the hourly time segment using (9).
Note that the average value of the solar irradiance of a specific state (I ap ) was used to extract the state output power from the module power curve.
where T cp is the temperature of the PV cells ( • C) at state p; T A is the ambient temperature of the site ( • C); K c is the current temperature coefficient (A/ • C); K v is the voltage temperature coefficient (V/ • C); NOT pv is the nominal operating temperature of PV cell ( • C); I sc is the short circuit current (A); V oc is the open circuit voltage (V); P pv,h (I ap ) is the output power of the PV module at state p; I ap is the average solar irradiance of state p; FF is the fill factor (-); I mppt is the current at the maximum power point (A); V mppt is the voltage at the maximum power point (V); and P rs is the solar panel rated power (MW). Availability factor P pv,h (I ap )/P rs (parameter AFp p,h ) is the output power of all PV modules at the ith bus. The parameter AFp p,h was computed using the hourly irradiance data of each state. The calculation of the AFp p,h parameter was global, as the parameter was time dependent and affected all the PV modules equally.
The number (n p,i ) of rated PV generators (P rs ) per bus is one of the decision variables for the problem. Multiplying the number of PV generators (n p,i ) and the rated PV generator (P rs ) with the parameter AFp p,h in the power balance equation gave the actual power output of the solar based DG units.

Calculation of Wind Generator Output Power
The wind turbine output power corresponding to each state depended on the wind speed of the site under study. After generating the Weibull PDF for a specific 24 h period of a typical day, the output power of the wind turbine during various states was determined from the wind turbine power curve using (10). The average value of the wind speed V aw of a specific state was used to extract that state output power.
where V ci , V rw , and V co are the wind turbine cut-in, rated, and cut-off speeds, respectively. P w,h (WS aw ) is the wind turbine output power at state w. V aw is the average wind speed of state w, and P rw is the wind turbine rated power (MW). Availability factor P w,h (WS aw )/P rw (parameter AFw w,h ) is the wind turbines' output power at the ith bus. The parameter AFw w,h was computed using the hourly wind speed data of each state. This calculation was global, as the parameter was time dependent and affected all the wind turbines equally.
The number (n w,i ) of rated wind turbines (P rw ) per bus is one of the decision variables for the problem. Multiplying the number of wind turbines (n w,i ) and the rated wind turbines (P rw ) with the parameter AFw w,h in the power balance equation gave the actual power output of the wind based DG units.

Mathematical Formulation of the Planning Problem
This section presents the proposed multi-stage optimization model that finds the optimal sizes, number, time, and locations of REHDGs, simultaneously focusing on solar PV, wind, and biomass (sugar-cane waste) DGs and capacitor banks. The major objective of the proposed model was to improve the small-signal stability and maximize the REHDG power generated and absorbed into DNS at a minimized cost. The model was formulated as a stochastic mixed integer linear programming (MILP) optimization problem in a MATLAB environment whose syntax is presented in Section 1. Furthermore, a linearized alternating current (AC) network model based on the fast decoupled power flow (FDPF) was used in the formulation of the problem in order to capture the characteristics of the network system effectively.

The Objective Function
The objective of this planning formulation was to maximize renewable power integration into the DNS from the distribution system operator's (DSO) perspective by optimally allocating solar PV, wind, and biomass DGs and capacitor banks and constraining the optimization variables related to small-signal stability to within the required level at the least cost. The objective function of the formulated MILP planning stage optimization problem was to minimize NPV of the total cost (11), subject to the linear constraints stated in Section 4.2. Minimize, The first term in (11), the cost term C I t , is the total investment cost amortized in annual installments throughout the lifespan of the installed components, considering there is reinvestment in an identical piece of component after the component's lifespan expiration, as done in [36,37]. This cost valuation was done using the principle of an infinite or perpetual planning horizon [42]. In this study, the overall investment cost was taken as the summation of investment costs on existing and new DGs and capacitor banks, as given in (12).
The capital recovery factor d(1+d) LT (1+d) LT −1 was used to weight all the investment costs to return interest on capital invested for all the components. The capital recovery factors for the components were where LT is the life time of each component (generators and capacitor banks) and d is the interest rate on investment. The formulations of investment variables (x −,i,t − x −,i,t−1 ) included in (12) enforced that the addition of the investment cost of a component in the summation was done just one time. This implies that such component(s) can only be put to utilization at the beginning of a planning year and not midway into the planning year, even if it is bought mid-year.
The second term in (11) corresponds to the operation and welfare costs throughout the time stages. This term consists of three cost terms vis-à-vis: total maintenance cost (C M t ), total energy cost (C E t ), and total emission cost (C X t ). Equation (13) models the maintenance cost, C M t , which is the sum of the respective maintenance costs of existing and new DGs (MC E g and MC N g ) and capacitor banks (MC E cb and MC N cb ) at each planning stage under the principle of a perpetual planning horizon. The cost term C E t denotes the total energy cost in the system based also on a perpetual planning horizon. C E t is the sum of the costs of power generated by the existing and new DGs and the purchase from the utility for each stage as modeled in (14). The costs of power generated by the existing and new DGs are the multiplication of the cost of unit energy produced (i.e., expected cost of operation) and the amount of power produced. The cost term C X t sums up the emission costs associated with existing and new DGs and the power from the utility feeders, as characterized in (15). The emission cost is the expected costs of emissions based on the power produced from the existing and new DGs and purchased from the utility, respectively. Here, the cost function of emission was assumed to be linear for the sake of simplicity. The cost function, in reality, was highly non-linear and non-convex [43].
where γ CO 2 e s,t , ER g , and OC g,i,s,h,t are the penalty for emissions ($/tCO 2 e), emission rate (tCO 2 e/MWh), and operation costs ($) for the respective components.
The third term in (11) presents the net present value of the costs incurred for production (maintenance and energy costs) and emissions (welfare) after the last planning stage. This cost is also referred to as the end effect, accounting for the residue values of the invested components. It should be noted that this term is also estimated based on the principle of a perpetual planning horizon and depends on the operation and emission costs of the last time stage.

The Constraints
Constraints were applied on the REHDGs allocation problem to exert restrictions on the optimization of the objective function(s) with respect to the decision variables. The constraints used in the optimal REHDG allocation problem formulations were as below: 1. Complex power flow: The network AC power flow commonly known as Kirchhoff's law of voltage is highly non-linear and non-convex. The equations representing it are given below.
The principle of the fast decoupled power flow (FDPF) model, postulated in [44], was applied to linearize (16) and (17) to yield (18) and (19). The expression for voltage magnitude at bus i is given as in (20).
where ∆V min Voltage magnitude and angle constraints: The voltage magnitude and angle must be constrained to maintain the voltage and small-signal stabilities of the system.
3. Power flow limits: On any line k, the power flow must be within the specified limits for that line as in (23) and (24).
4. Active and reactive power limits of power from transmission feeders: The power supply from the feeder should have minimum and maximum limits for technical reasons. Inequalities (25) and (26) enforce the limits. where Active and reactive power limits of REHDGs: The output of a unit generator should not exceed its capacity multiplied by the generation binary variables. This ensures that the power generation variable of a generator is zero when it is either unused or un-invested. The capacity at a given year y and state s is considered to be constrained between maximum and minimum values. The active and reactive capacity limits of exiting REHDGs are given in (27) and (28), respectively. Equations (29) and (30) give the corresponding limits for the new REHDGs.
Q Nmin 6. Reactive power limits of the capacitor bank: The reactive power supply from capacitor banks should be bounded. Inequalities (31) and (32) limit the reactive power generated to between zero and maximum capacity.
7. REHDGs and capacitor banks' penetration limits in the DNS: The total penetration limits in the system are determined by the government policy on the average penetration of the renewable DG units and reactive compensators.
where is the maximum penetration limit as a percentage of the load profile and AF g,h,t is the capacity utilization factor of the generators.
8. Active and reactive power balance (Kirchhoff current law): When the distribution network expansion planning includes new renewable generations with supply from the utility substation, then the new network is modeled as (35) and (36). At each node, power balance must be observed. This is achieved with the constraints:

Small-Signal Stability
The small-signal stability (SSS) of a power system is its capability to maintain or regain synchronism during small disturbances. Small-signal instabilities occur in power systems basically as a result of oscillations from the synchronous generators' rotor angle speed, the intermittent nature of renewable generators' power, and/or the market demand (sudden changes in load demand) that shifts the power system operating point into an overstretched state [45,46]. The dynamic small-signal oscillations due to the variability of REHDGs' (solar PV, wind, and biomass) output power became important issues in the present setup of DNS. The SSS issue in today's DN is mostly of poor damping of the system oscillations when variable renewable generations are injected into the network. SSS analysis gives important details on the intrinsic dynamic characteristics of DNSs. Therefore, the dynamic network model incorporating the suitable time-varying load model is a subject of priority any time the dynamic stability of the power network is being evaluated.

Power System Network Model
There are two network models that are commonly used in studying and analyzing power system dynamics and network power flows: the static and dynamic models. Each of these network models has specific areas of implementation, which depends on the kind of problem being considered. Studies related to small-signal stability are dynamic power studies that employ the dynamic model of the power network. Thus, the dynamic network model is described by some non-linear differential algebraic equations (DAEs), whose semi-explicit form is (37). Some of the assumptions made to simplify the mathematics of the network without loss of generality are that: there are no losses and bulk energy storage in the network, so that the total power injected equals the total load demand.
where x and y are the vectors of state and algebraic variables, respectively, and z 1 and z 2 are control inputs and exogenous parameters like load demand.

Dynamic Models of Power Generators
This section develops dynamic models of conventional or synchronous generator (SG) and renewable energy distributed generators (REDGs) such as PV, wind, and biomass generators. These generators are referred to as the components of the network. The first dynamics of physics is assumed in each of these models. In the allocation of DGs, many REDGs can be installed on a bus, while a single SG represents the supply from the transmission feeders and normally assigned to Bus 1 of the distribution network. The aggregation of individual REDGs allocated to a bus can be referred to as a solar PV farm, wind farm, and biomass farm (system). Different aggregated units are allowed to be connected to a single bus. There may be some unit-less buses where no load or generator is connected. Therefore, a general form of dynamic model incorporating differential algebraic equations (DAEs) of various components of a distribution system can be written as in (38). The network component (r) may be an SG, PV farm, wind farm, or biomass system.
x r = f r (x r , y r , z 1r ; γ r ) P r + jQ r = g r (x r , y r , z 2r ; γ r ) (38) where γ r is the rth component model parameter that depends on the operating point, which is zero for the initial operating point. The details of the functions in (38) are presented subsequently. A simplified dynamic model called the classical model of these generators (SG and REDGs) is considered in the analysis of this work. The simulation results from several works [47,48] proved that both full and simplified models of these renewable generators provide relatively the same result in the analysis of network small-signal stability; hence the justification for the use of the simplified models of the generators. The state space descriptions of the SG and REDGS models complying with the expression in (38) are presented in the subsequent subsections.

Dynamic Model of the Synchronous Generator
A synchronous generator is comprised of an exciter system, a prime-mover, and a synchronous machine. The exciter causes magnetization of the rotor when currents are induced in the exciter winding. The turbine converts the kinetic energy of a moving fluid to mechanical power, which rotates the rotor to generate amagnetic field. A synchronous machine transforms mechanical to electrical power, which is injected into the network. There are various models of SG that are being used in power system analysis. However, the classical model or the constant voltage behind the transient reactance (X d ) is adopted in this paper. This model is comprised of differential equations (DEs) of the motion of the generator known as electro-mechanical swing dynamics (39) and (40) and the network algebraic equations called electro-magnetic dynamics in (41)-(57). The mechanical power P m was assumed to be constant for the sake of simplicity. A simplified model of SG in state space is presented as differential algebraic equations (DEAs) below: For generator buses i = 1, . . . , m:α 2H i where the input P mi is the mechanical power applied to the ith generator; E sg,i > 0 is the constant voltage magnitude behind the transient reactance of the generator; V i is the voltage magnitude; θ i is the voltage angle at bus i; M i = 2H i ω s > 0 is the inertia; D i > 0 is the damping coefficient; X di > 0 is the transient reactance; α i > 0 is the rotor angle; and ω i is angular frequency of the ith generator.
The network algebraic equations' part of DAE including power flow due to the SGs' injections into the buses is given as: For generator buses i = 1, . . . , m, where G i,j ≥ 0 is the conductance and B i,j ≥ 0 the susceptance. The terms cos(α i − θ i ) represent the active and reactive power generated by the SG and absorbed into the network.

Dynamic Model of the Photovoltaic Generator
A PV generator of a solar farm is comprised of a semiconductor device, a DC/DC converter (buck and boost), a DC link, and a solid state synchronous DC/AC three phase voltage source inverter (VSI), which is similar to a synchronous machine excluding the revolving components. A PV array is comprised of a parallel connection of n ll circuits with each containing a number (n s ) of PV cells connected in series. A buck and boost converter ensures that the PV cell operates around the maximum power point (MPP) to maximize the output power of the cell. A DC link is a component with a capacitor (and sometimes with inductance) and used to filter ripples that result from the rectification stage. The DC/AC VSI is used for converting DC output voltage into AC output voltage and regulating active and reactive power flowing in the system. The PV generator supplies balanced sinusoidal voltages at fundamental frequency with several switching elements under various controllers. The I-V curve of the PV cell and switching functions are non-linear; thereby, it becomes excessively complex to model the system fully to evaluate its stability. Consequently, a linear function approximation, based on energy conversion principles and ideal switching device assumptions, was used. Hence, the PV generator was modeled using the DAEs in a simplified constant voltage behind the transient reactance model as follows [47][48][49].
For generator buses i = 1, . . . , m, where E pv = n pv V pveq and D swl is the conductance due to switching losses of DC/AC inverters. The network algebraic equations' part of DAE including power flow due to the PV generators' injections into the buses is given as: For generator buses i = 1, . . . , m, where V pv and I pv are the dynamic state variables of the PV generator; V and θ are the algebraic variables that satisfy the power balance equation of the network; n pv is number of PV generators in the farm. The terms cos(α i − θ i ) represent the power generated or absorbed by the PV generators into or from the network.

Dynamic Model of the Wind Generator
A typical wind farm consists of a wind turbine (variable speed) with a doubly fed induction generator (DFIG) and a set of inverters in back-to-back (B-to-B) configuration. A wind turbine transforms the kinetic energy of wind into mechanical power, P m , (47) available for use by the DFIG.
where ρ is the air density (kg/m 2 ); A is the wind turbine swept area (m 2 ); C p is the performance coefficient that depends on the blade tip speed ratio (λ) and pitch angle (θ p ); and V 3 m is wind speed in m/s. A DFIG converts the turbine mechanical power into electrical power. A three phase DFIG is comprised of a rotor and stator in which the rotor is connected to a B-to-B inverter configuration to regulate rotor winding voltage and the stator is connected to the turbine bus to send converted electric power to the grid. The back-to-back inverter is used to regulate generator rotor voltages and reactive power flow between the stator and inverter. The inverter configuration is comprised of two inverters vis-à-vis the rotor side (RSI) and grid side (GSI) inverters. Therefore, the dynamic model of DFIG is split into mechanical and electrical parts. The DAEs that describe the dynamic behavior of the DFIG simplified model in state space representation are as follows: For generator buses i = 1, . . . , m, where E wi is the equivalent voltage magnitude behind the transient reactance of DFIG; θ p is the pitch angle. The network algebraic equations' part of DAE including power flow due to the DFIGs injections into the buses is given as: For generator buses i = 1, . . . , m, where n w is number of wind generators in the farm. The terms cos(α i − θ i ) represent the power generated or absorbed by the DFIGs into the network.

Dynamic Model of the Gas-Turbine Generator
A biomass system is comprised of a gas turbine and a synchronous generator with their control systems. A gas turbine converts combusted hot flue gas into rotatory energy, which is transmitted to the electric generator. The synchronous generator converts the rotatory mechanical energy into electrical energy, which is transmitted to the grid. The DAEs that describe a simplified dynamic model of a biomass turbo-generator in state space representation consist of the power angle equation of the rotor and synchronous speeds, the equation of motion of the turbo-generator, and the network equation as below: For generator buses i = 1, . . . , m:σ where P m = E te + E mag − E cmp − E f ta * t, P e = E gnd * t = 3/2(E tg I q + E tg I d ) and E tg is the biomass generator' voltage behind the transient reactance.
The network algebraic equations' part of DAE including power flow due to the biomass turbo-generators injections into the buses is given as: For generator buses i = 1, . . . , m: where n b is number of biomass generators in the system. The terms cos(α ti − θ i ) represent the power (active and reactive) generated by the biomass generators and absorbed into the network.
The active and reactive power flow in the load buses (i = 1, . . . , m + n) due to all the network generators is given as:

Linearized Network Model
The linearization of a set of non-linear DAEs is the first step in the evaluation of a power system's small-signal stability. It can be accomplished either by performing a power flow analysis or solving the system in (39)-(57) using a numerical method [50]. The linearized compact semi-explicit DAE equation is thereby given as: where A and B are the coefficient matrices of system states and control variables, respectively. Some of the SSS analysis approaches being employed by researchers and utilities are the Fourier transform based Prony and eigenvalue analysis [47]. However, the eigenvalue analysis approach is superior and more common in terms of giving more detailed information and having less computational demand.
The generalized eigenvalues of matrices E and A can now be used for the evaluation of the small-signal stability of the linearized dynamic network.

Eigenvalue Analysis Approach
The eigenvalue is used to analyze and evaluate the small-signal stability of a power system. The computation of all the state matrix's (Matrix A) eigenvalues is a definite way of determining the oscillation performance of a power system operating point. If all the eigenvalues of Matrix A have negative real parts in the S-plane, then the system is said to be completely stable [45,50]. In a stable condition, the oscillation(s) that occurs in state variables owing to small perturbations or imbalances of the system operating point will die out gradually after some time. However, the system is unstable if at least one of the Matrix A' eigenvalues has a positive real part [45,51]. This condition shows that the oscillations that occur in the system increase gradually or remain the same in magnitude over a period of time. Consequently, some stability margin must be observed to limit oscillation occurrence by ensuring that the eigenvalues of the state Matrix A do lie within, but not too close to, the imaginary axis of an S-plane.

Case Studies
This section presents and discusses the results of two different case studies to validate the proposed model. A 14 bus system was first studied over a three year planning horizon to test the developed model. Furthermore, a larger network of a 118 bus system and a ten year planning horizon were subsequently used to validate the scalability of the proposed model. United States dollars ($) represent the currency used in the optimization of this study.
The following data and assumptions were common to both case studies: • The planning horizon was split into yearly decision stages, where 7.0% was considered as the interest rate.

•
For the sake of simplicity, 2% of investment costs was considered for the maintenance costs of the corresponding components.

•
The slack bus (Bus 1) was considered as the substation node whose voltage magnitude and angle were taken to be 1.0 pu and 0 • , respectively.

•
The voltage upper and lower bounds for other buses were 1.5 pu and 0.95 pu, respectively.

•
The yearly demand growth projection of 5% was considered throughout the planning horizon. For the sake of simplicity, the variable power resources (solar and wind) were assumed to be available in every bus based on the fact that the distribution networks covered a small geographical region and that these resources were equally distributed throughout the region. The dispatchable power resource, sugar-cane biomass, was equally available with the widespread cultivation in this study area (KwaZulu-Natal). • The REHDGs' penetration limit, , of 30% was considered, which was well above the target level for South Africa's 2030 renewable projections [53].

•
The capacitor bank with 0.1 MVAr of installed capacity was considered to be the minimum deployable reactive compensator in the system. The unit investment cost for the capacitor bank was $25/kVAr.

•
The emission rate of power supplied from the coal fired utility was taken as 0.4 tCO 2 e/MWh.

•
The emission cost for the first stage was taken as $25/tCO 2 e and $20/tCO 2 e additions for the subsequent years up until the last planning stage.

•
The price of electricity supplied by the substation (utility) was taken to be $0.11/kWh for South Africa.

•
The power factor, pf, of renewable DGs was set to 0.95 lagging, meaning that the DGs was Type II and always absorbed reactive power.
The optimizations, in this study, were implemented on an Acer Veriton with two Intel Core (TM) i5 650 processors at 3.20 GHz and 8 GB of RAM using MATLAB R2019a. An optimality gap of 0.1% was considered as the stopping criterion.

Numerical Results I: IEEE-14 Bus System
The single line diagram and network data of the IEEE-14 bus system can be found in many sources, for example [54]. This system had 9 load buses, 5 substation buses, and 20 branches. The time varying load demand (hourly load demand) of the test system as scaled with the actual South African load profile [41] was used in this study.

Results of the Optimal REHDGs' Allocation Problem
Intermittent renewable generations such as PV and wind based DGs usually operate with fixed lagging power factors [37], meaning such DGs always absorb reactive power from the network. Hence, the power factor of PV and wind DGs in this study was set to 0.95 lagging. However, recently, there have been photovoltaic and wind generator systems emerging that provide reactive power support. Yet, the optimization of renewable energy DGs' integration into DNS presented here was lacking reactive power support from these generations, and the optimization results are discussed below.
Based on the stopping criterion considered, the computation time required to obtain optimal solutions was 18 s. The optimal solutions for REHDGs (solar PV, wind, and biomass) and the capacitor banks are depicted in Tables 1 and 2, respectively. In the first stage, a larger percentage of the investments was made. This was due to the higher NPV of costs for maintenance, energy, and emission at Stage 1 (T1) than in the successive stages. This showed that investing more into REHDGs in the first stage was economically viable since the costs were reduced gradually over the planning horizon.  Biomass  2  3  3  3  Biomass  9  2  3  3  Solar  3  17 25  35  Biomass  10  2  2  2  Wind  3  74 94 115  Biomass  12  1  1  1  Biomass  3  6  6  6  Biomass  13  2  2  2  Biomass  4  5  6  6  Biomass  14  1  1  2   Table 2. Optimal investment solution of capacitor banks for distribution system planning. Stages  T1 T2 T3  T1 T2 T3  x cb,i,t  x cb,i,t   2  6  6  9  9  30  30 30  3  9  11  14  13  5  6  6  6  13 23  33  14  2 2 5 Table 1 shows that more wind DG units were installed than solar PV units despite being given equal parameters of integration. This was because wind generators had a higher capacity utilization factor than photovoltaic generators. Furthermore, the power from biomass (sugar-cane) generation provided supportive services such as spinning reserves to make up the inadequate supply from the intermittent DGs. Biomass power was a fast-response generation that could be maintained constant at any requested amount and time, without uncertainties and within the operating limits [55]. Biomass power presents a viable alternative to power storage devices since an efficient and economical electrical storage system is still being searched. The total capacities of REHDGs power (MW) installed in the three planning stages are shown in Figure 1. The totals of 265.2 MW, 71 MW, and 74.8 MW of renewable power were installed in the network for the first, second, and third stages, respectively.

Stages Bus
From Table 2, it can be deduced that the optimal locations of capacitor banks were mainly on the most loaded buses and the buses close to the end of the network. This is a normal power system operation where capacitor banks are used to compensate for the reactive power deficit in the system, thus helping to enhance network voltage stability by keeping the voltage within the limits. Figure 2 shows the total capacity of the capacitor banks invested in and installed throughout the planning horizon to be 9.7 MVAr, while 6.5 MVAr, 1.3 MVAr, and 1.9 MVAr were installed at each planning stage, respectively. The results from Tables 1 and 2 and captured in Figure 3 demonstrate the complementarity of intermittent renewable generations and biomass generation. Consequent upon these results, the hybrid renewable generators were optimally allocated close to one another and achieved a high optimality gap of 0.1%. The inclusion of reactive compensators greatly increased the capacity of renewable DG units that were integrated into the system by helping to maintain active and reactive power balance when reactive power consuming generators were installed. Ordinarily, the optimal capacity of REHDG units that could have been integrated would have been around 130 MW.  For this particular study, the total load demand of the South African load profile scaled IEEE-14 bus system was 223.1 MWh with 5% yearly demand growth projections throughout the planning horizon. Adding the 84 MW PV farms, 264.5 MW wind farms, and 62.5 MW biomass systems reduced the electricity production from coal fired plants by 30%, amounting to the total NPV investment costs of $20 B, $5.36 B, and $5.64 B for the three planning stages, respectively. This brought the total investment costs to $31 B. The NPV costs of maintenance, energy, and emission throughout the planning horizon were correspondingly equal to $620 M, $163.63 M, and $2.4 B, respectively. The overall total NPV cost for the whole planning stages was $39.74 B. With the integration of 411 MW of REHDGs' power, a total sum of $30.71 M would be saved if the whole power to satisfy the load demand had come only from coal fired conventional generations.

Assessment of Small-Signal Stability
Another important aspect of this REHDGs' allocation optimization analysis is the assessment of the impact of their integration on the long term dynamic voltage and small-signal stabilities of distribution networks. Figures 3-5 show samples of the voltage profiles and eigenvalue plots of the network without renewable DGs (base case), as well as with solar PV DGs, wind DGs, biomass DGs, and all the REHDG units during every operational period throughout the planning horizon.   Considerable improvements in voltage profile are seen in Figure 3. This figure indicates dynamic voltage profiles at any selected bus (Bus 1 in this case) during the operational periods of DNS without and with renewable generations, respectively. From Figure 3 (base case), it can be seen that the voltage levels were very close to the minimum permissible limits. Meanwhile, the voltages (color green) during REHDGs' integration were largely close to the 1.0 pu value with an average deviation of 0.13%. This allowed significant operational margins for the permissible voltage magnitude limits. Invariably, the substantial improvements in voltage profiles were due to the effect of the combined integration of renewable DGs and capacitor banks. Figure 4 shows typical eigenvalue plot results for the base case system. It is observed that the system's margin of oscillatory stability was very low. The system's three state variables were: rotor angle (α), rotor speed (w s ), and transient voltage (E ). The eigenvalues of the base case system lied away from the origin to the left half of the S-plane, except the quadrature transient internal voltage eigenvalues of Generators 1, 2, and 3 that were very near oscillatory instability. The oscillatory frequency and damping ratio were 3.17 Hz and 0.0377, respectively. From Figure 5, it is shown that all the eigenvalues of the component generators were situated away to the left half of the S-plane during the integration of all the renewable DGs. Figure 5 displays the eigenvalue plot of the whole network with the integration of REHDGs. Furthermore, Table 3 tabulates the eigenvalue analysis results of the critical modes of the base case system, with solar PV, wind, and biomass DGs, all REHDG units combined. The damping ratio of the critical mode improved from 0.0377 to 0.8568. It can be observed from Figure 3 (base case) that the network was voltage stable prior to the integration of REHDGs. Meanwhile, Figure 4 shows that SSS-wise, the system was marginally stable. The optimization results also indicate that as the voltage angle limits increased up until a maximum allowable value, the capacity of the renewable generators installed increased, indicating increased renewable power injections into the network. That means more power flowed in the network, and the network was more robust and could consume (contain) the effect of power variations (any small disturbances) from the intermittent renewable generations. This agrees with the power system property that the changes in active power flow were dependent on changes in voltage angles, but not on voltage magnitude, while changes in reactive power were governed by voltage magnitude changes. It is also deduced that setting limits or constraints on the voltage angle helped in constraining and enhancing the small-signal stability of the distribution system.

Numerical Results II: IEEE-118 Bus System
The proposed MILP model was applied to IEEE-118 bus test system whose network data and single line diagram were described in [56]. The system had 91 load buses, 54 substation buses, and 186 branches. The time varying load demand (hourly load demand) of the test system as scaled with the actual South African load profile [41] was used in this study.

Results of the Optimal REHDGs' Allocation Problem
As mentioned earlier, the optimization of renewable energy DGs' integration into DNS presented in this case study was lacking reactive power support from these generations, and the optimization results for the IEEE-118 test system are discussed below.
Based on the stopping criterion considered, the computation time required to obtain optimal solutions on the IEEE-118 system was 179 s. Tables 4 and 5 present the optimal solutions for REHDGs and the capacitor banks for a ten year planning stage, respectively. As also observed in Section 6.1.1, the larger percentage of the investment was made in the first year. This shows that more REHDG investments in the first stage were more economical than in the subsequent stages, because the costs were reducing gradually throughout the planning stages. Table 4 indicates that more of the wind DG units were installed than solar PV units despite being given equal parameters of integration. This is because wind generators had a higher capacity utilization factor than photovoltaic generators. Biomass generation was allowed to contribute 47.6% of renewable penetrations to make up the inadequate supply from the intermittent DGs. The biomass power generators were distributed throughout the length of the network to take economic advantage of the widespread sugar-cane plantations by reducing the cost accruable due to the transportation of biomass fuel. The total capacities of REHDG power (MW) installed throughout the planning stages are shown in Figure 6. The totals of 1117. 7  Similar to the observation in Section 6.1.1, it can be deduced from Table 5 that the optimal locations of capacitor banks were mainly on the most loaded buses and the buses close to the end of the network; thus helping to enhance network voltage stability by keeping the voltage within the limits. Figure 7 shows the total capacity of capacitor banks invested in and installed throughout the planning horizon to be 442. 7 Type   Stages  T1  T2  T3  T4  T5  T6  T7  T8  T9 T10  x g,i,t  Table 4. Cont. Stages  T1 T2 T3 T4 T5 T6 T7 T8 T9 T10 x g,i,t  In addition, the results from Tables 4 and 5 as captured in Figure 6 demonstrate the complementarity of intermittent renewable generations and biomass generation. Furthermore, the inclusion of reactive compensators greatly increased the capacity of renewable DG units that were integrated into the system by helping to maintain active and reactive power balance when reactive power consuming generators were installed. Without capacitor banks, the optimal capacity of REHDG units that could have been integrated would have been around 600 MW.

Type
For this case study, the peak load demand was 4.1963 GWh with 5% yearly demand growth projections throughout the planning horizon. The totals of 350.4 MW PV farms, 572.7 GW wind farms, and 855 MW biomass systems totaling 1.7788 GW were added to reduce the electricity generation from coal fired plants by 30%, amounting to the total NPV investment costs of $151.39 B, $7.57 B, $7.95 B, $8.35 B, $8.76 B, $9.20 B, $9.66 B, $10.14 B, $10.65 B, and $11.18 B for the ten planning stages, respectively. This brought the total investment costs to $234.85 B. The overall total NPV costs of maintenance, energy, and emission throughout the planning horizon are correspondingly listed in Table 6. The overall total NPV cost for the whole planning stages was $550.15 B. With the integration of 1.7788 GW of REHDG power, a total sum of $8.05 B would be saved if the whole power to satisfy the load demand had come only from coal fired conventional generations. Table 6. IEEE-118 bus system: Optimized cost information for the ten year planning horizon (×10 9 $). Overall  T1  T2  T3  T4  T5  T6  T7  T8  T9

Assessment of Small-Signal Stability
The assessment of the impact of REHDGs' integration on the long term dynamic voltage and small-signal stabilities of distribution network is presented in this section. Figures 8-10 show samples of the voltage profiles and eigenvalue plots of the network without renewable DGs (base case), as well as with solar PV DGs, wind DGs, biomass DGs, and all the REHDG units during every operational period throughout the planning horizon. In the same vein, considerable improvement in voltage profile is noticed in Figure 8. The two plots (base case and REHDG) in Figure 8 indicate dynamic voltage profiles at any selected bus (Bus 1 in this case) during the operational periods of DNS without and with renewable generations, respectively. From Figure 8 (base case), it can be seen that the voltage levels were very close to the minimum permissible limits. Meanwhile, the voltages in Figure 8 (REHDG) were largely close to the 1.0 pu value with an average deviation of 0.13%. This allowed significant operational margins for the permissible voltage magnitude limits. Invariably, the substantial improvements in voltage profiles were due to the effect of the combined integration of renewable DGs and capacitor banks.   Figure 9 shows typical eigenvalue plot results for the base case system. It can be seen that the system was critically damped with a small margin to oscillatory instability. The system's three state variables were: rotor angle (α), rotor speed (w s ), and transient voltage (E ). The eigenvalues of the base case system lied away from the origin to the left half of the S-plane except the quadrature transient internal voltage eigenvalues of six generators. The oscillatory frequency and damping ratio were calculated to be 3.2738 Hz and 0.001, respectively. From Figure 10, it is shown that all the eigenvalues of the component generators were situated away to the left half of S-plane during the integration of all these renewable DGs. Figure 10 shows the eigenvalue plot of the whole network with the integration of REHDGs. Furthermore, Table 7 tabulates the critical modes of the base case system, with solar PV, wind, and biomass DGs and all the REHDGs units combined. The damping ratio of the critical mode improved from 0.001 to 0.9894 for the overall network during the integration of REHDGs. It can be observed from Figure 8 (base case) that the network was voltage stable prior to the integration of REHDGs. Meanwhile, Figure 9 shows that SSS-wise, the system was marginally stable. The optimization results also indicated that as the voltage angle limits increased up until a maximum allowable value, the capacity of renewable generators installed increased, indicating increased renewable power injections into the network. That means more power flowed in the network, and the network was more robust and could consume (contain) the effect of power variations (any small disturbances) from the intermittent renewable generations. This agrees with the power system property that the changes in active power flow were dependent on changes in voltage angles, but not on voltage magnitude, while changes in reactive power were governed by voltage magnitude changes. It is also deduced that setting limits or constraints on the voltage angle helped in constraining and enhancing the small-signal stability of the distribution system.

Conclusions
This study presented a new multi-stage distribution expansion planning mathematical model to integrate and allocate large scale hybrid renewable DGs such as solar PV, wind, and biomass (sugar cane) and capacitor banks into distribution systems optimally. The scenario based probabilistic modeling approaches, beta and Weibull distributions, were applied to model the random behavior of solar irradiance and wind speed, respectively. Biomass DG was taken as a firm generation whose capacity could be determined as and when needed. The proposed planning model decided the optimal time of integration and the numbers, sizes and location of REHDGS and reactive compensators in the distribution networks simultaneously. The main objective of this optimization work, to maximize the REHDG power generated and absorbed into the distribution networks while the long term dynamic voltage and small-signal stabilities are maintained at the required levels and a least possible NPV of the total cost, was achieved. The model was formulated as a stochastic mixed integer linear programming (MILP) problem, while the non-linear AC network was made linear with the principle of fast decoupled power flow in order to characterize the network without the loss of generality, maintain accuracy, and reduce the computational complexity. Two standard test distribution systems, the IEEE-14 bus and IEEE-118 bus, were used to validate the proposed model successfully and conduct the required assessment based on the objectives of this work. The results of both case studies indicated that integrating biomass DGs and reactive compensators with variable renewable generations (PV and wind DGs) significantly increased the amount of renewable power injected into the networks. This eventually brought about a monumental decrease in the total cost compared to meeting the load demands with the conventional generations. For the IEEE-14 bus system, 411 MW of renewable power and 9.7 MVAr of reactive power were added to the system, while the IEEE-118 bus added 1.7788 GW REHDG power and 442.7 MVAr compensators to the network. In both case studies, the dynamic voltage and small-signal stabilities were highly enhanced.
The planning model for REHDGs and capacitor banks' integration proposed in this study demonstrated a significant improvement to the system in terms of dynamic stability enhancement, electricity and emission cost reductions, welfare, and environmental enhancement, and many other benefits accrued from it.
The formulation model proposed in this work was, therefore, a major step towards developing reliable and stable networks/grids that support the integration of large scale renewable generations.
Research in progress is focused on the application of intelligent search approaches for solving REHDGs' allocation problem, enhancing long term dynamic stability, and estimating the economic "end effect" of the hybrid distributed generation components. Future work will be devoted to the combined investigation of harmonic losses together with the long term dynamic stability of the distribution network during the integration of large scale hybrid renewable DGs. Another future research interest is the investigation of other distributed generation technologies in relation to system power quality and dynamic stability.