Stochastic Modelling and Stability Analysis of Large-Scale Wind Power Generation System with Dynamic Loads

Stochastic Modelling and Stability Analysis of Large-Scale Wind Power Generation System with Dynamic Loads

Joseph C. Attachie Christian K. Amuzuvi Godwin Diamenu*

Electrical and Electronic Engineering Department, University of Mines and Technology, (UMaT), P.O. Box 237, Tarkwa, Ghana

Renewable Energy Engineering Department, UMaT, P.O. Box 237, Tarkwa, Ghana

Corresponding Author Email:
8 June 2021
5 September 2021
18 September 2021
Available online: 
31 October 2021
| Citation

© 2021 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (



As a result of dwindling fossil fuel reserves and the negative impact of greenhouse gases (GHGs) on the environment, it is important that the search for a power grid that will comprise of only variable renewable energy (VRE) generation sources such as wind or solar energy which is available everywhere for free is given much attention. The main challenge associated with these sources of energy is their variability and random nature. It is as a result of instability introduced by the VRE generation sources, that is why there is a strict control measures put in place by the regulators as to how much VRE generation sources can feed into the power grid in the case of its integration into an existing power grid. It becomes imperative to consider implications for grid stability and reliability when considering a power grid that will consist of VRE generation only. Eigenvalue approach was used to analysed the performance of anticipated large-scale VRE grid to ascertain its behaviour. Eigenvalue approach is one of the methods to examine the stability of a power system’s dynamic performance. The results indicated that it is possible to attain the necessary stability provided the assumptions made during the modelling stage is revised to improve upon the model.


chapman-kolmogorov equation, eigenvalue, master equation, stability, stochastic, transition rate, variable renewable energy

1. Introduction

The purpose of a power system is to generate power, transmit this power and to distribute it to users at voltage levels and reliability that are appropriate to various users. The energy resources that are used to generate electric power can be categorised as fossil fuel, nuclear fuel and renewable resources. Fossil fuels include crude oil, coal, and natural gas. Fossil fuels are formed from dead plants and animals buried in the earth’s crust for millions of years under pressure and heat. The formation of fossil fuels takes millions of years and are considered non-renewable. Burning fossil fuels creates a wide range of pollution that includes the release of carbon dioxide (CO2), sulphur oxides, and the formation of nitrogen oxides. These are harmful gases that causes health and environmental hazards. Renewable energy resources (RERs) include hydropower, wind, solar, hydrogen, biomass, tidal, and geothermal. Renewable energy technologies (RETs) can produce sustainable, clean energy from renewable resources [1]. Carbon intensity is a measure of the amount of CO2 that is released into the atmosphere for every unit of energy produced and is dependent on the type of fossil fuel involved. Nuclear power plants do not emit CO2 however, carbon intensity of coal-fired power plants is higher. RERs such as wind, solar and hydropower plants emit no CO2.

Addressing climate change and the excessive greenhouse gas (GHG) emissions are among the top urging issues at a global level. Approximately 90% of our energy consumption comes from fossil fuels which is causing serious damage to our environment in various forms [2]. The main naturally occurring GHGs in the Earth’s atmosphere are water vapour and CO2. Of these two, it is water vapour which has the greatest GHG action. While CO2 concentrations are strongly influenced by human activity, atmospheric water vapour is almost entirely determined by climatic conditions and not human action [3]. Human activity is responsible for production of a number of other potent greenhouse gases, including methane, nitrous oxide, chlorofluorocarbons (CFCs) and hydrochlorofluorocarbons (HCFCs).

Climate change is caused by GHGs, and it is now widely accepted as a real condition that has potential consequences for human and other living organisms’ survival [4, 5]. Concerns about global warming and sustainability have recently spurned interest in renewable energy generation and usage. Generally, there are three main ways that CO2 emission can be reduced according to [6]:

(i) By moving from the traditional fossil fuel based such as coal, natural gas, and crude oil generation to renewable energy generation such as wind, solar, tidal, etc.;

(ii) By moving towards increased nuclear generation which is largely CO2-free; and

(iii) By removing CO2 from exhaust gases of traditional thermal generation using for example, carbon capture and storage technology.

For the past 200 years, the energy system is based solely on fossil fuels which has greatly promoted the development of human society. However, the awareness of serious consequences usually brought about by the large-scale use of fossil fuels is increasing, and at the same time, their fast rate of depletion, deterioration of the environment, as well as their roles in political and economic disputes of a number of nations and regions resulting in conflicts and wars has assumed an alarming proportion [2, 7].

According to the Ref. [8], there has been severe changes in climatic conditions in the past five years, including dramatic changes in temperature, sea level rise, and extreme weather events. GHG concentrations in the atmosphere have also increased to record levels, which is likely to worsen the warming trend for generations to come. Rising levels of GHGs concentration in the atmosphere is increasing the global average temperature, causing ocean acidification, rapid decrease in glaciation, and perturbing the climate system through, for instance, changed rainfall patterns. Storms are also more frequently occurring outside the traditional storm season. Rainfall patterns are not only getting more intense with climate change but are also swinging from wet to dry more erratically. Climate change is leading to observed increases in economic losses and human suffering, by displacing vulnerable populations and deepening existing inequalities.

About 90% of natural disasters experienced between the years 2015 - 2019 are weather related, with developing countries being the most affected due to limited adaptive capacities. The weather events include storms, flooding, heat waves and drought that continue to result in human and economic losses. Global average temperature has increased by 1.1℃ since the pre-industrial period, and by 0.2℃ compared to the year 2011 - 2015. The increase in temperature results in a rise in evaporation loses from reservoirs that support urban, rural, industrial and agricultural water needs. Such climatic conditions disrupt the agricultural activities, hydroelectric production and human settlements, and also affect critical communications and social infrastructure, reversing socio-economic development gains. Crude oil production is on the decline; resulting in exploration of crude oil from more difficult geological and environmental conditions. Not only has the cost of oil drilling increased, but also the energy consumed to produce the crude oil has also increased as well [9, 10].

For those countries where growing dependence on imported fossil fuel is a significant energy security issue, renewable energy resources (RERs) can provide alternative. Hence, most countries are introducing regulations and policies to minimise the use of fossil fuels in order to reduce the GHGs effects and to preserve the little remaining for the future generation. For that reason, the use of green and renewable (RERs) like wind, solar, geothermal, tidal waves, etc. that will have to provide for most energy in the future are being promoted by those countries [4].

Model-based and measurement-based analyses are the approaches for determining a power system’s behaviour. Measurement-based analysis is a real-time monitoring approach for power systems. It is different from model-based analysis. It uses phasor measurement units (PMUs) which are high-speed sensors that measure voltage and current synchro-phasors of the power system with the accuracy in the order of one microsecond to obtain the necessary information from the power system. Measurement-based analysis provides accurate data and the prediction of the stability behaviour of a power system based on synchro-phasor technology. When model-based analysis approach is used to determine the stability of a power system, the power system is described using mathematical models. Among the analysis approaches that can be applied to the mathematical model include Lyapunov method, eigenvalue analysis and Routh-Hurwitz stability criterion. However, among all these analyses approaches, eigenvalue analysis approach seems to be superior due to the fact it has a number of indices that can be used to determine the stability of the power system or otherwise. These properties or indices include the polarity, controllability, observability, participation factor, mode shape and sensitivity [11].

2. Materials and Method

In the literature, there are a lot on models related to RERs integrated into an existing power grid. However, much was not found on models pertaining to entirely RERs left alone on a large-scale. This section describes the formulation of one such model based on large-scale wind turbine (WT) power generation system and examined the stability of such a model. An entirely variable renewable energy (VRE) grid must ensure an increased power quality and system stability as being envisaged for the future.

2.1 Model formulation

Wind power generation system modelling based on Chapman-Kolmogorov (C-K) equation.

In stochastic model, the knowledge of the state variables is subject to uncertainties due to some inherent randomness. In other words, the variable states are not described by unique values, but rather by probability distributions. This type of model is suitable for the applications where the state information cannot be perfectly predictable. Markov processes are often used to model randomness. Markov process has no memory. A consequence of the Markov property is the Chapman-Kolmogorov equation. The differential form of the Chapman-Kolmogorov Eq. is the Master Eq. (ME) for jump processes or diffusive processes.

The ME describes the evolution of the probability distribution of states of a continuous-time Markov process and is commonly used to describe stochastic physical, chemical, or biological systems. Jump processes are characterised by discontinuous motion; thus, there is a bounded and non-vanishing transition probability per unit time. The ME is a stochastic model and is derived from the fundamental assumption on the dynamic properties of the underlying stochastic process referred to as the Markov property [12].

The state variables are subject to random phenomena. A large-scale power grid consisting of a number of off-shore and on-shore wind farms supplying power to dynamic loads is as given in Figure 1. As a result of wind speed intermittency or variations, the power generation into the power grid is stochastic or random in character impacting the stability and reliability of the power grid. The fluctuation of the renewable power generation can be regarded as random variables which are real-valued functions defined on a sample or state space with the assignment of possible probabilities to the possible values of the random variables [10, 13]. Due to the large-scale nature of the renewable power grid being considered, continuous random variable is used which can be easily handled analytically as compared to discrete random variables. Again, for large-scale modelling, continuous random variables are preferred [11, 14].

A typical WT power curve is given in Figure 2 depicting the three types of wind speeds usually exhibited by wind turbines operating on wind farms [12, 15]. Based on the wind turbine power curve, three transition states are assumed for the wind turbines for large-scale wind farm power production.

Figure 1. Power grid with stochastic generation from wind turbines

Figure 2. Typical wind turbine power curve

In this study, the intermittency nature and the changing pattern of the wind speed is transformed into a transition diagram as given in Figure 3 having three states. One can represent the transitions among the various states as random events that do occur at some rates. Each transition state is assigned a probability, p and in between the states is the transition rates represented by μ and β. Transitions between the different states occur from time to time.

Figure 3. State transition diagram for wind turbine systems

From Figure 3, the six possible transition states are (Table 1):

Table 1. Wind speed variation transition states


Current state

Next state

Transition rate



x $\rightarrow$ x+1




x $\rightarrow$ x-1




x-1 $\rightarrow$ x




x-1$\rightarrow$ x+1




x +1$\rightarrow$ x




x +1$\rightarrow$ x-1



x = rated wind speed (required wind speed);

x + 1 = cut-out wind speed (increase in wind speed);

x - 1 = cut-in wind speed (decrease in wind speed).

A stochastic process is a process that evolves probabilistically through various states attained at various times from a well-defined initial state x0 at t0. The probability density function (pdf) of the process depends on their respective states and times at particular moment. The probability that the process will be at state x1 at time, t1; then progress to state x2, at time, t2 then goes through all other subsequent configurations or states to be at state xn at time tn, given that it started from state x0 at time t0 is given by Eq. (1):

$\begin{align}  & {{P}_{\,n}}{{|}_{\,1}}=(({{x}_{\,n}},\,\,{{t}_{n}}),\,({{x}_{\,n-1}},\,\,{{t}_{\,n-1}}),(x{{\,}_{n-2}},\,\,t{{\,}_{n-2}}) \\  & .........({{x}_{1}},\,{{t}_{1}})|({{x}_{0}},\,\,{{t}_{0}})) \\ \end{align}$         (1)

The probability of the system being at state, x1 at time t1, multiplied by the probability of it being at state, x2 at time, t2 given that it has been at state, x1 at t1, and at state x0, at time, t0 respectively can be expressed as [16, 17]:

${{P}_{2\,|\,1}}(({{x}_{2}},\,{{t}_{2}}),\,({{x}_{1}},\,{{t}_{1}})|({{x}_{0}},\,{{t}_{0}}))=P(({{x}_{1}},\,\,{{t}_{1}})|({{x}_{0}},\,\,{{t}_{0}})){{P}_{1|2}}({{x}_{2}},\,\,{{t}_{2}})|({{x}_{1}},\,\,{{t}_{1}}),\,({{x}_{0}},\,\,{{t}_{0}}))$               (2)

If all such probabilities are known, the stochastic process is fully specified [16, 17].

The conditional probability density given by [18]:

$f\left(x_{2}, t_{2} / x_{1}, t_{1}\right)=\frac{f\left(x_{1}, x_{2} ; t_{1}, t_{2}\right)}{f\left(x_{1}, t_{1}\right)}$

$f\left(x_{k+1}, \ldots \ldots, x_{n} ; t_{k+1}, \ldots, t_{n} / x_{1}, \ldots \ldots x_{k} ; t_{1}, \ldots \ldots, t_{k}\right)=\frac{f\left(x_{1}, \ldots \ldots \ldots, x_{n} ; t_{1}, \ldots \ldots, t_{n}\right)}{f\left(x_{1}, \ldots \ldots \ldots x_{k} ; t_{1}, \ldots \ldots ., t_{k}\right)}$           (3)

where the relationship between the conditional and joint probabilities is given by:

$\begin{align}  & f({{x}_{1}},.......,{{x}_{n}};\,{{t}_{1}},.....,{{t}_{n}})\,\,=\,\, \\  & f({{x}_{n}},\,\,{{t}_{n}}/{{x}_{n-1}},....,{{x}_{1}};\,{{t}_{n-1}},........{{t}_{1}})\,\,x\,......... \\  & .....x\,f({{x}_{2}},\,\,{{t}_{2}}/{{x}_{1}},\,\,{{t}_{1}})\,x\,f({{x}_{1}},\,\,{{t}_{1}}) \\ \end{align}$

The ME for jump processes is given by Eq. (3):

$\frac{\partial p(x,\,t)}{\partial t}\,\,=\,\,\int{(w(x|{x}')p({x}',\,\,t)\,\,-\,w({x}'|x)p(x,\,\,t))d{x}'}$                (4)

where, p(x, t)=probability of being in state x at time, t;

$w\left(x \mid x^{\prime}\right)$=transition probability of moving from state x' to state x per unit time.

One-step stochastic process is are continuous-time Markov process whose range consists of the integer n, and whose transition probability per unit time, wnm permits only jumps between adjacent sites. Thus,

$\begin{align}  & {{\omega }_{nm}}\,=\,{{r}_{m}}{{\delta }_{n,m-1}}\,+\,{{g}_{m}}{{\delta }_{n,m+1}},\,\,(m\ne n) \\  & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\omega }_{nm}}\,=\,1\,-\,({{r}_{n}}\,+\,{{g}_{n}}) \\ \end{align}$                (5)

where, rn = probability per unit time that, being in state n, a jump occurs to site n - 1.

gn = probability per unit time that, being in state n, a jump occurs to state n + 1 [18].

A wind turbine (WT) transitioning from one state to another can be regarded as a jump process. For this modelling, every transition state is taken to be positive since the WT is not supposed to be stationary hence, a current state depends on all other previous states which can be treated as a union of disjoint intervals based on the theory and concept of continuous random variable. For a continuous random variable X with probability density function (PDF) fX(x), the following holds [13]:

$(\text{i})\,\,{{f}_{X}}(x)\,\,\ge \,\,0\,\,\text{for}\,\,\text{all}\,\,x\in \,\,\mathbb{R}$        (6)

$(\text{ii})\,\,\int_{-\infty }^{\infty }{{{f}_{X}}(u)\,}du\,\,=\,\,1$                    (7)

$(\text{iii})\,\,\,P(a\,\,<\,\,X\,\,\le \,\,\,b)\,\,=\,{{F}_{X}}(b)\,\,-\,\,{{F}_{X}}(a)\,\,=\,\,\int\limits_{a}^{b}{{{f}_{X}}}(u)\,du$                  (8)

$\begin{align}  & \text{(iv)}\,\,\text{More}\,\text{generally,}\, \\  & \text{for}\,\,\text{a}\,\text{set}\,A,\,\,P(X\in \,A)\,\,=\,\int\limits_{A}{{{f}_{X}}}(u)\,du\, \\ \end{align}$                (9)

To determine the probability of a continuous random variable X, for instance, P([1, 2] U [5, 6]); one can write [13]:

$P(X\,\in \,\,[1,\,2]\,\,\cup \,\,[5,\,6]\,\,=\,\,\int\limits_{1}^{2}{{{f}_{X}}(u)\,\,du\,\,+\,\,}\int\limits_{5}^{6}{{{f}_{X}}(u)\,\,du\,\,\,}$               (10)

Also, the probability density corresponding to the cumulative distribution F(x; t) is given by [18]:

$f(x;t)\,=\,\frac{\partial F(x;t)}{\partial x}$.                   (11)

For two instants t1 and t2 with their corresponding random variables ξ(t1) and ξ(t2), their joint cumulative distribution is given by:

$F({{x}_{1}},\,{{x}_{2}};\,{{t}_{1}},\,{{t}_{2}})\,=\,P\{\xi ({{t}_{1}})\le {{x}_{1}},\,\xi ({{t}_{2}})\le {{x}_{2}}\}$                 (12)

With the joint probability density defined as:

$f({{x}_{1}},\,{{x}_{2}};\,{{t}_{1}},\,{{t}_{2}})\,\frac{{{\partial }^{2}}F({{x}_{1}},\,{{x}_{2}};\,{{t}_{1}},\,{{t}_{2}})}{\partial {{x}_{1}}\partial {{x}_{2}}}$         (13)

In general,

$\begin{align}  & F({{x}_{1}},......,{{x}_{n}};\,{{t}_{1}},......,{{t}_{n}})\,\, \\  & =\,\,P\{\xi ({{t}_{1}})\le {{x}_{1}},.....,\xi ({{t}_{n}})\le {{x}_{n}}\} \\ \end{align}$             (14)

$f({{x}_{1}},........,{{x}_{n}};\,{{t}_{1}},......,{{t}_{n}}\,=\,\frac{{{\partial }^{n}}F(\mathop{x}_{_{1}},......,{{x}_{n}};\,{{t}_{_{1}}},.......,{{t}_{n}})}{\partial \,{{x}_{_{1}}}............\partial \,{{x}_{n}}}$               (15)

The nth-order distribution function completely determines the stochastic process.

Eq. (14) must satisfy the symmetry condition given for every permutation:

$({{j}_{1}},\,\,{{j}_{2}},........{{j}_{n}})\,\,\text{of}\,\,(1,\,\,2,\,.......,n)\text{ then}$

$\begin{align}  & F({{x}_{j1}},\,{{x}_{j2}},.......,{{x}_{jn}};{{t}_{j1}},\,{{t}_{j2}},......{{t}_{jn}})\,\,=  F({{x}_{1}},\,{{x}_{2}},.......,{{x}_{n}};\,\,{{t}_{1}},\,{{t}_{2}},.........{{t}_{n}}) \\\end{align}$                       (16)

This condition follows from the AND logic of the cumulative distribution function. For instance, the statements ‘A AND B AND C’ and ‘A AND C AND B’ are equivalent.

Equation (A) as well must satisfy the compatibility condition given which among other things states that for m < n; then:

$\begin{align}  & F({{x}_{1}},......{{x}_{m}},\,\infty ,......\infty ;\,{{t}_{1}},.....{{t}_{m}},{{t}_{m+1}},.....{{t}_{n}})\,\,=\,\,  F({{x}_{1}},.....,{{x}_{m}};{{t}_{1}},.......,{{t}_{m}}) \\ \end{align}$                  (17)

This statement again follows from the AND logic: the statements ‘A AND B AND C’ and ‘A AND TRUE AND C’ are equivalent [18].

Based on Eq. (3), the ME for jump processes is modified and applied to wind turbine (WT) power generation system which is given by Eq. (9):

$\frac{\partial p(x,\,t)}{\partial t}\,\,=\,\,\int{(w(x|{x}')p({x}',\,\,t)\,\,+\,\,w({x}'|x)p(x,\,\,t))d{x}'}$              (18)

Assuming the system is at state x (state 1 or initial state) at time, t moving to state x – 1 (state 2 or next state) at time, t + Δt. The probability of moving from state x (state 1 or initial state) to state x – 1 (state 2 or next state) is given by ΩΔt. The transition must occur during the time, Δt. For the transition to occur:

(i) The WT power generation system must be at state x (state 1 or initial state) with the probability, px(t); and

(ii) Transition must occur from state x (state 1 or initial state) to state x - 1 (state 2 or next state) with the probability, ΩΔt, where, Ω = transition rate.

The probability of transitioning from state x (state 1 or initial state) to state x – 1 (state 2 or next state) between the times t, t + Δt is the product of the two probabilities given by px(t)ΩΔt.

During that small time interval, Δt, only one change in state must occur, either x – 1 (decrement in wind speed) or x + 1 (increment in wind speed). Transitioning to the next state depends on the previous state and the transition rate. Based on the above, dynamical equations can be developed for state x (state 1 or initial state), state x – 1 (state 2 or next state) or x + 1 state (another possible state 2 or next state) as follows based on the rules of conditional probability [13, 17].

${{P}_{(x)}}(t\,+\,\Delta t)\,\,=\,\,{{P}_{(x)}}(t)\,\,+\,{{P}_{(x)}}(t)\Omega \Delta t\,$    (19)

That of state (x-1) transition is given by:

${{P}_{(x\,-1)}}(t\,+\,\Delta t)\,\,=\,\,{{P}_{(x-1)}}(t)\,\,+\,{{P}_{(x-1)}}(t)\Omega \Delta t\,$                (20)

And that of state (x + 1) will also be given by:

${{P}_{(x\,+1)}}(t\,+\,\Delta t)\,\,=\,\,{{P}_{(x+1)}}(t)\,\,+\,{{P}_{(x+1)}}(t)\Omega \Delta t\,$              (21)

Dividing through Eq. (19), Eq. (20) and Eq. (21) respectively by Δt and taking the limit as Δt → 0, Eq. (10), Eq. (11) and Eq. (12) becomes:

$\frac{d{{p}_{x}}}{dt}\,\,\,=\,\,\Omega \,{{p}_{x}}$       (22)

$\,\,\,\,\,\,\,\,\,\,\,\,\frac{d{{p}_{(x-1)}}}{d\,t}\,\,\,=\,\,\Omega \,\,\,{{p}_{\,(x-1)}}$              (23)

 $\frac{d{{p}_{(x+1)}}}{dt}\,\,\,=\,\,\Omega \,{{p}_{(x+1)}}$                     (24)

Generalizing and applying the principle to the entire three states with six transition diagrams of Figure 3 (with the transitions given in Table 1) which represents the movement pattern of WTs due to fluctuations in wind speed with 1, 2, 3…n wind farms comprising a number of WTs; the dynamical ME in matrix notation are:

$\frac{d}{d\,t}\left( \begin{smallmatrix}  {{p}_{1}}_{_{(x-1)}} \\  {{p}_{1}}_{_{x}}  \\  {{p}_{1}}_{_{(x+1)}} \end{smallmatrix} \right)\,\,\,\,\,=\,\,\,\,\left( \begin{smallmatrix}  {{\beta }_{1}}\,{{p}_{1}}{{\,}_{(x-1)}}\,\,\,+{{\beta }_{2}}\,{{p}_{1}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,{{\mu }_{2}}{{p}_{1}}_{x}\,\,\,\,\,\,\,\,\,{{\mu }_{1}}{{p}_{1}}_{\,(x+1)} \\  {{\beta }_{2}}\,{{p}_{1}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,\,({{\mu }_{2}}\,+\,\,{{\beta }_{3}})\,{{p}_{1}}{{\,}_{(x+1)}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\mu }_{3}}{{p}_{1}}_{x}  \\  {{\beta }_{1}}\,{{p}_{1}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,{{\beta }_{3}}\,{{p}_{1}}{{\,}_{(x)}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({{\mu }_{1}}\,+\,{{\mu }_{3}})\,\,{{p}_{1}}_{\,(x+1)} \end{smallmatrix} \right)$

$\frac{d}{d\,t}\left( \begin{smallmatrix}  {{p}_{2}}_{_{(x-1)}} \\  {{p}_{2}}_{_{x}}  \\  {{p}_{2}}_{_{(x+1)}} \end{smallmatrix} \right)\,\,\,\,\,=\,\,\,\,\,\left( \begin{smallmatrix}  {{\beta }_{\,1}}{{p}_{2}}{{\,}_{(x-1)}}\,\,\,+{{\beta }_{2}}{{p}_{2}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,{{\mu }_{2}}{{p}_{1}}_{x}\,\,\,\,\,\,\,\,\,{{\mu }_{1}}{{p}_{2}}{{\,}_{(x+1)}} \\  {{\beta }_{2}}{{p}_{2}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,\,({{\mu }_{2}}\,+\,\,{{\beta }_{3}})\,\,{{p}_{2}}_{(x+1)}\,\,\,\,\,\,\,\,\,\,{{\mu }_{3}}{{p}_{2}}_{x}  \\  {{\beta }_{\,1}}{{p}_{\,2}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,{{\beta }_{3}}{{p}_{1}}{{\,}_{(x)}}\,\,\,\,\,({{\mu }_{1}}\,+\,{{\mu }_{3}}){{p}_{2}}{{\,}_{(x+1)}} \end{smallmatrix} \right)$

$\begin{align}  & \vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots  \\  & \vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots  \\  & \vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\,\vdots \,\,\,\vdots \,\,\vdots \,\,\vdots \,\,\vdots  \\ \end{align}$

$\frac{d}{d\,t}\left( \begin{smallmatrix}  {{p}_{n}}_{{{\,}_{(x-1)}}} \\  {{p}_{n}}_{_{\,x}}  \\  {{p}_{n}}{{\,}_{_{(x+1)}}} \end{smallmatrix} \right)\,\,\,\,=\,\,\,\left( \begin{smallmatrix}  \beta {{\,}_{\,1}}\,{{p}_{n}}_{\,\,(x-1)}\,\,\,+{{\beta }_{\,2}}\,{{p}_{n}}_{\,\,\,(x-1)}\,\,\,\,\,\,\,\,\mu {{\,}_{2}}\,{{p}_{n}}_{x}\,\,\,\,\,\,\,\,\,\mu {{\,}_{1}}\,{{p}_{n}}\,{{\,}_{(x+1)}} \\  \beta {{\,}_{2}}\,{{p}_{n}}{{\,}_{(x-1)}}\,\,\,\,\,\,\,\,\,({{\mu }_{2}}\,+\,\,\beta {{\,}_{3}})\,{{p}_{n}}_{\,\,(x+1)}\,\,\,\,\,\,\,\,\,\,{{\mu }_{\,3}}\,\,{{p}_{n}}_{x}  \\  \beta {{\,}_{1}}\,{{p}_{n}}_{\,\,(x-1)}\,\,\,\,\,\,\,\,{{\beta }_{\,3}}\,{{p}_{n}}_{\,\,(x)}\,\,\,\,\,({{\mu }_{\,1}}\,+\,\mu {{\,}_{3}})\,{{p}_{n}}\,{{\,}_{(x+1)}} \end{smallmatrix} \right)$             (25)

The matrix-vector form of Eq. (25) is given by Eq. (26):

$\frac{d}{d t}\left(\begin{array}{l}p \\ 1(x-1) \\ p_{1 x} \\ p_{1}(x+1)\end{array}\right)=\left(\begin{array}{ccc}\left(\beta_{1}+\beta_{2}\right) & \mu_{2} & \mu_{1} \\ \beta_{2} & \left(\mu_{2}+\beta_{3}\right) & \mu_{3} \\ \beta_{1} & \beta_{3} & \left(\mu_{1}+\mu_{3}\right)\end{array}\right)\left(\begin{array}{l}p_{1}(x-1) \\ p_{1} x \\ p_{1}(x+1)\end{array}\right)$

$\frac{d}{d t}\left(\begin{array}{l}p_{2(x-1)} \\ p_{2 x} \\ p_{2}(x+1)\end{array}\right)=\left(\begin{array}{lcc}\left(\beta_{1}+\beta_{2}\right) & \mu_{2} & \mu_{1} \\ \beta_{2} & \left(\mu_{2}+\beta_{3}\right) & \mu_{3} \\ \beta_{1} & \beta_{3} & \left(\mu_{1}+\mu_{3}\right)\end{array}\right)\left(\begin{array}{l}p_{2}(x-1) \\ p_{2} x \\ p_{2}(x+1)\end{array}\right)$

$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$

$\,\frac{d}{d\,\,t}\left( \begin{smallmatrix}  {{p}_{_{n}\,{{\,}_{\,(x\,-\,1)}}}} \\  {{p}_{n}}_{_{x}}  \\  {{p}_{n}}_{_{\,\,(x\,+\,1)}} \end{smallmatrix} \right)\,\,\,=\,\,\,\left( \begin{smallmatrix} ({{\beta }_{\,\,\,1}}\,\,+\,{{\beta }_{\,\,\,2}}\,)\,\,\,\,\,\,\,\mu {{\,}_{\,\,2}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\mu }_{\,\,\,1}}\, \\  _{\,{{\beta }_{\,\,\,2}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({{\mu }_{\,\,2}}\,\,+\,{{\beta }_{\,\,3}})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\mu }_{\,\,3}}}  \\  \,{{\beta }_{\,\,\,1}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\beta }_{\,\,3}}\,\,\,\,\,\,\,\,\,\,\,\,\,({{\mu }_{\,\,\,1}}\,\,+\,{{\mu }_{\,\,\,3}}\,) \end{smallmatrix} \right)\left( \begin{smallmatrix} {{p}_{_{n}}}_{_{_{(x\,-\,1)}}} \\  {{p}_{_{n}}}_{x}  \\  {{p}_{_{n}}}_{_{_{(x+1)}}} \end{smallmatrix} \right)$             (26)

The state matrix for the first wind farm is given by:

$\,U\,\,=\,\,\,\left( \begin{align}  & ({{\beta }_{\,1}}\,+\,{{\beta }_{\,2}})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\mu }_{\,2}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\mu }_{\,1}} \\  & {{\beta }_{\,2}}_{\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({{\mu }_{\,2}}\,\,+\,{{\beta }_{\,3}})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\mu {{\,}_{3}}\,\,} \\  & \beta {{\,}_{1}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\beta }_{\,3}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({{\mu }_{\,1}}\,\,+\,{{\mu }_{\,3}}) \\ \end{align} \right)$                 (27)

The transition rate for the wind farms, which varies between 0.1 to 1.0, the least being 0.1 and the maximum being 1.0, the maximum being 1.0 is assumed for the wind farm model. Eq. (18) now becomes:

$\,U\,\,=\,\,\left( \begin{align}  & 2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,1.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,1.0 \\  & {{1.0}_{\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,1.0}} \\  & 1.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,1.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0 \\ \end{align} \right)$                    (28)

Even though large-scale WT power generation system with its associated size of state matrix, [U] should be in order of thousand or more, depending on the number of wind farms being considered that goes with its many eigenvalues, it is not necessary to compute all these eigenvalues [19].

2.2 Model assumptions

The following assumptions were made during the development of the model:

  1. The wind speed variation is a function of the weather conditions only;

  2. A transition to the next state may occur at any point in time;

  3. Only one transition can occur at a time (1-step transition);

  4. The transition to another state is gradual over time (not immediate);

  5. A transition from one state to the next state for all WTs constituting a wind farm occurs at the same time;

  6. The transition rates vary between 0.1 – 1.0;

  7. The wind power plants are spread geographically to deal with the issue of short-term variability;

  8. The WTs operate in variable-speed mode, and the control system regulates the rotor speed to obtain peak efficiency in fluctuating winds;

  9. The wind farms are sited in areas with high average wind speeds; and

  10. The electrical load model is based on Markov property (current state does not depend on previous state).

2.3 The electrical load modelling

The variation of real power, (P) and reactive power, (Q) taken by a load with the variation of the voltage is important when representing the load for power flow and stability studies. For system studies, usually the load on a substation is considered, and this load is of a composite nature consisting of industrial, commercial, and domestic consumers [20]. The dynamic load being considered here can vary between 1% to 100%, however, for this particular modelling, the variations is assumed to be between 30% to 90%. The least electrical load variation being 30% is represented by state A, that of 60% electrical load variation point is represented by state B and that of the 90% is represented by state C. Consequently, the transition among the three electrical load variation points or states is as indicated in Figure 4. The arrow indicates the transition and the rate at which the transition takes place being represented by ω and δ. All transition towards the right-hand side, the transition rate is indicated by δ, and that towards the left-hand side is indicated by ω.

Figure 4. State transition diagram for the electrical loads

From Figure 4, the six possible transition states for the electrical loads are given in Table 2:

Table 2. Electrical load transition states


Current state

Next state

Transition rate

























The matrix-vector forms of the MEs for the electrical load is given by Eq. (20):

$\frac{d}{d t}\left(\begin{array}{l}p_{1(B)} \\ p_{1(A)} \\ p_{1(C)}\end{array}\right)=\left(\begin{array}{ccc}\left(\delta_{1}+\delta_{2}\right) & \omega_{2} & \omega_{1} \\ \delta_{2} & \left(\omega_{2}+\delta_{3}\right) & \omega_{3} \\ \delta_{1} & \delta_{3} & \left(\omega_{1}+\omega_{3}\right)\end{array}\right)\left(\begin{array}{l}P_{1(B)} \\ P_{1(A)} \\ P_{1(C)}\end{array}\right)$

$\frac{d}{d t}\left(\begin{array}{l}p_{2(B)} \\ p_{2(A)} \\ p_{2(C)}\end{array}\right)=\left(\begin{array}{cccc}\left(\delta_{1}+\delta_{2}\right) & \omega_{2} & \omega_{1} \\ \delta_{2} & \left(\omega_{2}+\delta_{3}\right) & \omega_{3} \\ \delta_{1} & \delta_{3} & \left(\omega_{1}+\omega_{3}\right)\end{array}\right)\left(\begin{array}{l}P_{2(B)} \\ P_{2(A)} \\ P_{2(C)}\end{array}\right)$

$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$$\begin{align}  & \vdots  \\  & \vdots  \\  & \vdots  \\ \end{align}$

$\frac{d}{d t}\left(\begin{array}{l}p_{n(B)} \\ p_{n(A)} \\ p_{n(C)}\end{array}\right)=\left(\begin{array}{ccc}\left(\delta_{1}+\delta_{2}\right) & \omega_{2} & \omega_{1} \\ \delta_{2} & \left(\omega_{2}+\delta_{3}\right) & \omega_{3} \\ \delta_{1} & \delta_{3} & \left(\omega_{1}+\omega_{3}\right)\end{array}\right)\left(\begin{array}{c}P_{n}(B) \\ P_{n}(A) \\ P_{n}(C)\end{array}\right)$                 (29)

The state matrix for the electrical load is therefore given by:

$V\,\,=\,\,\left( \begin{align}  & ({{\delta }_{1}}\,\,\,+\,{{\delta }_{2}})\,\,\,\,\,\,\,\,\,{{\omega }_{2}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\omega }_{1}}\,\,\,\,\,\, \\  & {{\delta }_{2}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(\,{{\omega }_{2}}\,+\,\,\,{{\delta }_{3}})\,\,\,\,\,\,\,\,\,\,\,\,{{\omega }_{2}}\,\,\, \\  & {{\delta }_{1}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\delta }_{3}}\,\,\,\,\,\,\,\,\,\,\,\,(\,{{\omega }_{1}}\,+\,\,{{\omega }_{2}}\,)\,\, \\ \end{align} \right)$                   (30)

Assuming a transition rate of 1.0; Eq. (21) becomes:

$V\,\,\,=\,\,\,\,\left( \begin{align}  & \,2.0\,\,\,\,\,\,\,\,\,\,\,1.0\,\,\,\,\,\,\,\,\,\,\,1.0\,\,\,\,\,\, \\  & \,1.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0\,\,\,\,\,\,\,\,\,\,\,1.0\,\,\, \\  & \,1.0\,\,\,\,\,\,\,\,\,\,1.0\,\,\,\,\,\,\,\,\,\,\,\,2.0\,\, \\ \end{align} \right)$              (31)

2.4 The eigenvalue stability approach

The Eigen function usually represents stationary states of a system. Such a system can achieve that state under certain conditions and eigenvalues represent the values of that property of the system in that stationary state. The mathematical sign (positive or negative sign) pattern of the eigenvalues indicates the behaviour of the equilibrium state of the dynamical system under investigation. The eigenvalue approach is used to investigate the stability of the wind turbine power system supplying power to dynamic electrical loads.

The state matrix of both systems is combined after which the stability of the entire system is determined. Eq. (19) and Eq. (18) are combined resulting in Eq. (22).

$M=\left(\begin{array}{llll}2.0 & 1.0 & 1.0 \\ 1.0 & 2.0 & 1.0 \\ 1.0 & 1.0 & 2.0\end{array}\right)+\left(\begin{array}{cccc}-2.0 & 1.0 & 1.0 \\ 1.0 & -2.0 & 1.0 \\ 1.0 & 1.0 & -2.0\end{array}\right)$

$\,M\,=\,\,\left( \begin{align}  & 0.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0 \\  & {{2.0}_{\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,0.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0}} \\  & 2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,0.0 \\ \end{align} \right)\,\,\,\,$                (32)

The matrix, M has non-trivial solution given by Eq. (24) which says that the values of λ which satisfies:

$\det \left| M\,-\,\lambda I \right|\,\,=\,\,0$              (33)

where, λ = unknown eigenvalue(s) of the matrix, M, I = n x n identity matrix.

In general, for any n x n matrix M,

$M\,\,\,-\,\,\,\lambda I\,\,=\,\,\left[ \begin{align}  & {{m}_{11}}\,-\,\,\lambda \,\,\,\,\,\,\,\,\,\,\,{{m}_{12}}\,\,\,\,\,\,\,\,\,\,....\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{m}_{1n}} \\  & {{m}_{21}}\,\,\,\,\,\,\,\,\,\,\,\,\,{{m}_{22}}\,-\,\lambda \,\,\,\,\,\,\,\,....\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{m}_{2n}} \\  & ...\,.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,.....\,\,\,\,\,\,\,\,\,\,\,\,\,\,....\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,.... \\  & {{m}_{n1}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{m}_{n2}}\,.................{{m}_{nn}}\,-\lambda  \\ \end{align} \right]$                  (34)

The determinant of Eq. (34) is expanded along the first row as follows [21].

$\operatorname{Det}(M-\lambda I)=\left(m_{11}-\lambda\right)$ 

              $\operatorname{det}\left[\begin{array}{llll}m & 22^{-\lambda} & \ldots \ldots & m & 2 n \\ \ldots \ldots & & \ldots & \\ m & n 2 & \cdots \cdot & & m \\ & & & n n\end{array}\right]+\ldots$              (35)

Based on Eq. (34), Eq. (32) is solved as follows for the eigenvalues or characteristic roots of the wind farm 1 (WF-1) made up of a number of wind turbines and the electrical loads as follows:

$\left(\left[\begin{array}{ccc}\lambda & 0 & 0 \\ 0 & \lambda & 0 \\ 0 & 0 & \lambda\end{array}\right]-\left[\begin{array}{ccc}0.0 & 2.0 & 2.0 \\ 2.0 & 0.0 & 2.0 \\ 2.0 & 2.0 & 0.0\end{array}\right]\right)=0$              (36)

$\left[ \begin{align}  & 0\,-\,\,\lambda \,\,\,\,\,\,\,\,\,\,\,\,2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0 \\  & {{2.0}_{\,\,\,\,\,\,\,\,\,\,\,\,}}\,\,\,0\,-\lambda \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0 \\  & 2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,2.0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,0\,-\lambda  \\ \end{align} \right]\,\,\,\,\,=\,\,\,\,\,0$                 (37)

Expanding the determinants of Eq. (28) along the first row:

$-\lambda\left|\begin{array}{lr}-\lambda & 2 \\ 2 & -\lambda\end{array}\right|-2\left|\begin{array}{cc}2 & 2 \\ 2 & -\lambda\end{array}\right|+2\left|\begin{array}{rr}2 & -\lambda \\ 2 & 2\end{array}\right|=0$                  (38)

$-\lambda \,\,\left[ {{(-\lambda )}^{2}}\,\,-\,\,4 \right]\,\,-\,2\,\,\left[ -2\,\lambda -4 \right]\,\,\,+\,\,2\,\,\left[ 4\,\,+\,\,2\lambda \,\, \right]\,\,=\,\,0$               (39)

${{\lambda }^{3}}\,\,-12{{\lambda }^{2}}\,\,-16\,\,=\,\,\,0$             (40)

Eq. (40) is the eigenvalue Equation or characteristic Equation of the WF-1 and its loads. Based on Eq. (31), the eigenvalue equations or characteristic polynomial equations of the entire wind power system made of 1, 2…... n sizes of wind farms and electrical loads is given by Eq. (41).

$\begin{align}  & \,W{{F}_{1}}:\mathop{\lambda }_{1}^{3}\,\,-12{{\lambda }_{1}}\,\,-\,{{16}_{1}}\,\,\,\,=\,\,0 \\  & W{{F}_{2}}:\,\,\mathop{\lambda }_{2}^{3}\,\,-12{{\lambda }_{2}}\,\,-\,{{16}_{2}}\,\,=\,\,0 \\  & W{{F}_{3}}:\,\mathop{\lambda }_{3}^{3}\,\,\,-12{{\lambda }_{3}}\,\,-\,\,{{16}_{3}}\,\,=\,\,0 \\  & \,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\, \\  & \,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\, \\  & \,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\, \\  & W{{F}_{n}}:\,\mathop{\lambda }_{n}^{3}\,\,\,\,-12{{\lambda }_{n}}\,\,-\,{{16}_{n}}\,\,=\,0\,\,\,\,\,\,\, \\ \end{align}$               (41)

Determining the eigenvalues or characteristic roots of Eq. (40), thus:

${{\lambda }^{3}}\,\,-12\lambda \,\,-\,\,16\,\,=\,\,0$

By the method of inspection, (λ + 2) is a factor of Eq. (40), and can be written as:

$\left( \lambda +2 \right)\left( {{\lambda }^{2}}\,-\,\,2\lambda \,\,-\,\,8\, \right)\,\,=\,\,0$           (42)

$\left( \lambda +2 \right)\left( \lambda \,+\,\,2)(\,\lambda \,-\,\,4\, \right)\,\,=\,\,0$

Therefore, the three eigenvalues or characteristic roots of the first row of Eq. (41) are:

${{\lambda }_{\,1}}\,=\,-2,\,\,{{\lambda }_{\,2}}\,=\,\,-2,\,\,\,\text{and}\,\,{{\lambda }_{\,3}}\,=\,\,4\,\,\,\,\,$

Therefore, the eigenvalues or characteristic roots of the entire WT power generation system comprising 1, 2, 3……n wind farms can be stated as:

$\begin{align}  & W{{F}_{1}}:\,\,\,{{\lambda }_{\,1}}\,=\,-2,\,\,{{\lambda }_{\,2}}\,=\,\,-2,\,\,\,\text{and}\,\,{{\lambda }_{3}}\,=\,\,4 \\  & W{{F}_{2}}:\,\,{{\lambda }_{\,1}}\,=\,-2,\,\,{{\lambda }_{\,2}}\,=\,\,-2,\,\,\,\text{and}\,\,{{\lambda }_{3}}\,=\,\,4\, \\  & W{{F}_{3}}:\,\,\,{{\lambda }_{\,1}}\,=\,-2,\,\,{{\lambda }_{\,2}}\,=\,\,-2,\,\,\,\text{and}\,\,{{\lambda }_{3}}\,=\,\,4\,\, \\  & \,\,\vdots \,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \\  & \,\,\vdots \,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \\  & \,\,\vdots \,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\vdots \,\,\,\,\,\vdots \,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \\  & W{{F}_{n}}:\,\,\,{{\lambda }_{\,1}}\,=\,-2,\,\,{{\lambda }_{\,2}}\,=\,\,-2,\,\,\,\text{and}\,\,{{\lambda }_{3}}\,=\,\,4\, \\ \end{align}$                 (43)

2.5 Hurwitz stability approach

Hurwitz criterion is an approach used to determine the stability of a system by using the coefficients of the characteristic polynomial Eq. without counting the roots involved [22-24]. Given the characteristic polynomial:

$p\,\,(\lambda )\,\,\,=\,\,{{a}_{\,n}}\,\,\lambda {{\,}^{n}}\,\,+\,{{a}_{n-1}}\lambda {{\,}^{n-1}}\,+...........+\,\,{{a}_{1}}\lambda \,\,\,+\,\,\,{{a}_{\,0}}$          (44)


${{a}_{0}},...........,{{a}_{n}}\in \mathbb{R},\,\,\,\,\,\,\,{{a}_{0}},\,\,{{a}_{n}}>\,\,0.$

The characteristic polynomial, p(λ) with real coefficients is stable or a Hurwitz polynomial if all its zeros have negative real parts [25, 26]. Hurwitz matrix, H for Eq. (44) is given by Eq. (45):

$H{{\,}_{n}}\,\,\,=\,\,\,\left[ \begin{align}  & {{a}_{\,1}}\,\,\,\,\,1\,\,\,\,\,\,\,0\,\,\,\,\,\,\,\,0\,........\,\,...............\,0\, \\  & {{a}_{\,3}}\,\,\,{{a}_{\,2}}\,\,\,\,a{{\,}_{1}}\,\,\,\,1\,\,\,\,........\,.........\,...\,\,\,\,\,0 \\  & a{{\,}_{5}}\,\,\,{{a}_{\,4}}\,\,\,\,a{{\,}_{3}}\,\,\,a{{\,}_{2}}\,\,\,.......\,\,........\,..\,\,\,\,\,0 \\  & a{{\,}_{7}}\,\,\,a{{\,}_{6}}\,\,\,{{a}_{\,\,5}}\,\,\,a{{\,}_{4\,\,\,\,\,\,}}......\,\,.......\,\,..\,\,\,\,\,0 \\  & \vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots \,\,\,\,\,\,\, \\  & {{a}_{2n-1}}\,\,{{a}_{2n-2}}\,\,{{a}_{2n-3}}\,\,\,\,\cdots \,\,\,\,\cdots \,\,\,a\,{{\,}_{n}} \\ \end{align} \right]$                      (45)

The stability is determined using the sub-determinants of the characteristic polynomial, p(λ) as follows [22-24]:

${{D}_{1}}\,=\,\,{{a}_{1}}>0$                (46)

${{D}_{2}}\,=\,\,\left| \begin{align}  & {{a}_{1}}\,\,\,\,\,\,\,{{a}_{0}} \\  & {{a}_{3}}\,\,\,\,\,\,\,{{a}_{2}} \\ \end{align} \right|\,\,\,=\,\,{{a}_{1}}{{a}_{2}}\,\,-\,\,{{a}_{0}}{{a}_{3}}>0$          (47)

${{D}_{3}}\,=\,\,\left| \begin{align}  & {{a}_{1}}\,\,\,\,\,\,\,{{a}_{0}}\,\,\,\,\,\,\,\,\,\,0 \\  & {{a}_{3}}\,\,\,\,\,\,\,{{a}_{2}}\,\,\,\,\,\,\,\,{{a}_{1}} \\  & 0\,\,\,\,\,\,\,\,\,\,\,0\,\,\,\,\,\,\,\,\,\,{{a}_{3}} \\ \end{align} \right|\,\,\,=\,\,{{a}_{3}}{{D}_{2}}>0$                  (48)

The combined characteristic polynomial Eq. of the WT power generation system and the dynamic loads obtained is given by Eq. (40):

${{\lambda }^{3}}\,\,-12\lambda \,\,-\,\,16\,\,=\,\,0$

Comparing Eq. (44) with Eq. (40), the coefficients are:

$a_{0}=1, a_{1}=0, a_{2}=-12, a_{3}=-16$.

The sub-determinants are therefore, given by:


${{D}_{2}}\,=\,\,\left| \begin{align}  & {{a}_{1}}\,\,\,\,\,\,\,{{a}_{0}} \\  & {{a}_{3}}\,\,\,\,\,\,\,{{a}_{2}} \\ \end{align} \right|\,\,\,=\,\,{{a}_{1}}{{a}_{2}}\,\,-\,\,{{a}_{0}}{{a}_{3}}$

${{D}_{\,2}}\,=\,\,\left| \begin{align}  & {{a}_{1}}\,\,\,\,\,\,\,{{a}_{0}} \\  & {{a}_{3}}\,\,\,\,\,\,\,{{a}_{2}} \\ \end{align} \right|\,\,\,=\,\,(0)\,(-12)\,\,-\,\,(1)\,(-16)\,\,=\,\,\,16$

${{D}_{\,3}}\,=\,\,\left| \begin{align}  & {{a}_{\,1}}\,\,\,\,\,\,\,{{a}_{\,0}}\,\,\,\,\,\,\,\,0 \\  & {{a}_{\,3}}\,\,\,\,\,\,\,a{{\,}_{2}}\,\,\,\,\,\,\,\,{{a}_{\,1}} \\  & 0\,\,\,\,\,\,\,\,\,\,\,0\,\,\,\,\,\,\,\,a{{\,}_{3}} \\ \end{align} \right|\,\,\,=\,\,{{a}_{\,3}}{{D}_{\,2}}\,\,=\,\,(-16)\,(-16)\,\,=\,\,256$

3. Results and Discussions

Natural response of a system yields either pure exponential decay or damped sinusoidal natural responses. These natural responses decay to zero as time approaches infinity. The eigenvalues of a particular dynamic system are used to determine whether that dynamic system is stable or unstable. The particular stability behaviour depends upon the nature of the eigenvalues. If the eigenvalues are two or more and are all negative, then the oscillation clearly decays to zero exponentially and the system is stable. On the other hand, if at least one of the eigenvalues is positive, the origin is unstable.

A real, positive, and distinct eigenvalue simply indicates that the system under consideration is unstable. Such real and positive eigenvalues show an exponential plot when graphed against time. If the set of eigenvalues for a system under consideration has both positive and negative eigenvalues components, is an indication of instability. If the set of eigenvalues for the system under consideration has repeated real eigenvalues, then the stability of such a system depends on whether the eigenvectors associated with the eigenvalues are linearly independent, or orthogonal. This is the case of degeneracy, where more than one eigenvector is associated with an eigenvalue. The determination of the system’s behaviour under such a circumstance requires further analysis.

The three eigenvalues or characteristic roots of the WT power generation system and the dynamic loads obtained are:

${{\lambda }_{1}}\,=\,\,-2,\,\,{{\lambda }_{2}}\,=\,\,-2,\,\,\,\text{and}\,\,{{\lambda }_{3}}\,=\,\,4\,\,\,\,\,$

The real negative eigenvalue value describes a decaying mode. The larger its magnitude, the faster the decay. A positive real eigenvalue represents an unstable system. So far as one of the eigenvalues obtained is real positive, is an indication that the proposed future WT power generation system shall be unstable despite the fact that two other eigenvalues obtained are negative.

For a system to be Hurwitz stable, the determinants must meet the following criteria: D1 > 0; D2 > 0; and D3 > 0. The results obtained in the Hurwitz analysis are:


From the results obtained; determinants D2 and D3 respectively meets the stipulated criteria with the exception of determinant, D1. The results obtained in the eigenvalue stability and the Hurwitz analysis suggest that the large-large WT power generation system might not be stable.

In summary, the findings of this work include:

A stochastic model of a large-scale wind farm based on power curve for a wind turbine. The power curve for a WT has three speeds namely, cut-in speed, rated output speed and cut-out speed.

Based on the stated model assumptions, a large -scale wind power system or grid might not be stable. This was proven through the model analysis approaches adopted.

4. Conclusion

A complex pair of eigenvalues or a real eigenvalue crossing over to the right half of the ‘S’ plane (real positive value) results in instability of the system being analysed. For the WT power generation system to be stable, all the eigenvalues must be real and negative. The WT power generation system analysed has both real negative and positive components; based on that more analysis is required to establish the actual stability of the system or otherwise. The assumptions made during the modelling stage can be revised to improve upon the model’s stability. Ascertaining the actual behaviour pattern will ensure that the appropriate control systems are designed for the stabilization of the WT power system based on the model and its eigenvalues obtained.



real numbers


greenhouse gases


carbon dioxide


sulphur dioxide


nitrogen oxide 


stochastic differential equations


noise-to-state stability


stochastic differential algebraic equations


renewable energy resources


master equation


wind turbine


wind farm






variable renewable energy


nitrogen oxide




real power


reactive power

Greek symbols




dimensionless transition rate for WT


dimensionless transition rate for WT


dimensionless transition rate for electrical load


dimensionless transition rate for electrical load


dimensionless transition rate


[1] Rashid, M.H., Kumar, N., Kulkarni, A.R. (2014). Power Electronics: Devices, Circuits, and Applications, 4th Edition. Pearson Education Ltd., Edinburgh Gate Harlow, England, UK.

[2] Tianze, L., Hengwei, L., Chuan, J., Luan, H., Xia, Z. (2011). Application and design of solar photovoltaic system. In Journal of Physics: Conference Series, 276(1): 012175.

[3] Beggs, C. (2009), Energy: Management, Supply and Conservation, 2nd Edition. Butterworth-Heinemann, Elsevier Ltd., UK.

[4] Kandasamy, K., Perumal, R., Velu, S.K. (2018). NSGA-II algorithm for solving combined environmental economic dispatch problem incorporating electric vehicles and wind energy. Journal of Electrical Engineering, 18(2): 12. 

[5] Tubiello, F.N., Salvatore, M., Rossi, S., Ferrara, A., Fitton, N., Smith, P. (2013). The FAOSTAT database of greenhouse gas emissions from agriculture. Environmental Research Letters, 8(1): 015009.

[6] Machowski, J., Lubosny, Z., Bialek, J.W., Bumby, J.R. (2020). Power System Dynamics: Stability and Control. John Wiley & Sons. 

[7] Alrikabi, N.K.M.A. (2014). Renewable energy types. Journal of Clean Energy Technologies, 2(1): 61-64. 

[8] Anon. (2019). The global climate in 2015-2019., accessed on November 05, 2019.

[9] Chen, C.J. (2011). Physics of Solar Energy. John Wiley and Sons, Inc., New Jersey, U.S.A.

[10] Moriarty, P., Honnery, D. (2012). What is the global potential for renewable energy? Renewable and Sustainable Energy Reviews, 16(1): 244-252.

[11] Lin, D. (2015). Methods for analyzing power system small signal stability (Doctoral dissertation, Memorial University of Newfoundland). 84-87. 

[12] Smith, S., Shahrezaei, V. (2015). General transient solution of the one-step master Eq. in one dimension. Physical Review E, 91(6): 062119.

[13] Ross, S.M. (2014). Introduction to Probability Models, 11th Edition. Academic Press, Oxford, UK.

[14] Rodriguez-Gonzalez, P.T., Rico-Ramirez, V., Rico- Ramirez, R., Diwekar, U.M. (2019). A new approach to solving stochastic optimal control problems. Mathematics 2019, 7(12): 1207.

[15] Gitano-Briggs, H. (2012). Low speed wind turbine design. advances in wind power., accessed on November 20, 2020. 

[16] Haag, G. (2017). Modelling with the master Eq. - solution methods and applications in social and natural sciences. Springer International Publishing AG, Cham, Switzerland.

[17] Oliveira, L.R. (2014). Master equation: biological applications and thermodynamic description. PhD Dessertation. University of Bologna. Bologna, Italy.

[18] Scott, M. (2013). Applied Stochastic Processes in Science and Engineering. University of Waterloo, Ontario, Canada.

[19] Padiyar, K.R. (2007). FACTS controllers in power transmission and distribution. New Age International (P) Ltd. Publishers, New Dehli, India.

[20] Hase, Y., Khandelwal, T., Kameda, K. (2020). Power System Dynamics with Computer-based Modeling and Analysis. John Wiley & Sons. 

[21] James, G. (2011). Advanced Modern Engineering Mathematics, 4th Edition, Pearson Education Ltd., Edinburgh Gate, England, UK.

[22] Erawaty, N., Amir, A. K. (2019). Stability analysis for Routh-Hurwitz conditions using partial pivot. In Journal of Physics: Conference Series, 1341(6): 062017.

[23] Mahardika, R., Sumanto, Y.D. (2019). Routh-hurwitz criterion and bifurcation method for stability analysis of tuberculosis transmission model. In Journal of Physics: Conference Series, 1217(1): 012056.

[24] Nise, N.S. (2015). Control systems engineering, 7th Edition. John Wiley and Sons, Inc., River Street, Hoboken, USA., accessed on April 04, 2021.

[25] Adm, M., Garloff, J., Tyaglov, M. (2018). Total nonnegativity of finite Hurwitz matrices and root location of polynomials. Journal of Mathematical Analysis and Applications, 467(1): 148-170.

[26] Kim, J.H., Su, W., Song, Y.J. (2018). On stability of a polynomial. Journal of Applied Mathematics and Informatics. 36(3-4): 231-236.