RomAniuk On simulatiOn Of maintenance cOsts fOr water distributiOn system with fuzzy parameters

The main aim of each water distribution system (WDS, water pipes network) is to deliver water of desirable quality and assumed quantity for customers. In order to meet these requirements, maintenance services, like repairs and replacements of broken or malfunctioned parts, are required. The literature devoted to the problems of simulation of the WDS behaviour and to analysis of its reliability is very rich (see, e.g, [3, 4, 7, 10, 11, 15, 17, 19, 25, 29] for a description of various approaches, and [16, 28] for a more detailed review of methods and literature). Because of a necessary long-time horizon planning (e.g. 20 or even 50 years), an influence of a value of the money in the future on maintenance costs should be taken into account. Also such a problem is addressed in many papers. For example, in [20], the main aim is to estimate and validate various cost functions for different types of assets of a WDS, if their hydraulic (e.g., flow, pump head, pump power) and physical characteristics (e.g., volume, material, nominal diameter) are stored in a specially prepared data base. A method of a linear regression is used to model a dependency between various types of the costs (like an equipment cost) and the mentioned characteristics of a pipe (like a nominal diameter) based on data from several Portuguese urban water utilities. However, a method for a calculation of a present value of the total costs is not developed there in a more detailed way. Contrary to the previous paper, in [12], let’s say, a “macro-management” approach for a WDS rehabilitation problem with a real long-time horizon is widely discussed. In this paper, both an economics and hydraulic capacity of a WDS is analysed, if a deterioration of a pipe follows the Hazen-Williams equation. The total costs are related RomAniuk m. on simulation of maintenance costs for water distribution system with fuzzy parameters. Eksploatacja i niezawodnosc – maintenance and Reliability 2016; 18 (4): 514–527, http://dx.doi.org/10.17531/ein.2016.4.6.


Introduction
The main aim of each water distribution system (WDS, water pipes network) is to deliver water of desirable quality and assumed quantity for customers.In order to meet these requirements, maintenance services, like repairs and replacements of broken or malfunctioned parts, are required.
The literature devoted to the problems of simulation of the WDS behaviour and to analysis of its reliability is very rich (see, e.g, [3,4,7,10,11,15,17,19,25,29] for a description of various approaches, and [16,28] for a more detailed review of methods and literature).Because of a necessary long-time horizon planning (e.g.20 or even 50 years), an influence of a value of the money in the future on maintenance costs should be taken into account.Also such a problem is addressed in many papers.
For example, in [20], the main aim is to estimate and validate various cost functions for different types of assets of a WDS, if their hydraulic (e.g., flow, pump head, pump power) and physical characteristics (e.g., volume, material, nominal diameter) are stored in a specially prepared data base.A method of a linear regression is used to model a dependency between various types of the costs (like an equipment cost) and the mentioned characteristics of a pipe (like a nominal diameter) based on data from several Portuguese urban water utilities.However, a method for a calculation of a present value of the total costs is not developed there in a more detailed way.
Contrary to the previous paper, in [12], let's say, a "macro-management" approach for a WDS rehabilitation problem with a real long-time horizon is widely discussed.In this paper, both an economics and hydraulic capacity of a WDS is analysed, if a deterioration of a pipe follows the Hazen-Williams equation.The total costs are related RomAniuk m. on simulation of maintenance costs for water distribution system with fuzzy parameters.Eksploatacja i niezawodnoscmaintenance and Reliability 2016; 18 (4): 514-527, http://dx.doi.org/10.17531/ein.2016.4.6.

maciej RomAniuk
On simulatiOn Of maintenance cOsts fOr water distributiOn system with fuzzy parameters O symulOwaniu kOsztów utrzymania dla systemu dystrybucji wOdy O parametrach rOzmytych* In this paper we propose a model for evaluation of maintenance costs of a water distribution system (WDS).The set of possible states of each connection (i.e. a pipeline in the WDS) is related to various possible degrees of quality of the pipe and types of its malfunctions.The process of transitions between these states forms a semi-Markov process.Using Monte Carlo simulations, the length of services and the times of necessary replacements and repairs of the connections are obtained.These values are then used as an input for estimation of the maintenance costs of the whole WDS.During this step we take into account the concept of present value of money.Contrary to other approaches, instead of a constant yield, a stochastic process (the one-factor Vasicek model) of an interest rate is assumed.Then various simulated measures of reliability and the maintenance costs can be analysed, like an influence of various parameters of the pipes (e.g.intensities of damages) on the final costs of the performed services.They can be crucial in the analysis of risk for various possible decisions.Apart from the crisp approach, the Monte Carlo simulations are also applied, if some of the parameters of the WDS are fuzzified.Therefore uncertainty and experts' knowledge can be easily incorporated into the proposed procedure of the estimation of the maintenance costs.Observed differences between the crisp and the fuzzy output are highlighted.Simulation algorithms, necessary for both of these approaches, are also provided.
Keywords: water distribution system, maintenance costs, semi-Markov process, Monte Carlo simulations, present value.

sciENcE aNd tEchNology
to a rehabilitation of a pipe and some additional maintenance costs (mainly breakage repair costs).To model the breakage rate, the exponential dependency on an age of a pipe is assumed.The present value of the total costs is then calculated using a constant discount rate.
The main aim of the whole procedure is to minimize the total costs associated with a rehabilitation, if some hydraulic constraints, necessary for a quality and a quantity of the supplied water, are preserved.As mentioned in this paper, an analysis period for such a procedure is typically equal to 30-60 years.
There are also more specialized approaches to the problem of the costs estimation.For example, in [25], an optimisation procedure for a problem of minimising the costs is also described.In this case, the authors express the total cost of a renewal, risk and an unavailability of a WDS as one function of time.But then this cost equation is minimized using the idea, that groups of pipes may be renewed together, if the benefits of grouping are balanced with the costs of shifting renewal investments in time.A special "greedy algorithm" is developed in order to find a solution for such a setup, and this numerical method is then applied for a WDS in a medium-sized Norwegian municipality.The authors informed about a high level of possible savings, if their "grouping" method is used.
Apart from planning a renewal schedule, the estimation of the costs can have impact on other types of decisions, even regarding physical characteristics of pipes.In [21], the costs related to a simple WDS are analysed for up to 50 years period.Two kinds of scenarios (with and without a replacement) are discussed, for which minimization procedures are applied using the installation and repair costs with, so called, a damage and inconvenience cost multiplier (which is related to the claims from damages and other inconveniences caused by breaks).In this case a genetic algorithm is used to solve the optimization problem.Then practical conclusions about the diameter of pipes, which should be installed, and an age, when instead of repairing, a replacement should be done, are drawn.Once again, the authors emphasizes possible savings, if an appropriate analysis of the problem is conducted.
The authors of [3] also focus on physical characteristics of a pipe.For various types of materials used for pipes, different models of break rates, based on an exponential hazard function, are incorporated.Such an approach leads to different functions of the costs, in which the costs and possible benefits of inspection technologies are also taken into account.It should be noted that in this paper, only a nominal value of the money is considered.
In some papers, time, which is necessary for a whole repair procedure, is completely neglected.Such an assumption may be accepted, even in practical situations (see, e.g., [28] for a more detailed discussion), but in many papers the required time is modelled using some random variable.For example, in [11], the generalized exponential density is applied.The authors also describe various types of the costs related to repairs, like the costs of a procedure itself, the costs of an extra energy (which is necessary to satisfy the public demands), the costs of water losses and a loss of revenues (which takes place, when the water demands are not satisfied).
A completely different scale (let's say, a "micro-management") of the problem of the costs is considered in [1].The authors, using USA-CREL condition index (CI) system, optimize a maintenance time of a wastewater plumbing system for a single building -Esteghlal Hotel in Iran.The costs of repairs are given by a virtual variable related to a condition index of a component, and the costs of replacements are based on some physical characteristics of a pipe.In order to optimize the costs, a saving to investment method is applied.It allows the authors to analyse the maintenance management costs and it leads them to formulate practical advices concerning the time when the whole water system in the building should be replaced instead of extending its maintenance or additional repairs.
Then other approaches to a function of the costs are also possible.In [32], a reliability of sand filters, as a part of a whole WDS, is considered.For the costs function, a concept of Life Cycle Costing is adopted.Apart from various sources of possible payments (during a construction phase, a usage phase and, finally, a disposal), some other factors are taken into account: an economical operational readiness rate, a target operational readiness rate and a critical operational readiness rate.They form a whole concept of a risk analysis related to the possible costs.
As it is seen from these exemplary papers, the problem of an analysis and an estimation of the costs is very important during a construction and a further maintenance of a WDS.The approaches to such a problem are various, starting from a form of the cost function (which can depend on physical characteristics of a pipe, time of a repair etc.), through a scale of a WDS (a single building or a whole system), considered parts of a WDS (pipes, filters, intakes etc.), a scope of an analysis undertaken (which can include an additional optimization procedure, a calculation of some statistical measures etc.) and so on.However, we can point out some ideas, which are very common in these papers.Firstly, usually a long time period is considered, even 50-60 years.Secondly, some numerical algorithm is applied, which generates possible times of occurrences of the faults, varied values of the costs etc. Thirdly, this algorithm is related to some random process, which is assumed for the considered model.Fourthly, because of the mentioned length of the analysed period, the value of the money in time is usually taken into account.
Such an approach (i.e.introduction of the relation between the time and the value of the cash flows) is widely studied in financial mathematics, because it is known, that one unit of money, which is paid now, has different value than the same unit, which will be paid in the future.Therefore, we also adopt this idea and analyse the estimated present value (i.e. the value of money for now) of the costs of future services.The maintenance costs are one of the key factors during selection from a set of various possible decisions, if the financial risk is taken into account.However, in various papers, the related discounting factor is only a minor, almost negligible, part of the whole simulated model.Usually, the yield is given by a factor, which is constant in time.In our considerations, the stochastic model of the interest rate is directly embedded into the Monte Carlo simulations, as important part of a reliability analysis.
Additionally, it is given by the one-factor Vasicek model.Therefore, such a yield can be better adjusted to real life data (e.g.estimated using statistical methods), especially if a long time horizon is considered.Also other, more complex models of the interest rate process can be directly used in the general approach presented in this paper.To the best knowledge of the author, an incorporation of a stochastic model of the yield is not considered in other papers concerning the analysis of the WDS costs at all.Further, using the Monte Carlo simulations, we show that there are significant differences in the estimated output, if, instead of the stochastic model, more classical approaches (like a constant yield) are taken into account.It means that, if we ignore a stochastic nature of the yield (which, of course, can not be constant all the time), some significant error is introduced into our analysis of the costs.
Additionally, in many of the papers, only one source of an uncertainty is modelled, i.e., some random process describes a behaviour of a WDS, e.g. the times of its defaults.However, there are other possible sources of uncertainties, which can be related to an incorporation of the experts' knowledge.Therefore, in our setting, the parameters of the model are described by fuzzy numbers.It means, that they are not completely precise (i.e., "crisp") but they are, in some way, uncertain ("near to / about").Such an idea is close to the real life, because, e.g., an unconditional time of a replacement of a pipe is rather "about 5 years", than strictly "always and only 5 years and not one day more / less".
In this paper, the costs of the services are related to a set of possible states of the single connection.These states reflect various types of quality and degradation of the pipeline and its current behaviour (i.e. if it operates).We assume that the random transitions between the states of the single connection form a semi-Markov process.These transitions are also influenced by a deterministic, unconditional replacement age, which can be part of a maintenance routine.
The connections in the WDS can be grouped by their types, depending on various technical properties (like quality, size etc.), which leads to different sets of statistical and reliability parameters.Then, the present value of the future possible costs of the maintenance services is estimated using Monte Carlo simulations for the whole WDS, which consists of the mentioned groups.Apart from the estimator of the present value, other statistical measures related to the repairs and the replacements, as well as their costs, can be also directly obtained using simulations.These measures can be an important input for the reliability analysis of the WDS.
This paper can be seen as an important step in generalization of the approach proposed in [26].The currently presented contribution is fourfold.Firstly, as it was previously mentioned, instead of the constant yield, the interest rate is modelled using a stochastic process.The Monte Carlo simulations of trajectories of this process are directly embedded into a procedure of generation of the costs of the services of the whole WDS.The stochastic model allows us for a better adjustment of the yield to the real life data and to minimise the possible error related to an assumption of the constant yield.This is very important for the practitioners in a proper estimation of the future maintenance costs.Secondly, instead of a classical Markov process with a discrete state space, the semi-Markov process of the transitions between the states, with the introduction of the deterministic, unconditional replacement age, is proposed.It allows us to apply other random distributions for the transitions, not only the exponential one.Thirdly, we enrich our model introducing fuzzification of some parameters of each connection in the WDS, like the costs of the services and the intensities of the transitions.This fuzzification reflects the uncertainty observed in real life data, which can be overcome with the experts' knowledge.It allows us, to some extent, to measure "an error" produced by these uncertain parameters, because a simulated output forms also a fuzzy number.This is important for a maintenance analysis, because it leads to a better estimation of the future, uncertain costs.Fourthly, apart from the Monte Carlo simulations of the WDS for the crisp case, some examples of the analysis of the WDS, if the parameters are given by fuzzy triangular numbers, are presented.It should be mentioned, that other types of L-R fuzzy numbers can be also directly applied for the considered approach.
This paper is organized as follows.In Section 2, the model of the WDS is proposed.In this part, the set of connections, the model of the maintenance costs related to repairs and replacements, as well as the one-factor Vasicek model are presented.In Section 3, a Monte Carlo approach is used to simulate a behaviour of the WDS, if all of the considered parameters are crisp.Apart from the description of the relevant algorithm, some examples of various settings are analysed.Section 4 is devoted to the problem of fuzzification of the parameters of the model.Firstly, some definitions and necessary symbols are introduced.Then, various examples of the output for the fuzzy environment, which are based also on the Monte Carlo simulations, are discussed.

Model of WDS
Let us suppose that the considered WDS is modelled by a graph of connections G.In this graph, each connection (i.e. a pipeline which is a part of the whole WDS) is represented as an edge, and possible sources or outflows are denoted by nodes.In practice, various types of nodes are possible, e.g., intake points (sources), supply points (sinks), branching nodes (see, e.g., [18] for a more detailed discussion).

States of the connection
In the following, we focus only on the edges of the graph G, i.e. the connections of the considered WDS.We assume that the state of each connection for a given time t can be described by one of the possible states from a set S={0,1,2,3,4}.These numbered states are related to the various types of quality of this pipeline (which depends on its degradation and the probability of possible malfunction) and its current behaviour (i.e., if this connection operates or not).In this way we have: State 0 - The considered section of WDS is under replacement now, because it is completely broken, so it does not work.

State 1 -
The section is under repair, because of its temporary malfunction, and it can not be used now.

State 2 -
The section works and can be described as to be in a burn-in phase, when some initial defects just after starting this part (as a new one) are possible.

State 3 -
The section operates and it is in its normal operations state (or a standard state), because the previous period of occurrence of some initial failures has ended.

State 4 -
The section works, but a wear-out period after the normal operations state has been reached, so this part of the WDS may fail with higher probability due to deteriorating stresses from its previous lifetime.Some of the considered states, namely burn-in, normal operations and wear-out phases, are described in a more detailed way in [4].The quality of a pipeline can be also modelled using other approaches, see, e.g., [7,17,25,28,29].However, the presented setting is more flexible, because it allows for an introduction of more than the five considered states.Therefore, the quality of a pipeline can be modelled even in a more complex way, if it is necessary for some practical purposes.Additionally, the parameters of each state can be statistically estimated independently of other states.
Our main aim is to analyse the state of the j-th connection at the time t, which is defined by a process of the state of the connection S t j ( ) ( ).From our previous assumption, we get S t j ( ) ( )∈{ } 0 1 2 3 4 , , , , , i.e. for any time t, the state of the j-th connection is described by some state from the set S.
Then, S j ( ) ( ) 0 is a starting state for the j-th connection, i.e., the state of the pipe during initialization of the whole algorithm (or estimation procedure, equivalently).The behaviour of the process S t j ( ) ( ) for the fixed j is described by random sequences L L j j 1 2 and S S j j 1 2 ( ) ( ) , , .Variables in the first sequence are periods between the moments when the j-th connection enters into some state and then leaves this state in order to proceed to another state.The second sequence describes the consecutive states for the j-th connection at the We assume that the process formed by S L L , 1 is a Markov renewal process.Then, straightforwardly, is a Markov renewal counting process (i.e., the number of the transitions between the states in our case) and S t S Additionally, we assume that for every i,j the processes S t i ( ) ( ) and S t j ( ) ( ) are independent of each other in probabilistic sense.From sciENcE aNd tEchNology the practical point of view, it means that there is no "information flow" between the connections and that the state of one part does not depend on the condition of other pipeline.
If for any j, L L j j 1 2 ( ) ( ) … , , are iid random variables distributed according to the exponential distribution, then the process S t j ( ) ( ) directly reduces to a Markov chain with continuous time t and the discrete finite state space S (see, e.g., [22]).However, due to next assumptions and incorporation of the fuzzy approach, such a simplification is not relevant for the case discussed in this paper.
Let R j ( ) denotes the deterministic and unconditional replacement age (i.e. the planned replacement).It means that when the length of period after the last replacement of the j-th connection surpasses R j ( ) , then this pipeline is always replaced regardless of its previous history (see also [19]).Mathematically, we calculate: and when R R the process S t j ( ) ( ) immediately changes its value to the state 0 (i.e.under replacement).
It is easily seen, that if the deterministic transitions given by R j ( ) are not taken into account, the behaviour of the process S t j ( ) ( ) is completely described by the Markov renewal process , . 1 In order to conduct the simulations, the transition probabilities: i.e. the distribution of the time t of the transition from the state k to the state l for each connection j, should be set.In the following we assume that: where X X .It means that the distribution of Q kl j ( ) is given by the minimum of random times of transitions from the state k to each other possible state.In the case of the exponential distribution, the formula (3) is consistent with a classical stochastic approach, i.e. with the previously mentioned Markov chain.However, in practical applications other kinds of distributions could be also used to model , like, e.g., some lognormal distribution.Especially, the formula (3) leads directly to the method which can be applied in numerical simulations of the transitions between the states even in the fuzzy case considered in Section 4.
If we restrict ourselves to the exponential distributions for , for the fixed connection j, then we have an intensity matrix: In practice, the values in ( ) j Λ should be estimated from real-life observations.In order to do this, we should note, that a mean sojourn time for the state k is given by: and a one-step probability of the transition between the states k and l is equal to: Then, using a set of equations based on ( 4) and ( 5) for all of the states, it is possible to set values of the intensities in the intensity matrix calibrated to our observations.For example, τ 1 j ( ) is a mean time for a repairing procedure, τ 2 j ( ) is a mean time, when a pipe is in its burn-in phase, Pr j 12 ( ) is a probability, that a pipe is repaired in such a way, that it can be treated to be in the burn-in phase afterwards etc.Of course, this procedure may require some additional experts' knowledge (e.g. about a mean time for the burn-in phase).However, this requirement directly coincidences with an introduction of a fuzzy approach described in Section 4.
If, instead of the exponential distributions, other kinds of densities are considered, the equations similar to (4) and ( 5) can be also used.We have (see, e.g., [13]) where ( ) () t is a cdf of the time spent in the state k and associated with the transition from k to l, and for a relevant kernel matrix.These integrals can be analytically intractable for more complex distributions, but it is still possible to calibrate the necessary parameters for the assumed densities even in practical applications (see, e.g., [13] for a more detailed discussion and some examples).
As it was previously mentioned, in the case of the set S considered in this paper, only some transitions between the states are possible.
Therefore, only some intensities λ kl j ( ) , related to the exponential dis- tributions are not equal to zero.So, the intensity matrix has the form where X denotes a non-zero entry.Rows and columns of the matrix (6) are numbered in the same way as the states from the set S, i.e. from 0 up to 4 These inequalities reflect the fact that the burn-in and wear-out states have greater intensities of both replacements and repairs than the normal operations phase, as stated previously.
During the Monte Carlo simulations considered in this paper, the process of the state of the connection S t j ( ) ( ) is, in addition, modified by the instantaneous transitions caused by the relevant value of R j ( ) .
In the following we assume that all of the connections could be grouped in K groups by their types which depend on various technical data for the pipes, e.g., their dimensions, a kind of material which is used etc. (see, e.g., [3] for additional details).Straightforwardly, these types are related to different intensity matrices and, e.g., values of the starting states and unconditional replacement ages.

Model of maintenance costs
Based on the model described in Section 2.1, the Monte Carlo approach can be applied to simulate a trajectory of the process S t j ( ) ( )for each connection j, and then, in a similar way, to generate a behaviour of the whole WDS.The lengths of the periods between the transitions are part of such an output and the times of entries to the states 0 or 1 (i.e. the necessity of replacement / repair of some part) can be directly calculated.In the following, these times are denoted by t 1 ,t 2 ,….
In this paper we focus only on the maintenance costs related to replacements and repairs.Of course, other types of costs (like costs of water losses, loss of revenues etc. -see, e.g., [4]) can be also incorporated into the considered numerical model.
We assume that the mentioned costs depend on the type of maintenance service (i.e., if replacement or repair of some part is necessary), the length of this service and the type of the connection.Therefore we have: where c j i t ( ) ( ) denotes a total sum of costs for the given j-th con- nection and the moment t i of starting the necessary service, c const j ( ) is some constant value independent of the length of the maintenance, i.e.
it is a fixed cost, and c var j ( ) denotes a variable cost, i.e. a value which depends on the length of this service.It is possible to apply various approaches for calculation of the variable costs using the Monte Carlo simulations.For example, these costs may be directly related to the length of period of the labour (and then strictly deterministic), or modelled by an additional random distribution which depends on the mentioned length, the type of connection and other parameters (which directly leads to a random value of the payment).

Model of interest rates
In the considered setting we assume that the value of money depends on time, i.e. we apply the concept of the present value, which is widely known in financial mathematics (see, e.g.[5,27]).Such a setup is especially useful if we are interested in the long time horizon T (e.g., 10-20 years) for which the estimated costs of the maintenance services should be found.Then, the values evaluated for different scenarios can be compared for now (i.e. when t=0).These scenarios can be directly related to various possible decisions (e.g., different values of R j ( ) ), so they lead to a selection of the best decision if the financial risk is taken into account.
In the following, we evaluate the present value of the total sum of the costs of repairs and replacements PV(c), which is given by: In order to calculate ( 8), the relevant model of the interest rate should be used to find the discounting factor PV(.) for each c There are many such models known in financial mathematics (see, e.g., [5]), but in this paper we focus on the one-factor Vasicek model described by the formula: where r t is a value of the interest rate at time t, W t is the standard Brownian motion, and a, b, σ are parameters of this model.Moreover, b characterizes a long term mean level (i.e. the trajectory of r t is directed to this value during its long run), a reflects a speed of reversion towards b, and σ is an instantaneous volatility (variability) of the trajectory introduced by the random component of the Brownian motion.
In our setting, the model of the interest rate is directly embedded into the Monte Carlo simulations, as explained in Section 3. Therefore, it is necessary to apply the relevant iterative scheme for generation of increments Δr t and the discounting factor.In the case of the model ( 9), such a scheme (see, e.g., [5] for more details) is based on evaluation of r t at the fixed moments 0 0 1 = < < … < t t t n , using the formula: where , , , … are iid samples from N(0,1).Then, the cdf of for the given value of r t i is equal to: In the same way, the factor: The conditional covariance inside the pair r fv for the fixed value of r fv ) is equal to: .

Parameters of the model
Taking into account our previous considerations, the parameters of the whole model can be divided into two groups: Parameters of the given type 1.
k=1,…,K of the connection related to its reliability and its maintenance costs: S setup: a -the speed of reversion, b -the long term mean level, σ -the instantaneous volatility.

Simulations -the crisp case
In this part we assume that all of the parameters mentioned in Section 2.4 are given by crisp values (i.e.real numbers).Because of this assumption, the relevant algorithm for simulation of the behaviour of the WDS and the interest rates is more straightforward than in the fuzzy case considered in Section 4.

Algorithm
In order to calculate the present value of the maintenance costs and other output which is useful for the reliability and maintenance analysis (see Section 3.2 for some examples), the three main phases of the algorithm should be repeated n times, where n is an overall number of simulations (see Appendix, Algorithm 1).
During the first phase, the relevant Markov renewal process is simulated using the random transition probabilities given by (2) (in the general case) or (3) (for the exponential distribution considered further on in this paper).Additionally, the condition (1) should be checked and, if R R , the state of the connection is deterministically set to 0. Then, all of the transitions to the states 0 or 1 for each trajectory are found and the relevant times and periods of the maintenance services are calculated for these events.It leads us to an evaluation of the nominal total sums of the costs for each generated time of the necessary service using the formula (7).
In the second phase, the iterative schemes (10) for r t and (11) for fv t 0, ( ) are used to generate the trajectory of the in- terest rate process (9), and, simultaneously, the discounting factor.During this Monte Carlo step, the ordered times of the transitions to the states 0 or 1 of all of the trajectories simulated in the first phase are taken into account.
In the third phase, the estimator of the total discounted costs of the maintenance services is calculated.

Examples of analysis
In each of the following examples n = 1 000 000 simulations are conducted.

Example I
In the following, for illustration purposes, i.e. to validate a correctness of the introduced algorithms and to present possible outcomes, we use artificial parameters for the considered WDS.However, some of them were obtained from an oral communication with the experts.Let us assume that only one type of the pipeline is considered and it is described by the parameters: const rpl const rpr var rpl , , , c n var rpr (12) Additionally, the intensity matrix of the transitions is given by: ( ) 1 0 0 12 0 0 0 0 24 26 0 0.5 1 0 1 0 0.4 0.9 0 0 0.3 0.6 1.1 0 0 0 and the unit of time can be identified with one year.Therefore an average time which is necessary for complete replacement of the broken connection is about one month, an average time of repair (if then the burn-in phase is achieved) is half of this time, the unconditional replacement age is equal to five years etc.Of course, in practical applications the relevant parameters of the considered connections should be estimated from real data (or based on the experts' opinions), as noted in Section 2.1.
For the interest rate model it is assumed that: a=0.1,b=0.05,r0 =0.04,σ=0.001(14) and we are interested in a rather long twenty years time horizon of the financial analysis of the maintenance costs (i.e., T=20).After using the Monte Carlo simulations, the following estimators are found: x rpr = 360 597 .
(an average number of replacements),  .(the discounted value of the maintenance costs) for the whole considered period.Some other useful estimators of times and costs can be also directly found from the same sample (see Table 1 for some examples).
The obtained estimator of the discounted value of the maintenance costs for the Vasicek model PV(c) can be directly compared with the similar values for more classical approaches -a nominal value of the cash flow (denoted further by PV nomin (c)) and the model with the constant yield r (denoted by PV const (c)).Using the respectively modified Monte Carlo approach, we get PV nomin (c)=2454.9and PV const (c)=1547.42.In this second model, r = 0.05 is set, which is equal to the long term mean level b for the Vasicek model given by (14).The relative differences between the outputs from these models are rather significant -about 50.6813%, if we compare the nominal approach with the Vasicek model, and about −5.01964%, if the constant yield is taken into account.Therefore, a selection of the appropriate interest rate model seems to be an important step in a decision making process.

Example II
The considered Monte Carlo approach could be also useful for decision makers to examine an influence of various parameters or scenarios on the estimated output like the costs of maintenance services or the reliability statistics.For example, from the previous analysis the average number of the unconditional replacements may be seen as too high, so one may be tempted to increase the unconditional replacement age and set the higher value of R (1) .Of course, such an action could have undesirable effect in the long-time horizon and, if R (1) is set too high, then also the number of "usual" (i.e.not planned) replacements or repairs caused by deterioration processes may increase.But the simulations, similar to the previously described ones, can be directly used to support the process of taking appropriate decisions.
For example, for the parameters considered in our setting, such a course is desirable, because for R (1) =10 we have x rpr = 361 715 ., which is similar to the value from the previous case (but now this average is also higher, the relative difference is about 0.310041%), x rpl = 173 847 . is significantly lower value, x rpl R , .= 1 2931 is only 7.14% of the previous average of number of the unconditional replacements, x rpl NR , .= 172 554 (which is similar to the previous case, but once again the relative difference is positive and equal to 0.286468%), and the total discounted value of the cash flow is also reduced to PV(c)=1578.87.
A similar analysis can be done also for a whole interval of the possible values of R (1) .Figure 1 presents the relation between the unconditional replacement age and PV(c), in the similar way as in Example I, for different models of the interest rate: the Vasicek model (circles), the constant rate (squares) and the nominal value (rhombuses).All of these present values are exponentially decreasing functions, but the differences between them are, once again, significant.
The similar analysis can be done also for other values, which are very important for the practitioners and the decision makers, like the average number of repairs x rpr (see Figure 2) and the average number of not planned replacements x rpl NR , (see Figure 3).In the both cases, the average number of services is a slowly increasing function of R (1)  for the given set of the parameters.

Fuzzification of the parameters
As it is known, some sources of uncertainty may be easily modelled by the fuzzy approach (see, e.g., [6]).In such a case the value of some uncertain parameter is based on expert's knowledge.This approach is especially very important when data is sparse and various data analysis methods, like statistics, are not usable or even not Fig. 1.Graph of relation between R (1) and PV(c) from Example II Fig. 2. Graph of relation between R (1) and x rpr = 361 715 .from Example II Fig. 3. Graph of relation between R (1) and x rpl NR , .= 172 554 from Example II sciENcE aNd tEchNology possible.Then, based on opinions of the experts, the necessary parameters of the model can be evaluated.However, usually these opinions have not completely precise form (e.g.like real numbers), they are rather stated as "the value of this parameter is about 5", so they can be transformed into special form, known as fuzzy numbers.The fuzzy numbers are also used during making statistical decisions (see, e.g., [9] for an example of application to some reliability data).
In the following, we assume that some of the parameters describing the considered type of the connection (enlisted in Section 2.4) are fuzzified.Such an assumption requires different Monte Carlo approach than in the crisp case, therefore new, relevant algorithm is discussed in Section 4.2.A comparison of the Monte Carlo simulations of crisp and fuzzy random numbers can be found, e.g., in [8].
For example, the planned replacement age R (j) is treated further as a fuzzy number.Theoretically, such a parameter is completely deterministic and crisp value, e.g.we may say "R (1) is equal to 5 years".However, in the practical cases, this value is not completely exact and crisp and it is rather "about 5 years plus / minus a few days / months", which is caused by a temporary lack of funds for the given period, urgent and more important repairs in other areas just in the given moment (hence, lack of the experienced staff at that time), problems with preparing plans for the traffic exactly in the given date (like busy working day) etc.Therefore, the introduced fuzziness has practical basis.
However, some quantitative information about the planned replacement age or other parameter is still required.Let us assume, that the experts express their opinions.These sentences can have, e.g., a form of linguistic variables (like "a time of a planned replacement is rather later than sooner").But in the setting considered in this paper, it is necessary to "translate" them into relevant sets of real numbers with membership functions, i.e. fuzzy numbers (see, e.g., [30]).Then, these numbers, which have to satisfy some additional requirements about a monotonicity of their membership functions (like L-R numbers), are used during simulations in order to evaluate a fuzzy output, e.g., fuzzy maintenance costs.Applying such a procedure, we gain an additional knowledge about a dependency between the uncertain (fuzzy) parameters (i.e. the experts' opinions) and a simulated (uncertain) outcome (i.e.conclusions, which are important for our management decisions).It means, that a level of the uncertainty in the input is directly reflected in the uncertainty of the output and can be also easily measured.For example, a difference between the fuzzy maintenance costs for two scenarios -"the planned replacement age is equal to 5 years plus / minus a half of the year" and "the planned replacement age is equal to 4.5 years plus / minus a half of the year" -can be exactly seen and measured, and a relevant decision can be taken.
It should be also noted, that an incorporation of the fuzzy variables is related to a different type of uncertainty, than one modelled by probabilistic approaches.Then, we enrich our model in a significant way.

Preliminaries
Now we present basic definitions and notation concerning the fuzzy approach, which will be used in the further part of the paper.Additional details can be found in, e.g., [14].
For a fuzzy subset  A of the set of real numbers R we denote by µ ,⋅,/ is an operator for the fuzzy numbers (related to the equivalent operator +,−,⋅,/ in the crisp case), then for the fuzzy numbers a ,b , their outcome is also a fuzzy number.Using α-cuts and the interval arithmetic we have: ( ) = 0 for x≤0.
A triangular fuzzy number, denoted further by [a,b,c], is a L-R number with the membership function of the form: In the following, to approximate a desired fuzzy output f , we check monotonicity of the underlying function f(x), when all of the parameters of the whole model, except of the single argument x, are held fixed.
The same idea, related to the extension principle, applies also for λ >0 given by a L-R fuzzy number, where λ ̃ is the fuzzy intensity parameter of the exponential probability density function.Therefore, if f(λ) is an increasing function, then for the given α, the left end point f ̃L [α] is approximated using the crisp value λ ̃L [α] as the intensity of the exponential random variables generated in Monte Carlo simulations.In the same way λ ̃U [α] is used to approximate f ̃U [α] (see also, e.g., [2,23,24]).

Simulations in the fuzzified environment
We assume that the parameters of each type of the connection, = ( ) ∈ λ , i.e. such a matrix consists of five rows and five columns of positive fuzzy numbers and the pattern of strictly positive entries is the same as for (6).
Using the Monte Carlo simulations, the output similar to the one discussed in Section 3.2 is then obtained.But now these values, like the discounted value of the maintenance costs, are given by fuzzy numbers, so notation like PV c  ( ) will be used further.The fuzzy output is approximated by α-level sets PV c  ( )[ ] α eval- uated using Algorithm 2 (see Appendix).Firstly, we should decide if the left or right end points of PV c  ( ) should be estimated during the simulations.Then a binary variable left=1 or left=0 is set, respectively.The other parameters are: a starting value α ∈ ( , ] 0 1 , an upper bound α α 1 0 1 ∈ ( , ] and an increment Δα>0.In the first phase, for the given α, the relevant α-level cuts of the fuzzy parameters for all types of connections are found.Based on the monotonicity of PV(c), for each possible fuzzified argument and the fixed value of left, the left or right end point of the α-level cut of each parameter is then selected.
During the second step, the simulations of the whole WDS are conducted using Algorithm 1 (see Appendix).The generated crisp output To obtain the mentioned approximation of the whole fuzzy output, the value of α in the above procedure is gradually set to subsequent values, starting from α 0 up to α 1 with the given increment Δα.After evaluation of the left end points of PV(c) (for left=1), the right end points are found in the same manner (i.e.left=0 is applied).Then the obtained intervals are putting on one another to form the approximated fuzzy number.
The above method is based on the extension principle introduced by Zadeh (see [31]) and a similar approach using α-level cuts is also applied in financial mathematics for pricing of the derivatives (see, e.g., [23,24]), in optimization of queuing systems (see, e.g., [2]) etc.
Because the considered function for the fuzzy result (e.g.PV c  ( ) is not given by an exact, analytical formulae, the simulations are required for its evaluation.As previously mentioned, it should be determined for each of the fuzzy parameter, if the left or the right end point of its α-level cut is used to calculate the left or right end point of the output interval during the first step of the above procedure.This choice is strictly related to the monotonicity of a function f(p) of the considered outcome, i.e. if the output (like the discounted maintenance costs) is an increasing or decreasing function of the given parameter p (like the constant costs of repair).In Algorithm 2 (see Appendix) the considered approach is applied for PV c  ( ) , but the same method can be used for other output functions f(.), which are important from the practical point of view.

Example III
For simplicity, in the following, we assume that the fuzzy parameters of the considered single type of the connection are given by some triangular numbers.However, other types of L-R numbers can be also directly used in our approach.Let us start from the set of the parameters: , , , ,, , ,  c var rpr (15) Other parameters are the same as in Section 3.2.1.It is easily seen from (15), that only the constant and variable costs of repairs and replacements are strictly fuzzy triangular numbers in the considered case.Therefore, they are the only source of uncertainty in the following example.
Using the Monte Carlo approach described by Algorithm 2, the fuzzy discounted maintenance costs PV c  ( ) for the Vasicek model (see Figure 4, circles) and the fuzzy average of the nominal costs of replacements  c rpl (see Figure 5) are approximated.The obtained fuzzy values are almost triangular symmetrical fuzzy numbers for which the 1-level cut sets are equal to the output estimated in Section 3.2.2.Some approximated values of the intervals for different α can be found in Table 2.
Both PV c  ( ) and  c rpl are increasing functions of the fuzzy parameters (15).Therefore, to calculate the left end point of the interval PV c L ( )[ ] α (or the right end point PV c U ( )[ ] α , respectively) for the given value of α, only the left end points (or the right points) of the α-level cuts for ( 15) should be considered.The same applies for  c rpl .Using the Monte Carlo simulations, the present values of the maintenance costs can be also found for the more classical models of the interest rates, i.e. the nominal costs (see Figure 4, squares) and the constant yield (see Figure 4, rhombus).The obtained approximations of the fuzzy numbers are similar in shape, but significantly different.

Example IV
Of course, not only the maintenance costs can be considered as uncertain values in the practical situations.Therefore we analyse situation if the unconditional replacement age is also given by a triangular, fuzzy number, denoted further by ( 1 )  R  .Other fuzzy parameters ) ) are described by (15) and the crisp parameters are the same as in Section 3.2.1.As previously, to approximate a desired fuzzy outcome, we should check the monotonicity of the considered function for the introduced fuzzy parameters.As noted in Section 4.2.1, PV(c) is an increasing  The fuzzy values, given by ( 17), describe only the fuzzy intensities of the exponential distributions of the transitions between the states from the set S, as discussed in Section 4.1.It is easily seen, that the fuzzy numbers (17) are "close" to the crisp values (13) assumed in Example I. Therefore, the fuzzy output obtained during current analysis can be easily compared with the values from the previously considered examples.
As previously mentioned, in order to evaluate the fuzzy output, it is necessary to check the monotonicity of the function for the considered parameter.And PV(c) is a decreasing function of ( ) λ , ( ) λ , ( ) λ , because for the higher values of these intensities, the expected time of replacement or repair is lower, so the final cost is also lower.All of the other parameters can be examined in the same way.The approximated fuzzy discounted maintenance costs, which are evaluated using the Monte Carlo simulations are also compared with the similar output obtained in Example III (see Figure 9).As it is seen, in the considered case PV c  ( ) is a L-R fuzzy number, almost a triangular one, with a wider support comparing to the fuzzy number from Example III.
In a similar way, the fuzzy average (nominal) costs of replacements is approximated.In Figure 10, the obtained values are compared with the similar fuzzy output from Example III.Some of the evaluated intervals for different α can be found in Table 3.
As previously, the fuzzy approximations of the average number of repairs  x rpr (see Figure 11) and the average number of not planned replacements , rpl NR x  (see Figure 12) are also obtained using the introduced approach.These fuzzy L-R numbers are almost triangular, but their supports are wider comparing to the outputs from the examples, where the intensities are strictly crisp values.

Conclusions
In this paper the model of the transitions between the states of the single connection of the WDS is proposed.This model is related to a semi-Markov process and the deterministic jumps caused by the unconditional replacement age for each pipeline.Using the Monte Carlo approach, the behaviour of the whole WDS is simulated and the costs of the maintenance services (i.e.repairs and replacements) are evaluated.Then, based on a stochastic process (the one-factor Vasicek model of interest rates), the present value of the future cash flows and other important statistical measures of the mentioned costs are calculated.The estimated differences between the classical approach (i.e., a constant yield) and the proposed model are emphasized.Apart from the strictly crisp setup, the fuzzification of some parameters of the connection is considered.An introduction of the fuzzy numbers leads to a better incorporation of the experts' knowledge and more real-life modelling of the uncertain parameters.The necessary numerical algorithms and some relevant examples of the simulated output for both the crisp and the fuzzy settings are also provided.
There are various possible fields for future works based on the approach introduced in this paper.Firstly, other kinds of distributions, instead of the exponential ones, can be assumed -e.g., the repairing time can be described by some lognormal random variable.Secondly, if a different cdf is applied, a relevant case can be calibrated using some real-life data, so statistical tests concerning a validity of this cdf (versus, e.g., the exponential one) can be proposed.It also leads to a possibility of measuring a level of a difference in the estimated maintenance costs, which can be important for the practitioners.
independent random variables and each X kl j ( ) is given by the exponential distribution with the intensity λ kl j repair), n k -the number of the connections of this type.Parameters of the interest rate model related to the financial 2.
(an average number of the unconditional, planned replacements), x rpl NR , .= 172 061 (an average number of the replacements without the planned ones, i.e. the not planned replacements) and PV c ( ) = 1629 2  a is a fuzzy subset of R for which µ  A is a normal, upper-semicontinuous, fuzzy convex function with a compact support.Then for each α∈[0,1], the α-level set  a[α] is a closed interval of the form  a k=1,…,K, are given by L-R numbers.Therefore in the following they are denoted by  R k ( ) ,  Λ ( ) k etc.Then for the intensity matrix of the transitions we have   Λ

Fig. 4 .
Fig. 4. Fuzzy discounted and nominal maintenance costs from Example III

Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .
Fig. 9. Fuzzy discounted maintenance costs from Example V (circles) and Example III (squares) sciENcE aNd tEchNology ment of the broken part by a new one, which begins its work in the burn-in phase afterwards, and λ 30 j ( ) describes the intensity of break- ing down with a necessary replacement as its consequence, if this connection was in its normal operations state previously.Additionally, from the description of the states, we assume that: . For example, λ 02 j ( ) corresponds to the intensity of replace-

Table 1 .
Exemplary estimators from Example I does not contain zero for all α∈[0,1] in the last case.A fuzzy number a ĩs called positive (a ≥0) if µ

Table 2 .
Exemplary α-level cuts of the fuzzy discounted maintenance costs

Table 3 .
Exemplary α-level cuts of the fuzzy discounted maintenance costs PV c  ( ) and the fuzzy average of the nominal costs of replacements  c rpl from Example V