text:
Nam Hoang

Nguyen

Abdel Lisser
email: abdel.lisser@centralesupelec.fr

Vikas Vikram Singh
email: vikassingh@maths.iitd.ac.in

Distributionally robust chance-constrained Markov decision processes with random transition probabilities

Keywords: Markov decision processes, Distributionally robust chance-constrained optimization, Random transition probabilities, Machine replacement problem MSC Classification: 90C15, 90C25, 90C40, 90C59

Markov decision process (MDP) is a decision making framework where a decision maker is interested in maximizing the expected discounted value of a stream of rewards received at future stages at various states which are visited according to a controlled Markov chain. Many algorithms including linear programming methods are available in the literature to compute an optimal policy when the rewards and transition probabilities are deterministic. In this paper, we consider an MDP problem where the reward vector is known and the transition probability vector is a random vector which follow a discrete distribution whose information is not completely known. We formulate the MDP problem using distributionally robust chance-constrained optimization framework under various types of moments based uncertainty sets, and statistical-distance based uncertainty sets defined using ϕ-divergence and Wasserstein distance metric. For each uncertainty set, we propose an equivalent mix-integer bilinear programming problem or a mix-integer semidefinite programming problem with bilinear constraints. As an application, we study a machine replacement problem and perform numerical experiments on randomly generated instances.

Introduction

An MDP is a decision making framework to model the performance of a stochastic system which evolves over time according to a controlled Markov chain. We consider the case where the system has a finite number of states. At time t = 0, the system is at some initial state s 0 ∈ S, where S is a finite state space, according to an initial distribution γ, and a decision maker chooses an action a 0 ∈ A(s 0 ), where A(s 0 ) denotes the set of finite number of actions available to the decision maker at state s 0 . As a consequence a reward R(s 0 , a 0 ) is earned and at time t = 1, the system moves to a new state s 1 with probability p(s 0 , a 0 , s 1 ). The same thing repeats at time t = 1 and it continues for the infinite horizon. We assume that the reward and transition probabilities are stationary, i.e., R(X t = s, A t = a) = R(s, a) and P (X t+1 = s ′ |X t = s, A t = a) = p(s, a, s ′ ) for all t; X t and A t denote the state and action at time t, respectively. The decision taken at time t, which could be deterministic or randomized, may depend on the history h t at time t, where h t = (s 0 , a 0 , s 1 , . . . , s t-1 , a t-1 , s t ). Let H t be the set of all possible histories at time t. A history dependent decision rule f t at time t is defined as f t (h t ) ∈ ℘(A(s t )) for every history h t with final state s t , where ℘(A(s t )) denotes the set of probability distributions on the action set A(s t ). A sequence of history dependent decision rules f h = (f t ) ∞ t=0 is called a history dependent policy. A history dependent policy (f t ) ∞ t=0 is called a stationary policy if there exists a decision rule f such that f t = f for all t and it depends only on the current state. We denote a stationary policy, with some abuse of notations, by f and define f = (f (s)) s∈S such that f (s) ∈ ℘(A(s)) for every s ∈ S. As per stationary policy f , whenever the Markov chain visits state s, the decision maker chooses an action a with probability f (s, a). We denote the set of all history dependent and stationary policies by P O HD and P O S , respectively. A history dependent policy f h ∈ P O HD and initial distribution γ defines a probability measure P f h γ over the state and action trajectories, and E f h γ denotes the expectation operator corresponding to the probability measure P f h γ . For a given policy f h and an initial distribution γ, the expected discounted reward at a discount factor α ∈ (0, 1) is defined as [START_REF] Altman | Constrained Markov Decision Processes[END_REF]; [START_REF] Puterman | Markov Decision Processes[END_REF]

V α (f h , p) = (1 -α)E f h ∞ t=0 α t R(X t , A t ) = s∈S a∈A(s)
m(f h , p; s, a)R(s, a).

(1)

The set { m(f h , p; s, a)} (s,a) is the occupation measure defined by m(f h , p; s, a) = (1 -α) ∞ t=0 α t P f h γ (X t = s, A t = a), ∀ s ∈ S, a ∈ A(s).

(2)

It is well known that for a discounted MDP problem there exists an stationary optimal policy. For a given transition probabilities p(s, a, s ′ ) for all s, s ′ ∈ S, a ∈ A(s), it follows from Theorem 3.2 of [START_REF] Altman | Constrained Markov Decision Processes[END_REF] that the set of occupation measures m(f, p), for f ∈ P O S is equivalent to the following set

Q(p) = ρ ∈ R |K| (s,a)∈K
ρ(s, a) δ(s ′ , s) -αp(s, a, s ′ ) = (1 -α)γ(s ′ ), ∀ s ′ ∈ S, ρ(s, a) ≥ 0, ∀ (s, a) ∈ K , where δ(s ′ , s) is the Kronecker delta and K = {(s, a) | s ∈ S, a ∈ A(s)}. For each ρ ∈ Q(p), the stationary optimal policy f can be defined as f (s, a) = ρ(s, a)

a∈A(s) ρ(s, a)

, ∀ (s, a) ∈ K, whenever the denominator is nonzero if it is zero, we choose f (s) arbitrarily from ℘(A(s)) [START_REF] Altman | Constrained Markov Decision Processes[END_REF]. We restrict our attention to the stationary policies in the rest of the paper.

The rewards and transition probabilities represented by |K|-dimensional vector R = (R(s, a)) (s,a)∈K and |K|.|S|-dimensional vector p = (p(s, a, s ′ )) (s,a,s ′ )∈K×S , respectively, are considered as the parameters of an MDP model and are assumed to be exactly known. However, in practice R and p are not known in advance and are estimated from historical data. This leads to errors in the optimal policies [START_REF] Mannor | Bias and variance approximation in value function estimates[END_REF]. Most efforts to take into account this uncertainty focused on the study of robust MDPs where the rewards or the transition probabilities are known to belong to a prespecified uncertainty set [START_REF] Iyengar | Robust dynamic programming[END_REF]; [START_REF] Nilim | Robust control of Markov decision processes with uncertain transition matrices[END_REF]; Varagapriya et al (2022); [START_REF] White | Markov decision processes with imprecise transition probabilities[END_REF]; [START_REF] Wiesemann | Robust Markov decision processes[END_REF]; [START_REF] Ho | Robust ϕ-divergence MDPs[END_REF]; [START_REF] Goyal | Robust Markov decision processes: Beyond rectangularity[END_REF]. However, [START_REF] Delage | Percentile optimization for Markov decision processes with parameter uncertainty[END_REF] showed that the robust MDP approach usually leads to conservative policies. For this reason, a chance-constrained Markov decision process (CCMDP) was introduced in Delage and [START_REF] Delage | Percentile optimization for Markov decision processes with parameter uncertainty[END_REF], where the controller obtains the expected discounted reward with certain confidence. In [START_REF] Delage | Percentile optimization for Markov decision processes with parameter uncertainty[END_REF], the case of random rewards and random transition probabilities are considered separately and it is shown that a CCMDP is equivalent to a second-order cone programming (SOCP) problem when the running reward vector follows a multivariate normal distribution and the transition probabilities are exactly known. When the transition probabilities follow Dirichlet distribution and the running rewards are exactly known, the CCMDP problem becomes intractable and the optimal policies can be computed using approximation methods. [START_REF] Varagapriya | Joint chance-constrained Markov decision processes[END_REF] considered a constrained MDP problem where the running cost vectors are random vectors and the transition probabilities are known. They formulated it as a joint chance-constrained MDP problem and proposed two SOCP based approximations which give upper and lower bounds to the optimal value of joint chance-constrained MDP problem if the cost vectors follow multivariate elliptical distributions and the dependence among the constraints is driven by a Gumbel-Hougaard copula.

In many practical situations, it is often the case that only a partial information about the underlying distribution is available based on the historical data. In that case, a distributionally robust approach, is used to model the uncertainties, which assumes that the true distribution belongs to an uncertainty set based on its partially available information. Such an approach has been used in modelling the uncertainties of many optimization and game problems [START_REF] Jiang | Data-driven chance constrained stochastic program[END_REF]; [START_REF] Liu | Distributionally robust chance constrained geometric optimization[END_REF]; [START_REF] Singh | Distributionally robust chance-constrained games: Existence and characterization of Nash equilibrium[END_REF]. There are at least two popular ways to construct an uncertainty set for the distribution of the uncertain parameters. The first one is based on the partial information on moments of the true distribution and the second one is based on the statistical distance between the true distribution and a reference distribution. The moments-based uncertainty sets assume certain conditions on the first two moments [START_REF] Cheng | Distributionally robust stochastic knapsack problem[END_REF]; [START_REF] Delage | Distributionally robust optimization under moment uncertainty with application to data-driven problems[END_REF]; [START_REF] Popescu | Robust mean-covariance solutions for stochastic optimization[END_REF]. The statistical distancebased uncertainty sets contain all the distributions which lie inside a ball of small radius and center at a reference distribution which is usually considered to be an empirical distribution or a normal distribution [START_REF] Esfahani | Data-driven distributionally robust optimization using the Wasserstein metric: Performance guarantees and tractable reformulations[END_REF]; [START_REF] Jiang | Data-driven chance constrained stochastic program[END_REF]. To define a distance between the distributions, either a ϕ-divergence Ben-Tal et al ( 2013); [START_REF] Jiang | Data-driven chance constrained stochastic program[END_REF] or Wasserstein distance metric is used [START_REF] Esfahani | Data-driven distributionally robust optimization using the Wasserstein metric: Performance guarantees and tractable reformulations[END_REF]; [START_REF] Gao | Distributionally robust stochastic optimization with Wasserstein distance[END_REF]; [START_REF] Zhao | Data-driven risk-averse stochastic optimization with Wasserstein metric[END_REF]. A discounted MDP problem with uncertain transition probabilities under distributionally robust optimization framework is considered in the literature [START_REF] Xu | Distributionally robust Markov decision processes[END_REF]; Zhi [START_REF] Chen | Distributionally robust optimization for sequential decision-making[END_REF]. They considered the case where the decision maker aims to find an optimal policy which maximizes worst-case expected discounted reward. To the best of our knowledge, an MDP problem with uncertain transition probabilities under distributionally robust chance constraint based payoff criterion is not considered in the literature.

In this paper, we consider an infinite horizon MDP with discounted payoff criterion where the reward vector is known and transition probabilities are defined by a random vector. The transition probability vector is assumed to follow a discrete distribution which is not completely known and belong to a given uncertainty set. We formulate the random transition probability vector with a distributionally robust chance constraint which guarantees the maximum reward for a given policy with at least a given level of confidence. We call this class of MDP a distributionally robust chance-constrained Markov decision process (DRCCMDP), where we consider both moments and statistical distance based uncertainty sets. The main contributions of the paper are as follows.

1. We consider three different types of uncertainty sets based on the moments of the random transition probabilities. We show that the DRCCMDP problem can be reformulated as a mixed-integer bilinear programming (MIBP) problem or a mixed-integer semidefinite programming (MISDP) problem with additional bilinear constraints. The MIBP problems can be solved efficiently using the GUROBI solver, while the MISDP problems with bilinear constraints can be handled by the CUTSDP solver available in the YALMIP toolbox of Matlab, which is time consuming, without any guarantee of the running time. 2. We consider two types of uncertainty sets based on statistical distance between the true distribution of transition probabilities and a reference distribution. The uncertainty sets are constructed using either ϕ-divergences distance or Wasserstein distance. For ϕ-divergences distance, we explore four distinct types of ϕ-divergences (Kullbach-Leibler, Variation, Modified χ 2 , Hellinger) to construct uncertainty sets. For each uncertainty set, we show that the DRCCMDP problem is equivalent to MIBP problem. 3. We illustrate our theoretical results on a machine replacement problem.

The paper is organized as follows. In Section 2, we define a DRCCMDP under a discounted payoff criterion with random transition probabilities and known reward vector. In Section 2.1 and Section 2.2, we propose equvalent reformulations of DRC-CMDP under different moments based and statistical distance based uncertainty sets. We present how to solve reformulations using existing solvers in Section 2.3. The numerical results on a machine replacement problem is given in Section 3. We conclude the paper in Section 4.

Distributionally robust chance-constrained Markov decision process

In this section, we consider an MDP model defined in Section 1, where the running reward vector R is exactly known, nonnegative and the transition probability vector p are random variables. For each triplet (s, a, s ′ ) ∈ K×S, we assume that the p(s, a, s ′ ) is an 1-dimensional random variable defined on a probability space (Ω, F, P). Therefore, for each realization ω ∈ Ω, the term p(s, a, s ′ )(w) ∈ [0, 1] and s ′ ∈S p(s, a, s ′ )(ω) = 1. It represents the probability of moving to a new state s ′ , when an action a in taken at state s. Assume that p follows a discrete distribution F p , whose support is taken by the set of historical data on the transition probabilities. Denote this set by E p = {p 1 , p 2 , . . . , p J }. For a given policy f ∈ P O S and a realization ω ∈ Ω, the expected discounted reward V α (f, p)(ω) can be written as

V α (f, p)(ω) = s∈S a∈A(s) m(f, p, s, a)(ω)R(s, a), (3)
where m(f, p)(ω) is an occupation measure corresponding to transition probability vector p(ω). Since the transition probabilities are random variables, it is clear that m(f, p) is a |K|-dimensional random vector and V α (f, p) is an 1-dimensional random variable defined on same probability space (Ω, F, P). We consider the case where decision maker is interested in maximizing the expected discounted reward which can be obtained with at least a given confidence level (1 -ϵ). This leads to the following chance-constrained optimization problem

(CCMDP) sup y∈R, f ∈P O S y s.t. P p (V α (f, p) ≥ y) ≥ 1 -ϵ, (4)
In most of the practical situations, we only have partial information about the underlying probability distributions of p based on historical data of the transition probabilities. Such situations can be modelled with the distributionally robust optimization approach, where the decision marker believes that the distribution of p belongs to some uncertainty set D p . To ensure that the chance constraint P(V α (f, p) ≥ y) ≥ 1-ϵ holds, we assume that it holds for any distribution which belongs to the uncertainty set. This leads to the following distributionally robust chance-constrained optimization problem

(DRCCMDP) sup y∈R, f ∈P O S y s.t. (i) inf Fp∈Dp P p (V α (f, p) ≥ y) ≥ 1 -ϵ. (5)
In the following sections, we consider two types of uncertainty sets with moments based uncertainty sets and statistical distance based uncertainty sets.

Moments based uncertainty sets

We consider the uncertainty sets which are constructed based on available information about the first two moments of the true distribution of p. For each s ′ ∈ S, consider a |K|-dimensional sub-vector p(s ′ ) = (p(s, a, s ′ )) (s,a)∈K of p. We estimate sample mean vector µ ∈ R |S|×|K| of p and positive definite sample covariance matrix Σ(s ′ ) of p(s ′ ) for all s ′ ∈ S by observing sufficiently large number of data. We consider three most popular uncertainty sets, based on estimates µ and Σ(s ′ ), which have been studied in the literature Delage and Ye (2010); [START_REF] Cheng | Distributionally robust stochastic knapsack problem[END_REF]; [START_REF] Calafiore | On distributionally robust chance-constrained linear programs[END_REF]. They are defined as follows:

1. Uncertainty set with known mean and known covariance matrix:

D 1 =    F p ∈ M + Ep (i) E(1 {p∈Ep} ) = 1, (ii) E(p) = µ, (iii) E[(p(s ′ ) -µ(s ′ ))(p(s ′ ) -µ(s ′ )) T ] = Σ(s ′ ), s ′ ∈ S.    , (6
) 2. Uncertainty set with known mean and unknown covariance matrix:

D 2 =    F p ∈ M + Ep (i) E(1 {p∈Ep} ) = 1, (ii) E(p) = µ, (iii) E[(p(s ′ ) -µ(s ′ ))(p(s ′ ) -µ(s ′ )) T ] ⪯ δ 0 Σ(s ′ ), s ′ ∈ S.    , (7)
3. Uncertainty set with unknown mean and unknown covariance matrix:

D 3 =    F p ∈ M + Ep (i) E(1 {p∈Ep} ) = 1, (ii) [E (p(s ′ )) -µ(s ′ )] T Σ(s ′ ) (-1) [E (p(s ′ )) -µ(s ′ )] ≤ δ 1 , s ′ ∈ S, (iii) E[(p(s ′ ) -µ(s ′ ))(p(s ′ ) -µ(s ′ )) T ] ⪯ δ 2 Σ(s ′ ), s ′ ∈ S.    , (8)
where M + Ep is the set of all probability distributions on E p , and δ 1 ≥ 0, δ 2 , δ 0 ≥ 1. The notation A ⪯ B implies that B -A is a positive semidefinite matrix. We denote the set of n × n (positive semidefinite) symmetric matrices by (S n + ) S n and • denotes the Frobenius product. Using these notations, we present the deterministic reformulation of DRCCMDP problem (5) for each type of moments based uncertainty set defined above. Theorem 1. For DRCCMDP problem (5), the following results hold.

(i) If the true distribution of p belongs to the uncertainty set D 1 , then (5) is equivalent to the following deterministic problem

sup y s.t. (i) -v -w T µ - s ′ ∈S z(s ′ ) • Σ(s ′ ) ≥ 1 -ϵ, (ii) 1 {Vα(f,pj )≥y} + v + w T p j + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (iii) z(s ′ ) ∈ S |K| , ∀ s ′ ∈ S (iv) a∈A(s) f (s, a) = 1, ∀ s ∈ S, f (s, a) ≥ 0, ∀ (s, a) ∈ K. (9)
(ii) If the true distribution of p belongs to the uncertainty set D 2 , then the optimization problem (5) is equivalent to the following deterministic problem sup y

s.t. (i) -v -w T µ - s ′ ∈S z(s ′ ) • δ 0 Σ(s ′ ) ≥ 1 -ϵ, (ii) 1 {Vα(f,pj )≥y} + v + w T p j + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (iii) z(s ′ ) ∈ S |K| + , ∀ s ′ ∈ S, (iv) a∈A(s) f (s, a) = 1, ∀ s ∈ S, f (s, a) ≥ 0, ∀ (s, a) ∈ K. ( 10
)
(iii) If the true distribution of p belongs to the uncertainty set D 3 , then the optimization problem (5) is equivalent to the following deterministic problem sup y

s.t. (i) -v - s ′ ∈S z(s ′ ) • δ 2 Σ(s ′ ) ≥ 1 -ϵ, (ii) 1 {Vα(f,pj )≥y} + v - s ′ ∈S Q(s ′ ) • N j (s ′ ) + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (iii) z(s ′ ) ∈ S |K| + , Q(s ′ ) ∈ S |K|+1 + , ∀ s ′ ∈ S, (iv)
a∈A(s) f (s, a) = 1, ∀ s ∈ S, f (s, a) ≥ 0, ∀ (s, a) ∈ K. ( 11
)
where

N j (s ′ ) = Σ(s ′ ) p j (s ′ ) -µ(s ′ ) (p j (s ′ ) -µ(s ′ )) T δ 1 . ( 12
)
Proof. For any f ∈ P O S and y ∈ R, we consider the following optimization problem inf

Fp∈Dp P p (V α (f, p) ≥ y) . ( 13
)
Note that p is a discrete distribution with finite support E p = {p 1 , . . . , p J }, then we can represent its true distribution F p by its probability mass function which is a J-dimensional vector q = (q 1 , . . . , q J ) such that J j=1 q j = 1, q j ≥ 0 for all j = 1, 2, . . . , J. By representing optimization problem (13) in terms of variable q, we show the deterministic reformulation of (5) for each uncertainty set.

(i) For uncertainty set D 1 , the optimization problem (13) can be rewritten as follows

inf q≥0 J j=1 q j 1 {Vα(f,pj )≥y} s.t. (i) J j=1 q j = 1, (ii) J j=1 q j p j = µ, (iii) J j=1 q j [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T = Σ(s ′ ), ∀ s ′ ∈ S. (14)
The dual problem of ( 14) can be written as follows sup

(v,w,z) -v -w T µ - s ′ ∈S z(s ′ ) • Σ(s ′ ), s.t. (i) 1 {Vα(f,pj )≥y} + v + w T p j + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (ii) z(s ′ ) ∈ S |K| , ∀ s ′ ∈ S.
where v ∈ R, w ∈ R |K|.|S| and z(s ′ ), s ′ ∈ S, are the dual variables of (i), (ii), and (iii) of ( 14). Note that ( 14) is a linear program(LP), and therefore the strong duality holds which in turn implies that the DRCCMDP problem ( 5) can be reformulated as ( 9). (ii) For the case of uncertainty set D 2 , the optimization problem (13) can be rewritten as follows

inf q≥0 J j=1 q j 1 {Vα(f,pj )≥y} s.t. (i) J j=1 q j = 1, (ii) J j=1 q j p j = µ, (iii)
J j=1 q j [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ⪯ δ 0 Σ(s ′ ), ∀ s ′ ∈ S. (15)
The dual problem of ( 15) is given by sup -v -w T µ -

s ′ ∈S z(s ′ ) • δ 0 Σ(s ′ ), s.t. (i) 1 {Vα(f,pj )≥y} + v + w T p j + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (ii) z(s ′ ) ∈ S |K| + .
Since ( 15) is a semidefinite programming (SDP) problem, the strong duality holds, which ensures that the DRCCMDP problem (5) can be reformulated as ( 10). (iii) For the case of uncertainty set D 3 , it follows from Schur complement that the constraints (ii) in ( 8) are equivalent to N (s ′ ) ∈ S |K|+1 + , for any s ′ ∈ S, where

N (s ′ ) = Σ(s ′ ) E(p(s ′ )) -µ(s ′ ) (E(p(s ′ )) -µ(s ′ )) T δ 1 ,
Moreover, for any s ′ ∈ S, we can write N (s ′ ) = J j=1 q j N j (s ′ ), where N j (s ′ ) is defined in (12). The rest of the proof follows from the same arguments used in the case of uncertainty set D 2 .

The deterministic reformulations ( 9), ( 10) and ( 11) have infinite number of constraints due to indicator function 1 {Vα(f,pj )≥y} . Using the fact that 1 {Vα(f,pj )≥y} can take only binary variables, we further reformulate ( 9), ( 10) and ( 11) as mixed integer optimization problems. We present the proof for (9) and for ( 10) and ( 11) the proof follows using the same arguments. Let M = 1 T |K| R, where 1 |K| denotes |K|-dimensional vector with all components equal to 1, and η > 0 be an infinitesimal number. The following theorem gives a mixed-integer reformulation of (9). Theorem 2. Suppose initial distribution γ = (γ(s)) s∈S satisfy γ(s) > 0 for all s ∈ S. Then, the optimization problem (9) can be reformulated as the following mixed-integer problem

(Bilinear-K-K) sup y,v,w,z,β,(ρj ) J j=1 y s.t. (i) -v -w T µ - s ′ ∈S z(s ′ ) • Σ(s ′ ) ≥ 1 -ϵ, (ii) β j + v + w T p j + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (iii) z(s ′ ) ∈ S |K| , ∀ s ′ ∈ S, (iv) ρ T j R -y ≤ M β j -η(1 -β j ), y -ρ T j R ≤ M (1 -β j ), j = 1, . . . , J, (v) β j ∈ {0, 1} , j = 1, . . . , J, (vi) (s,a)∈K ρ j (s, a) (δ(s ′ , s) -αp j (s, a, s ′ )) = (1 -α)γ(s ′ ), ∀ s ′ ∈ S, (vii) ρ j (s, a) ≥ 0, ∀ (s, a) ∈ K, j = 1, . . . , J, (viii) ρ 1 (s, a) a ′ ∈A(s) ρ j (s, a ′ ) = ρ j (s, a) a ′ ∈A(s) ρ 1 (s, a ′ ), ∀ (s, a) ∈ K, j = 2, . . . , J. (16)
Proof. For an f ∈ P O S , we can write V α (f, p j ) = m(f, p j ) T R, where m(f, p j ) is the occupation measure vector defined by (2) corresponding to policy f and transition probabilities vector p j . Therefore, V α (f, p j ) ∈ [0, M ] which implies that any optimal solution y of (5) belongs to [0, M ]. Then, we can add the constraint y ∈ [0, M ] without loss of optimality. We show that β j = 1 {Vα(f,pj )≥y} is equivalent to following constraints

(i) V α (f, p j ) -y ≤ M β j -η(1 -β j ), y -V α (f, p j ) ≤ M (1 -β j ), j = 1, . . . , J, (ii) β j ∈ {0, 1}. (17) 1. Let β j = 1 {Vα(f,pj )≥y} . If V α (f, p j ) < y, then β j = 0, y -V α (f, p j ) ≤ M and V α (f, p j ) -y ≤ -η for sufficiently small positive η. If V α (f, p j ) ≥ y, then β j = 1, V α (f, p j ) -y ≤ M and y -V α (f, p j ) ≤ 0.
Therefore, the constraints given by ( 17) are satisfied. 2. Suppose the constraints given by ( 17) are satisfied. If β j = 0, then V α (f, p j ) < y.

If

β j = 1, then V α (f, p j ) ≥ y. Therefore, β j = 1 {Vα(f,pj )≥y} .
It follows from the discussion in Section 1 that for all j = 1, . . . , J, the set of occupation measures { m(f, p j ) | f ∈ P O S } is equivalent to the set Q(p j ). Since γ(s) > 0 for all s ∈ S, a ′ ∈A(s) ρ j (s, a ′ ) > 0 for all j = 1, 2, . . . , J. This implies that for every f ∈ P O S , there exists ρ j ∈ Q(p j ) such that

V α (f, p j ) = ρ T j R, ∀ j = 1, . . . , J, f (s, a) = ρ j (s, a) a ′ ∈A(s) ρ j (s, a ′ ) , ∀ (s, a) ∈ K, j = 1, . . . , J.
Hence, optimization problem ( 9) is equivalent to the optimization problem ( 16). The constraint (viii) ensures that all ρ j , j = 1, . . . , J corresponds to the same policy f .

Remark 1. If (y * , v * , w * , z * , β * , (ρ * j ) J j=1 is an optimal solution of (16), the optimal policy of DRCCMDP problem (5) is defined by

f * (s, a) = ρ * 1 (s, a) a ′ ∈A(s) ρ * 1 (s, a ′ ) , ∀ (s, a) ∈ K.
Similarly, the reformulations of DRCCMDP problem (5) for the case of uncertainty sets D 2 and D 3 are given as below (Bilinear-K-U) sup y,v,w,z,β,(ρj 16).

) J j=1 y s.t. (i) -v -w T µ - s ′ ∈S z(s ′ ) • δ 0 Σ(s ′ ) ≥ 1 -ϵ, (ii) β j + v + w T p j + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (iii) z(s ′ ) ∈ S |K| + , ∀ s ′ ∈ S, (iv) -(viii) of (
(

) (Bilinear-U-U) sup y,v,Q,z,β,(ρj ) J j=1 y s.t. (i) -v - s ′ ∈S z(s ′ ) • δ 2 Σ(s ′ ) ≥ 1 -ϵ, (ii) β j + v - s ′ ∈S Q(s ′ ) • N j (s ′ ) + s ′ ∈S z(s ′ ) • [p j (s ′ ) -µ(s ′ )][p j (s ′ ) -µ(s ′ )] T ≥ 0, ∀ j = 1, . . . , J, (iii) z(s ′ ) ∈ S |K| + , Q(s ′ ) ∈ S |K|+1 + , ∀ s ′ ∈ S, 18
(iv) -(viii) of ( 16). ( 19)

Statistical distance based uncertainty sets

In this section, we consider uncertainty sets defined by ϕ-divergence and Wasserstein distance metrics. In such uncertainty sets, a reference distribution ν is known to the decision maker based on the available estimated data of transition probabilities. The decision maker believes that the true distribution of transition probabilities vector p, denoted by F p belongs to a ball centered at the reference distribution ν. We assume that both ν and F p are discrete distributions on same support E p . Let (q 0 j ) J j=1 be the probability mass function of the reference distribution ν, i.e., q 0 j is the weight of jth atom p j . We assume that q 0 j > 0 for all j = 1, . . . , J.

Uncertainty set with ϕ -divergence distance

The ϕ-divergence distance between two discrete probability distributions ν 1 and ν 2 with support E p is given by

I ϕ (ν 1 , ν 2 ) = J j=1 ϕ ν 1 (p j ) ν 2 (p j ) ν 2 (p j ),
where ν 1 (p j ) (resp. ν 2 (p j )) is the weight of ν 1 (resp. ν 2 ) on the jth atom p j of E p and ϕ is a convex function on R + . For general ϕ-divergence with the choices of function ϕ and its properties, we refer to Ben-Tal et al (2013). The uncertainty set of the distribution of p based on ϕ-divergence is defined by

D 4 = F p ∈ M + Ep | I ϕ (F p , ν) ≤ θ ϕ , (20)
where θ ϕ > 0 denotes the radius. The following definition of the conjugate of a function is useful for our subsequent analysis. Definition 1. The conjugate of ϕ is a function ϕ

* : R → R ∪ ∞ such that ϕ * (r) = sup t≥0 {rt -ϕ(t)} , ∀ r ∈ R.
Table 1 presents four special types of ϕ-divergences with their conjugate. -

1, r ≤ -1, r, -1 ≤ r ≤ 1, ∞, r > 1. Modified χ 2 -distance (t -1) 2 . -1, r ≤ -2, r + r 2 4 , r > -2.
Hellinger distance

( √ t -1) 2 . r 1-r , r < 1, ∞, r ≥ 1.
Lemma 1. The dual of an optimization problem

inf Fp∈D4 P Fp ((V α (f, p) ≥ y), (21)
is given by

sup λ>0,β∈R β -λθ ϕ -λP ν (V α (f, p) ≥ y) ϕ * β -1 λ -λ [1 -P ν (V α (f, p) ≥ y)] ϕ * β λ , (22)
and the strong duality between (21) and ( 22) holds.

Proof. The primal problem ( 21) can be written as a following optimization problem

v P = inf q≥0 J j=1 q j 1 {Vα(f,pj )≥y} s.t. (i) J j=1 q 0 j ϕ q j q 0 j ≤ θ ϕ , (ii) J j=1 q j = 1. ( 23
)
The dual problem of ( 23) is given by

v D = sup λ≥0,β∈R β -λθ ϕ + inf q≥0 J j=1
q j 1 {Vα(f,pj )≥y} -βq j + λq 0 j ϕ q j q 0 j , where λ is the dual variable of the constraint (i) of ( 23) and β is the dual variable of the constraint (ii) of ( 23). Note that ( 23) is a convex optimization problem because ϕ is a convex function. Since θ ϕ > 0, it is clear that for q = q 0 , (i) of ( 23) is strictly feasible. Hence, the Slater's condition is satisfied which guarantees the strong duality, i.e., v P = v D . It remains to show that v D is same as ( 22). For a given β ∈ R, consider a function

F (λ, q) = J j=1 q j 1 {Vα(f,pj )≥y} -βq j + λq 0 j ϕ qj q 0 j . It is clear that F (λ, q)
is a linear function of λ. Let G(λ) = -λθ ϕ + inf q≥0 F (λ, q). For any t ∈ [0, 1] and

λ 1 ≥ 0, λ 2 ≥ 0, we have inf q≥0 F [tλ 1 + (1 -t)λ 2 , q] = inf q≥0 (tF (λ 1 , q) + (1 -t)F (λ 2 , q)) ≥ t inf q≥0 F (λ 1 , q) + (1 -t) inf q≥0 F (λ 2 , q).
This implies that G(λ) is a concave function of λ on [0, ∞). Consider the optimization problem sup λ≥0 G(λ). We prove that G(0) ≤ sup λ>0 G(λ). In fact, if G(0) > sup λ>0 G(λ). Then, there exists a t sufficiently close to 1 such that tG( 0

)+(1-t)G(1) > sup λ>0 G(λ). Since G(λ) is a concave function, G(1 -t) ≥ tG(0) + (1 -t)G(1) which in turn implies that G(1 -t) > sup λ>0 G(λ)
. This gives a contradiction. Therefore, G(0) ≤ sup λ>0 G(λ) and we can restrict on λ > 0 without loss of optimality in the dual problem v D . Hence, the dual problem of ( 21) can be rewritten as follows

v D = sup λ>0,β∈R β -λθ ϕ + inf q≥0 J j=1 q j 1 {Vα(f,pj )≥y} -βq j + λq 0 j ϕ q j q 0 j , = sup λ>0,β∈R β -λθ ϕ + J j=1 inf qj ≥0 q j 1 {Vα(f,pj )≥y} -βq j + λq 0 j ϕ q j q 0 j , = sup λ>0,β∈R β -λθ ϕ + J j=1 inf qj ≥0 λq 0 j ϕ q j q 0 j - q j q 0 j β -1 {Vα(f,pj )≥y} λ , = sup λ>0,β∈R β -λθ ϕ - J j=1 sup qj ≥0 λq 0 j q j q 0 j β -1 {Vα(f,pj )≥y} λ -ϕ q j q 0 j . Let t j = qj q 0 j
, using the definition of the conjugate of a function, the dual problem v D can be rewritten as follows

v D = sup λ>0,β∈R β -λθ ϕ - J j=1 λq 0 j ϕ * β -1 {Vα(f,pj )≥y} λ , = sup λ>0,β∈R β -λθ ϕ -λP ν (V α (f, p) ≥ y) ϕ * β -1 λ -λ [1 -P ν (V α (f, p) ≥ y)] ϕ * β λ ,
which completes the proof.

Lemma 2. Consider the DRCCMDP problem (5) under the uncertainty set defined by (20) for the ϕ-divergences listed in Table 2. Then, the DRCCMDP problem (5) can be rewritten as follows

sup y∈R, f ∈P O S y s.t. (i) P ν (V α (f, p) ≥ y) ≥ g * (θ ϕ , ϵ), ( 24
)
where the values of θ ϕ , ϵ and g * (θ ϕ , ϵ) for different ϕ-divergences are given in Table 2.

Table 2 The function f for selected ϕ-divergences

Divergence g * (θ ϕ , ϵ) θ ϕ , ϵ K-L inf x∈(0,1) e -θ ϕ x 1-ϵ -1 x-1 θ ϕ > 0 0 < ϵ < 1 Variation 1 -ϵ + θ ϕ 2 θ ϕ > 0 0 < ϵ < 1 Modified χ 2 1 -ϵ + θ 2 ϕ +4θ ϕ (ϵ-ϵ 2 )-(1-2ϵ)θ ϕ 2θ ϕ +2 θ ϕ > 0 0 < ϵ < 1 2 Hellinger -B+ √ ∆ 2 , where B = -(2 -(2 -θ ϕ ) 2 )ϵ - (2-θ ϕ ) 2 2 , C = (2-θ ϕ ) 2 4 -ϵ 2 , ∆ = B 2 -4C = (2 -θ ϕ ) 2 4 -(2 -θ ϕ ) 2 ϵ(1 -ϵ). 0 < θ ϕ < 2 - √ 2 0 < ϵ < 1
Proof. The details of the proof for the Hellinger distance case is given in Appendix A. The proofs for Kullback-Leibler, Variation distance and Modified χ 2 -distance follow from Propositions 2, 3 and 4 of [START_REF] Jiang | Data-driven chance constrained stochastic program[END_REF], respectively. Theorem 3. Suppose initial distribution γ = (γ(s)) s∈S satisfy γ(s) > 0 for all s ∈ S. Then, the optimization problem (24) can be reformulated as the following mixed-integer optimization problem

(Bilinear-Phi) sup y,β,(ρj ) J j=1 y s.t. (i) J j=1 q 0 j β j ≥ g * (θ ϕ , ϵ), ∀ j = 1, . . . , J,
(iv) -(viii) of ( 16). ( 25)

Proof. Using Lemma 2, the optimization problem (24) can be written as sup y

s.t. (i) J j=1 q 0 j 1 {Vα(f,pj )≥y} ≥ g * (θ ϕ , ϵ), (ii) a∈A(s) f (s, a) = 1, ∀ s ∈ S, f (s, a) ≥ 0, ∀ (s, a) ∈ K. ( 26
)
By taking β j = 1 {Vα(f,pj )≥y} for all j = 1, . . . , J, and using the same arguments used in Theorem 2 the optimization problem ( 26) is equivalent to the mixed integer optimization problem (25)

Uncertainty set with Wasserstein distance

The Wasserstein distance W d (F p , ν) between true distribution F p and reference distribution ν with support on E p , is given by

W d (F p , ν) = inf J i,j=1 ω ij ||p i -p j || d 2 | J j=1 ω ij = q 0 i , J i=1 ω ij = q j , ω ij ≥ 0, ∀ i, j 1 d ,
where d ≥ 1 is some constant. For each i and j, w ij is the weight transported from the ith atom of the reference distribution ν to the jth atom of the true distribution F p and q 0 i is the weight of ith atom of reference distribution ν and q j is the weight of jth atom of F p . For more details on the Wasserstein distance metric, we refer to [START_REF] Villani | Optimal transport: old and new[END_REF]; [START_REF] Villani | Topics in optimal transportation[END_REF]. The uncertainty set of p using Wasserstein distance is defined by

D 5 = F p ∈ M + Ep | W d (F p , ν) ≤ θ W , (27)
where θ W > 0. We have the following lemma. Lemma 3. The DRCCMDP problem (5) under the uncertainty set defined by ( 27) is equivalent to the following optimization problem

sup y,v,h,f y s.t. (i) - J i=1 q 0 i v i -hθ d W ≥ 1 -ϵ, (ii) h ≥ 0, (iii) 1 {Vα(f,pj )≥y} + v i + hc d ij ≥ 0, ∀ i, j = 1, . . . , J, (iv) a∈A(s) f (s, a) = 1, ∀ s ∈ S, f (s, a) ≥ 0, ∀ (s, a) ∈ K, ( 28
)
where

c d ij = ||p i -p j || d 2 .
Proof. The optimization problem inf Fp∈D5 P p (V α (f, p) ≥ y) can be equivalently written as follows inf ω≥0,q≥0 J j=1

q j 1 {Vα(f,pj )≥y} s.t. (i) J j=1 ω ij = q 0 i , ∀ i = 1, . . . , J, (ii) J i=1 ω ij = q j , ∀ j = 1, . . . , J, (iii)
J i,j=1 ω ij c d ij ≤ θ d W . ( 29
)
By replacing q j = J j=1 ω ij in (29), we have

inf ω≥0 J i,j=1 ω ij 1 {Vα(f,pj )≥y} s.t. (i) J j=1 ω ij = q 0 i , ∀ i = 1, . . . , J (ii)
J i,j=1 ω ij c d ij ≤ θ d W . (30)
The dual problem of ( 30) is given by sup v∈R J , h≥0

- J i=1 q 0 i v i -hθ d W s.t. (i) 1 {Vα(f,pj )≥y} + v i + hc d ij ≥ 0, ∀ i, j = 1, . . . , J. (31)
Note that ( 30) is an LP, which ensures that the strong duality holds. Therefore, the DRCCMDP problem ( 5) is equivalent to (28).

Theorem 4. Suppose initial distribution γ = (γ(s)) s∈S satisfy γ(s) > 0 for all s ∈ S.

Then, the optimization problem (28) can be reformulated as the following mixed-integer optimization problem

(Bilinear-Wasserstein) sup y,v,h,β,(ρj ) J j=1 y s.t. (i), (ii) of (28) (ii) β j + v i + hc d ij ≥ 0, ∀ i, j = 1, . . . , J, (iv) -(viii) of (16). ( 32
)
Proof. By taking β j = 1 {Vα(f,pj )≥y} for all j = 1, . . . , J, and using the same arguments used in Theorem 2 the optimization problem ( 28) is equivalent to the mixed integer optimization problem (32)

Solving mixed integer programming problems with bilinear and positive semidefinite cone constraints

The equivalent mixed-integer programming problems proposed in Sections 2.1 and 2.2 to solve DRCCMDP problem (5) for different types of moments based and statistical distance based uncertainty sets are quite similar. For example, they have the same objective function, and the same set of integer and bilinear equality constraints.

Additionally, these problems have (i) linear constraints corresponding to moments based uncertainty set defined by ( 6) and statistical distance based uncertainty sets defined by ( 20) and ( 27), (ii) linear constraints as well as positive semidefinite cone constraints corresponding to moments based uncertainty sets defined by ( 7) and ( 8).

The equivalent optimization problems corresponding to the uncertainty sets defined by ( 6), ( 20) and ( 27) are called MIBP problems. The equivalent optimization problems corresponding to the uncertainty sets defined by ( 7) and ( 8) are called MISDP problems with bilinear constraints.

The MIBP problems can be solved efficiently by GUROBI 9.0 or higher version. It deals bilinear constraints by relaxing them using linear constraints based on McCormick lower and upper envelopes [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF], which depend on the local bounds of the variables present in the bilinear terms. In our case, all the mixed integer programming problems have following bilinear constraints in common

ρ 1 (s, a) a ′ ∈A(s) ρ j (s, a ′ ) = ρ j (s, a) a ′ ∈A(s) ρ 1 (s, a ′ ), ∀ (s, a) ∈ K, j = 2, . . . , J. (33) Let a ′ ∈A(s) ρ j (s, a ′ ) = M j (s), ∀ j = 1, . . . , J, s ∈ S, ρ 1 (s, a)M j (s) = C j (s, a), ∀ j = 2, . . . , J, s ∈ S, a ∈ A(s), ρ j (s, a)M 1 (s) = D j (s, a), ∀ j = 2, . . . , J, s ∈ S, a ∈ A(s).
By introducing auxiliary variables M j (s), C j (s, a), D j (s, a), the constraints (33) are equivalent to the following set of constraints (i) a ′ ∈A(s) ρ j (s, a ′ ) = M j (s), ∀ j = 1, . . . , J, s ∈ S, (ii) C j (s, a) = D j (s, a), ∀ j = 2, . . . , J, s ∈ S, a ∈ A(s), (iii) ρ 1 (s, a)M j (s) = C j (s, a), ∀ j = 2, . . . , J, s ∈ S, a ∈ A(s), (iv) ρ j (s, a)M 1 (s) = D j (s, a), ∀ j = 2, . . . , J, s ∈ S, a ∈ A(s).

(34)

We propose tight bounds for ρ j (s, a) and M j (s) which are used in generating the McCormick envelopes of bilinear terms present in (iii) and (iv) of (34). It is well known that for the occupation measure ρ j there exists a policy f such that (see [START_REF] Altman | Constrained Markov Decision Processes[END_REF])

ρ j (s, a) = f (s, a)M j (s), ∀ s ∈ S, a ∈ A(s), j = 1, . . . , J,
and

M j = (M j (s)) s∈S = (1 -α)γ T [I -αP j (f )] (-1) , = (1 -α)γ T I + ∞ i=1 (αP j (f )) i , ∀ j = 1, . . . , J,
where P j (f ) is a transition probability matrix induced by transition probabilities p j and stationary policy f , which is given by

P j (f )(s, s ′ ) = a∈A(s) f (s, a)p j (s, a, s ′ ), ∀ s, s ′ ∈ S.
Since a∈A(s) f (s, a) = 1, we have P j (f )(s, s ′ ) ≥ min a∈A(s) p j (s, a, s ′ ), for any s, s ′ ∈ S. Let P j,min (s, s ′ ) = min a∈A(s) p j (s, a, s ′ ) and P j,min = (P j,min (s, s ′ )) s,s ′ ∈S be a |S| × |S| matrix. We derive lower bound for M j (s) as follows

M j (s) = (1 -α)γ T I + ∞ i=1 (αP j (f )) i s ≥ (1 -α)γ T I + ∞ i=1 (αP j,min ) i s = (1 -α)γ T [I -αP j,min ] (-1) s . Denote M l j (s) = (1 -α)γ T [I -αP j,min ]
(-1) s and M u j (s) = 1s ′ ̸ =s M l j (s ′ ). Since s∈S M j (s) = 1, a lower bound and an upper bound of M j (s) are given by

M l j (s) ≤ M j (s) ≤ M u j (s), ∀ j = 1, . . . , J, s ∈ S. (35)
Since 0 ≤ ρ j (s, a) ≤ M j (s), for all s ∈ S, a ∈ A(s), j = 1, . . . , J, the lower and upper bounds of ρ j (s, a) are given by 0 ≤ ρ j (s, a) ≤ M u j (s), ∀ j = 1, . . . , J, s ∈ S, a ∈ A(s).

The MISDP problems with bilinear constraints are difficult to solve in general. To the best of our knowledge, there is no commercial solver which can solve efficiently this type of optimization problem. We propose using CUTSDP solver, which is an internal solver in YALMIP toolbox of MATLAB. One of the core ideas in YALMIP is to rely on external solvers for time-consuming tasks. A positive semidefinite constraint A ⪰ 0, where A is an n × n matrix, is equivalent to an infinite number of linear constraints x T Ax ≥ 0, for all x ∈ R n . CUTSDP solver uses a cutting plane (outer approximation) approach, which relaxes positive semidefinite constraint by finite number of linear constraints, solves the relaxation problem and add violated linear cuts to the model and repeat the same procedure. The main idea of CUTSDP is to appropriately generate several vectors x to construct a finite number of linear constraints, that leads to an outer approximation of the original problem. The approximation problem is a MIBP problem which is solved by an external solver (such as GUROBI). If the solution of approximation problem does not satisfy the original positive semidefinite constraint, the solver will add a cutting plane based on a negative eigenvalue, and repeats the process, with a hope that after some iterations, it will eventually satisfy the semidefinite constraint.

Machine replacement problem

In this section, we present a series of numerical experiments performed on a machine replacement problem to compare the approaches discussed earlier. These comparisons aim to evaluate the performance and efficiency of the different reformulations and solvers in solving the respective problems. By conducting these comparisons, we can gain insights into the strengths and limitations of each approach and make informed decisions based on the specific problem characteristics. All the numerical results below are performed using Matlab R2023a on an Intel Core i5-1135G7, Processor 2.4 GHz (8M Cache, up to 4.2 GHz), RAM 16G, 512G SSD.

We consider a machine replacement problem where a machine in a factory has a life-time of N years. At every stage a maintenance of the machine is scheduled but a factory owner can decide whether to repair the machine or not. There is a high probability that the machine behaves like a new one if it is being repaired and its life gets reduced by a year if it is not being repaired. The factory owner incurs maintenance cost if he decides to repair the machine. It can be modelled as an MDP problem where the life of a machine represents the state of underlying Markov chain, i.e., there are N + 1 states. The first state represents a brand new machine. At each state there are two actions: i) "repair", ii) "do not repair". Figure 1 presents a case of fixed transition probabilities of the Markov chain with respect to each action. The maintenance cost corresponding to every state-action pair is not exactly known and is realised after the decision is made. Therefore, it is modelled with a random variable. We assume that for every state action pair (s, a), the maintenance cost is defined as ĉ(s, a) = K + Ẑ(s, a), where K represents the fixed cost and Ẑ(s, a) represents a variable cost which is a random variable. The machine generates a revenue L(s, a) at state-action pair (s, a) and the profit for each (s, a) ∈ K is given by R(s, a) = L(s, a) -K -Ẑ(s, a).

(37)

We define the reward vector R by taking the expected value of R, i.e., R(s, a) = E( R(s, a)), ∀ s ∈ S, a ∈ A(s). The reward vector is given by

R(s, a) = L(s, a) -K -µ Ẑ (s, a), (38)
where K = 10 and L, µ Ẑ are given in Table 3, µ Ẑ is the mean vector of Ẑ. We randomly simulate 100 transition probabilities p j , j = 1, . . . , 100 as follows. For any s, s ′ = 1, . . . , N + 1 and a = 1, 2, we consider the following cases.

• Case 1: If a = 1, i.e., the decision is to repair the machine. Assume that we are actually at state s. If 2 ≤ s ≤ N + 1, we come back to state 1 with high probability q 1,s and stay at state s with probability 1 -q 1,s . Then, p(s, a, s ′ ) = 1 -q 1,s if s ′ = s, p(s, a, s ′ ) = q 1,s if s ′ = 1 and p(s, a, s ′ ) = 0, otherwise. We randomly generate q 1,s on [0.7, 1]. If s = 1, then p(s, a, s ′ ) = 1 if s ′ = 1 and p(s, a, s ′ ) = 0, otherwise. • Case 2: If a = 2, i.e., the decision is to do not repair the machine. If 1 ≤ s ≤ N , we move to the next state s + 1 with high probability q 2,s and stay at actual state s with probability 1 -q 2,s . Then, p(s, a, s ′ ) = 1 -q 2,s if s ′ = s, p(s, a, s ′ ) = q 2,s if s ′ = s + 1 and p(s, a, s ′ ) = 0, otherwise. We randomly generate q 2,s on [0.7, 1]. If s = N + 1, then p(s, a, s ′ ) = 1 if s ′ = N + 1 and p(s, a, s ′ ) = 0, otherwise.

The transition probability vector p is |S|.|K|-dimensional random vector with support on E p = {p 1 , . . . , p J }. Let µ be the sample mean of p and Σ(s ′ ) be the sample covariance matrix of p(s ′ ), for any s ′ ∈ S, where µ and Σ(s ′ ) are given as follows

µ = 1 J J j=1 p j , Σ(s ′ ) = 1 J -1 J j=1 (p j (s ′ ) -µ(s ′ ))(p j -µ(s ′ )) T .
The factory owner is interested in maximizing the expected discounted profit. We assume that the factory owner has a finite number of the same machines which are modelled using the same Markov chain. Therefore, we compute the optimal repair policy with respect to a single machine and the same repair policy can be applied for all other machines. In our numerical experiments, we set the number of states to 10, the threshold value ϵ = 0.1, the discount parameter α = 0.85 and the initial distribution of states γ to be uniformly distributed. For the above instance, |K| = 20. We use the above µ and Σ(s ′ ), s ′ ∈ S for all the moments based uncertainty sets. The reference distribution in all statistical distance based uncertainty sets is assumed to be uniformly distributed on E p . We summarize the other parameters related to all the uncertainty sets in Table 4. state 9, the factory owner decides to repair for moments based uncertainty sets with high probability whereas for statistical distance based uncertainty sets repair action is taken with small probability. At last state repair action is taken with high probability for all types of uncertainty sets. The optimal value of the DRCCDMP problem is more for statistical distance based uncertainty sets as compared to moments based uncertainty sets. We present the time analysis by considering the number of states for all the uncertainty sets between 10 and 5000. All the parameters are taken similar to the case of 10 states. The results are presented in Figure 2 which shows that DRCCMDP problem with statistical distance based uncertainty sets can be solved efficiently up to 5000 states while it takes significantly longer computation time to solve the model with moments based uncertainty sets. For uncertainty set D 1 , the Gurobi solver perform efficiently up to 1000 states whereas for uncertainty sets D 2 and D 3 CUTSDP solver manages only up to 100 states.

Conclusions

We study a DRCCMDP problem with random transition probabilities under various moments and statistical distance based uncertainty sets defined using ϕ-divergence and Wasserstein distance metric. We propose equivalent MIBP problems and MISDP problem with bilinear constraints for the DRCCMDP problem depending on the choice of the uncertainty set. All these optimization problems can be solved efficiently using commercial solver GUROBI, except ones with both bilinear and positive semidefinite constraints. Using randomly generated data, the numerical experiments are performed on a machine replacement problem up to 5000 states which shows that a DRCCMDP problem with statistical distance based uncertainty sets can be solved very efficiently, while it takes more time in the case of moments based uncertainty sets.

Acknowledgement

This research was supported by the French government under the France 2030 program, reference ANR-11-IDEX-0003 within the OI H-Code.

A Proof of Lemma 2 -Case Hellinger distance

Given a policy f ∈ P O S and y ∈ R, let O = {p ∈ E p | V α (f, p) ≥ y}. From Table 1, the conjugate of ϕ has the following form

ϕ * (r) = r 1-r , if r < 1, ∞, if r ≥ 1. ( 39
) Let L = sup λ>0,β∈R β -λθ ϕ -λϕ * -1 + β λ P ν (O) -λϕ * β λ (1 -P ν (O)) . (40)
The constraint (i) of ( 5) is equivalent to

L ≥ 1 -ϵ. (41)
We consider two cases as follows:

Case 1: Let β λ < 1. Since λ > 0, the following inequality holds

β -1 λ < β λ < 1.
From (39), we have: ϕ * β λ = β λ-β , ϕ * β-1 λ = β-1 λ+1-β . Consequently, it follows from (40) that: L = sup λ>0,β<λ P ν (O) If u * ≤ 1 -P ν (O), we deduce that (1 -P ν (O), ∞) ⊂ (u * , ∞). Since G is decreasing on (u * , ∞), it implies that G is decreasing on (1 -P ν (O), ∞). Hence, the optimal value of G is attained when u = 1 -P ν (O), i.e, η = 0. From (42), L = 0 which violates the constraint (41). Therefore, u * > 1 -P ν (O) > 0. Since, G is decreasing on (u * , ∞) and increasing on (1 -P ν (O), u * ), then u = u * is the optimal solution of G(u) and L = -2 √ a 1 b 1 + c 1 . Therefore,

L =-2 (2 -θ ϕ ) 2 4 1 - (2 -θ ϕ ) 2 4 P ν (O)(1 -P ν (O)) + 1 -P ν (O) + (2 -θ ϕ ) 2 (2P ν (O) -1) 4 .

Fig. 1

 1 Fig.1Machine replacement MDP with two actions: "repair" (with solid lines) and "do not repair" (with dashed lines)

Fig. 2

 2 Fig. 2 CPU time (in seconds) vs number of states.

  β)(λ-β+1) -β 2 λ-β -λθ ϕ . Let η = λ -β. 2 -θ ϕ ) -η . Let g(λ, η) = λ 2 Pν (O) η(η+1) -1 η + λ(2 -θ ϕ ) -η.It is a second-order polynomial of λ and the coefficient of λ 2 is negative because 0 ≤ P ν (O) ≤ 1 and η > 0. It is well known that the maximum value of a second order polynomial f (x) = ax 2 + bx + c with a < 0 is c -b 2 4a and it holds at x = -b 2a . Hence, the maximum value of g(λ, η)holds at λ * = η(η+1)(2-θ ϕ ) 2(1+η-Pν (O)) . Since θ ϕ < 2, λ * > 0.Therefore, for a given η > 0, the optimal value L holds at λ * and L= c -b 2 4a , where c = -η, b = 2 -θ ϕ , a = Pν (O) u = η + 1 -P ν (O), then η > 0 is equivalent to u > 1 -P ν (O) and we can write θ ϕ ) 2 P ν (O)(P ν (O)where G(u) = a 1 u + b1 u + c 1 such that a 1 = (2 -θ ϕ ) 2 4 -1, b 1 = (2 -θ ϕ ) 2 P ν (O)(P ν (O) -1) 4 , c 1 = 1 -P ν (O) + (2 -θ ϕ ) 2 (2P ν (O) -1) 4 . Since 0 < θ ϕ < 2 and 0 ≤ P ν (O) ≤ 1, a 1 < 0 and b 1 ≤ 0. It is clear that G is decreasing on (u * , ∞),increasing on (-u * , u * ) and decreasing on (-∞, -u * ), where θ ϕ ) 2 4 -(2 -θ ϕ ) 2 P ν (O)(1 -P ν (O)), (43) G(u * ) = a 1 u * + b 1 u * + c 1 = -2 a 1 b 1 + c 1 .

Table 1

 1 List of selected ϕ-divergences with their conjugate

	Divergence	ϕ(t), t ≥ 0	ϕ * (r)
	Kullback-Leibler	t log(t) -t + 1.	e r -1
	Variation distance	|t -1|.

Known mean unknown covariance

δ 0 = 0.9

Unknown mean unknown covariance

We compute the optimal policies of the DRCCMDP problem (5) for each uncertainty Wasserstein (p,1-p) 1 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 2 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 3 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 4 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 5 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 6 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 7 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 8 (0, 1) (0, 1) (0, 1) (0, 1) (0, 1) 9 (0. set by solving corresponding equivalent mixed-integer optimization problem using the existing solvers mentioned in Section 2.3. The optimal policies corresponding to all the uncertainty sets are summarized in Table 5, where p is the probability of "repair" action and 1-p is the probability of "do not repair" action. It is clear that the optimal repair policy corresponding to all the uncertainty sets for first eight states is same. At Then, ( 41) is rewritten equivalently as follows

By taking the square on both side of (44), we get

By rewriting (45), we get the following second-order inequality in P ν (O)

and B, C, ∆ are given in Table 2. It is clear that ( 44) is equivalent to either P ν (O) ≥ x max or P ν (O) ≤ x min . Moreover, x max and x min are solutions of the following two equalities

which implies that x max is a solution of (46) and x min is a solution of (47). Hence, the condition P

which violates the constraint (44). Then, ( 44) is equivalent to P ν (O) ≥ x max , i.e., the constraint (i) of ( 5) is equivalent to: P ν (ρ T R ≥ y) ≥ -B+ √ ∆ 2 .

Case 2: Let 1 ≤ β λ . From (39), ϕ * β λ = ∞, which in turn implies that L = -∞ and it violates the constraint (41).


Felipe Bordeu

Pierre-Alain Boucard

Pierre Gosselet
email: gosselet@lmt.ens-cachan.fr

Balancing Domain Decomposition with Nonlinear Relocalization: Parallel Implementation for Laminates

Keywords: Composite Material, Partitioning, Parallelization, Domain decomposition, Newton-Krylov-Schur methods, Damage

Over the past thirty years, composite materials have been used increasingly in industry, especially in the aeronautical and spatial industries. Therefore, there is a great interest in the prediction of their degradation. The objective of this work is to develop a program capable of simulating structures on an industrial level using the latest models developed. In the first section, we present a review of the damage mesomodel used for the simulations. The second section presents the partition technique, and the parallel resolution technique utilized for the simulation. A new Newton loop is added to the classic algorithm to improve the performance in the case of localized nonlinearities. The third section present the implementation on a C++ home made finite element code. And the technique utilized to parallelized the problem. Finally, the forth section gives illustrations showing the level of performance which can be expected from such an approach.

Introduction

Today, the simulation of the degradation mechanisms in composite structures until final failure is still an industrial challenge, specially in the aeronautical and the spatial industry. This is mainly due to: (i) the complexity of the degradation mechanisms involved in the failure of composite structures and (ii) the size and the geometric complexity of the parts. Also structural phenomena like buckling can appear in these cases when stability is lost.

Degradation in composites is a complex phenomenon involving both damage and inelastic residual deformation. We observe fiber-matrix interface debonding, plasticity, progressive transverse cracking, brittle fracture of fibers (in plies) and delamination (between plies). These mechanisms are highly nonlinear, fiber-matrix debonding, transverse cracking and delamination present a strong unilateral feature which depends on whether the cracks are open or close. In order to take all these phenomena into account, a model proposed in [START_REF] Ladevèze | A damage computational method for composite structures[END_REF] was chosen. In our work an enhanced version of this model was used [START_REF] Lubineau | Construction fo a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard[END_REF]. This model is defined at the mesoscale characterized by the thickness of the plies [START_REF] Ladevèze | Sur la mécanique de l'endommagement des composites[END_REF]. Then, the laminated structure is described as a stacking sequence of homogeneous layers and interlaminar interfaces. The model of each mesoconstituent (i.e. the elementary layer and the interface) are introduced using the internal variable framework. Damage is quantified by mesodamage indicators, which can be connected directly to the loss of stiffness. An important point is that the state of damage is assumed to be uniform throughout the thickness of the elementary layer, but not throughout the whole thickness of the laminate.

On can observe that the scale of the model regarding to the scale of the composite structure are very different. Indeed, this difference results in a high computational cost in the case of a simulation of an industrial size structure. Today the increasing power of computers and the efficiency of parallel resolution algorithms offer the opportunity to perform simulation of very complex geometries.

Damage mesomodel

The first description of the mesomodel for laminates can be found in [START_REF] Ladevèze | Sur la mécanique de l'endommagement des composites[END_REF]. Two main assumptions lead to such a mesomodel. The first one is that the behavior of any laminated structure can be reconstructed starting from two elementary mesoconstituents: the elementary layer [START_REF] Ladevèze | Damage modelling of the elementary ply for laminated composites[END_REF] and the interface [START_REF] Allix | Damage analysis of interlaminar fracture specimens[END_REF] (fig. 1). The second assumption is that the damage state is uniform throughout the whole thickness of the elementary layer (but not throughout the thickness of the laminated composite).

In the earlier version of the model, the behavior law of each mesoconstituent was considered to be intrinsic, but recent works showed that at the mesoscale the behavior of an interface is coupled with the internal variables of the adjacent layers.

The elementary layer

The elementary layer is defined in classical energy form: 2e = σ t S layer σ [START_REF] Lubineau | Construction fo a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard[END_REF], where σ is the stress and S layer the compliance operator. This equation is written in a local reference frame defined by the fiber's direction (N 1 ) and the transverse direction (N 2 ). Operator S layer takes into account the damage state of the elementary layer [START_REF] Lubineau | Construction fo a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard[END_REF]. The damage indicators, which are constant throughout the thickness, are associated with the following degradation mechanisms:

• fiber breakage,

• diffuse intralaminar degradation,

• transverse microcracking For the sake of simplicity, we present the general form of the evolution laws (for more details, one can refer to [START_REF] Lubineau | Construction fo a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard[END_REF]). The damage evolution laws use classical thermodynamic forces. The calculation of these forces involves a "mean value" operator whose role is to calculate the mean value of the thermodynamic forces throughout the thickness of the elementary layer. This operator is very important because it ensures that the damage indicator remains constant throughout the thickness of the elementary layer. Then, the damage indicators for fiber breakage and diffuse intralaminar degradation are updated using the evolution laws. These evolution laws depend on the material and can be modeled as progressive or brittle. For a complete description of these functions see [START_REF] Ladevèze | A mesomodel for localisaton and damage computation in laminates[END_REF].

The case of transverse microcracking is a little different. Transverse microcracking is taken into account by using an equivalent microcracking rate ρ for the elementary layer [START_REF] Lubineau | Construction fo a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard[END_REF]. Then, with this microcracking rate and the evolution law, one can update the damage indicators for this phenomenon.

Other phenomena, like fiber hiperelasticity, matrix plasticity or fiber rotation can be easily taken into account.

The interface

The interface, which is introduced in order to model the debonding of two adjacent plies, can be interpreted as the thin layer of matrix between two plies and is a two-dimensional entity which ensures the transfer of stresses and displacements from one elementary layer to another.

The displacement jump between the upper and lower surfaces of the interface is expressed in the local reference frame defined by the orthotropic directions N 1 , N 2 and N 3 ; N 3 is normal to the interface, and N 1 and N 2 are the bisectors of the angle formed by the fibers of the upper and lower plies (fig. 2).

Figure 2: Orthotropic directions of the interface

As in the case of the elementary layer, the elastic strain energy of the interface (2e = σ t S interf ace σ ) is defined by a compliance operator S interf ace . This operator take into account the damage indicators as well as the unilateral behavior of the interface in mode I.

In the current version of the mesomodel, the interface's damage evolution laws are strongly coupled with the state variables of the adjacent layers [START_REF] Lubineau | Construction fo a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard[END_REF]. This coupling is achieved through the mean value of the microcracking rate ρ of the upper and lower adjacent layers [START_REF] Ladevèze | A damage computational method for composite structures[END_REF].

ρ = ρ + + ρ - 2 (1)
The mean values operator used for the layer and the coupled evolution law of the interface make the model highly nonlocal. Special care must be taken to ensure the correct implementation and integration of the evolution laws.

Parallel Nonlinear Solving Method

The previous section allow us to describe all the main mechanisms involved in the ruin of a composite structures. All this makes the problem highly nonlinear and non local, and also because the model is written at the mesoscale, a large number of d.o.f. is expected.

To treat these kind of problems, parallel solvers seems to be the best choice. The Balancing Domain Decomposition method (BDD) [START_REF] Mandel | Balancing domain decomposition[END_REF] was chosen to parallelized the problem. As we will further explain, this method rely on Schur non-overlapping decomposition of the structure and Krylov iterative solvers. This method is used as a linear solver inside a Newton-Raphson scheme to solve the nonlinear problem. Like every parallel solver a partition of the domain (the mesh) is required.

Preprocessing and partitioning

The discretization and partitioning of a composite structure prior to a parallel simulation with the mesomodel becomes a very delicate task. Because the behavior law is nonlocal and the integration of the evolution law is done at Gauss points, the mesh must be regular in the thickness of the composite (fig. 3). To achieve this, the mesh is created only with hexahedral and wedges elements, this ensure that the integration points are always aligned in the thickness of the mesh. This point add some difficulty in the meshing process but is not a real restriction because virtually all composite structures are made from a stack of single laminas.

It is known that the performance of parallel solvers are very dependent of the quality of the partitioning [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]. On one hand, one would seek partitions with the same number of elements in order to achieve a good balance of the load among the processors during the parallel resolution. On the other hand, one also wants to minimize the interaction zones between subdomains in order to minimize communications among processors. Finally, one must take into account the fact that the behavior law is nonlocal and requires information of the neighboring elements in order to be integrated correctly.

In the case of the use of a classical algorithm, two elements that are aligned in the thickness of the composite could end in different subdomains. In this case, one will be forced to transfer information across subdomains during the behavior law integration. One way to prevent this is to force the partition algorithm not to cut between element that are aligned in the thickness of the composite structure (fig. 4).

In order to restrict the way the partitioning is performed, we calculate the dual graph of the mesh (fig. 5.a and 5.b). Then this graph is reduced in such a way that elements that are aligned in the thickness of the composite are condensed into one single vertex. The reduce graph is then weighted to take into account the number of elements In our case we used Metis [START_REF] Kasypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF], an open source software capable of performing fast and high quality partitions for irregular graphs. After the partition, the graph is expanded (fig.

5.e) to finally calculate the partition of the mesh (fig. 5.f). At the same time that the partition is calculated, a file with the connectivity of the subdomains is created. This information will be used to optimize the reading of the meshes once the parallel simulation begins.

This technique enable us to perform fast, hight quality and automatic partitions of the mesh.

The resolution algorithm

As we said before, a BDD method was used to parallelize the resolution. The principle of the BDD is to condense the internal unknowns of each subdomain on the boundary, and then solve the condensed problem in parallel using an iterative method. First, for each substructure s we construct the Schur complement of the tangent operator and the condensed residue through the elimination of inner degrees of freedom ( 2) and (3). In equations ( 2) and (3) the i subscript indicates the internal unknowns and b the unknowns on the boundaries. These Schur complements and condensed residues computed for each subdomain can be assembled to construct a global linear problem over all boundary unknowns (4). Solving the global linear problem enables to compute the displacements on the boundaries. Finally a localization step is needed to determine the inner displacement in each subdomain [START_REF] Allix | Damage analysis of interlaminar fracture specimens[END_REF].

S s T = K s T bb -K s T bi K s -1 T ii K s T ib (2) b s = r s b -K s T bi K s -1 T ii r s i (3) S T ∆u b = b (4) ∆u s i = K s -1 T ii (r s i -K s T ib ∆u s b ) (5)
The equations ( 2), ( 3) and ( 5) are defined within subdomains, and only the equation ( 4) is defined over all the boundary unknowns. To avoid the heavy transfer (in the case of a direct resolution) the equation ( 4) is solved using a Krylov solver. The conjugate gradient, used to solve the global problem, only requires local matrix-vector products and vectors assemblies, avoiding the assemblage of the global system. More detail about this method can be found in [START_REF] Mandel | Balancing domain decomposition[END_REF].

Then, this linear solver is used inside a Newton loop which treats the non linear problem. Algorithm 1 summarized the structure of the Newton-Krylov-Schur method.

Algorithm 1: Classic Newton-Krylov-Shur algorithm repeat -Assemble local tangent operators K s T -Factorize K s T ii (eq. ( 2)) -Assemble local residue r s (eq. ( 3)) -Solve condensed global problem (eq. ( 4)) (Krylov solver) -Compute inner displacement in each subdomain (eq. ( 5)) until error< ε global In general the evolution of the nonlinearity in a composite structure is not uniform. Indeed, in most of the cases the nonlinearity is localized in a reduce zone (holes, crack tips, junctions), having sometime little or no influence on the rest of the structure. When the algorithm 1 is used, the convergence of the Newton is directly related to the convergence of the localized nonlinearity. To correctly treat this type of nonlinearity at the correct scale, an enhanced version of the NKS solver was proposed [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: Application to post-buckling analyses[END_REF]. In this version, a new Newton loop is introduced instead of the equation ( 5). This local Newton forces the subdomains to obey the nonlinear behavior and the equilibrium. Algorithm 2 summarizes the structure of the Newton-Krylov-Schur method with a nonlinear localization.

Algorithm 2: Newton-Krylov-Shur with nonlinear localization repeat -Assemble local tangent operators K s T -Factorize K s T ii (eq. ( 2)) -Assemble local residue r s (eq. ( 3)) -Solve condensed global problem (eq. ( 4))

repeat -Compute inner displacement ∆u s i = K s -1 T ii (r s i -K s T ib u s b ) -Update u i , K s T ii , K s T ib and r s i until error< ε local until error< ε global

Implementation

The strategy presented in the previous section was implemented in a home made 3D finite element software. The software was written in C++ using a preprocessor [START_REF] Leclerc | Platform metil: Optimizations and facilities linked to the generation of code[END_REF] written in Python that was in charge of generating the C++ code for elements and behavior laws. This makes the code very scalable for new functions, behavior laws or element types. Because the code generated is already optimized, the performance of the final code is guaranteed.

Most of todays supercomputer are SMP (Symmetric Multi-Processors) Clusters, this means that the cluster is composed of multi-core (or multi-processor) nodes. Figure 6 show a simple schema of a SMP Cluster. One can see that each node contains more than one core (or processor). All cores belonging to the same node share some resources, like RAM, local hard drive and network connections. Two level of synchronization and communication can be distinguished, the node level (between threads) and the cluster level (between nodes). To take advantage of this configuration, the idea is to run only one process per node, and make use of multithreading to take advantage of all the available cores.

Communication between threads is almost instantaneous (because they share the same RAM) only the synchronization must be taken care of. The synchronization at this level was ensured by the use of an open source library named ZThread. The second level of communication and synchronization was implemented using MPI.

Simulation Results

To verify the implementation of the domain decomposition method in parallel an elastic problem was treated. Figure 7.a shows the mesh and one of the partition (50 subdomains) used for the test. The different partitions were done with a classical algorithm, because the behavior law is elastic, there is no need to used the algorithm presented before. The boundary condition of the problem were blocked displacement in all directions in one side (face o in fig. 7.b), and an unitary displacement in the three direction on the other side (face p in fig. 7.b). The material is linear elastic, the piece has 150.000 degrees of freedom. Figure 7.b show the deformation and the strain on the piece. First, we solved the problem for a different number of subdomains. Figure 8 show that the number of iteration of the iterative method rest almost constant, this is impor- tant to ensure the extensibility of the method. The number of iteration may vary from decomposition to decomposition, due mainly to the partition algorithm [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]. Figure 8 show the memory used for the resolution.

Then, using a fix number of subdomains (50), the simulation was executed for a different number of nodes. Each cluster node has 2 cores at 2.4 GHz and 8 GB of RAM. In Figure 9 one can see the time spend and the speed up of the resolution in function of the number of node used.

The second case is holed carbon fiber reinforced laminate plate, figure 10 show the geometry and the partition used for the simulation, fiber orientation of the layers are [45 4 90 4 -45 4 0 4 ] s . The partitioning of the mesh was done with the proposed technique.

In this simulation all the degradation mechanisms of the damage model were active. The loading (0.075% of deformation) was apply in 6 steps. Table 1 show, the number of iteration executed for a classic Newton algorithm and for the algorithm with nonlinear localization.

Each global Newton computations of the classic algorithm implies one global resolution with the Krylov solver plus one local computation per subdomain. In the case of the algorithm with nonlinear localization the number of local resolution may varies from subdomain to subdomain, according to the local nonlinearities. Even though the

Conclusion

In the first section, we reviewed briefly one of the possible model for the simulation of composites and pinpointed the difficulties of implementing such a model. In the second section we proposed a technique for partitioning a mesh correctly so no additional effort would be necessary in a parallel calculation. Also we proposed a calculation strategy for the simulation of large composite structures in parallel. The strategy enable us to reduce the iteration number executed by global newton.

Third section showed the techniques used for implementation of the strategy. These techniques enable us to run the software in parallel in SMP cluster.

Finally, we showed some preliminary results of the partitioning tool and calculation strategy.

Figure 1 :

 1 Figure 1: The model of a laminate

Figure 3 :

 3 Figure 3: Regular (left) and unregular (right) mesh (grey area: Integration point aligned in the thickness)

Figure 4 :

 4 Figure 4: Poor partitioning (left), Good partitioning (right)

Figure 5 :

 5 Figure 5: (a) Original mesh and its graph, (c) dual graph of the mesh (the bold lines represent element aligned in the thickness), (c) weighted reduce graph, (d) partitioned weighted graph, (e) expanded graph, (f) partitioned final mesh

Figure 6 :

 6 Figure 6: Schema of a SMP Cluster Super-calculator

Figure 7 :

 7 Figure 7: The mesh (each color represent a subdomain) (left), deformation and strain (right)

Figure 8 :Figure 9 :

 89 Figure 8: Iterations of the Krylov solver vs subdomain number (left), memory usage vs subdomain number (right)

Figure 10 :Table 1 :

 101 Figure 10: Holed carbon fiber reinforced laminate plate geometry and partition


Natasha Bahrani
email: natasha.bahrani@ifsttar.fr

Eyal Levenberg

Juliette Blanc

Pierre Hornych

Inverse analysis of pavement layer moduli based on data collected by buried accelerometers and geophones

Keywords: Pavement analysis, Sensor technology, Back calculation of pavement layer moduli

The aim of this study is to develop a method for identifying asphalt pavement layer properties based on readings from accelerometers and geophones that are embedded near the ride surface. These sensors are relatively small in size and easily embeddable, making them an ideal choice for wide-area applications in the live transportation network. As a first step, a section within the IFSTTAR accelerated pavement testing (APT) facility was instrumented with accelerometers and geophones; also installed was an anchored displacement sensor to serve as a reference/validation device. The APT facility offers the ability to control the loading configuration and intensity, travel speed, and wander (i.e., lateral offset) position relative to the sensor locations. Thus, it becomes possible to isolate the task of property identification through inverse analysis from other real-world complications. The paper commences by describing the experimental setup, and presenting some raw sensor measurements during a single pass of the APT's wheel carriage. Then, assuming a layered-elastic model, a method is proposed and demonstrated for estimating the pavement moduli. The method is based on best-matching measured velocities and accelerations for the geophones and accelerometers (respectively), with the model predictions-without integrating the signals to convert them into deflections. Very good match is obtained for the sensor readings, and the inferred moduli closely agree with reference values. This outcome means that there is great potential in building a pavement condition monitoring system with nearsurface accelerometers and geophones.

Introduction

Embedded sensing technologies are being developed and evaluated in recent years for wide-area instrumentation of pavement systems [START_REF] Lajnef | Smart pavement monitoring system[END_REF][START_REF] Levenberg | Wireless pavement sensors for wide-area instrumentation[END_REF]. In general terms, the aim is to carryout mechanical condition monitoring as well as usage monitoring. The current study investigates the application of near-surface single-axis geophones and accelerometers to determine pavement layer properties. The considered geophones and accelerometers measure, respectively, the vertical velocity and vertical acceleration of the pavement ride surface resulting from a passing vehicle.

A classical approach for interpreting the sensor measurements consists in integrating the signal traces with respect to time (once for the geophones and twice for the accelerometers), in order to arrive at deflections. Then after, the deflection basin can then be used in to analyse the pavement behaviour and back-calculate pavement layer properties [START_REF] Clark | Calculating pavement deflections with velocity transducers[END_REF][START_REF] Liu | Measurement and evaluation on deterioration of asphalt pavements by geophones[END_REF]Arraigada et al. 2012;Ryynanen et al. 2012;[START_REF] Duong | Monitoring of pavements using the geophones[END_REF]. This is typically done using methods similar to those applied for the Falling Weight Deflectometer (FWD) [START_REF] Tholen | Comparison of FWD with other deflection testing devices[END_REF][START_REF] Nazarian | Determination of deflection of pavement systems using velocity transducers[END_REF]). In the work of [START_REF] Duong | Monitoring of pavements using the geophones[END_REF], geophone signals have also been used to identify vehicle types (number of axles) and travel speeds. An approach based on calculation of deflections from accelerometer and geophone signals is also included in Bahrani and al. (2020).

A different approach for interpreting the sensor measurements consists of matching directly the geophone or accelerometer measurements (or both) with the predictions of a pavement model [START_REF] Levenberg | Inferring pavement properties using an embedded accelerometer[END_REF][START_REF] Levenberg | Backcalculation with an implanted inertial sensor[END_REF]. This approach completely avoids the need for integration, and is therefore well suited for noisy signals. An iterative or optimization procedure can be used to guide the matching and determine the sought layer properties. In the current paper the method is applied to analyze both geophone and accelerometer measurements made on the IFSTTAR accelerated pavement test (APT) facility as a means of back-calculating the elastic moduli of the pavement layers. The quality of the results is checked against measurements made with an anchored displacement sensor, and also with backcalculated layer moduli from FWD measurements.

Experimental setup

Three types of sensors were used for the tests: one single-axis geophone, and two single-axis accelerometers. An anchored deflectometer (displacement transducer) was also installed in the pavement, to serve as reference, because it measures the displacements with a high accuracy; however this type of sensor is less appropriate for measurements on live road network, due to its high cost and difficult installation. The characteristics of the geophone and accelerometers are presented on The pavement section investigated herein was a 22 m long sector; the structure consisted of two asphalt concrete (AC) layers, with a total thickness of 11 cm and a 30 cm thick 0/20 mm granular subbase, overlaying a clayey-sandy soil subgrade extending to a depth of 3 m and resting on a 20 cm concrete slab which is supposedly fixed. The structure was built in September 2017. Laboratory complex modulus tests have been performed on the AC, characterizing the material across different temperature and frequencies (e.g., the norm of the complex modulus at 15 °C and 10 Hz was 11543 MPa). FWD tests and Benkelman beam measurements have also been performed, giving an average elastic modulus of 110 MPa for the subgrade, and 145 MPa for the granular layer.

The applied loads consisted of dual tires, with a contact surface of approximately 258 mm×220 mm (for each tire) and center-to-center spacing of 352 mm (see Figure 3). Tests were conducted under different load levels, between 45 KN and 65 KN, different loading speeds, between 1 rpm and 12 rpm (i.e. between 7 km/h and 87 km/h) and for 11 different lateral wheel positions.

Instrumentation

The section was instrumented in 2018 with the geophone and accelerometers, and also with an anchored vertical displacement sensor (based on LVDT S-series displacement sensor). The geophone and accelerometers were installed on the centerline of the pavement, one centimeter below the surface. They were glued with resin inside blind-holes drilled into the upper AC lift. They were oriented to measure vertical responses in the z direction with positive values in the downward direction. The effective measurement depth is from the mid height of the sensors to the top of the concrete slab (3 m deep). The geophone was ION LF-24, which uses the principle of a spring-mounted magnetic mass moving within a wire coil. The accelerometers were Silicon Design 22106-002 and MEMSIC CXL04GP1, which use MEMS capacitive technology.

Data Analysis methodology

Pavement response calculation methodology

The analysis of the deflection sensor, geophone and accelerometer measurements is based on the use of the multilayer linear elastic pavement calculation software ELLEA1(ver0.96), developed by [START_REF] Levenberg | Validation of NCAT Structural Test Track Experiment using INDOT APT Facility[END_REF] 1 . The approach consists of using the software to calculate the vertical displacement, vertical velocity and vertical acceleration at the position of each sensor, and matching the model predictions with the sensor measurements. For that purpose, an iterative procedure is used to adjust the model parameters (the pavement layer moduli), until the best match between the modelling results and the measurements is obtained. The final set of moduli corresponds to the estimated moduli of the pavement layers. This modelling methodology is successively applied to the deflection sensor, the geophone and the two accelerometers, and the pavement layer moduli obtained with each sensor will be compared.

Optimization of layer moduli

For the pavement modelling, the initial characteristics considered for the pavement layers (elastic modulus and thickness of each layer) are those given in part 2 except for AC modulus which was adjusted according to the real test temperature and frequency. The wheel loading is 65 kN and the dual wheel load distribution described in Fig. 3 is replicated for the simulation assuming circular contact areas. The Poisson ratio for all layers was set to 0.35, and the tire-pavement contact stress to 0.6 MPa. The vertical displacements are calculated over a horizontal distance of approximately 5m from the center of the load.

The optimization method, to determine the pavement layer moduli, is similar for all the sensors, and is based on the following steps:

1) Computing the pavement response (vertical displacement) with the modelling software for a passing wheel; 2) Calculating vertical velocities and accelerations by numerical derivation (finite difference formulas); 3) Comparing the computed response traces to the measured traces; 4) Finding the best match between the two response by adjusting the modulus parameters of the pavement;

The procedure allows obtaining a "best estimate" of the pavement layer moduli.

Only the modulus values of the pavement model are adjusted in the optimization process, the other parameters are fixed.

Comparison of measurements and model prediction and estimation of pavement layer moduli

The optimization procedure described above was applied to the measurements made in the APT experiment using the different sensors, for different loading conditions (different load levels, loading speeds and load positions). The data presented in this paper correspond only to measurements made at an average temperature of 18°C in the AC layers, and speeds of 16 m/s and 20 m/s (corresponding to AC layer moduli of 9038 MPa at 16m/s and 9441 MPa at 20 m/s). The sensors were placed in the center of the wheel path, between the dual wheels represented in Figure 3. For each load condition, measurements were made during several successive load applications, but for the modelling, only a single wheel load signal was considered. The duration of the signal was about 1 s to 1.5 s, corresponding to about 500 to 600 recorded measurement points.

A global error indicator is defined, R, characterizing the relative difference (in percent) between the model response and the reference values (measured sensor response). This indicator calculates the pointwise difference between the two signals, divided by the absolute maximum value of the measured response as defined in Eq. 1, where N is the total number of reference points.

(1) This difference R is then minimized by optimizing the modulus values of the model to obtain the best match of the calculated and measured values. The minimization is carried out using Solver which is the gradient optimizer in Excel (Fylstra, 1998). In the solver, maximum and minimum attainable modulus values are fixed for each layer, as represented in Table 2, to avoid getting unrealistic solutions. The choice for these bounds can be considered inconsequential, as long as they encapsulate the optimal solution. The different signals selected for the comparison with the ELLEA1 model predictions are presented in Figure 5. Each graph represent the response measured by one sensor (displacement sensor, geophone or accelerometer), under the passage of a moving wheel at 16 m/s or 20 m/s, and the corresponding predicted values, before and after optimization. Positive values represent downwards displacements, velocities and accelerations.

The displacement signals indicate a positive deflection when the wheel is approaching, and then a return to zero when the wheel drives away. The measured displacement signal is not completely symmetrical, most probably due to the viscoelastic behaviour of the AC layers (the return to zero is slower). The vertical velocity signals are positive when the wheel approaches, and then negative. The vertical acceleration signal is first positive, when the wheel is approaching, then negative under the wheel, and then positive again when wheel drives away. The acceleration signals are also slightly asymmetrical, and present more noise than the other signals, due to the low acceleration levels.

As can be seen in Figure 5 for all sensors, the modelling results are in good agreement with the measurements, except that the calculated signals are symmet-rical, due to the elastic model used. After optimization, the match with the measurement is only slightly improved. The calculated final differences between the modelling results and the measurements (error indicator R) are summarized in table 3. The average modulus values obtained after the optimization are given in tables 4 and 5. The back-calculated moduli of the pavement layers are generally close to the reference values obtained from laboratory tests and FWD measurements. However, in some cases, the calculated modulus is equal to the limit defined in the optimization procedure, which indicates some difficulties to converge. This occurs in particular for the asphalt layer modulus, for accelerometer CX (for the two speeds), and also for the geophone and deflectometer, for the speed of 16 m/s. This seems to indicate that the asphalt layer modulus has a relatively low influence on the calculated solution, for this pavement structure.

It is planned to perform other tests, on different pavement structures, to further evaluate and improve this back calculation procedure. However, these first results are very encouraging, and indicate that instrumentation with accelerometers and geophones can be used to determine pavement layer properties.

Conclusion

The overall approach of this work is to use embeddable sensors in the asphalt layer to measure the pavement response (vertical displacement, velocity and acceleration). For this purpose, an experimental pavement was tested on the IFSTTAR APT facility, instrumented with an anchored deflectometer, a geophone and 2 types of accelerometers, and submitted to 65 KN dual wheel loads, and different loading speeds. The data collected for speeds of 20m/s and 16m/s is considered in this work. The measured responses have been compared with calculations performed with a layer-elastic model, and a procedure for back calculation of pavement layer moduli has been proposed. The first results show a good match between the modelled and measured responses, and realistic back-calculated moduli. An error indicator R has been used to estimate the relative difference between the responses. The deflectometer gives the best match among the different sensors, however this type of sensor is not very suitable for real road monitoring, due to its high cost and difficulty of installation. The results obtained with the two accelerometers are very similar, with maximum relative differences of 2.78% and 3.1 % respectively. The geophone gives slightly less accurate results.

Following the APT test, it is planned to install the sensors on a real road section, to perform measurements under real heavy vehicle traffic, and also to improve the modelling approach, by taking into account a visco-elastic model for the asphalt materials.

Figure 1 .

 1 Figure 1. They were installed in an experimental pavement tested on the IFSTTAR Fatigue Carrousel, a circular outdoor Accelerated Pavement Testing facility (APT). The carrousel consists of a circular test track, of 40 m diameter, and a central loading system. A picture of the facility is shown in Figure 2. The test track is generally split into several sectors, each accommodating a different pavement structure. There are four identical loading arms at right angles, each housing a wheel carriage at the end. These wheel carriages can accommodate different tire arrangements, different load levels (up to 130 kN per arm) and different lateral wander positions. The mean travel radius of the loading carriage is 19 m, i.e. a circumference of about 120 m. The width of the test pavement is between 3 and 6 meters. The maximum loading speed is 100 km/h, and the typical loading rate is 500 000 carriage passes per month.

Fig. 1 Fig. 2

 12 Fig. 1 Geophones and Accelerometers selected for the study

Fig. 5

 5 Fig. 5 Signal comparison between measured, calculated and best match responses with respect to two speeds

Table 2

 2 Limits for the optimization of the back calculation

	Pavement Layer	Min Modulus (MPa)	Max Modulus (MPa)
	Bituminous layer	8000	12000
	Granular base	100	200
	Subgrade	50	150

Table 3

 3 Relative Differences between the measured and calculated responses

		Anchor Deflec-	Geophone	Accelerometer	Accelerometer
		tometer	Ion	CX	SD
	Difference 16m/s	0.144 mm	0.693mm/s	29.81mm/s²	29.513mm/s²
	Percentage %	2.7%	9.5%	3.1%	2.78%
	Difference 20m/s	0.009 mm	0.758 mm/s	32.693 mm/s²	33.53 mm/s²
	Percentage %	1.9%	8.4%	2.6%	2.36%

Table 4

 4 Modulus values obtained after the optimization compared to the reference, at the speed of 16m/s

	Pavement Layer Reference	Deflectometer	Geophone	Accelerometer	Accelerometer
		(MPa)	(MPa)	Ion (MPa)	CX (MPa)	SD (MPa)
	Bituminous layer	9038	8000	8000	12000	9967
	Granular base	145	155	154	100	100
	Subgrade	110	120	150	107	122

Table 5

 5 Modulus values obtained after the optimization compared to the reference, at the speed of 20m/s

	Pavement Layer	Reference	Deflectometer	Geophone	Accelerometer	Accelerometer
		(MPa)	(MPa)	Ion (MPa)	CX (MPa)	SD (MPa)
	Bituminous layer	9441	9500	8569	12000	9998
	Granular base	145	100	132	175	100
	Subgrade	110	137	150	67	101

https://orbit.dtu.dk/files/125076641/ELLEA1_ver0.96_.xlsb

Acknowledgements

"The research presented in this paper was carried out as part of the H2020-MSCA-ETN-2016. This project has received funding from the European Union's H2020 Program for research, technological development and demonstration under grant agreement number 721493"