A novel approach for solving fuzzy differential equations using Runge-Kutta and Biogeography-based optimization

Abstract

The objective of this work is to present a methodology for solving the Kolmogorov’s differential equations in fuzzy environment using Runga-Kutta and Biogeography-based optimization (RKBBO) algorithm. As most of the data are usually uncertain and hence the fuzzy set theory has been used for handling the uncertainties and therefore the corresponding differential equations becomes fuzzy Kolmogorov’s differential equations. The obtained differential equations has been solved by Runga-Kutta method and improvement of their solution has been found by formulating a nonlinear optimization problem and solved by BBO for finding the global or near to global solutions. Performance index has been used for finding the most critical component of the system for increasing the performance of system by considering reliability, availability and maintainability index. Final results are compared with the existing results. To show the application of the proposed method, a case of the thermal power plant, a repairable industrial system, has been taken for evaluating the fuzzy reliability of the system.

Keywords

Biogeography-based optimization Runge-Kutta method differential equation fuzzy numbers reliability repairable system

1 Introduction

In the field of computational mechanics, many engineering problems can be viewed as a continuous domain and their physics are represented by using differential equations. These equations can be solved analytically; however, in most practical cases, the domains are too complicated for an analytical approach and hence analytical methods are not applicable or at least are very complicated to use it in solving the important problems in engineering and science, especially nonlinear ones. This leads to the creation of numerical approximation of the differential equations which are enabling up to find the approximate solutions of the differential equations. Several numerical techniques such as finite element, Euler, Runga-Kutta, etc. are used for finding the approximate solution of these types of differential equations. But today, the system becomes complex and complex these days and then the parameters which are associated with these differential equations becomes so complicated and massive and hence the traditional analysis does not give the real trend of the system due to unable to handle the uncertainties in the data. For handling the uncertainties issue, fuzzy set theory [33] has been used and hence the corresponding differential equations become fuzzy differential equations (FDEs) which was first formulated by Kaleva [22] and Seikala [30]. Based on these FDEs, the fuzzy initial value problem has been occurring in most of the science and engineering which are satisfying the fuzzy initial conditions. Nowadays, the topic of the FDE has been rapidly growing in recent years and an extensive number of papers have been published on the numerical solutions for various types of differential equations [1 , 26]. For instance, Garg [17] presented a method for solving the fuzzy Kolomorgov’s differential equation using Runga-Kutta method. Chakraverty and Tapaswini [11] proposed a method to solve nth order fuzzy differential equations using collocation method by using Legendre polynomials. Arqub et al. [7] presented a method for solving fuzzy differential equations based on the reproducing kernel theory under strongly generalized differentiability. Allahviranloo et al. [4] presented a fuzzy variational iteration method for solving the nonlinear fuzzy differential equations. Dhayabaran and Kingston [12] presented a method for finding the numerical solution of fuzzy differential equations using Runga-Kutta with harmonic mean of three quantities. Ahmadian et al. [3] presented a numerical algorithm for solving the first-order hybrid fuzzy differential equation based on high order Runga-Kutta method. Dong et al. [14] presented a solution of the linear-order random fuzzy differential equations to some particular initial value problems. However, each of these numerical techniques has its own operational restrictions, which strictly confine their functioning domain. In addition, most of them are based on classical mathematical tools and hence these approaches do not guarantee exact optimal solutions, but they achieve reasonably good solutions for hard problems in relatively short time periods. However, the heuristic techniques require derivatives for all non-linear constraint functions, that are not derived easily because of the high computational complexity. To overcome this difficulty meta-heuristics have been selected and successfully applied to handle a number of optimization problems. These heuristics include Genetic Algorithm (GA), particle swarm optimization (PSO), Biogeography-based optimization (BBO) and hence provides us a powerful tools for solving the most sophisticated problems and have strong tools to reinforce their problem-solving insight and expand their horizons [16 , 20]. Mastorakis [28] provided several new ideas on the formulation of ODEs as optimization problems. He utilized finite element concepts to develop the fitness functions. Many others also applied evolutionary algorithms to differential equations [10, 29]. As most of them have analyzed the performance by solving the system of equations in analytical form. As we know that the complexity of the system increases day-by-day and hence formulating and solving their corresponding differential equation up to a desired degree of accuracy is an challenging task for the system analyst or plant personal for maintaining the performance of the system. Thus, this paper provide an alternative way to achieve the solution of fuzzy differential equation by formulating non-linear optimization problem.

Biogeography-based optimization (BBO), proposed by Simon [31], is a new entrant in the domain of global optimization based on the theory of biogeography. The basic idea of BBO is based on the biogeography theory, which is the study of the geographical distribution of biological organisms. Different from other population based algorithms, in BBO, poor solutions can improve the qualities by accepting new features from good ones. BBO is developed through simulating the emigration and immigration of species between habitats in the multidimensional solution space, where each habitat represents a candidate solution. Like other evolutionary algorithms, BBO probabilistically shares information between candidate solutions. Biogeography not only gives a description of species distributions, but also a geographical explanation [21, 32]. Biogeography is modeled in terms of such factors as habitat area and immigration rate and emigration rate, and describes the evolution, extinction and migration of species. BBO has certain features in common with other population-based optimization methods like GA and PSO, BBO can share information between solutions. This makes BBO applicable to many of the same types of problems that GAs and PSO are used for, including unimodal, multimodal and deceptive functions. In BBO and PSO, each solution stays survive to the end of optimization procedure, but in most of evolutionary based algorithms, solutions die at the end of each generation. However, BBO also has some unique features that clearly differ from other population-based optimization methods. Also, in some of evolutionary due to crossover step, good solutions lose their efficiency, but in BBO do not have crossover step. One more characteristic of BBO is that its set of solutions is maintained and improved from one generation to the next by migrating. Also, for each generation, BBO uses the fitness of each solution to determine its emigration and immigration rate [18]. Moreover, the structure of BBO is similar to PSO in terms of maintaining the solutions from one iteration to others, but each solution can learn from its neighbors and adapt itself as long as the algorithm progresses. In PSO, the solution is updated indirectly through velocity vector while the solutions of BBO are changed directly via migration from other solutions, i.e., BBO solutions directly share their attributes with other solutions. After tests on many benchmarks, and comparisons with many other widely used heuristic algorithms like GAs, PSO and others, BBO outperformed most of the other algorithms on most of the benchmarks [18 , 32].

As far as availability/reliability field is concerned, almost all the authors have studied the behavior of the industrial system under the steady state conditions only. However, the system analyst always wants to desire the solution of the system under transient state also. Moreover, in order to achieve the desired goals of production, it is expected that production system should remain operative for maximum possible duration. But unfortunately, failure is an unavoidable phenomenon associated with technological products and systems. Over time, however, a given system suffers failures and has to be brought back into the serviceable state through appropriate maintenance and repairs [15, 17]. Aggarwal et al. [2] presented a method for solving Chapman-Kolmogorov differential equations by using Runga-Kutta fourth order method. Ding et al. [13] proposed a Petri net and ordinary differential equations (ODE)-based method for the system performance analysis in which system is modeled by a Petri net and then converted it into a set of ODE. Finally, the performance metrics can be derived from the ODE solutions that can be obtained in polynomial complexity. Malinowski [27] presented a fuzzy stochastic differential equations with solutions of decreasing fuzziness. For it, they used two kinds of sequences of approximate solutions and then showed that each sequence of approximate solutions can be used to prove existence and uniqueness of solution to fuzzy stochastic differential equations of decreasing fuzziness.

To this effect, the knowledge of the behavior of system, their component(s), is customary in order to plan and adapt suitable maintenance strategies. Therefore, the objective of this paper is to transform the fuzzy Kolomogorv differential equations problem into an unconstrained optimization problem while meta-heuristic or an evolutionary algorithm combined with Runga - Kutta method is employed to search for the optimal solution. In the proposed approach, firstly differential equations has been modeled in fuzzy environment for handling the uncertainties in the data. The solution of these differential equations has been computed, corresponding to each α- cut, by formulating a nonlinear optimization problem. Results obtained from proposed technique are compared with the existing Markov and Runga-Kutta method results. The approach has been illustrated by taking the case study of the thermal power plant, a complex repairable industrial system situated in the northern part of India. The effect of failure and repair rates on the system reliability has also been analyzed through the performance index and hence critical components of the system has been examined.

The rest of the manuscript has been described as follows. Section 2 describes the basic notations and assumptions that have used in the entire manuscript. Section 3 presents the proposed approach for solving the differential solution while Biogeography-based optimization algorithm is discussed in 4. In Section 5, approach has been illustrated through a case study of condensed system of the thermal power plant. The final results by considering the case of transient as well as steady state are discussed in Section 5.3 and the performance analysis for finding the critical component of the system has been analyzed in Section 6. Finally the concrete conclusion drawn are discussed in Section 7.

2 Assumptions and notations

2.1 Assumption

The following assumptions have been used in developing the probabilistic model.

The states of all components are mutually independent (statistically independent).

A repaired system is as good as new, performance wise, for a specified duration and standby unit are of same nature as that of active ones.

Failure rates and repair rates are represented by trapezoidal fuzzy numbers and are independent with each other.

At any given time, the system is either in operating state or in the failed state.

2.2 Notations

○

indicates the system is in full working state.

□

indicates the system is in failed state.

A, B, C, D, E, F

represent full working states of subsystems.

a, b, c, d, e, f

represent failed states of subsystems.

P₀ (t)

Probability of the system working with full capacity at time ‘t’.

P₁ (t)

Probability of the system in standby state at time ‘t’.

P₂ (t) to P₁₂ (t)

Probability of the system in failed state.

l_i, i = 1, 2, …, 6

failure rates of A,B,C,D,E,F subsystems respectively.

m_i, i = 1, 2, …, 6

repair rates of A,B,C,D,E,F subsystems respectively.

system availability

3 Proposed approach

The general nth order linear differential equation is represented as $\begin{matrix} a_{n} y^{(n)} \oplus a_{n - 1} y^{(n - 1)} \oplus \dots \oplus a_{1} y^{(1)} \oplus a_{0} y = g (x) \\ y^{(1)} (0) = γ_{0}, y^{(2)} (0) = γ_{1}, \dots, y^{(n - 1)} (0) = γ_{n - 1} \end{matrix}$ where the coefficients a_i and $y^{(i)} = \frac{d^{i} y}{{dx}^{i}}$ are crisp in nature. But in real life situations, their values are depending on the information available from various resources and hence they are imprecise in nature because the information collected from the various resources are generally either out of date or represent the past behavior of the data. Hence by using as such data in the analysis does not tell the exact behavior of the system. Therefore, to handle these issues, the above problem has been converted into fuzzy differential equation by representing each terms in terms of fuzzy numbers. Thus, the nth order fuzzy linear differential equation is represented as $\begin{matrix} {\tilde{a}}_{n} {\tilde{y}}^{(n)} \oplus {\tilde{a}}_{n - 1} {\tilde{y}}^{(n - 1)} \oplus \dots \oplus {\tilde{a}}_{1} {\tilde{y}}^{(1)} \oplus {\tilde{a}}_{0} \tilde{y} = \tilde{g} (x) \\ {\tilde{y}}^{(1)} (0) = \tilde{γ_{0}}, {\tilde{y}}^{(2)} (0) = \tilde{γ_{1}}, \dots, {\tilde{y}}^{(n - 1)} (0) = {\tilde{γ}}_{n - 1} \end{matrix}$ where ${\tilde{y}}^{(i)} = \frac{d^{i} \tilde{y}}{{dx}^{i}}$ for i = n, n - 1, …, 1 and ${\tilde{a}}_{n}$ is a non-zero fuzzy number and ${\tilde{a}}_{n - 1}$ , ${\tilde{a}}_{n - 2}, \dots, {\tilde{a}}_{1}$ , ${\tilde{a}}_{0}$ are any types of fuzzy numbers. In order to find the solution of these types of fuzzy differential equation, the following procedure may be adopted.

Find the α- cuts corresponding to each fuzzy parameters ${\tilde{a}}_{n}, {\tilde{a}}_{n - 1}, \dots, {\tilde{a}}_{1}, {\tilde{a}}_{0}$ , ${\tilde{y}}^{(n)}, {\tilde{y}}^{(n - 1)}, \dots, {\tilde{y}}^{(1)}, \tilde{y}$ and ${\tilde{γ}}_{0}, {\tilde{γ}}_{1}, \dots, {\tilde{γ}}_{n - 1}$ respectively, as follow. $\begin{matrix} {\tilde{a}}_{n} (x, α) = [a_{n (L)} (x, α), a_{n (R)} (x, α)]; \\ {\tilde{a}}_{n - 1} (x, α) = [a_{(n - 1) (L)} (x, α), a_{(n - 1) (R)} (x, α)]; \\ \dots, {\tilde{a}}_{0} (x, α) = [a_{0 (L)} (x, α), a_{0 (R)} (x, α)]; \\ {\tilde{y}}^{(n)} (x, α) = [y_{(L)}^{(n} (x, α), y_{(R)}^{(n)} (x, α)]; \\ {\tilde{y}}^{(n - 1)} (x, α) = [y_{(L)}^{(n - 1)} (x, α), y_{(R)}^{(n - 1)} (x, α)]; \\ \dots, \tilde{y} (x, α) = [y_{(L)} (x, α), y_{(R)} (x, α)] and \\ {\tilde{γ}}_{0} (0, α) = [γ_{0 (L)} (0, α), γ_{0 (R)} (0, α)]; \\ {\tilde{γ}}_{1} (0, α) = [γ_{1 (L)} (0, α), γ_{1 (R)} (0, α)]; \\ {\tilde{γ}}_{n - 1} (0, α) = [γ_{(n - 1) (L)} (0, α), γ_{(n - 1) (R)} (0, α)]; \end{matrix}$

Based on these α- cuts, the fuzzy initial value problem for the nth order fuzzy linear differential equation is converted into the following nth order differential equation: $\begin{matrix} [a_{n (L)} (x, α), a_{n (R)} (x, α)] [y_{(L)}^{(n)} (x, α), y_{(R)}^{(n)} (x, α)] \\ + [a_{n - 1 (L)} (x, α), a_{n - 1 (R)} (x, α)] [y_{(L)}^{(n - 1)} (x, α), \\ y_{(R)}^{(n - 1)} (x, α)] + \dots \dots + [a_{0 (L)} (x, α), a_{0 (R)} \\ (x, α)] [y_{(L)} (x, α), y_{(R)} (x, α)] \\ = [g (x), g (x)], \\ Initial Conditions \\ [y_{(L)} (0, α), y_{(R)} (0, α)] \\ = [γ_{0 (L)} (0, α), γ_{0 (R)} (0, α)], \dots, \\ [y_{(L)}^{n - 1} (0, α), y_{(R)}^{n - 1} (0, α)] \\ = [γ_{n - 1 (L)} (0, α), γ_{n - 1 (R)} (0, α)] \end{matrix}$ which can be rewritten in its equivalent form as $\begin{matrix} \sum_{j = 0}^{n} [a_{j (L)} (x, α), a_{j (R)} (x, α)] \\ [y_{(L)}^{(j)} (x, α), y_{(R)}^{(j)} (x, α)] = [g (x), g (x)] \\ [y_{(L)}^{(j)} (0, α), y_{(R)}^{(j)} (0, α)] = [γ_{j (L)} (0, α), γ_{j (R)} (0, α)] \end{matrix}$

The above formulated nth order differential equation is converted into the following ordinary differential equations.

$\begin{matrix} \begin{matrix} \sum_{j = 0}^{n} b_{j} y^{(j)} = g (x) \\ y_{(L)}^{(j)} (0, α) = γ_{j (L)} (0, α) \end{matrix}} and \\ \begin{matrix} \sum_{j = 0}^{n} c_{j} y^{(j)} = g (x) \\ y_{(R)}^{(j)} (0, α) = γ_{j (R)} (0, α) \end{matrix}} \end{matrix}$ (1) where $\begin{matrix} b_{j} y^{(j)} = min (a_{j (L)} (x, α) y_{(L)}^{(j)} (x, α), \\ a_{j (L)} (x, α) y_{(R)}^{(j)} (x, α), a_{j (R)} (x, α) y_{(L)}^{(j)} (x, α), \\ a_{j (R)} (x, α) y_{(R)}^{(j)} (x, α)) \\ and & c_{j} y^{(j)} = max (a_{j (L)} (x, α) y_{(L)}^{(j)} (x, α), \\ a_{j (L)} (x, α) y_{(R)}^{(j)} (x, α), a_{j (R)} (x, α) y_{(L)}^{(j)} (x, α), \\ a_{j (R)} (x, α) y_{(R)}^{(j)} (x, α)) \\ for j = 0, 1, \dots, n - 1, n \end{matrix}$

For a complex or large structural systems, the expression of these b_jy^(j) and c_jy^(j) are highly nonlinear in nature and hence if we solve the respective problem by using existing Laplace transformation methods then the corresponding results obtained do not tell the exact nature of the system. To overcome this problem, a nonlinear optimization problem has been constructed by utilizing the quantified fuzzy a_j’s and y^(j) at cut level α in the form of bounded interval as a decision variable. Once, the quantified input data at cut level α in the form of bounded interval is substituted in the expression of each b_j and c_j, the finally computed values at cut level α has wide range of solutions and it becomes smaller and smaller as the analysis progresses further i.e. cut level α increases from 0 to 1. Thus, the lower and upper boundary values of these systems of equations are computed at cut level α by solving the optimization problem (2) $\begin{matrix} Minimize / Maximize & F (a_{j}, y^{(j)}) \\ Subject to & μ_{a_{j}} \geq α \end{matrix}$ (2) $\begin{matrix} μ_{y^{(j)}} \geq α \\ 0 \leq α \leq 1 \end{matrix}$ where F is the fitness function which are obtained by solving the differential equation (1) using Runga-Kutta method. The obtained minimum and maximum value of F are denoted by F_min and F_max respectively and are both α that is, $\begin{matrix} μ_{\tilde{F}} (F_{min}) = μ_{\tilde{F}} (F_{max}) = α \end{matrix}$ Since the problem is nonlinear in nature, so it requires an efficient technique for finding the global or near to global solutions. A variety of methods and algorithms exist for optimization of such problems and applied in various technological fields during the last decades. BBO is one of the most popular and successful evolutionary algorithms which are used for finding the most compressible solution of the nonlinear optimization problem [18 , 32]. Thus, in the light of applicability, BBO is used as a tool to solve the optimization problem (2) in the process of determining the fuzzy membership function which has optimized spread. The objective function for maximization problem and the reciprocal of the objective function for minimization problem is taken as the fitness function. To stop the optimization process maximum number of generations and change in population fitness value are used as terminationcriteria.

The obtained values of y_(L) (x₀, α) and y_(R) (x₀, α), after solving the optimization problem (2), corresponding to real number x = x₀ represents the α- cut for the solution if the following conditions are met.

y_(L) (x₀, α) is a monotonically increasing function as α increases from 0 to 1.

y_(R) (x₀, α) is a monotonically decreasing function as α increases from 0 to 1.

If [y_(L) (x₀, α) , y_(R) (x₀, α)] defines the α- cut of a fuzzy number then the fuzzy solution of the fuzzy differential equation exits and [y_(L) (x₀, α) , y_(R) (x₀, α)] represents the α-cut corresponding to fuzzy solution

\tilde{y} (x_{0})

, otherwise fuzzy solution

\tilde{y} (x_{0})

of fuzzy differential equation does not exit.

4 Biogeography based optimization

Biogeography-based optimization (BBO), proposed by Simon [31], is a stochastic optimization technique for solving multimodal optimization problems. It is based on the concept of biogeography, which deals with the distribution of species that depend on different factors such as rain fall, diversity of topographic features, temperature, land area, etc. Biogeography describes how species migrate one island or habitat to another, how new species arise and how species become extinct. A habitat is an island which is geographically isolated from other habitats. The habitat having high suitability index (HSI) means that the habitat is geographically well suited for species living. The species have more possibility to migrate into the neighboring habitat when (i) the high HSI habitats have more species (ii) habitat has low HSI. This process is named as emigration. In immigration, the species move towards the high HSI habitat having few species. The immigration and emigration of species in a habitat are known as migration. The migration of species in a single habitat is explained in Fig. 1. From this figure, it has been seen that when there is no species in a single habitat, then the immigration rate (λ) of the habitat will be maximized (I). On the other hand, when the number of species increases, then only a few can successfully survive. Therefore, the immigration rate is decreased and become zero at S = S_max. Similarly, considering the emigration curve in the figure, the emigration rate (μ) increases with the increases of the species and hence the emigration rate is higher for the more crowded habitat than the less populated habitat. The optimization steps of the BBO algorithm are described as follows.

4.1 Initializing the BBO parameters and habitats

The optimization problem, in general, is expressed as $\begin{matrix} Max & f (x) \\ s . t . & g (x) \leq 0 \\ x_{i}^{min} \leq x_{i} \leq x_{i}^{max} i = 1, 2, \dots, n \end{matrix}$ where f (x) is the objective function, g (x) is the set of all constraints, x = (x₁, x₂, …, x_n) be the decision variables; $x_{i}^{min}$ and $x_{i}^{max}$ be the lower and upper limit of the variable x_i for i = 1, 2, …, n. The initialization phase of the BBO algorithm is to set the parameter of BBO algorithm. The parameters are maximum species count (S_max), maximum migration rate (E), immigration rate (I), the maximum mutation rate (m_max), habitat modification probability (p^mod), the number of habitats and the number of iterations.

In BBO algorithm, each individual or decision variable is considered as a “habitat” with a habitat suitability index (HSI) to measure the individual where the length of the habitat is depending on the number of decision variables. This HSI measures the goodness of the solution and is same as that of fitness value in population-based optimization algorithms. The corresponding variables of individuals that characterize habitability feature is called as suitability index variables (SIVs). These initialized randomly within the maximum and minimum limits as $x_{i}^{min} + Ud (x_{i}^{\max} - x_{i}^{min})$ where Ud is the uniformly distributed random number lies between (0, 1). Each habitat, in a population size NP, is represented by N- dimensional vector as H = [SIV₁, SIV₂, …, SIV_N] where N is the number of the SIVs (features) to be evolved for optimal HSI. HSI is the degree of acceptability that is determined by evaluating the objective function i.e.HSI = f (H).

4.2 Perform BBO operators

There are mainly two main operators: migration (which include both emigration and immigration) and mutation. Migration is the adaptive process which is used to change the existing habitat’s SIV. Each individual has its own immigration rate λ and emigration rate μ for improving the solution and thus evolve a solution to the optimization problem and is modified based on its probability P_mod, that is user defined parameter. Based on this probability, its immigration rate (λ) of each habitat is decided whether to modify each SIV’s or not. Another habitat is selected based on emigration rate μ, its SIV will migrate randomly to the selected habitat’s SIV. A good solution is a habitat with relatively high μ and low λ, while the converse is true for a poor solution. High HSI habitat resists change more the low HSI habitat. Poor solutions accept more useful information from good solution, which improve the exploitation ability of the algorithm. At initially, the immigration and emigration rates are the functions of the number of species in the habitat and they are evaluated as $\begin{matrix} λ_{k} = I (1 - \frac{k}{S_{max}}); μ_{k} = E (\frac{k}{S_{max}}) \end{matrix}$ where I is the maximum possible immigration rate; E is the maximum possible emigration rate; k is the number of the species of the kth individual; and S_max is the maximum number of species. I and E are often set equal to 1, or slightly less than 1. The pseudo code of the migration operator is depicted in Algorithm 1.

On the other hand, in BBO, mutation is used to increase the diversity of the population to get the best solution. The low and high HSI valued habitats have less possibility to mutate while average HSI valued habitats have more possibility to mutate. For instance, consider there are a total of S_max habitat and the Sth habitat, there are S species. For the Sth island, the immigration rate is λ_s and emigration rate is μ_s. By considering P_S as the probability that there are exactly S species live in a habitat so P_S changes from time t to t + ▵ t is given as

$\begin{matrix} P_{S} (t + ▵ t) = P_{S} (t) (1 - λ_{s} ▵ t - μ_{s} ▵ t) + \\ P_{S - 1} (t) λ_{s - 1} ▵ t + P_{S + 1} (t) μ_{k + 1} ▵ t \end{matrix}$ (3)

By taking the limit of Equation (3) as t → 0, we have $\begin{matrix} \dot{P_{S}} = {\begin{matrix} - (λ_{s} + μ_{s}) P_{s} + μ_{S + 1} P_{S + 1} & ; S = 0 \\ - (λ_{s} + μ_{s}) P_{s} + μ_{S - 1} P_{S - 1} + μ_{S + 1} P_{S + 1} & ; 1 \leq S \leq S_{max} - 1 \\ - (λ_{s} + μ_{s}) P_{s} + μ_{S - 1} P_{S - 1} & ; S = S_{max} \end{matrix} \end{matrix}$ where P_S-1, P_S and P_S+1 are the species count probability; λ_S-1, λ_S, λ_S+1 are immigration rate and μ_S-1, μ_S, μ_S+1 are emigration rate of the habitat with S - 1, S and S + 1 species respectively, S_max is the maximum species count in the habitat. The steady state probability of the quantity of each species is given as below $P_{0} = {[1 + \sum_{j = 1}^{S_{max}} \prod_{i = 1}^{l} \frac{λ_{i - 1}}{μ_{i}}]}^{- 1}$ (4) $P_{k} = P_{0} \prod_{i = 1}^{k} (\frac{λ_{i - 1}}{μ_{i}})$ (5) $\begin{matrix} = & \frac{1}{1 + \sum_{j = 1}^{S_{max}} \prod_{i = j}^{l} \frac{μ_{i}}{λ_{i - 1}} + \sum_{l = k + 1}^{S_{max}} \prod_{i = k + 1}^{l} \frac{λ_{i - 1}}{μ_{i}}}, \\ k \in {1, \dots, S_{max}} \end{matrix}$

Thus, based on that the randomly selected SIV habitat modified it based on mutation rate m, in the case of E = I, and is expressed as $\begin{matrix} m_{i} = m_{max} \times (1 - \frac{P_{i}}{p_{max}}) \end{matrix}$ where m_max ∈ [0, 1] is the predefined parameter called as maximum mutation probability and p_max = max(P₁, P₂, …, P_NP), NP is the number of habitats, P_i is the probability of habitat i which can be calculated by Equations (4) and (5). This mutation scheme tends to enhance the diversity among the population, which helps to decrease the chance of getting trapped in local optima. The pseudo code of the mutation operator is depicted in Algorithm 2.

4.3 Update the habitat

The generation process is repeated until a number of generations has been met, as a termination criteria.

5 Case study

The above mentioned technique for analyzing the system reliability from the fuzzy kolmogorov’s differential equation is illustrated through a case study of thermal power plant, a repairable industrial system situated in the Northern part of India. Thermal plants are large capital-oriented engineering systems in which each of them comprises subsystems/units namely coal handling, steam generation, cooling water, crushing, ash handling, power generation, feed water system etc., arranged in predetermined configurations. The present analysis is based on the study of one of the important unit i.e., condensed system [17, 23]. The brief description of this unit is as follows.

5.1 System description

A thermal power plant is a complex engineering system comprising of various systems: Coal handling, Steam Generation, Cooling Water, Crushing, Ash handling, Power Generation and Feed water system. The components of the condensed system in thermal power plant is described as below [17, 23]:

Condenser (A): It consists of single unit arranged in series configuration such that failure of it will cause the entire failure of the system.

Gland steam condenser (B): It consists of single unit connected in series with condenser. Failure of this unit causes the complete failure of the system.

Drain Cooler (C): It consists of one unit arranged in series. Failure of this unit causes the complete failure of the system.

Heaters (D): In this subsystem three units of low pressure heaters are arranged in series configuration. Failure of any one will cause the complete failure of the system.

Deaerator (E): It consists of single unit arranged in series. Failure of this unit causes the complete failure of the system.

Extraction pumps (F): Herein two units of condensate extraction pumps are arranged in parallel configuration with one operative and other in cold standby. Failure of any one will reduce the efficiency of the system. Complete failure of pumps occurs when both the components will fail.

The failure and repair rates corresponding to each main component of the system, are summarized in Table 1 [17, 23], in the form of trapezoidal fuzzy numbers for evaluating the fuzzy availability of the system. On the other hand, the transition diagram for representing the system component behavior is described in Fig. 2.

5.2 Fuzzy Kolmogorov’s differential equations for the system

Let the probability of n occurrences in time t is denoted by P_n (t) i.e. P (x = n, t) = P_n (t), n = 0, 1, 2, …. Then P₀ (t) represents the probability of zero occurrences in time t. The probability of zero occurrences in time (t + ▵ t) is given by

$P_{0} (t + ▵ t) = (1 - l ▵ t) P_{0} (t) + (m ▵ t) P_{1} (t)$ (6)

Similarly,

$P_{1} (t + ▵ t) = (l ▵ t) P_{0} (t) + (1 - m ▵ t) P_{1} (t)$ (7) where l and m represents the failure and repair rates of the system. The above equations shows the probability of one occurrence in time (t + ▵ t) and is composed of two parts, namely, (i) probability of zero occurrences in time t multiplied by the probability of one occurrence in the interval ▵t and (ii) the probability of one occurrence in time t multiplied by the probability of no occurrences in the interval ▵t. Thus, by taking ▵t → 0, we get $\begin{matrix} (\frac{d}{dt} + l) P_{0} (t) = m P_{1} (t) \end{matrix}$

Using this concept the following differential equations associated with the transition diagram (Fig. 2) of condensed system are formed in fuzzy environment.

$\begin{matrix} \frac{d {\tilde{P}}_{0} (t)}{dt} \oplus {\tilde{δ}}_{1} {\tilde{P}}_{0} (t) = \sum_{j = 1}^{5} {\tilde{m}}_{j} {\tilde{P}}_{j + 1} (t) \oplus {\tilde{m}}_{6} {\tilde{P}}_{1} (t) \\ \frac{d {\tilde{P}}_{1} (t)}{dt} \oplus {\tilde{δ}}_{2} {\tilde{P}}_{1} (t) = \sum_{j = 1}^{6} {\tilde{m}}_{j} {\tilde{P}}_{j + 6} (t) \oplus {\tilde{l}}_{6} {\tilde{P}}_{0} (t) \\ \frac{d {\tilde{P}}_{2} (t)}{dt} \oplus {\tilde{m}}_{1} {\tilde{P}}_{2} (t) = {\tilde{l}}_{1} {\tilde{P}}_{0} (t), \\ \frac{d {\tilde{P}}_{3} (t)}{dt} \oplus {\tilde{m}}_{2} {\tilde{P}}_{3} (t) = {\tilde{l}}_{2} {\tilde{P}}_{0} (t), \\ \frac{d {\tilde{P}}_{4} (t)}{dt} \oplus {\tilde{m}}_{3} {\tilde{P}}_{4} (t) = {\tilde{l}}_{3} {\tilde{P}}_{0} (t), \\ \frac{d {\tilde{P}}_{5} (t)}{dt} \oplus {\tilde{m}}_{4} {\tilde{P}}_{5} (t) = {\tilde{l}}_{4} {\tilde{P}}_{0} (t), \\ \frac{d {\tilde{P}}_{6} (t)}{dt} \oplus {\tilde{m}}_{5} {\tilde{P}}_{6} (t) = {\tilde{l}}_{5} {\tilde{P}}_{0} (t), \\ \frac{d {\tilde{P}}_{7} (t)}{dt} \oplus {\tilde{m}}_{1} {\tilde{P}}_{7} (t) = {\tilde{l}}_{1} {\tilde{P}}_{1} (t), \\ \frac{d {\tilde{P}}_{8} (t)}{dt} \oplus {\tilde{m}}_{2} {\tilde{P}}_{8} (t) = {\tilde{l}}_{2} {\tilde{P}}_{1} (t), \\ \frac{d {\tilde{P}}_{9} (t)}{dt} \oplus {\tilde{m}}_{3} {\tilde{P}}_{9} (t) = {\tilde{l}}_{3} {\tilde{P}}_{1} (t), \\ \frac{d {\tilde{P}}_{10} (t)}{dt} \oplus {\tilde{m}}_{4} {\tilde{P}}_{10} (t) = {\tilde{l}}_{4} {\tilde{P}}_{1} (t), \\ \frac{d {\tilde{P}}_{11} (t)}{dt} \oplus {\tilde{m}}_{5} {\tilde{P}}_{11} (t) = {\tilde{l}}_{5} {\tilde{P}}_{1} (t), \\ \frac{d {\tilde{P}}_{12} (t)}{dt} \oplus {\tilde{m}}_{6} {\tilde{P}}_{12} (t) = {\tilde{l}}_{6} {\tilde{P}}_{1} (t), \end{matrix}}$ (8) where ${\tilde{δ}}_{1} = \sum_{i = 1}^{6} {\tilde{l}}_{i}$ and ${\tilde{δ}}_{2} = {\tilde{δ}}_{1} \oplus {\tilde{m}}_{6}$ with fuzzy initial conditions ${\tilde{P}}_{0} (0) = (0.95, 0.955, 0.965, 0.97)$ , ${\tilde{P}}_{1} (0) = (0.004, 0.0045, 0.0055, 0.006)$ and ${\tilde{P}}_{j} (0) = (0, 0, 0, 0), j = 2 to 12$ .

As the system of differential equation is very complex, it is difficult to find its analytical solution, particularly using Laplace transformation method. Thus, in order to get the desired availability of the system so that it will work for a longer period of time so as to maximize profit under the various combinations of the failure and repair rates of the systems’ component. For this, initially the differential Equation (8) has been solved by using Runga-Kutta method and hence based on that an optimization model has been formulated for improving the availability of the system by taking an objective function $\tilde{Av} (t) = {\tilde{P}}_{0} (t) \oplus {\tilde{P}}_{1} (t)$ . The optimal values of these availability function has been computed by using BBO algorithm.

5.3 Results and discussion

In the present analysis, an algorithm is implemented in Matlab (MathWorks) and the program has been run on a T6400 @ 2GHz Intel Core(TM) 2 Duo processor with 2GB of Random Access Memory. The following BBO parameters have been used for evaluating the performance of the system: habitat modification probability p^mod= 0.8, maximum immigration and migration rates of each island is 1 and mutation probability m_max is 0.002. Habitat size is set to be randomly as 20 × D where D is the dimension of the problem and the maximum number of the iterations are set to be 1000. In order to eliminate stochastic discrepancy, in each case study, 30 independent runs are made for each of the optimization methods involving 30 different initial trial solutions, and their best values are reported. The termination criterion has been set either limited to a maximum number of generations (k_max = 1000) or to the order to a relative error equal to 10^-6, whichever is achieved first. Under the information extraction phase, the data related to failure rates (l_i) and repair rates (m_i) are collected from historical records related to the main components of the system. As most of the collected data are usually uncertain or represent the past behavior of the system so to handle these uncertainties in the data, the data are represented in the form of trapezoidal fuzzy numbers and are summarized in Table 1. Using these data, the system behavior has been analyzed by following the steps of the proposed approach for mission time t = 48 hrs with left and right spreads with an increment of 0.1 confidence level. The computed results by proposing approach RKBBO are tabulated in Table 2 along with the existing Markov and Runga - Kutta results and from these results it has been concluded that

The results computed by the traditional method (crisp) [24] do not give an exact idea about the behavior of the system. As these theories are independent of the degree of the confidence level and thus the result obtained by these methodologies do not consider the vagueness/uncertainties in the data. Therefore, these methods are suitable only for a system whose data are precise.

The results computed by Kumar and Lata [23] and Garg [17] by simultaneously handling the uncertainties in the data in the form of trapezoidal fuzzy numbers using Runga-Kutta approach. From their corresponding results, it has been observed that the level of uncertainties are decreases in [17] as compared to [23]. For instance, the system availability corresponding to α = 0.5 are lies in the interval [0.8227397, 0.8369285] by the Garg [17] approach while by Kumar and Lata [23] approach their interval are [0.822748, 0.836866]. But their approach are suitable only for a components whose system structure is small or precisely known. As the components of the system increases then the analysis have contain a large amount of uncertainties during the analysis.

On the other hand, the results computed by the proposed approach at different levels of significance has been computed for the system availability at t = 48 hrs and found that the proposed results gives a better range than other methodologies [17 , 24] ranges at any cut level of significance. For instance, the system availability corresponding to α = 0.5 are lies in the interval [0.8233522, 0.8359435] by the proposed approach while by Garg [17], Kumar and Lata [23] and Kumar and Aggarwal [24] approaches this interval is [0.8227397, 0.8369285], [0.822748, 0.836866] and [0.824961, 0.834418] respectively. Thus, their corresponding average value is 0.8296478 by proposed approach while 0.8298341 by Garg [17] approach, 0.829807 by Kumar and Lata [23] approach and 0.8296895 by traditional approach [24]. Therefore, if the system analyst or plant personnel want to improve the performance of the system, then the new target of the system at 0.5 level of significance would be greater than 0.8296478 rather than other values. This means that uncertainties existing during the analysis are reduced up to the desired degree and hence decision makers/ system analysts may use these results for predicting the behavior of the system in a more realistic manner which leads to a more sound and effective decision for future course of actions in lesser time. Similarly the analysis has been done for different level of significance. The reason behind is that BBO gives near to the optimal solution. In order to analyze the decrease in spread (in % ) by proposed approach from the existing results, an analysis has been done, which computes the support of the availability parameters and found that there is 16.939% and 17.567% decrease in the spread from proposed results as compared to existing Markov and Runga-Kutta methods results. Moreover, it is clear from that results that the sides of membership functions of availability parameters are parabolic, not linear as were taken as initially. The defuzzified value of the parameters computed from their Markov, Runga-Kutta and proposed results are 0.82980750, 0.829844625 and 0.829852375. Thus, based on these defuzzified values plant personnel may have changed their target goals from traditional analysis. For instance, in order to improve the availability of the system, then the new target of system availability should be greater than 0.829852375 rather than other values comes from existing approaches. Due to this and their reduced region of prediction, the value obtained through proposed technique are conservative in nature which may be beneficial for a system expert/ analyst for future course of action, i.e. now the maintenance will be based on the defuzzified values rather than crisp values. Thus, we can also say that the proposed approach is beneficial for the system analyst for increasing the performance of the system. The complete results are summarized in Table 2 which indicates that the results obtained by the proposed approach acts as a bridge between the Markov approach (crisp values) and other existing technique.

6 Performance analysis

For depicting the most critical component of the system for increasing the performance and to save money, manpower, time, an analysis has been conducted by defining the composite measure of reliability, availability and maintainability parameters and named as the Performance - Index. The major advantages of this index are that it simultaneously shows the impact of the various parameters on the performance of the system by varying their corresponding failure and repair rates. The mathematical expression of this index is given as $\begin{matrix} Performance Index = w_{1} R_{s} + w_{2} {Av}_{s} + w_{3} M_{s} \end{matrix}$ where R_s, Av_s and M_s represents the reliability, availability and maintainability of the system while w_i, i = 1, 2, 3 are their corresponding weight such that $\sum_{i = 1}^{3} w_{i} = 1, w_{i} \in [0, 1]$ . Based on this index, the system analyst/plant personnel may have predicted the impact of components’ parameters on its performance. The effect of the each component parameter on the system performance index has been analyzed for different values of their weights W₁, W₂ and W₃. Here, W₁ = [0.33, 0.34, 0.33] signifies that DM does not want the preference towards the objectives; that is, they pay equal attention towards each objective, weight set W₂ = [0.25, 0.50, 0.25] signifies that DM wants to pay two times more attention towards the availability parameters than other two while the set W₃ = [0, 1, 0] shows that the decision maker give full weightage to the system availability parameters. Performance analysis has been conducted with initial preference of the weights towards the objectives, as W₁ and results corresponding to these sets are summarized in Table 3. Not satisfied with this outcome or willing to know other possible solutions, keeping this results in view, let DM be given more importance on availability parameter i.e. weight set W₂. Similarly, process has been repeated for different combinations of weight for different satisfaction levels. The complete summary of the results along with components ranges of their failure and repair rates are shown in Table 3. Based on this analysis, it has been concluded that for improving the performance of the system, more attention should be given to the components as per the preferential order; extraction pumps, Gland steam condenser, Heaters, condenser, Deaerator and drain cooler. Thus the system manager can analyze the critical behavior of the system and hence can plan for suitable maintenance strategy for achieving the highest goals. The maintenance engineer may therefore determine the repair policy, according to the optimal design by considerations and the company strategies.

7 Conclusion

An attempt has been made in this article for analyzing the behavior and performance of the system by using the concept of Runga-Kutta and Biogeography-based optimization. The data corresponding to the system components are measured in terms of failure and repair rates while the uncertainties corresponding to them are handled with the help of defining the trapezoidal fuzzy number. A condensed system from the thermal power plant, a repairable industrial system has been taken for illustration. The functioning/non-functional states of the system are modeled by a Markov process with constant failure and repair rates. The obtained fuzzy Markov model has been solved by formulated a nonlinear optimization model in which Runga-Kutta method has been used for finding the initial solution and that solutions have been improved by using the BBO algorithm. A promisable solution has been obtained by the proposed approach and compared their results with the existing techniques results. From these results, it has been concluded that the obtained results have less range of uncertainties as compared to others and hence beneficial for the system analyst for improving the performance of the system. Moreover, the most critical component of the system has been examined by performing an experiment by using the performance index. Different preferences on the system reliability, availability and maintainability index has been given and hence based on that the performance of the system has been analyzed. Based on this analysis, system analyst pay attention towards the components as per the preferential order; extraction pumps, Gland steam condenser, Heaters, condenser, Deaerator and drain cooler. Thus the optimal design parameters can facilitate the system design analyst to select the suitable design policy so that the system can achieve high production goals while maintaining its performance.

References

Abbasbandy

and Allahviranloo

, Numerical solutions of fuzzy differential equations by taylor method, Journal of Computational Methods in Applied Mathematics2 (2002), 113–124.

Aggarwal

A.K.

, Kumar

, Singh

and Garg

T.K.

, Markov modeling and reliability analysis of urea synthesis system of a fertilizer plant, Journal of Industrial Engineering International11(1) (2015), 1–14.

Ahmadian

, Salahshour

and Chan

C.S.

, A runga-kutta method with reduced number of function evaluations to solve hybrid fuzzy differential equations, Soft Computing19 (2015), 1051–1062.

Allahviranloo

, Abbasbandy

and Behzadi

, Solving nonlinear fuzzy differential equations by using fuzzy variational iteration method, Soft Computing18(11) (2014), 2191–2200.

Allahviranloo

, Kiani

N.A.

and Motamedi

, Solving fuzzy differential equatios by differential transformation method, Information Sciences179(7) (2009), 956–966.

Allahviranloo

and Salahshour

A new approach for solving first order fuzzy differential equation, in: Hullermeier

, Kruse

, Hoffmann

(eds.), Information Processing and Management of Uncertainty in Knowledge-Based Systems. Applications. vol. 81 of Communications in Computer and Information Science, Springer, Berlin Heidelberg, 2010, pp. 522–531.

Arqub

O.A.

, AL-Smadi

, Momani

and Hayat

, Numerical solutions of fuzzy differential equations using reproducing kernel hilbert space method, Soft Computing (2015) 1–20. http://dx.doi.org/10.1007/s00500-015-1707-4

Buckley

J.J.

and Feuring

, Introduction to fuzzy partial differential equations, Fuzzy sets and Systems105 (1999), 241–248.

Buckley

J.J.

and Feuring

, Fuzzy differential equations, Fuzzy sets and Systems110 (2000), 43–54.

10.

Cao

, Kang

and Chen

, Evolutionary modeling of systems of ordinary differential equations with genetic programming, Genetic Programming and Evolvable Machines1 (2000), 309–337.

11.

Chakraverty

and Tapaswini

, Standard orthogonal polynomialsbased solution of fuzzy differential equations, International Journal of Fuzzy Computation and Modelling1(1) (2014), 51–71.

12.

Dhayabaran

D.P.

and Kingston

J.C.

, A method for solving fuzzy differential equations using runge-kutta method with harmonic mean of three quantities, International Journal of Engineering Science and Innovative Technology4 (2015), 90–96.

13.

Ding

, Zhou

and Zhou

, A polynomial algorithm to performance analysis of concurrent systems via petri nets and ordinary differential equations, IEEE Transactions on Automation Science and Engineering12(1) (2015), 295–308.

14.

Dong

S.L.

, Vu

and Nguyen

, The formulas of the solution for linear-order random fuzzy differential equations, Journal of Intelligent and Fuzzy Systems28(2) (2015), 795–807.

15.

Garg

, An approach for analyzing fuzzy system reliability using particle swarm optimization and intuitionistic fuzzy set theory, Journal of Multiple-Valued Logic and Soft Computing21 (2013), 335–354.

16.

Garg

, Performance analysis of reparable industrial systems using PSO and fuzzy confidence interval based methodology, ISA Transactions52(2) (2013), 171–183.

17.

Garg

, An approach for analyzing the reliability of industrial system using fuzzy kolmogorov’s differential equations, The Arabian Journal for Science and Engineering40(3) (2015), 975–987.

18.

Garg

, An efficient biogeography based optimization algorithm for solving reliability optimization problems, Swarm and Evolutionary Computation24 (2015), 1–10.

19.

Garg

and Sharma

S.P.

, Stochastic behavior analysis of industrial systems utilizing uncertain data, ISA Transactions51(6) (2012), 752–762.

20.

Garg

and Sharma

S.P.

, Multi-objective reliability-redundancy allocation problem using particle swarm optimization, Computers & Industrial Engineering64(1) (2013), 247–255.

21.

Gong

, Cai

, Ling

and Li

, A real-coded biogeographybased optimization with mutation, Applied Mathematics and Computation216(9) (2010), 2749–2758.

22.

Kaleva

, Fuzzy differential equations, Fuzzy Sets and Systems24(3) (1987), 301–317.

23.

Kumar

and Lata

, Reliability evaluation of condensate system using fuzzy markov model, Annals of Fuzzy Mathematics and Informatics4(2) (2012), 281–291.

24.

Kumar

and Aggarwal

K.K.

, Petri net modeling and reliability evaluation of distributed processing systems, Reliability Engineering and System Safety41(2) (1993), 167–176.

25.

and Simon

, Blended biogeography-based optimization for constrained optimization, Engineering Applications of Artificial Intelligence24(3) (2011), 517–525.

26.

, Friedman

and Kandel

, Numerical solutions of fuzzy differential equations, Fuzzy sets and Systems105 (1999), 133–138.

27.

Malinowski

M.T.

, Fuzzy stochastic differential equations of decreasing fuzziness: Approximate solutions, Journal of Intelligent and Fuzzy Systems29(3) (2015), 1087–1107.

28.

Mastorakis

N.E.

, Unstable ordinary differential equations: Solution via genetic algorithms and the method of nelder-mead, Stevens Point, in: Proceedings of the 6th WSEAS International Conference on Systems Theory & Scientific Computation, World Scientific and Engineering Academy and Society (WSEAS), Wisconsin, USA, 2007.

29.

Reich

, Simulation of imprecise ordinary differential equations using evolutionary algorithms, in: Proceedings of the 2000 ACM Symposium on Applied Computing1, 2000.

30.

Seikala

, On the fuzzy intial value problem, Fuzzy Sets and Systems24 (1987), 319–330.

31.

Simon

, Biogeography - based optimization, IEEE Transactions on Evolutionary Computation12(6) (2008), 702–713.

32.

Simon

, Ergezer

, Du

and Rarick

, Markov models for biogeography-based optimization, IEEE Transactions on Systems, Man, and Cybernetics - Part B: Cybernetics41(1) (2011), 299–306.

33.

Zadeh

L.A.

, Fuzzy sets, Information and Control8 (1965), 338–353.