Age-usage semi-Markov models

Applied Mathematical Modelling 35 (2011) 4354–4366

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Age-usage semi-Markov models Guglielmo D’Amico ⇑ Dpt of Drug Sciences, University ‘‘G. D’Annunzio’’, 66013 Chieti, Italy

a r t i c l e

i n f o

Article history: Received 15 January 2010 Received in revised form 22 February 2011 Accepted 2 March 2011 Available online 10 March 2011 Keywords: Non-homogeneity Reward process System reliability Algorithm

a b s t r a c t This paper presents a non-homogeneous age-usage semi-Markov model with a measurable state space. Several probability functions useful to assess the system’s reliability are investigated. They satisfy the same family of equations we call indexed Markov renewal equations. Sufficient conditions to assure the existence and uniqueness of their solutions are provided. The numerical analysis of these equations is executed through the construction of a process discrete in time and space, which is shown to converge to the continuous one in the Skorohod topology. An algorithm useful for solving the discretized system of equations is presented by using a matrix representation.  2011 Elsevier Inc. All rights reserved.

1. Introduction Age-usage reliability models are drawing even more interest in the study of the performance of systems. The reason is that, frequently, to place the usage measure beside the age measure provides a more accurate description of the system’s reliability. For example, the reliability of an automobile depends on age and on usage (miles driven). Often, the higher the accumulated usage, the higher is the failure probability of a device. There are three main approaches to analyzing age-usage models. The first approach consists in the introduction of a correlation between the two measures; see for example [1] and the references therein. The second approach uses a bivariate renewal theory. This was proposed for the first time by Baggs and Nagagaja [2] and recently it has received increasing attention by researchers (see [3–6]). The third approach mixes the two scale measures into a single composite one (see [7,8]). In general, this is done by considering a functional dependence between the usage and the age measures. Many reliability studies are based on Markov processes. The Markov property has been frequently criticized as well as the use of a discrete even countable state space [9,10]. To overcome these drawbacks, semi-Markov processes have been frequently applied in reliability studies [11–16]. In spite of this, few works deal with semi-Markov processes with arbitrary state space. In [10,17], definitions and general results are established for the time homogeneous case and applications in the reliability modeling are discussed. In [18], the entropy for a time homogeneous semi-Markov process with Borel state space is investigated and asymptotic equirepartition properties and invariance principles are assessed. In [19], recurrence relations for generalized hitting times for semi-Markov processes are studied. To the best of the author’s knowledge, semi-Markov processes were never applied to construct age-usage models. In this paper we propose an age-usage model that is constructed specifying a functional relationship between the two measures. The evolution of the system is described by non-homogeneous semi-Markov processes with Borel state space. The kernel of the process is indexed by a stochastic process. The index represents the accumulated usage of the system. In this way, ⇑ Corresponding author. Tel.: +39 08713554609; fax: +39 08713554622. E-mail address: [email protected] 0307-904X/$ - see front matter  2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.03.006

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

4355

systems with different usages exhibit different future dynamics. The usage process is assumed to be a function of the states of the semi-Markov process, of the time and of the usage accumulated. The investigation of the age-usage model leads to the study of Markov renewal equations that are indexed by the usage stochastic process. We refer to these types of equations as an indexed Markov renewal equation (IMRE). A sufficient condition for the existence and uniqueness of the solution is provided. The solution of the IMRE is made possible after a discretization of the state space and of the time. We discuss the convergence of the discrete approximating process to the continuous one in the sense of the Skorohod topology. This provides a theoretical validation of the adopted discretization scheme. The solutions of the discretized equations are obtained by using a matrix representation of the equations that makes use of a matrix operator that we called the backslash matrix operator. The results established here generalize, at the same time, those obtained in a Markovian environment with finite state space [20], those concerning ordinary Markov renewal equations for an homogeneous semi-Markov process with an arbitrary state space [10,17] and those on the numerical solution of Markov renewal equations obtained in the finite state space case [21,22,14]. These extensions are of practical relevance as they permit the construction of a model that considers simultaneously relevant aspects of mathematical modeling such as the non-homogeneity, the usage accumulation, the general state space. The plan of the paper is as follows. In Section 2 we define the concepts related to the proposed age-usage semi-Markov model. Transition probabilities relevant to the system reliability analysis are expressed via indexed Markov renewal equations. Conditions assuring the existence and the uniqueness of their solutions are given. In Section 3, we present the numerical analysis of these equations by using quadrature formulas. The validity of this numerical scheme is established by proving the convergence of the discrete approximating process to the corresponding continuous one in the Skorohod topology. In Section 4 we present an algorithm to compute the solutions of the discrete equations. Finally, in Section 5 we present our conclusion. 2. Age-usage semi-Markov models In this section, following the line of research in [10,17], we define age-usage semi-Markov models and provide some results about them. Let consider a measurable space ðE; EÞ such that for all x 2 E it results that fxg 2 E. Definition 1 (Non-homogeneous semi-Markov kernel). A function Q ðx; s; A; tÞ; x 2 E; s 2 Rþ ; A 2 E; t 2 Rþ ; t P s, is called a nonhomogeneous semi-Markov kernel if: (i) Q ðx; s; A; tÞ 8x 2 E; s 2 Rþ ; A 2 E is a non-decreasing right continuous real function with Q(x, s; A, s) = 0; (ii) Q(, s; , t) 6 1 "s, t such that s 6 t; (iii) p(, s; ) = Q(, s; , 1) is a Markov transition probability function from ðE; EÞ to itself. For each ðx; sÞ 2 E  Rþ there exists a probability space ðX; F ; Pðx;sÞ Þ and two sequences of random variables (Jn, Tn) such that

Pðx;sÞ ½J 0 ¼ x; T 0 ¼ s ¼ 1; Pðx;sÞ ½J nþ1 2 A; T nþ1 6 tjF n  ¼ Pðx;sÞ ½J nþ1 2 A; T nþ1 6 tjJ n ; T n  ¼ Qðx; s; A; tÞ: In this way (Jn, Tn) is a Markov process with state space E  Rþ and transition probability function given by the semi-Markov kernel. We call this process a non-homogeneous Markov Renewal Process. The sequence (Jn, n P 0) gives the successive states of E in time and the sequence (Tn, n P 0) gives the times at which transitions occur. Set S(0) = 0 and define

SðtÞ ¼ supfn 2 N : T n 6 tg;

t > 0;

where S(t) is the number of transitions that occur by time t. Definition 2 (Non-homogeneous semi-Markov process (NHSMP)). The process {X(t); t P 0}, defined by

XðtÞ ¼ J SðtÞ is a non-homogeneous semi-Markov process (NHSMP) with arbitrary state space E and kernel Q(x, s, A, t). The semi-Markov process defined here is time inhomogeneous because its kernel is a function of times s and t and not only of the difference t  s. The non-homogeneity property gives the possibility to face the age problem of a system by assigning different transition probabilities depending on the starting time (age) s. Systems with equal age may experience different performance due to diverse accumulated usage. For this reason it is preferable to consider also the usage measure.

4356

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

Definition 3 (The usage process). Let U n : X ! Rþ denote the total usage accumulated by the system up to the nth transition:

U n ¼ U n1 þ

Z

Tn

f ðJ n1 ; s; U n1 Þds;

ð1Þ

T n1

where f : E  Rþ  Rþ ! Rþ is a Borel measurable bounded function and U 0 2 Rþ is known and non-random. The process Un can be interpreted as a reward process with the function f as a measure of the usage per unit time. The function f depends on the state of the system Jn1, on the time s and on the accumulated usage Un1. To construct an age-usage model we should specify a dependence structure between the age variable and the usage variable. Toward this end we adopt the following assumption: Assumption 1. The three-dimensional stochastic process {Jn, Tn, Un} meets the following Markov renewal property:

P½J nþ1 2 A; T nþ1 6 t; U nþ1 6 MjrðJ h ; T h ; U h Þ; h 6 n; J n ¼ x; T n ¼ s; U n ¼ m ¼ P½J nþ1 2 A; T nþ1 6 t; U nþ1 6 MjJn ¼ x; T n ¼ s; U n ¼ m;

ð2Þ

where r(Jh, Th, Uh), h 6 n is the natural filtration of the three-dimensional process. Assumption 1 asserts that, the knowledge of the values Jn, Tn, Un suffices to give the conditional distribution of the triplet Jn+1, Tn+1, Un+1 whatever the values of the other past variables might be. It should be noted that Un+1 is measurable with respect to (Tn+1, Tn, Un, Jn), consequently, all information about the probabilistic evolution of the future of the process are in the function:

Qðx; s; m; A; tÞ :¼ P½J nþ1 2 A; T nþ1 6 tjJ n ¼ x; T n ¼ s; U n ¼ m The function Q(x, s, m; A, t) has a fundamental role in the theory we are going to expose. In recognition of its importance, we call it indexed semi  Markov kernel. The formal definition of the indexed semi-Markov kernel is given now. Definition 4 (Indexed non-homogeneous semi-Markov kernel (ISMK)). A function Q ðx; s; m; A; tÞ; x 2 E; s 2 Rþ ; m 2 Rþ ; A 2 A; t 2 Rþ ; t P s, is called a ISMK if: (i) Q ðx; s; m; A; tÞ 8x 2 E; s 2 Rþ ; m 2 Rþ ; A 2 E is a non-decreasing right continuous real function with Q(x, s, m; A, s) = 0; (ii) Q(, s, m; , t) 6 1 "s, t such that s 6 t and m 2 Rþ ; (iii) p(, s, m; ) = Q(, s, m; , 1) is a Markov transition probability function from ðE; EÞ to itself.

Remark 1. The Markov Renewal Process (Jn, Tn), which is embedded in the ISMK, depends on the usage process Un, the latter acts as a stochastic index. Moreover, the usage process Un depends on (Jn, Tn) through the functional relationship (1). To describe the behavior of the age-usage model observed at some fixed time t we also need additional stochastic processes. Definition 5 (Indexed semi-Markov (ISMP) and usage processes). Given the three-dimensional process {Jn, Tn, Un} and the ISMK Q(x, s, m; A, t), denote by

NðtÞ ¼ supfn 2 N : T n 6 tg; ZðtÞ ¼ J NðtÞ ; UðtÞ ¼ U NðtÞ þ

Z

t

T NðtÞ

f ðJ NðtÞ ; s; U NðtÞ Þds:

They represent the number of transitions up to time t, the state of the system at time t and the usage accumulated up to t, respectively. We refer to Z(t) as an indexed semi-Markov process (ISMP). Remark 2. Special cases of ISMP with finite state space have found some applications. For example, in [23] Janssen and De Dominicis propose a non-homogeneous semi-Markov model indexed by the process counting the number of transitions. This cases is obtained simply by defining Un = n. A similar model has been proposed in [24] to have the possibility of considering the seniority effects in the modeling of stochastic pension funds. This model is recovered by setting Un = Un1 + Tn  Tn1 for a given U0. : Let denote by p(x, s, m; A) ¼ P[Jn+1 2 AjJn = x, Tn = s, Un = m]. It is simple to show that

pðx; s; m; AÞ ¼ lim Q ðx; s; m; A; tÞ: t!1

The survival probability in state x 2 E is

: Hðx; s; m; tÞ ¼ P½T nþ1 6 tjJ n ¼ x; T n ¼ s; U n ¼ m ¼ Q ðx; s; m; E; tÞ:

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

4357

Because 8x 2 E; s < t 2 Rþ ; m 2 Rþ ; A 2 E; Q ðx; s; m; A; tÞ 6 pðx; s; m; AÞ the measure Q(x, s, m; , t) is absolutely continuous with respect to p(x, s, m; ), then it exists a E-measurable function F(x, s, m, A; ) such that

Q ðx; s; m; A; tÞ ¼

Z

Fðx; s; m; y; tÞpðx; s; m; dyÞ:

A

F can be chosen as a distribution function expressing the probability:

Fðx; s; m; y; tÞ ¼ P½T nþ1 6 tjJ n ¼ x; J nþ1 ¼ y; T n ¼ s; U n ¼ m: Definition 6 (ISMK n-fold convolution). Given the ISMK Q(x, s, m; A, t), define its n-fold convolution as follows:

Q

ðnÞ

(R R t E

ðx; s; m; A; tÞ ¼

  Rs Q ðx; s; m; ds; dyÞQ ðn1Þ y; s; m þ s f ðx; n; mÞdn; A; t if n P 1; t > s;

s

0

if s P t:

ð3Þ

where Q(0)(x, s, m; A, t) = v{x 2 A}v{t > s} and v(L) is the indicator function of event L. It is simple to show that

Q ðnÞ ðx; s; m; A; tÞ ¼ P½J nþh 2 A; T nþh 6 tjJ h ¼ x; T h ¼ s; U h ¼ m;

8h 2 N:

The function

Rðx; s; m; A; tÞ :¼

1 X

Q ðnÞ ðx; s; m; A; tÞ ¼ E½NðtÞ  NðsÞjJNðsÞ ¼ x; T NðsÞ ¼ s; U NðsÞ ¼ m

n¼0

is called the indexed Markov renewal function. Definition 7. The ISMK is said to be normal if R(x, s, m; A, t) < 1 for all x, s, m, t. To properly assess the probabilistic behavior of the system, we introduce some functions. Definition 8. For all x 2 E; s < t 2 Rþ ; m 2 Rþ ; A 2 E; B 2 E define the probabilities:

/ðx; s; m; A; t; MÞ ¼ P½ZðtÞ 2 A; UðtÞ 6 MjZðsÞ ¼ x; T NðsÞ ¼ s; UðsÞ ¼ m; /ðnÞ ðx; s; m; A; t; MÞ ¼ P½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ njZðsÞ ¼ x; T NðsÞ ¼ s; UðsÞ ¼ m;

cðx; s; m; A; t; MÞ ¼ P½ZðuÞ 2 A; 8u 2 ½s; t; UðtÞ 6 MjZðsÞ ¼ x 2 A; T NðsÞ ¼ s; UðsÞ ¼ m; hB ðx; s; m; A; t; MÞ ¼ P½ZðtÞ 2 A; UðtÞ 6 M; ZðuÞ R B; 8u 2 ½s; tjZðsÞ ¼ x R B; T NðsÞ ¼ s; UðsÞ ¼ m Proposition 1. The probabilities /, /(n), c and h verify the following equations:

 f ðx; n; mÞdn 6 M  m ð1  Hðx; s; m; tÞÞ s   Z s Z Z t Q ðx; s; m; dy; dsÞ/ y; s; m þ f ðx; n; mÞdn; A; t; M ; þ

/ðx; s; m; A; t; MÞ ¼ vðx 2 AÞv

Z

t

s

E

ð4Þ

s

Z t  f ðx; n; mÞdn 6 M  m ð1  Hðx; s; m; tÞÞ /ðnÞ ðx; s; m; A; t; MÞ ¼ vðn ¼ 0Þvðx 2 AÞv s   Z s Z Z t Q ðx; s; m; dy; dsÞ/ðn1Þ y; s; m þ f ðx; n; mÞdn; A; t; M ; þ s

E

ð5Þ

s

with /(1)(x, s, m; A, t, M) = 0.

 f ðx; n; mÞdn 6 M  m ð1  Hðx; s; m; tÞÞ s   Z s Z Z t Q ðx; s; m; dy; dsÞc y; s; m þ f ðx; n; mÞdn; A; t; M ; þ

cðx; s; m; A; t; MÞ ¼ vðx 2 AÞv A

Z

s

 f ðx; n; mÞdn 6 M  m ð1  Hðx; s; m; tÞÞ s   Z t Z s Q ðx; s; m; dy; dsÞhB y; s; m þ f ðx; n; mÞdn; A; t; M :

Z EnB

ð6Þ

s

hB ðx; s; m; A; t; MÞ ¼ vðx 2 AÞv þ

t

s

Z

t

s

ð7Þ

4358

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

Proof. By reason of similarity we prove only Eq. (5) for /(n). We have

/ðnÞ ðx; s; m; A; t; MÞ ¼ Pðx;z;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n ¼ Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n; T NðsÞþ1 > t þ Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n; T NðsÞþ1 6 t:

ð8Þ

Observe that

Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n; T NðsÞþ1 > t ¼ Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ njT NðsÞþ1 Z t  > tP ðx;z;mÞ ½T NðsÞþ1 > t ¼ vðn ¼ 0Þvðx 2 AÞv f ðx; n; mÞdn 6 M  m ð1  Hðx; s; m; tÞÞ: s

On the other hand

Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n; T NðsÞþ1 6 t ¼ Eðx;z;mÞ ½vðT NðsÞþ1 6 t; J NðsÞþ1 2 EÞP½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ njFT NðsÞþ1  where FT NðsÞþ1 ¼ rfðJ n ; T n ; U n Þ; n 6 NðsÞ þ 1g. The use of the strong Markov property gives

Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n; T NðsÞþ1 > t ¼ Eðx;s;mÞ ½vðT NðsÞþ1 6 tÞPðJNðsÞþ1 ;T NðsÞþ1 ;UNðsÞþ1 Þ ½ZðtÞ 2 A; UðtÞ Z Z t Q ðx; s; m; dy; dsÞPðy;s;UðsÞÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n  1 6 M; NðtÞ  NðT NðsÞþ1 Þ ¼ n  1 ¼ E

s

Because JN(s) = x, TN(s) = s, UN(s) = m and TN(s)+1 = s

UðsÞ ¼ m þ

Z s

f ðx; n; mÞdn:

s

By substitution we obtain

Pðx;s;mÞ ½ZðtÞ 2 A; UðtÞ 6 M; NðtÞ  NðsÞ ¼ n; T NðsÞþ1 6 t   Z Z t Z s ¼ Q ðx; s; m; dy; dsÞ/ðn1Þ y; s; m þ f ðx; n; mÞdn; A; t; M ; E

s

s

which completes the proof.

h

Eqs. (4)–(7) can be seen as particular cases of the general equation

wðx; s; m; A; t; MÞ ¼ gðx; s; m; A; t; MÞ þ

Z Z E

t

  Z s Q ðx; s; m; dy; dsÞw y; s; m þ f ðx; n; mÞdn; A; t; M :

s

ð9Þ

s

Eq. (9), with the unknown function w, is called indexed Markov renewal equation. The next theorem establishes conditions that ensure existence and uniqueness of solution of IMRE (9). We extend some results by Çinlar [25] and Limnios and Opriçan [10] by considering simultaneously the non-homogeneity, an arbitrary state space and the influence of the usage process through an indexed semi-Markov kernel. Our proof is based on an opportune mixing of the argument by Iosifescu-Manu [26] and Limnios and Opriçan [10]. Let M be the Banach space of bounded real function defined on E equipped with norm kfk = supx2Ejf(x)j. If wð; s; m; A; t; MÞ 2 M, then kwk(s, m; A, t, M) = kw(, s, m; A, t, M)k. Let Lþ be the set of all non-negative functions w(, , ; A, , ) such that "A (i) "s P t or "m P M w(x, s, m; A, t, M) = 0; (ii) kwk(s, m; A, t, M) is bounded on compact sets [s, t]  [m, M]; (iii) "x 2 E w(x, , ; A, , ) is Borel-measurable. Theorem 1. Let assume the following (i) f : E  Rþ  Rþ ! Rþ is a Borel measurable bounded function; (ii) g 2 Lþ ; (iii) Q is a normal indexed non-homogeneous semi-Markov kernel; then the indexed Markov Renewal Eq. (9) admits a unique solution w 2 Lþ .

4359

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

Proof. For all A 2 E let us define the sequence of functions fwðnÞ ðx; s; m; A; t; MÞgn2N  as follows:

wð0Þ ðx; s; m; A; t; MÞ ¼ gðx; s; m; A; t; MÞ; and for all n P 1

wðnÞ ðx; s; m; A; t; MÞ ¼

Z Z E

Let wðnÞ ðx; s; m; A; t; MÞ ¼ and w(h) P 0 for all h.

t

  Z s Qðx; s; m; dy; dsÞwðn1Þ y; s; m þ f ðx; n; mÞdn; A; t; M : s

s

Pn

ðhÞ þ ðnÞ h¼0 w ðx; s; m; A; t; MÞ, then the sequence fw ðx; s; m; A; t; MÞg is non decreasing because g 2 L

Moreover it results that

wðnÞ ðx; s; m; t; A; MÞ ¼

n X

wðhÞ ðx; s; m; A; t; MÞ ¼ wð0Þ ðx; s; m; A; t; MÞ þ

h¼0

¼ gðx; s; m; A; t; MÞ þ

wðhÞ ðx; s; m; A; t; MÞ

h¼1

Z Z

t

Q ðx; s; m; dy; dsÞ

s

E

n X

n X

  Z s wðh1Þ y; s; m þ f ðx; n; mÞdn; t; A; M s

h¼1

Then

wðnÞ ðx; s; m; t; A; MÞ ¼ gðx; s; m; A; t; MÞ þ

Z Z

t

  Z s Q ðx; s; m; dy; dsÞwðh1Þ y; s; m þ f ðx; n; mÞdn; t; A; M :

s

E

ð10Þ

s

Assume that wðn1Þ ðs; x; m; A; t; MÞ 2 Lþ , fix A 2 E and set

rðs; m; A; t; MÞ ¼ sup sup kwðn1Þ kðs; m; A; u; hÞ: s6u6t m6h6M

Then, by induction on n we have that

Z Z  jðQ  wðn1Þ Þðx; s; m; A; u; hÞj ¼ 

  Z s  Q ðx; s; m; dy; dsÞwðn1Þ y; s; m þ f ðx; n; mÞdn; A; u; M  s E s    Z s Z Z u   Q ðx; s; m; dy; dsÞwðn1Þ y; s; m þ f ðx; n; mÞdn; A; u; M  6 E s s   Z Z u Z s 6 Q ðx; s; m; dy; dsÞr s; m þ f ðx; n; mÞdn; A; u; h E

s

Z Z 6 E

u

s u

Q ðx; s; m; dy; dsÞrðs; m; A; u; hÞ ¼ Hðx; s; m; tÞrðs; m; A; u; hÞ 6 rðs; m; A; u; hÞ:

s

Hence kQ  wðn1Þ kðs; m; A; t; MÞ 6 rðs; m; A; t; MÞ implies that Q  wðn1Þ 2 Lþ that is wðnÞ 2 Lþ . This means that the sequence wðnÞ is non-decreasing and bounded, then it admits a limit. Suppose to denote by

e s; m; A; t; MÞ ¼ lim wðnÞ ðx; s; m; A; t; MÞ: wðx; n!1

If we take the limit in Eq. (10) we obtain

lim wðnÞ ðx; s; m; t; A; MÞ ¼ lim gðx; s; m; A; t; MÞ þ lim

n!1

n!1

n!1

Z Z E

t

  Z s Qðx; s; m; dy; dsÞwðh1Þ y; s; m þ f ðx; n; mÞdn; t; A; M

s

s

and by Beppo Levi convergence monotone theorem we have

e s; m; t; A; MÞ ¼ gðx; s; m; A; t; MÞ þ wðx;

Z Z E

t

  Z s e y; s; m þ Q ðx; s; m; dy; dsÞ w f ðx; n; mÞdn; t; A; M :

s

ð11Þ

s

e satisfies the IMRE and consequently represents a solution of (10). To prove uniqueEq. (11) indicates that the function w ness, define the function

: e Dðx; s; m; A; t; MÞ ¼ wðx; s; m; A; t; MÞ þ Kðx; s; m; A; t; MÞ

4360

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

and suppose D(x, s, m; A, t, M) is a solution of the IMRE. Hence

  Z s Q ðx; s; m; dy; dsÞD y; s; m þ f ðx; n; mÞdn; t; A; M s E s   Z s Z Z t e Q ðx; s; m; dy; dsÞ w y; s; m þ f ðx; n; mÞdn; t; A; M ¼ gðx; s; m; A; t; MÞ þ E s s   Z s Z Z t Qðx; s; m; dy; dsÞK y; s; m þ f ðx; n; mÞdn; t; A; M : þ

Dðx; s; m; A; t; MÞ ¼ gðx; s; m; A; t; MÞ þ

E

Z Z

t

s

s

Because (11) holds true, we have

Kðx; s; m; A; t; MÞ ¼

Z Z E

t

  Z s Q ðx; s; m; dy; dsÞK y; s; m þ f ðx; n; mÞdn; t; A; M :

s

ð12Þ

s

Formula (12) can be represented in a compact form by using matrix convolution notation

Kðx; s; m; A; t; MÞ ¼ ðQ  KÞðx; s; m; A; t; MÞ: By induction we get

Kðx; s; m; A; t; MÞ ¼ ðQ ðnÞ  KÞðx; s; m; A; t; MÞ: The hypothesis of normality of the IMRK entails that the corresponding IMRF is finite, i.e.

Rðx; s; m; A; t; MÞ ¼

1 X

Q ðnÞ ðx; s; m; A; t; MÞ < 1:

n¼0

As a consequence limn?1Q(n)(x, s, m; A, t, M) = 0 and limn?1(Q(n)  K)(x, s, m; A, t, M) = 0; this implies K(x, s, m; A, t, M) = 0.

h

Remark 3. Existence and uniqueness of solution of Markov renewal equations were investigated by Çinlar [25] for homogeneous semi-Markov processes with finite state space, by Iosifescu-Manu [26] when the semi-Markov process is nonhomogeneous with a finite state space and finally by Limnios and Opriçan [10] when the semi-Markov process is time homogeneous with arbitrary state space. Theorem (1) generalizes all these cases.

3. Numerical analysis of indexed Markov renewal equations In this section we present the numerical solutions of the IMRE (9). Numerical methods for the solution of MRE have been proposed in the recent years by many researchers (see [21,22,14]. This approach uses a discrete semi-Markov process to approximate the continuous one. In [27] and subsequently in [28] two discrete time semi-Markov processes were constructed to bound solutions of the corresponding continuous time MRE. Here we propose a computation procedure for solving IMRE through a time and space discretization. The proposed approximation scheme is applicable to any IMRE where the function g is r-additive, i.e. "x, s, m, T " A1, A2, . . . , An, . . . with Ai Aj = ;

  X [ g x; s; m; Ai ; t; M ¼ gðx; s; m; Ai ; t; MÞ: iP1

Notice that this is the case of probability measures, then it works well for Eqs. (4)–(7). If g is not r-additive the matrix representation we present in next section is not valid. For this reason, henceforth, we shall assume g to be a r-additive function. Our approach extends that of [21] to the more general framework. Let G = E  [s, t] be a closed set with E  R. First, consider a time grid of discretizetion step h1:

 1 ¼ fx1 ðlÞ ¼ lh1 g x with l = u, u + 1, . . . , K, uh1 = s, h1K = t. Then, consider a state space grid with discretization step h2:

 2 ¼ fx2 ðiÞ ¼ ih2 g x with i = d, d + 1, . . . , N, dh2 = inf E, Nh2 = sup E.  1 and x  2 is: The grid generated by the topological product of x

 ¼x 2 x  2 ¼ fxli ¼ ðx1 ðlÞ; x2 ðiÞÞg; x where l = u, u + 1, . . . , K, i = d, d + 1, . . . , N, uh1 = s, h1K = t, dh2 = infE, Nh2 = supE.

4361

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

 : x2 ðiÞ 2 Ag. For all usage levels m and M, define the following grid Now, for each set A # E, denote by Ah2 ¼ fx2 ðiÞ 2 x function:

1 x 2 ! R D:x Dah2 ;Ah2 ðuh1 ; m; Kh1 ; MÞ ¼

X

ð13Þ

gðah2 ; uh1 ; m; x2 ðiÞ; kh1 ; MÞ

x2 ðiÞ2Ah2

: Now consider a generic quadrature formula to approximate the double integral on x

Z Z

t

  Z s Q ðx; s; m; dy; dsÞw y; s; m þ f ðx; n; mÞdn; A; t; M

s

E

s

N X K X



i¼d

wu;k;d;N;i;l ðQ ðah2 ; uh1 ; m; ði þ 1Þh2 ; ðl þ 1Þh1 Þ þ Q ðah2 ; uh1 ; m; ih2 ; lh1 Þ  Q ðah2 ; uh1 ; m; ih2 ; ðl þ 1Þh1 Þ

l¼u

 Q ðah2 ; uh1 ; m; ði þ 1Þh2 ; lh1 ÞÞw ih2 ; lh1 ; m þ

l X

! f ðah2 ; l2 h2 ; mÞwu;k;l ; Ah2 ; Kh1 ; M :

ð14Þ

l2 ¼u

Define the grid function

v ah ;ih 2

2

ðuh1 ; m; lh1 Þ ¼ wu;k;d;N;i;l ðQ ðah2 ; uh1 ; m; ði þ 1Þh2 ; ðl þ 1Þh1 Þ þ Q ðah2 ; uh1 ; m; ih2 ; lh1 Þ  Qðah2 ; uh1 ; m; ih2 ; ðl þ 1Þh1 Þ  Qðah2 ; uh1 ; m; ði þ 1Þh2 ; lh1 ÞÞ

ð15Þ

and set





wih2 ;Ah2 lh1 ; m þ F ah2 ðlh1 ; mÞ; Kh1 ; M ¼ w ih2 ; lh1 ; m þ

l X

! f ðah2 ; l2 h2 ; mÞwu;k;l ; Ah2 ; Kh1 ; M :

ð16Þ

l2 ¼u

Combining (13)–(16) and (19) yields

wah2 ;Ah2 ðuh1 ; m; Kh1 ; MÞ ¼ Dah2 ;Ah2 ðuh1 ; m; Kh1 ; MÞ þ

N X K X

v ah ;ih 2

i¼d

2

ðuh1 ; m; lh1 Þwih2 ;Ah2 ðlh1 ; m þ F ah2 ðlh1 ; mÞ; Kh1 ; MÞ:

l¼u

ð17Þ It is clear that Eq. (17) produces an approximation of Eq. (9). Next results provides a theoretical support to the proposed approach. Let bxc be the floor function. For any h1 > 0 and h2 > 0 define the following stochastic processes:

J hn2 ¼



Jn h2 ; h2 h

1;2 þ U hn1;2 ¼ U n1



Tn h1 ; h1   h1;2 2 f J hn1 ; s; U n1 ds:

T hn1 ¼ Z

h

T n1 h

1 T n1

If we consider the simplest quadrature method (rectangle formula), the following can be proved: h

Lemma 1. The processes J hn2 ; T hn1 ; U n1;2 have a ISMK with support measure h1 N  h2 N and 2 1 qh1;2 ðx; s; m; kh2 ; jh1 Þ :¼ P½J hnþ1 ¼ kh2 ; T hnþ1 ¼ jh1 jJ n ¼ x; T n ¼ s; U n ¼ m

¼ Q ðx; s; m; ðk þ 1Þh2 ; ðj þ 1Þh1 Þ þ Qðx; s; m; kh2 ; jh1 Þ  Q ðx; s; m; ðk þ 1Þh2 ; jh1 Þ  Q ðx; s; m; kh2 ; ðj þ 1Þh1 Þ Proof 2 1 qh1;2 ðx; s; m; kh2 ; jh1 Þ ¼ P½J hnþ1 6 kh2 ; T hnþ1 6 jh1 jJ n ¼ x; T n ¼ s; U n ¼ m 2 1 þ P½J hnþ1 6 ðk  1Þh2 ; T hnþ1 6 ðj  1Þh1 jJ n ¼ x; T n ¼ s; U n ¼ m 2 1 6 kh2 ; T hnþ1 6 ðj  1Þh1 jJ n ¼ x; T n ¼ s; U n ¼ m  P½J hnþ1 2 1  P½J hnþ1 6 ðk  1Þh2 ; T hnþ1 6 jh1 jJ n ¼ x; T n ¼ s; U n ¼ m

¼ P½J nþ1 6 ðk þ 1Þh2 ; T nþ1 6 ðj þ 1Þh1 jJ n ¼ x; T n ¼ s; U n ¼ m þ P½J nþ1 6 kh2 ; T nþ1 6 jh1 jJ n ¼ x; T n ¼ s; U n ¼ m  P½J nþ1 6 ðk þ 1Þh2 ; T nþ1 6 jh1 jJ n ¼ x; T n ¼ s; U n ¼ m  P½J nþ1 6 kh2 ; T nþ1 6 ðj þ 1Þh1 jJ n ¼ x; T n ¼ s; U n ¼ m ¼ Q ðx; s; m; ðk þ 1Þh2 ; ðj þ 1Þh1 Þ þ Qðx; s; m; kh2 ; jh1 Þ  Q ðx; s; m; ðk þ 1Þh2 ; jh1 Þ  Qðx; s; m; kh2 ; ðj þ 1Þh1 Þ: 

ð18Þ

4362

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

Theorem 2. Let Z(t) be a regular ISMP with arbitrary state space. Let U(t) be the corresponding usage process. Given x 2 X, if, 8s 2 Rþ ; f ð; s; Þ is a bounded continuous function, then (i) the process Z h1;2 ðtÞ converges to Z(t) for h1, h2 ? 0 in the Skorohod topology; ii) the process U h1;2 ðtÞ converges to U(t) for h1, h2 ? 0 in the Skorohod topology. Proof. It must be shown that "T > 0 it exists an increasing homomorphism fkh2 g from [0, 1) such that supt6T jkh2 ðtÞ  tj ! 0 for h2 ? 0 and supt6T jZ h1;2 ðkh2 Þ  ZðtÞj ! 0, supt6T jU h1;2 ðkh2 Þ  UðtÞj ! 0 as h1, h2 ? 0. Because Z(t) is regular, Tn ? 1 therefore we can consider T = Tn "n and prove the theorem. First we choose h2 such that 2 h2 < min jT hk2  T hk1 j;

16k6n

T h02 ¼ 0:

ð19Þ

Now let us consider, as in [21] the function

khn2 ðtÞ ¼ T hn2 þ t  T n ;

8t 6 T n :

ð20Þ

It is simple to verify that (20) is an increasing homomorphism mapping the intervals [Tk1, Tk] in notice that

2 ½T hk1 ; T hk2 .

Additionally

jkhn2 ðtÞ  tj ¼ jt þ ðT hn2  T n Þ  tj ¼ jðT hn2  T n Þj ! 0 as h2 ! 0: Condition (19) entails that the processes T hn2 and Tn make the same number of jumps in any time interval, in fact

Nh2 ðkh2 ðtÞÞ ¼ supfn 2 N : T hn2 6 khn2 ðtÞg ¼ supfn 2 N : T hn2 6 T hn2 þ t  T n g ¼ supfn 2 N : T n 6 tg ¼ NðtÞ: Consequently we have 1 jZ h1;2 ðkh2 ðtÞÞj ¼ jJ Nh2 ðkh2 ðtÞÞ  J NðtÞ j ¼ jJ hNðtÞ  J NðtÞ j n

1 and because supt6T n jJ hNðtÞ  J NðtÞ j ! 0 as h1 ? 0 the point (i) is true. As regards to the usage process U(t) we have that

  Z Tn     Z T hn2  h1;2    h1;2    h1;2 h1 f ðJ n1 ; s; U n1 Þds: U n  U n  6 U n1  U n1  þ  h f ðJ n1 ; s; U n1 Þds    T2 T n1 n1

h

h

1;2 Notice that U 01;2 ¼ U 0 , and if we assume that jU n1  U n1 j ! 0 as h1, h2 ? 0, then we have

 Z h2 Z Tn   Tn h1;2   h1 ; s; U n1 Þds  f ðJ n1 ; s; U n1 Þds  h f ðJ n1   T2 T n1 n1  Z Z h2  Z Tn Z Tn   Tn  Tn  h1;2 h1;2 h1;2   h1 h1 1 6  h f ðJn1 ; s; U n1 Þds  f ðJ n1 ; s; U n1 Þds þ  f ðJ hn1 ; s; U n1 Þds  f ðJ n1 ; s; U n1 Þds   T2 T n1 T n1 T n1 n1  Z Z T hn2 Z T hn2 Z Tn   T n1 h1;2 h1;2 h1;2 h1;2   1 1 1 1 ; s; U n1 Þds þ f ðJ hn1 ; s; U n1 Þds  f ðJ hn1 ; s; U n1 Þds  h f ðJ hn1 ; s; U n1 Þds ¼  h f ðJ hn1   T2 2 T T T n1 n1 n n1  Z T n  Z T n1 Z Tn     h1;2 h1;2 h1;2  1 1 1 þ  f ðJ hn1 ; s; U n1 Þ  f ðJ n1 ; s; U n1 Þds ¼  h f ðJ hn1 ; s; U n1 Þds  h f ðJ hn1 ; s; U n1 Þds   T2 2 T n1 T n n1  Z T n   h1;2 1 þ  f ðJ hn1 ; s; U n1 Þ  f ðJ n1 ; s; U n1 Þds: T n1

Now, since f is bounded, it exists mðn  1; h1 ; h2 Þ 2 R and Mðn  1; h1 ; h2 Þ 2 R such that h

1;2 1 mðn  1; h1 ; h2 Þ 6 f ðJ hn1 ; s; U n1 Þ 6 Mðn  1; h1 ; h2 Þ:

From mean value theorem it results that 2 mðn  1; h1 ; h2 ÞðT n1  T hn1 Þ6

Z

T n1 h2 T n1

h

1;2 1 2 f ðJ hn1 ; s; U n1 Þds 6 Mðn  1; h1 ; h2 ÞðT n1  T hn1 Þ:

4363

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

Because mðn  1; h1 ; h2 Þ 2 R and Mðn  1; h1 ; h2 Þ 2 R we have that: 2 mðn  1; h1 ; h2 Þ  ðT n1  T hn1 Þ ! 0 as h2 ! 0; 2 Mðn  1; h1 ; h2 Þ  ðT n1  T hn1 Þ ! 0 as h2 ! 0

R T n1

h

1;2 1 f ðJ hn1 ; s; U n1 Þds ¼ 0. R Tn h1;2 1 Similar arguments produces limh1;2 !0 h2 f ðJ hn1 ; s; U n1 Þds ¼ 0. Tn It remains to prove that

which implies that limh1;2 !0

lim

h1;2 !0

Z

Tn

T n1

h

2 T n1

h

1;2 1 f ðJ hn1 ; s; U n1 Þds ¼

Z

Tn

f ðJ n1 ; s; U n1 Þds:

T n1

This holds true because by integrability of f we have

Z

Tn

T n1

h

h

h

1;2 1;2 1;2 1 1 1 f ðJ hn1 ; s; U n1 Þds ¼ FðJ hn1 ; T n ; U n1 Þ  FðJ hn1 ; T n1 ; U n1 Þ

and because "sf(, s, ) is a continuous function then F(, s, ) is continuous too and then h

h

1;2 1;2 1 1 lim FðJ hn1 ; T n ; U n1 Þ  FðJ hn1 ; T n1 ; U n1 Þ ¼ FðJ n1 ; T n ; U n1 Þ  FðJ n1 ; T n1 ; U n1 Þ ¼

h1;2 !0

Z

Tn

f ðJ n1 ; s; U n1 Þds: 

T n1

Remark 4. The convergence in the Skorohod topology was discussed in [21,22] in the cases of non-homogeneous and homogeneous semi-Markov processes with finite state space respectively. Theorem (2) extends those results by considering the indexed semi-Markov process with arbitrary state space.

4. Matrix representation and algorithm In this section we propose an algorithm useful to solve IMRE by using a matrix representation of the discretized system of Eq. (17). The main difficulty to solve system (17) is due to the presence of the usage index which changes value to the passage of time. Let us denote by jEj the cardinality of E after that the discretization has been executed. Let introduce the following matrices:

Wðlh1 ; m þ F dh2 ðlh1 ; mÞ; Kh1 ; MÞ ¼ ½Wah2 ;Ah2 ðlh1 ; m þ F dh2 ðlh1 ; mÞ; Kh1 ; MÞ; Dðuh1 ; m; Kh1 ; MÞ ¼ ½Dah2 ;Ah2 ðuh1 ; m; Kh1 ; MÞ; Vðuh1 ; m; lh1 Þ ¼ ½v ah2 ;ih2 ðuh1 ; m; lh1 Þ; 3 2 wðlh1 ; m þ F dh2 ðlh1 ; mÞ; Kh1 ; MÞ 7 6 6 wðlh1 ; m þ F ðdþ1Þh2 ðlh1 ; mÞ; Kh1 ; MÞ 7 7 ~ 1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ ¼ 6 Wðlh 7: 6 .. 7 6 . 5 4

ð21Þ

wðlh1 ; m þ F Nh2 ðlh1 ; mÞ; Kh1 ; MÞ

Consequently we can represent system (17) in matrix form:

Wðuh1 ; m; Kh1 ; MÞ ¼ Dðuh1 ; m; Kh1 ; MÞ þ

K X

Vðuh1 ; m; lh1 Þ  Wðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ

ð22Þ

l¼uþ1

Let MðBi; CÞ be the set of (Bi  C) dimensional matrices. Definition 9. Define the i-backslash operator as the operator

ni : MðBi; CÞ ! MðBi; BCÞ: To each matrix A 2 MðBi; CÞ the i-backslash operator associates a diagonal block matrix ni A 2 MðBi; BCÞ having B blocks each of dimension i  C. The ath block has generic element (a, b) with (a  1)i < a 6 ai, (a  1)C < b 6 aC defined by

ni Aða; bÞ ¼ Aða; b  ða  1ÞCÞ:

ð23Þ

4364

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

The following examples make clarifies how the backslash operator acts. 2 3 2 3 a11 a12 0 0 0 a11 a12 0 4 5 4 If A ¼ a21 a22 then n1 A ¼ 0 0 a21 a22 0 0 5: 0 0 0 0 a31 a32 a31 a32 2

b11 6 b21 If B ¼ 6 4 b31 b41

3 2 b11 b12 7 6 b21 b22 7 then n2 B ¼ 6 4 0 b32 5 b42 0

b12 b22 0 0

0 0 b31 b41

3 0 0 7 7: b32 5 b42

The following proposition holds true: Proposition 2. The Eq. (22) has the following representation:

n1 Wðuh1 ; m; Kh1 ; MÞ ¼ n1 Dðuh1 ; m; Kh1 ; MÞ þ

K X

~ 1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ n1 Vðuh1 ; m; lh1 Þ  njEj Wðlh

ð24Þ

l¼uþ1

Proof. Let consider the (w, z) element of matrix Eq. (24):

n1 Ww;z ðuh1 ; m; Kh1 ; MÞ ¼ n1 Dw;z ðuh1 ; m; Kh1 ; MÞ þ

K X X l¼uþ1

~ k;z ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ:  n1 V w;k ðuh1 ; m; lh1 Þ  njEj W

k

ð25Þ The application of the n1 operator entails that each row constitutes a block, consequently to be n1Ww,z(uh1, m; Kh1, M) – 0 it should be true that (w  1)jEj < z 6 wjEj, 0 < w 6 jEj. Moreover

X

~ k;z ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ n1 V w;k ðuh1 ; m; lh1 Þ  njEj W

k

¼

wjEj X

~ k;z ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ n1 V w;k ðuh1 ; m; lh1 ÞnjEj W

ð26Þ

k¼ðw1ÞjEjþ1

because different values of the index k select elements outside the non zero blocks. From the definition of backslash operator we have:

n1 Ww;z ðuh1 ; m; Kh1 ; MÞ ¼ Ww;zðw1ÞjEj ðuh1 ; m; Kh1 ; MÞn1 Dw;z ðuh1 ; m; Kh1 ; MÞ ~ k;z ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ ¼ Dw;zðw1ÞjEj ðuh1 ; m; Kh1 ; MÞn1 V w;k ðuh1 ; m; lh1 ÞnjEj W ~ k;zðw1ÞjEj ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ ¼ V w;kðw1ÞjEj ðuh1 ; m; lh1 ÞW ¼ V w;kðw1ÞjEj ðuh1 ; m; lh1 ÞWkðw1ÞjEj;zðw1ÞjEj ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ

ð27Þ

Combining (25)–(27) yields

Ww;zðw1ÞjEj ðuh1 ; m; Kh1 ; MÞ ¼ Dw;zðw1ÞjEj ðuh1 ; m; Kh1 ; MÞ þ

K X

wjEj X

V w;kðw1ÞjEj ðuh1 ; m; lh1 ÞWkðw1ÞjEj;zðw1ÞjEj ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ;

l¼uþ1 k¼ðw1ÞjEjþ1

ð28Þ that is

Ww;zðw1ÞjEj ðuh1 ; m; Kh1 ; MÞ ¼ Dw;zðw1ÞjEj ðuh1 ; m; Kh1 ; MÞ þ

jEj K X X

V w;b ðuh1 ; m; lh1 ÞWb;zðw1ÞjEj ðlh1 ; m þ Fðlh1 ; mÞ; Kh1 ; MÞ:

ð29Þ

l¼uþ1 b¼1

Eq. (29) coincides with the equation for the (w, z  (w  1)jEj) element of system (22). Then a solution of system (24) is a solution of system (22). h To solve (22), it is possible to use the following algorithm: – as first step, "a = u, . . . , K by putting uh1 = Kh1 = ah1 and expressing (22) in matrix notation we get

n1 Wðah1 ; m; ah1 ; MÞ ¼ n1 Dðah1 ; m; ah1 ; MÞ ~ which are known terms for all values of m and M. By relation (21) we have n1 Wðah 1 ; m; ah1 ; MÞ.

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

4365

– as second step, "a = u, . . . , K  1, by putting uh1 = ah1, Kh1 = (a + 1)h1 we get

n1 Wðah1 ; m; ða þ 1Þh1 ; MÞ ¼ n1 Dðah1 ; m; ða þ 1Þh1 ; MÞ þ

aþ1 X

~ 1 ; m þ Fðlh1 ; mÞ; ða þ 1Þh1 ; MÞ n1 Vðah1 ; m; lh1 Þ  njEj Wðlh

l¼aþ1

¼ n1 Dðah1 ; m; ða þ 1Þh1 ; MÞ þ n1 Vðah1 ; m; ða þ 1Þh1 Þ ~  n Wðða þ 1Þh1 ; m þ Fðða þ 1Þh1 ; mÞ; ða þ 1Þh1 ; MÞ jEj

ð30Þ

~ ~ all terms on the right-hand side are known, then we get for all m and M n1 Wðah 1 ; m; ða þ 1Þh1 ; MÞ and Wðah1 ; m; ða þ 1Þh1 ; MÞ. – as third step, "a = u, . . . , (K  2), by putting uh1 = ah1 and Kh1 = (a + 2)h1 we get

n1 Wðah1 ; m; ða þ 2Þh1 ; MÞ ¼ n1 Dðah1 ; m; ða þ 2Þh1 ; MÞ þ

aþ2 X

~ 1 ; m þ Fðlh1 ; mÞ; ða þ 2Þh1 ; MÞ n1 Vðah1 ; m; lh1 Þ  njEj Wðlh

l¼aþ1

~ ¼ n1 Dðah1 ; m; ða þ 2Þh1 ; MÞ þ n1 Vðah1 ; m; ða þ 1Þh1 Þ  njEj Wðða þ 1Þh1 ; m þ Fðða þ 1Þh1 ; mÞ; ða þ 2Þh1 ; MÞ þ n1 Vðah1 ; m; ða þ 2Þh1 Þ ~  njEj Wðða þ 2Þh1 ; m þ Fðða þ 2Þh1 ; mÞ; ða þ 2Þh1 ; MÞ

ð31Þ

all terms on the right-hand side are known. In general, by iteration of the procedure, it is possible to solve Eq. (22) for all time and usage values. Remark 5. To reduce the computational effort, it is possible to adapt the approach by [29] which consists in the formulation of an initial value problem for computing the state probabilities. In this case the computational effort is less demanding. In spite of this, our results have been presented to extend the more general scheme proposed by Janssen and Manca [21]. 5. Conclusion This paper proposes a non-homogeneous age-usage semi-Markov model with a measurable state space. The system analysis leads to the investigation of indexed Markov renewal equations. Sufficient condition for existence and uniqueness of the solutions of indexed Markov renewal equations are provided. The numerical solution of the equations is made available by means of a discretization scheme. The discrete process is showed to converge to the continuous one in the Skorohod topology. The solution of the discretized system of equations is obtained by means of a particular matrix representation based on the backslash matrix operator. The model considers simultaneously relevant aspects of mathematical modeling such as the non-homogeneity, the usage accumulation, the general state space. References [1] G.E. Baggs, H.N. Nagagaja, Reliability properties of order statistics from bivariate exponential distributions, Commun. Stat. Stochastic Models 12 (1996) 611–631. [2] J. Hunter, Renewal theory in two dimensions: basic results, Adv. Appl. Probab. 6 (1974) 376–391. [3] N.D. Singpurwalla, S.P. Wilson, Failure models indexed by two scales, Adv. Appl. Probab. 30 (4) (1998) 1058–1072. [4] J.A. Sang-Chin Yang Nachlas, Bivariate reliability and availability modeling, IEEE Trans. Reliab. 50 (1) (2001) 26–35. [5] D.K. Manna, S. Pal, S. Sinha, A use-rate based failure model for two-dimensional warranty, Comput. Ind. Eng. 52 (2007) 229–240. [6] N. Jack, B.P. Iskandar, D.N.P. Murthy, A repair-replace strategy based on usage rate for items sold with a two-dimensional warranty, Reliab. Eng. Syst. Saf. 94 (2009) 611–617. [7] K.B. Kordonsky, I. Gertsbakh, System state monitoring and lifetime scales-I, Reliab. Eng. Syst. Saf. 47 (1995) 1–14. [8] K.B. Kordonsky, I. Gertsbakh, System state monitoring and lifetime scales-II, Reliab. Eng. Syst. Saf. 49 (1995) 145–154. [9] R.E. Barlow, L.C. Hunter, System efficiency and reliability, Technometrics 2 (1) (1960) 43–53. [10] N. Limnios, G. Opriçan, A unified approach for reliability and performability evaluation of semi-Markov systems, Appl. Stochastic Mod. Busin. Indus. 15 (1999) 353–368. [11] M. Carravetta, R. De Dominicis, M.R. D’Esposito, Semi-Markov processes in modelling and reliability: a new approach, Appl. Math. Model. 5 (4) (1981) 269–274. [12] G. Ciardo, R.A. Marie, B. Sericola, K.S. Trivedi, Performability analysis using semi-Markov reward processes, IEEE Trans. Comput. 39 (1990) 1252–1264. [13] V. Barbu, M. Boussemart, N. Limnios, Discrete time semi-Markov model for reliability and survival analysis, Commun. Stat. Theory Methods 33 (2004) 2833–2868. [14] A. Blasi, J. Janssen, R. Manca, Numerical treatment of homogeneous and non-homogeneous semi-Markov reliability models, Commun. Stat. Theory Methods 33 (2004) 697–714. [15] J. Janssen, R. Manca, Semi-Markov Risk Models for Finance, Insurance and Reliability, Springer, 2007. [16] G. D’Amico, The crossing barrier of a non-homogeneous semi-Markov chain, Stochastics Int. J. Probab. Stochastic Process. 81 (6) (2009) 589–600. [17] N. Limnios, G. Oprisßan, Semi-Markov Processes and Reliability, Birkhäuser, Boston, 2001. [18] V. Girardin, N. Limnios, Entropy for semi-Markov processes with Borel state spaces: asymptotic equirepartition properties and invariance principles, Bernoulli 12 (2006) 1–19. [19] D.S. Silvestrov, Recurrence relations for generalized hitting times for semi-Markov processes, Ann. Appl. Probab. 6 (2) (1996) 617–646. [20] M. Telek, A. Horvath, G. Horvath, Analysis of inhomogeneous Markov reward models, Linear Algebra Appl. 386 (2004) 383–405. [21] J. Janssen, R. Manca, Numerical solution of non-homogeneous semi-Markov processes in transient case, Methodol. Comput. Appl. Probab. 3 (2001) 271–293. [22] G. Corradi, J. Janssen, R. Manca, Numerical treatment of homogeneous semi-Markov processes in transient case: a straightforward approach, Methodol. Comput. Appl. Probab. 6 (2004) 233–246. [23] J. Janssen, R. De Dominicis, An algorithmic approach to non-homogeneous semi-Markov processes, Insur. Math. Eco. 3 (1984) 157–165.

4366 [24] [25] [26] [27]

G. D’Amico / Applied Mathematical Modelling 35 (2011) 4354–4366

J. Janssen, R. Manca, A realistic non-homogeneous stochastic pension funds model on scenario basis, Scand. Actuar. J. (1997) 113–137. E. Çinlar, Markov renewal theory, Adv. Appl. Probab. 1 (1969) 123–187. A. Iosifescu-Manu, Non-Homogeneous semi-Markov processes, Stud. Cercet Mat. 24 (1972) 529–533. D.A. Elkins, M.A. Wortman, On numerical solution of the markov renewal equation: tight upper and lower kernel bounds, Methodol. Comput. Appl. Probab. 3 (2001) 239–253. [28] S. Mercier, Numerical bounds for semi-Markovian quantities and application to reliability, Methodol. Comput. Appl. Probab. 10 (2008) 179–198. [29] M.C. Moura, E.L. Droguett, Mathematical formulation and numerical treatment based on transition frequency densities and quadrature methods for non-homogeneous semi-Markov processes, Reliab. Eng. Syst. Saf. 94 (2009) 342–349.