Carlos Casorrán

Bernard Fortz

Martine Labbé

Fernando Ordóñez

A study of general and security Stackelberg game formulations

Keywords: Integer programming, discrete optimization, game theory, bilevel optimization

In this paper, we analyze different mathematical formulations for general Stackelberg games (GSGs) and Stackelberg security games (SSGs). We consider GSGs in which a single leader commits to a utility maximizing strategy knowing that one of p possible followers optimizes its own utility taking this leader strategy into account. SSGs are a type of GSG that arise in security applications where the strategies of the leader consist in protecting subsets of targets and the strategies of the p followers consist in attacking a single target. We compare existing mixed integer linear programming (MILP) formulations for GSGs, sorting them according to the tightness of their linear programming (LP) relaxations. We show that SSG formulations are projections of GSG formulations and exploit this link to derive a new SSG MILP formulation that i) has the tightest LP relaxation known among SSG MILP formulations and ii) its LP relaxation coincides with the convex hull of feasible solutions in the case of a single follower. We present computational experiments empirically comparing the difficulty of solving the formulations in the general and security settings. The new SSG MILP formulation is computationally efficient, in particular as the problem size increases.

Introduction

Stackelberg games model situations where players strive to optimize their individual objectives in a single sequential encounter. These models assume a player, referred to as the leader, can commit to a strategy that optimizes its utility function and then players that respond to the leader's decision, referred to as followers, take this decision into account when deciding how to optimize their own utility functions. Stackelberg games were introduced to model market competition [von Stackelberg, 2011] and have been used in diverse applications since, such as traffic equilibrium [START_REF] Krichene | Stackelberg routing on parallel networks with horizontal queues[END_REF], network toll setting [START_REF] Labbé | A bilevel model of taxation and its application to optimal highway pricing[END_REF], and security [START_REF] Brown | Defending critical infrastructure[END_REF], Jain et al., 2010].

In this work we consider normal form Stackelberg games with finite sets of actions for the leader and followers. We refer to these as general Stackeblerg games (GSG). The utility functions of GSGs are described by matrices, where each combination of actions for the leader and follower gives a reward value for each participant. Selecting a single action corresponds to a pure strategy, while a mixed strategy corresponds to a probability distribution over the set of actions for the player. Therefore, for GSGs the utility functions are bilinear functions of the players' mixed strategies.

Stackelberg games can be expressed as bilevel optimization problems, where the top level represents the leader's decision problem and includes the followers' responses as the optimal solution to the second level problem [START_REF] Colson | An overview of bilevel optimization[END_REF]. Mixed integer formulations of GSGs have been introduced thanks to the bilinear objective functions and the linearization of the second level optimality conditions with the use of integer variables [Bard, 1998]. The manner in which the bilinear objectives and second level problem optimality conditions are linearized give rise to the different mixed integer linear programming (MILP) formulations considered in this work. For instance, using big M constraints to linearize both the leader objective and the second level optimality conditions give rise to the (D2) formulation [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF]. The (DOBSS) formulation considers a single big M constraint but introduces new variables representing the product of the leader and follower strategies, [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF]. Finally, (MIP-p-G) is a formulation without big M constraints [Yin and Tambe, 2012]. Which of these MILP formulations of the bilevel stackelberg game problem is more convenient for computational efficiency is an underlying question of this work. When the leader in a GSG faces a single follower the problem can be solved in polynomial time, see [START_REF] Conitzer | Computing the optimal strategy to commit to[END_REF]. The same reference shows that if there are multiple followers then the problem is NP-hard. A solution for the multiple followers problem can be obtained by using the algorithm for the single follower instance on a Harsanyi transformation of the problem, [START_REF] Harsanyi | A Generalized Nash Solution for Two-Person Bargaining Games with Incomplete Information[END_REF], which combines the multiple adversaries into a single adversary with exponentially many actions. Solution methods based on mixed integer formulations of the multiple follower problem have been presented, for example, by [START_REF] Jain | Quality-bounded solutions for finite bayesian stackelberg games: Scaling up[END_REF] and [START_REF] Yang | Scaling-up security games with boundedly rational adversaries: A cutting-plane approach[END_REF].

Recent work has applied Stackelberg games in security settings where a leader has a limited budget to protect a set of targets while a follower aims to attack a single target. In this domain, the payoff matrices are structured with only two payoff values for every participant depending on whether or not the defender strategy protects the target attacked. We refer to problems that have this structure as Stackelberg security games (SSGs), which are introduced in detail in Section 2. Some SSG applications have included assigning Federal Air Marshals to transatlantic flights [START_REF] Jain | Software assistants for randomized patrol planning for the lax airport police and the federal air marshal service[END_REF], determining randomized port and waterways patrols for the U.S. Coast Guard [START_REF] Shieh | Protect: A deployed game theoretic system to protect the ports of the united states[END_REF], preventing fare evasion in public transport systems [Yin et al., 2012], and protecting endangered wildlife [START_REF] Yang | Adaptive resource allocation for wildlife protection against illegal poachers[END_REF]. The SSG models considered are closely related to the Interdiction games literature, [START_REF] Mcmasters | Optimal interdiction of a supply network[END_REF], specially when there is a fortification step. Such fortification-interdiction problems are multi-level optimization problems where a defender decides a limited fortification of a network, so that an interdictor (attacker) blocks a number of edges in the network and an operator tries to maximize flow or minimize a path over the network. If the optimal operation response can be subsumed in the interdictor's decision problem, then the problem has the structure of a Stackelberg security game. There are many variants and extensions of such fortification-interdiction games that allow multiple/sequential interdictions and problem specific formulations and algorithms, see reviews in [START_REF] Smith | Algorithms for network interdiction and fortification games[END_REF], Snyder et al., 2016, Fischetti et al., 2018]. However, to the best of our knowledge there is no polyhedral study of different mixed integer optimization formulations that arise due to the bilevel nature of the interaction between the defender and the attacker.

In this paper we focus on the polyhedral analysis of different mixed integer formulations for GSGs and SSGs. In particular we provide the following four key contributions. First, we provide an exhaustive comparative study of existing MILP formulations for Stackelberg games. Starting from the natural bilevel representation of Stackelberg games, we use wellknown integer programming techniques such as Fourier-Motzkin elimination [START_REF] Dantzig | Fourier-Motzkin Elimination and Its Dual[END_REF] and Reformulation Linearization Technique [START_REF] Sherali | A hierarchy of relaxations and convex hull characterizations for mixed-integer zero-one programming problems[END_REF] to derive known MILP formulations. Our study leads to a ranking of these MILP formulations in terms of the strength of their linear programming (LP) relaxations. Second, we explicit a formal link through projections of variables between the polyhedra of the LP relaxation of the GSGs formulations and those of SSGs. This allows to extend our study of GSG formulations to the security setting, leading to a comparison of SSG MILP formulations.

Third, we derive (SDOBSS q,y,s ) and (MIP-p-S q,y ), two new SSG MILP formulations. We show that (MIP-p-S q,y ) is the MILP formulation with the tightest linear relaxation among SSG formulations. We further show that if we restrict (MIP-p-S q,y ) to a single attacker type, its LP relaxation provides a complete linear description of the convex hull of its feasible solutions. Fourth, we provide computational experiments that compare solution times of the MILP formulations in both settings. Our experiments show that the formulations with the tightest LP relaxations have faster solution times as the problem size increases. In particular (MIP-p-S q,y ) scales better than competing formulations, being able to tackle larger-sized instances.

The remainder of this paper is organized as follows. In Section 2, we define general and security Stackelberg games. In Section 3, we derive GSG formulations from the literature. We provide theoretical results comparing the formulations presented. In Section 4, we describe and analyze computational experiments for the formulations in Section 3.

In Section 5, we present SSG formulations using projections, in the appropriate space of variables, of the formulations in Section 3, and derive (SDOBSS q,y,s ) and (MIP-p-S q,y ), new MILP formulations for SSGs. We then extend our theoretical comparisons of the general formulations to the security formulations. In Section 6, we describe and analyze the computational experiments for the security formulations. We conclude with some closing remarks in Section 7.

Notation and definition of the problem

In this section, we provide a formal definition of the two types of problems we study.

General Stackelberg games-GSGs

Let K be the set of p followers. We denote by I the set of leader pure strategies and by J the set of follower pure strategies. The leader has a known probability of facing follower k ∈ K, denoted by π k ∈ [0, 1]. We denote the n-dimensional simplex by S n = {a ∈ [0, 1] n : n h=1 a h = 1}. A mixed strategy for the leader consists in a vector x ∈ S |I| such that for i ∈ I, x i is the probability with which the leader plays pure strategy i. Analogously, a mixed strategy for a follower k ∈ K is a vector q k ∈ S |J| such that, q k j is the probability with which follower k replies with pure strategy j ∈ J. The rewards or payoffs for the leader and each follower, resulting from their choice of strategy, are encoded in a different matrix for each follower. These payoff matrices are denoted by (R k , C k ), where R k ∈ R |I|×|J| is the leader's reward matrix when facing follower k ∈ K and C k ∈ R |I|×|J| is the reward matrix for follower k. The expected reward of the leader and follower k, respectively, can be expressed as follows:

i∈I j∈J k∈K π k R k ij x i q k j , (1)
i∈I j∈J C k ij x i q k j , ∀k ∈ K. (2)
For all k ∈ K, we define the function B k : S |I| -→ S |J| as the function that, given the leader's mixed strategy x, returns a best response q k for each follower k. The solution concept used in these games is the Strong Stackelberg Equilibrium (SSE), introduced in [Leitman, 1978] and defined below.

Definition 1. A profile of mixed strategies (x, {B k (x)} k∈K ) form an SSE if:

1. The leader always plays a payoff-maximizing strategy:

x T R k B k (x) ≥ x T R k B k (x ) ∀x ∈ S |I| , ∀k ∈ K.
2. Each follower always plays a best-response,

B k (x) ∈ F k (x), where ∀k ∈ K, F k (x) = arg max q k {x T C k q k : q k ∈ S |J| }
is the set of best responses for each follower.

3. Each follower breaks ties optimally in favor of the leader:

x T R k B k (x) ≥ x T R k q k ∀q k ∈ F k (x).
An SSE assumes that the follower breaks ties in favor of the leader by choosing, when indifferent between different follower strategies, the strategy that maximizes the payoff of the leader. An SSE is in practice always achievable as the leader can always induce one by selecting a sub-optimal mixed strategy arbitrarily close to the equilibrium, causing the follower to prefer the desired strategy [von Stackelberg, 2011].

Proposition 1. For any leader strategy x and any k ∈ K, there is a best response to the k-th follower's problem that is given by a vector q k ∈ {0, 1} |J| such that j∈J q k j = 1.

Proof. Assume that B k (x) = qk ∈ {0, 1} |J| . We show that any canonical vector e jk such that qk j > 0, is also a best response vector, i.e., e jk ∈ F k (x) and x T R k e jk ≥ x T R k q k for all q k ∈ F k (x). Since qk = j∈J qk j e jk , with e jk ∈ S |J| , and

x T C k e jk ≤ x T C k qk for all j ∈ J, we have that x T C k qk = j∈J qk j (x T C k e jk ) ≤ j∈J qk j (x T C k qk ) = x T C k qk .
This implies that for any qk j > 0 we have x T C k e jk = x T C k qk , giving e jk ∈ F k (x). A similar argument shows that for any j such that qk j > 0 we have x T R k e jk = x T R k qk ; Hence, e jk is a best response vector.

This result shows that we can restrict the follower's best response only to pure strategies without influencing the SSE solution concept, as done in [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF].

In mathematical optimization, Stackelberg games are formulated as bilevel programming (BP) problems [START_REF] Bracken | Mathematical programs with optimization problems in the constraints[END_REF]. In BP the optimization problems have two levels where the top level problem considers some variables that are the optimal solution to another, second level, optimization problem. Important BP surveys are those by [Kolstad, 1985, Savard, 1989, Anandalingam and Friesz, 1992, Labbé and Violin, 2016]. In our setting, the first level problem corresponds to the leader's decision problem and the nested problem corresponds to the follower's decision problem. The following model, (BIL-p-G x,q ), is a bilevel program for the general Stackelberg game problem:

(BIL-p-G x,q ) Max x,q i∈I j∈J k∈K π k R k ij x i q k j (3) s.t.
x ∈ S |I| (4)

q k ∈ arg max r k    i∈I j∈J C k ij x i r k j    ∀k ∈ K, (5)
r k j ∈ {0, 1} ∀j ∈ J, ∀k ∈ K, ( 6
) j∈J r k j = 1 ∀k ∈ K. ( 7
)
The objective function maximizes the leader's expected reward. Condition (4) characterizes the mixed strategies considered by the leader. The second level problem defined by ( 5)-( 7) indicates that the follower maximizes its own payoff by giving a best response with a pure strategy to the leader's mixed strategy. Recall that such a pure strategy always exists as shown in Proposition 1. If there are multiple optimal strategies for the follower, the main level problem selects the one that benefits the objective of the leader.

Stackelberg security games-SSGs

In a Stackelberg security game (SSG) the defender allocates security resources to protect a subset of targets. Let J be the set of n targets that could be attacked and assume there are security resources to protect up to m < n of these targets. The set I of defender pure strategies is composed by all m i=1 n i subsets of at most m targets of J that the defender can protect simultaneously. With a slight abuse of notation, we refer to i ∈ I in this context as both the index running through the set of defender pure strategies I and as i ⊂ J the corresponding subset of J with at most m targets that are protected by security resources.

Similar to GSGs, the elements j ∈ J constitute the pure strategies of each attacker, which for SSG represents the single target attacked by the follower. In SSGs, payoffs for the players only depend on whether the target attacked is protected or not. This means that many of the strategies have identical payoffs. The authors in [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF] use this fact to construct a compact representation of the payoffs.

We denote by D k the utility of the defender when facing an attacker k ∈ K and by A k the utility of attacker k. Associated with each target and each player there are two payoffs depending on whether or not the target is protected, see Table 1. [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF] Protected Unprotected

Defender D k (j|p) D k (j|u) Attacker A k (j|p) A k (j|u)
Table 1: Payoff structure in an SSG when target j is attacked by an attacker k take advantage of the aforementioned compact representation to define a protection vector c whose components, c j , represent the frequency with which target j is protected. The components of the vector c satisfy

c j = i∈I:j∈i x i ∀j ∈ J, (8)
i.e., the frequency with which target j is protected is expressed as the sum of all probabilities of the strategies that protect that target. Variables q k j indicate whether an attacker k strikes a target j.

The defender's and attacker k's expected rewards, are, respectively:

j∈J k∈K π k q k j {c j D k (j|p) + (1 -c j )D k (j|u)}, ( 9
) j∈J q k j {c j A k (j|p) + (1 -c j )A k (j|u)}, ∀k ∈ K. ( 10
)
As with GSGs, such a game can be modeled by means of bilevel programming.

(BIL-p-S x,c,q ) Max j∈J k∈K

π k q k j {c j D k (j|p) + (1 -c j )D k (j|u)} s.t. ( 4
), ( 8),

q k ∈ arg max r k    j∈J r k j (c j A k (j|p) + (1 -c j )A k (j|u)    ∀k ∈ K, r k j ∈ {0, 1} ∀j ∈ J, ∀k ∈ K, j∈J r k j = 1 ∀k ∈ K.
The objective function maximizes the defender's expected reward. Constraints (4) and ( 8) characterize the exponentially many mixed strategies considered by the defender and relate them to the frequencies with which targets are protected. The remaining constraints constitute the second level optimization problem which ensures that the attacker maximizes its profit by attacking a single target that is the best response to the defender's selected strategy. Notice that a more compact formulation-one involving a polynomial number of variables and constraints-can be obtained if projecting out the exponentially many x variables does not lead to exponentially many constraints. This would give a polynomial size formulation involving only the c and the q variables. Given an optimal solution to this compact formulation-an optimal protection vector c and an optimal attack vector q-a probability vector x, solution to this game in extensive form, can be obtained by solving the system of linear inequalities defined by conditions ( 4) and ( 8). As this system involves n + 1 equalities, there exists a solution in which the number of variables x i with a positive value is not larger than n + 1, i.e., the output size of an SSG, under extensive form, is polynomial in the input size. See Section 5 for more details.

General Stackelberg games-GSGs

In Section 3.1, we present equivalent MILP formulations for the p follower GSG. In Section 3.2 we compare the polyhedra of the LP relaxations for the different formulations.

General

Stackelberg games: single level formulations [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF] tackle the problem of solving the bilevel formulation presented earlier, (BIL-p-G x,q ) by using a MILP reformulation. They replace the second level nested optimization problem, described by ( 5)-( 7), by the following set of constraints:

j∈J q k j = 1 ∀k ∈ K, (11)
q k j ∈ {0, 1} ∀j ∈ J, ∀k ∈ K, ( 12
) 0 ≤ (s k - i∈I C k ij x i ) ≤ (1 -q k j ) • M ∀j ∈ J, ∀k ∈ K, (13)
where s k ∈ R for all k ∈ K and M is an arbitrarily large positive constant. The two inequalities in constraints ( 13) ensure that q k j = 1 only for a pure strategy that maximizes the follower's payoff. The problem defined by ( 3)-( 4) and ( 11)-( 13) is referred to as (QUAD x,q,s ). It is possible to eliminate the nonlinearity in the objective function of (BILp-G x,q ) by adding additional variables that represent the product between x and q. To be more precise, use z k ij = x i q k j for all i ∈ I, j ∈ J and k ∈ K. This gives rise to formulation (DOBSS q,z,s ) introduced in [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF]:

(DOBSS q,z,s ) Max i∈I j∈J k∈K π k R k ij z k ij s.t. (11), (12), j∈J z k ij = j∈J z 1 ij ∀i ∈ I, ∀k ∈ K, ( 14
) i∈I z k ij = q k j ∀j ∈ J, ∀k ∈ K, (15)
z k ij ≥ 0 ∀i ∈ I, ∀j ∈ J, ∀k ∈ K, ( 16
) 0 ≤ s k - i∈I j ∈J C k ij z k ij ≤ (1 -q k j ) • M ∀j ∈ J, ∀k ∈ K, ( 17
) s ∈ R |K| .
Alternatively the quadratic term in the objective of (BIL-p-G x,q ) can be addressed by adding

|K| new variables and introducing a second family of constraints involving a big M constant.

This gives rise to formulation (D2 x,q,s,f ) below (a DOBSS variant with 2 big M constraints that has appeared in [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF]):

(D2 x,q,s,f ) Max k∈K π k f k (18) s.t.
(4), ( 11) -(13),

f k ≤ i∈I R k ij x i + (1 -q k j ) • M ∀j ∈ J, ∀k ∈ K, (19)
s, f ∈ R |K| ∀k ∈ K.
Additionally, we project the real variables s k in constraints ( 13) and ( 17) out by using Fourier-Motzkin elimination [START_REF] Dantzig | Fourier-Motzkin Elimination and Its Dual[END_REF]. This gives rise to constraints:

i∈I (C k ij -C k i )x i ≤ (1 -q k ) • M ∀j, ∈ J, ∀k ∈ K, ( 20
) i∈I j ∈J (C k ij -C k i )z k ij ≤ (1 -q k ) • M ∀j, ∈ J, ∀k ∈ K. (21)
Replacing ( 13) by ( 20) in (D2 x,q,s,f ) and ( 17) by ( 21) in (DOBSS q,z,s ) yields (D2 x,q,f ) and

(DOBSS q,z ). We analyze the behavior of these last two new formulations compared to that of (D2 x,q,s,f ) and (DOBSS q,z,s ) to see if removing variables s at the expense of adding constraints is worthwile.

Another equivalent MILP formulation for the p-follower GSG can be obtained by replacing constraints (17) with the following set of constraints:

i∈I (C k ij -C k i )z k ij ≥ 0 ∀j, ∈ J, ∀k ∈ K. ( 22
)
These constraints are derived by multiplying constraints (20) by q k , reorganizing and replacing the nonlinear terms x i q k j by z k ij . This leads to (MIP-p-G q,z ): 11), ( 12), ( 14) -( 16), ( 22).

(MIP-p-G q,z ) Max i∈I j∈J k∈K π k R k ij z k ij s.t. (
The linear relaxation of (MIP-p-G q,z ) appears in [Yin and Tambe, 2012]. The MILP formulation is a p-follower extension to the single follower formulation (MIP-1-G q,z ), due to [START_REF] Conitzer | Commitment to correlated strategies[END_REF]. Formal proofs that the formulations seen thus far are equivalent MILP formulations, i.e., that they are valid for the p-follower GSG, appear in [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF] for (DOBSS q,z,s ) and [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF] and [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF] for (D2 x,q,s,f ). These proofs show that each of them is equivalent to (QUAD x,q,s ). The equivalence of (DOBSS z,q ) and (D2 x,q,f ) is obtained from the Fourier-Motzkin elimination procedure [START_REF] Dantzig | Fourier-Motzkin Elimination and Its Dual[END_REF]. The equivalence proof for (MIP-p-G q,z ) is analogous to the proof used to show the equivalence for (DOBSS q,z,s ) and is omitted here.

[ [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF] state that the big M constants used are arbitrarily large. To be as computationally competitive as possible, we provide the tightest value for each big M constant in the formulations discussed thus far.

Proposition 2. The tightest values for the positive constants M are:

1. In (19), M = max i∈I {max ∈J R k i -R k ij } ∀j ∈ J, ∀k ∈ K.
2. In ( 13) and ( 17)

, M = max i∈I {max ∈J C k i -C k ij } ∀j ∈ J, ∀k ∈ K.
3. In ( 20) and ( 21)

, M = max i∈I {C k ij -C k i }, ∀j, ∈ J, ∀k ∈ K.

Comparison of the formulations

Given a formulation F, we denote by F its linear (continuous) relaxation and by P(F) the polyhedral feasible region of F. Further, let

Q = {(x, z) ∈ R n × R m : Ax + Bz ≤ d}.
Then the projection of Q into the x-space, denoted P roj x Q, is the polyhedron given by [START_REF] Pochet | Production Planning by Mixed Integer Programming[END_REF].

P roj x Q = {x ∈ R n : ∃z ∈ R m for which (x, z) ∈ Q}, see
First, we introduce an additional formulation which we denote by (DOBSS x,q,z,s,f ). This formulation is equivalent to (DOBBS q,z,s ), in the sense that the values of their LP relaxations coincide. In this formulation, we introduce variables f k for all k ∈ K to rewrite the objective function so that it matches the objective function of (D2 x,q,s,f ). We also add variables x i for all i ∈ I by rewriting ( 14) as j∈J z k ij = x i for all i ∈ I and all k ∈ K. Using this last condition, we can simplify ( 17) to ( 13). The formulation (DOBSS x,q,z,s,f ) is as follows.

(DOBSS x,q,z,s,f ) Max 11) -( 13), ( 15), ( 16),

k∈K π k f k s.t. (
f k = i∈I j∈J R k ij z k ij ∀k ∈ K, ( 23
) j∈J z k ij = x i ∀i ∈ I, ∀k ∈ K, ( 24
) s ∈ R |K| .
Further, note that from the Fourier Motzkin elimination procedure we have that P(D2 x,q,f ) = P roj x,q,f P(D2 x,q,s,f ) and, P(DOBSS q,z ) = P roj q,z P(DOBSS q,z,s ).

Proposition 3. P roj x,q,s,f P(DOBSS x,q,z,s,f ) ⊆ P(D2 x,q,s,f ). Further, there exist instances for which the inclusion is strict.

Proof. Note that all the constraints of P(D2 x,q,s,f ) can be found in the description of P(DOBSS x,q,z,s,f ) except for constraints (4) and ( 19). Constraints (4) are implied by constraints ( 11), ( 15), ( 16) and (24).

Further, the projection of P(DOBSS x,q,z,s,f ) on the (x, q, s, f )-space can be obtained by applying Farkas' Lemma [Farkas, 1902]. Constraints (15), ( 16), ( 23) and ( 24) are the only ones involving variables z k ij and are separable by k ∈ K. For a fixed k ∈ K the projection is given by:

A k = {(x, q, f ) : αf k + i∈I β i x i + j∈J γ j q k j ≥ 0 ∀(α, γ, β) : αR k ij + β i + γ j ≥ 0 ∀i ∈ I, ∀j ∈ J} (25) For a fixed j ∈ J, define α = -1, β i = R k ij for all i ∈ I, γ j = 0 and γ = max i∈I (R k i -R k ij
) for all ∈ J with = j. This definition of the parameters satisfies αR k ij + β i + γ j ≥ 0 for all i ∈ I, j ∈ J. Substituting these parameters in the generic constraints of A k yields

f k ≤ i∈I R k ij x i + ∈J: =j max i∈I (R k i -R k ij )q k ∀j ∈ J, ∀k ∈ K. ( 26
)
Constraints (26) imply constraints (19) for the tight value of M provided in Proposition 2 since for all j ∈ J and k ∈ K,

∈J: =j max i∈I (R k i -R k ij )q k ≤ max i∈I max ∈J R k i -R k ij ∈J: =j q k = max i∈I max ∈J R k i -R k ij (1-q k j ).
This proves the inclusion. To show that the inclusion may be strict, consider the following example where |I| = |J| = 3 and |K| = 1. Let the payoff matrix for the game be

(R, C) =      (1, 0) (0, 0) (0, 0) (0, 0) (1, 0) (0, 0) (0, 0) (0, 0) (0, 0)     
and consider the point defined by x = (1, 0, 0) t , q = ( 1 3 , 1 3 , 1 3 ) t , s = 10 and f = 2/3. Such a point is feasible for (D2 x,q,s,f ) but violates constraints (26) for j = 2 and is therefore infeasible for P roj x,q,s,f P(DOBSS x,q,z,s,f ).

Next, we compare the polyhedra P(MIP-p-G q,z ) and P roj q,z P(DOBSS q,z,s ).

Theorem 1. P(MIP-p-G q,z )⊆ P(DOBSS q,z )=P roj q,z P(DOBSS q,z,s ). Further, there exist instances for which the inclusion is strict.

Proof. The description of P(DOBSS q,z ) differs from that of P(MIP-p-G q,z ) by only one set of constraints: (21) must hold instead of ( 22). Hence, the remainder of the proof consists in showing that (21) are implied by ( 11), ( 14)-( 16), ( 22) and the nonnegativity of the q variables. The LHS of ( 21) can be rewritten as:

i∈I (C k ij -C k i )z k i + i∈I j ∈J:j = (C k ij -C k i )z k ij ≤ i∈I j ∈J:j = (C k ij -C k i )z k ij , using (22), ≤ max i∈I {C k ij -C k i } j ∈J:j = i∈I z k ij ≤ M j ∈J:j = q k j ,
given Proposition 2 and (15)

= M (1 -q k ), by (11).

To show that the inclusion may be strict consider the p-follower GSG between a leader and a fixed follower k ∈ K where the payoff bimatrix is:

(R k , C k ) =   (0, 1) (1, 0) (0, 0) (0, 0)  
The point with coordinates x = (1/2, 1/2) t , q k = (1/2, 1/2) t and

z k =   1/4 1/4 1/4 1/4  
has an objective value of 1/4 and is feasible in P(DOBSS q,z ). However it is not a feasible point in P(MIP-p-G q,z ) as it doesn't verify constraints (22) when j = 2 and = 1.

From an interpretation point of view, (MIP-p-G q,z ) can be seen as the result of applying Reformulation Linearization Technique (RLT) [START_REF] Sherali | A hierarchy of relaxations and convex hull characterizations for mixed-integer zero-one programming problems[END_REF] to (DOBSS q,z ).

Indeed, by multiplying both sides of constraints (20) by variable q k and noticing that q k (1 -q k ) = 0 since q is binary, one obtains i∈I (C k ij -C k i )x i q k ≤ 0 which, once linearized by introducing variables z k i , yields (22). For a given formulation F, we denote its optimal value by v(F) and the optimal value of its LP relaxation by v(F). Since (D2 x,q,s,f ) and (DOBSS x,q,s,f ) and (DOBSS q,z ) and (MIP-p-G q,z ) have the same objective function, the following corollary holds.

Corollary 1. v(MIP-p-G q,z ) ≤ v(DOBSS q,z ) = v(DOBSS x,q,s,f ) ≤ v(D2 x,q,s,f ).
Finally, when (MIP-p-G) is restricted to a single follower type, [START_REF] Conitzer | Commitment to correlated strategies[END_REF] showed that the integrality costraints are redundant, i.e., the remaining constraints in (MIP-1-G) provide a complete linear description of the convex hull of feasible solutions.

Computational experiments for GSGs

Here, we present computational experiments for the formulations in Section 3. The machine used for these experiments is an Intel Core i7-4930K CPU, 3.40GHz, equipped with 64 GB of RAM, 6 cores, 12 threads and running the Ubuntu operating system release 12.10 (kernel Linux 3.5.0-41-generic). The experiments were coded in the programming language Python and GUROBI version 6.5.1 was the optimization solver used with a 3 hour solution time limit.

The instances solved in the computational experiments are randomly generated. We consider two different ways of randomly generating the payoff matrices for the leader and the different follower types. First, we consider matrices where all the elements are randomly generated between 0 and 10 and second, we consider matrices where 90% of the values are between 0 and 10 but we allow for 10% of the data to deviate between 0 and 100. In the first case we say that there is no variability in the payoff matrices, in the sense that all the data is uniformly distributed, whereas in the second case, we refer to the payoff matrices as matrices with variability.

A general Stackelberg game instance is defined by three parameters: |I|, the number of leader pure strategies, |J|, the number of follower pure strategies and |K|, the number of follower types. For the purpose of these experiments, we have considered instances where |I| ∈ {10, 20, 30}, |J| ∈ {10, 20, 30} and |K| ∈ {2, 4, 6}. For each instance size, 5 instances are generated without variability in the payoff matrices and 5 are generated with variability.

In total, we consider 135 instances without variability and 135 instances with variability.

Performance profiles summarize our results, with respect to the following 4 measures: total running time employed to solve the integer problem, running time employed to solve the linear relaxation of the integer problem, total number of nodes explored in the branch and bound (B&B) tree and percentage optimality gap at the root node. The percentage optimality gap at the root node is calculated by comparing the optimal values of the formulation and of its LP relaxation: v(F)-v(F ) v(F )

• 100. A performance profile graph plots the total percentage of problems solved for each value of these measures.

We study the behavior of (D2 x,q,s,f ), (D2 x,q,f ), (DOBSS q,z,s ), (DOBSS q,z ) and (MIP-p-G q,z ).

Figures 1 and2 compare the performance profiles when the payoff matrices are generated without variability and with variability, respectively.

We observe that the instances where variability is introduced in the payoff matrices solve faster than those where no variability is considered. When there is no variability, (DOBSS q,z,s ) and (MIP-p-G q,z ) are the two most competitive formulations. (D2 x,q,s,f ) can also be solved efficiently for the mid-range instances but slows down for the more difficult instances. Introducing variability in the payoff matrices, however, leads to a dominance of (MIP-p-G q,z ) with (DOBSS q,z,s ) coming in a close second and (D2 x,q,s ) becoming noncompetitive for these instances. Regarding the time spent solving the linear relaxation of the problems, formulation (MIP-p-G q,z ) is the hardest to solve due to the fact that is has the most variables and constraints, O(|K||J| 2 ). On the other hand, (D2 x,q,s,f ), with O(|K||J|) variables and constraints, is the fastest. With respect to the number of nodes and gap percentage, our theoretical findings are corroborated: (MIP-p-G q,z ) is the tightest formulation and therefore uses the fewest nodes. This is even more the case when variability is Table 2 summarizes the mean percentage optimality gap at the root node obtained across the instances solved. Finally, note that the formulations obtained through Fourier-Motzkin, (D2 x,q,f ) and (DOBSS q,z ), explore slightly less nodes in the B&B tree than their counterparts, (D2 x,q,s,f ) and (DOBSS q,z,s ), but because of the increase in the number of constraints, the time to solve each linear relaxation increases. This increases the overall solution time of the Fourier-Motzkin formulations.

Stackelberg security games-SSGs

In this section, we derive three SSG formulations: (ERASER c,q,s,f ), due to [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF], and (SDOBSS q,y,s ) and (MIP-p-S q,y ). We derive these formulations by exploring the inherent link between the general setting, considered up to now and the security setting, defined in Section 2.2. In this setting, the defender pure strategies i ∈ I correspond to the different ways in which up to m targets can be protected simultaneously. With a slight abuse of notation, i ∈ I refers both to the index running through the set of pure strategies I and to the subset of at most m targets protected by pure strategy i ∈ I. Recall that the payoff matrices of SSGs satisfy:

R k ij =    D k (j|p) if j ∈ i D k (j|u) if j / ∈ i (27) C k ij =    A k (j|p) if j ∈ i A k (j|u) if j / ∈ i (28)
The payoff for the leader that commits to a pure strategy i ∈ I and a follower of type k ∈ K responds by selecting strategy j ∈ J is either a reward if pure strategy i ∈ I protects attacked target j ∈ J, or, a penalty if strategy i does not protect target j. The same argument explains the link between payoffs for the attackers.

Stackelberg security games: single level formulations

The first formulation we derive is based on (D2 x,q,s,f ). Consider (D2 c,x,q,s,f ), an extended description of (D2 x,q,s,f ) where we introduce the c variables through constraints (8) (see Section 2.2). We further use relations ( 27) and ( 28) to adapt the payoff structure:

(D2 c,x,q,s,f ) Max k∈K π k f k s.t. ( 4
), ( 8), ( 11), ( 12),

0 ≤ s k -A k (j|p)c j -A k (j|u)(1 -c j ) ≤ (1 -q k j ) • M ∀j ∈ J, ∀k ∈ K, (29)
f k ≤ D k (j|p)c j + D k (j|u)(1 -c j ) + (1 -q k j ) • M ∀j ∈ J, ∀k ∈ K, (30) s, f ∈ R K .
This extended formulation is equivalent to (D2 x,q,s,f ), because, even though they are defined in different spaces of variables, the value of their LP relaxations coincide.

The formulation above has a large number of non-negative variables since in the security setting, the set I of all defender pure strategies is exponential in the number of targets as it contains all subsets of at most m targets of J that the defender can protect simultaneously.

In order to avoid having exponentially many non-negative variables in our formulation, we project out variables x i , i ∈ I, from the formulation. Note that only constraints (4) and ( 8) involve said variables.

Proposition 4. Consider the following two sets:

A = P roj c (x, c) ∈ R |I| × R |J| : (4), (8) B =    c ∈ R |J| : j∈J c j ≤ m, c j ∈ [0, 1] ∀j ∈ J    Then, A = B.
Proof. Observe first that using Farkas' Lemma [Farkas, 1902]:

A =    c ∈ R |J| : j∈J α j c j + α |J|+1 ≥ 0 ∀α ∈ R |J|+1 : j∈J:j∈i α j + α |J|+1 ≥ 0 ∀i ∈ I : |i| ≤ m and α |J|+1 ≥ 0    , Thus A ⊆ B. Indeed, the following 2|J| + 1 vectors in R |J|+1 : ∀j ∈ J, e j ∈ R |J|+1 : e j j = 1, e j k = 0 ∀k ∈ J : k = j and e j |J|+1 = 0, ∀j ∈ J, f j ∈ R |J|+1 : f j j = -1, f j k = 0 ∀k ∈ J : k = j and f j |J|+1 = 1 and g ∈ R |J|+1 : g j = -1 ∀j ∈ J and g |J|+1 = m,
satisfy j∈J:j∈i α j + α |J|+1 ≥ 0 and α |J|+1 ≥ 0. Additionally, when we substitute the above vectors into the generic constraints defining A, they yield all the constraints defining B.

To show that A = B, it remains to show that any other inequality

j∈J α j c j + α |J|+1 ≥ 0 (31)
such that α satisfies j∈J:j∈i

α j + α |J|+1 ≥ 0 ∀i ∈ I : |i| ≤ m and α |J|+1 ≥ 0, ( 32
)
is dominated by some nonnegative linear combination of the constraints defining B.

First, note that we can restrict our attention to constraints such that α j ≤ 0 for all j ∈ J.

If there exists ĵ ∈ J such that α ĵ > 0, since α must satisfy (32) and |i \ { ĵ}| ≤ |i| ≤ m, it follows that ᾱ with ᾱĵ = 0 and ᾱj = α j for all j ∈ J \ { ĵ} also satisfies (32) and since c ≥ 0, we have that

j∈J ᾱj c j + ᾱ|J|+1 ≤ j∈J α j c j + α |J|+1 .
Therefore, the constraint defined by α is dominated by the constraint defined by ᾱ. We thus distinguish two cases of α satisfying (32):

Case 1. |{j : α j < 0}| = k ≤ m, and Case 2. |{j :

α j < 0}| = k > m.
In Case 1, by considering a linear combination of inequalities c j ≤ 1 for 1 ≤ j ≤ k with respective weights -α j ≥ 0, we obtain that:

0 ≤ k j=1 α j c j - k j=1 α j ≤ j∈J α j c j + α |J|+1 ,
since α j = 0 for all j > k and α satisfies (32) for i = {1, . . . , k}.

For Case 2, assume w.l.o.g that α 1 ≤ α 2 ≤ . . . ≤ α k < 0 and α j = 0 for all j > k.

Then, build a linear combination of inequality j∈J c j ≤ m with weight -α m ≥ 0 and inequalities c j ≤ 1 for 1 ≤ j ≤ m with respective weights α m -α j ≥ 0. The valid inequality thus obtained is:

0 ≤ m j=1 α j c j + j>m α m c j - m j=1 α j ≤ j∈J α j c j - m j=1 α j , since α j ≥ α m for all j > m ≤ j∈J α j c j + α |J|+1 ,
since α satisfies (32) for i = {1, . . . , m}.

Proposition 4 leads to the following formulation based on (D2 c,x,q,s,f ): 11), ( 12), ( 29), (30),

(ERASER c,q,s,f ) Max k∈K π k f k s.t. (
j∈J c j ≤ m, 0 ≤ c j ≤ 1 ∀j ∈ J, s, f ∈ R K .
The above formulation involves a polynomial number of variables and constraints and was presented in [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF]. The next result is also an immediate consequence of Proposition 4.

Corollary 2. P roj c,q,s,f P(D2 c,x,q,s,f ) = P(ERASER c,q,s,f ).

We now derive new SSG formulations based on (DOBSS q,z,s ) and (MIP-p-G q,z ). We first present extended descriptions of both formulations by considering y k j variables satisfying:

y k j = i∈I: ∈i z k ij ∀j, ∈ J, ∀k ∈ K. ( 33
)
We use ( 27) and ( 28) to adapt the payoffs to the security setting leading to:

(DOBSS q,z,y,s )

Max 11), ( 12), ( 14) -( 16), ( 33), 11), ( 12), ( 14) -( 16), ( 33),

j∈J k∈K {π k (D k (j|p)y k jj + D k (j|u)(q k j -y k jj ))} (34) s.t. (
0 ≤ s k -A k (j|p) j ∈J y k jj - A k (j|u)(1 - j ∈J y k jj ) ≤ (1 -q k j ) • M ∀j ∈ J, ∀k ∈ K, ( 35
) s ∈ R |K| . ( 36
) (MIP-p-G q,z,y ) Max j∈J k∈K π k (D k (j|p)y k jj + D k (j|u)(q k j -y k jj )) s.t. (
A k (j|p)y k jj + A k (j|u)(q k j -y k jj )- A k ( |p)y k j -A k ( |u)(q k j -y k j ) ≥ 0 ∀j, ∈ J, ∀k ∈ K. (37)
Further, consider the following constraints:

j∈J y k j = j∈J y 1 j ∀ ∈ J, ∀k ∈ K, ( 38
)
and let us define the following polyhedra C and D:

C := {(q, z, y, s) ∈ [0, 1] |K||J| × [0, 1] |K||I||J| × [0, 1] |K||J| 2 × R |K| :
(11), ( 15), ( 16), ( 33), ( 35), ( 36), (38)} 15), ( 16), ( 33), ( 37), (38)} Lemma 1. C ⊇ P(DOBSS q,z,y,s ) and D ⊇ P(MIP-p-G q,z,y )

D := {(q, z, y) ∈ [0, 1] |K||J| × [0, 1] |K||I||J| × [0, 1] |K||J| 2 : (11), (
Proof. Consider constraints (14) and sum over all i ∈ I such that ∈ i:

i∈I: ∈i j∈J z k ij = i∈I: ∈i j∈J z 1 ij ∀ ∈ J, ∀k ∈ K. (39)
Applying ( 33) to (39) yields (38) and the result follows.

We now project the z variables from the larger polyhedra C and D. Said variables only appear in constraints ( 15), ( 16) and (33).

Lemma 2. Consider the following two sets; X = P roj q,y (q, z, y) ∈ R |K||J| 2 +|K||J|+|I||J||K| : (15), ( 16), ( 33)

Y = {(q, y) ∈ R |K||J| 2 +|K||J| : ∈J y k j ≤ mq k j ∀j ∈ J, ∀k ∈ K, 0 ≤ y k j ≤ q k j ∀j, ∈ J, ∀k ∈ K} Then, X = Y.
Proof. Note that constraints (15), ( 16) and ( 33) can be treated independently for each k ∈ K and each j ∈ J. First consider the case where q k ĵ = 0 for ĵ ∈ J and k ∈ K. Constraints (15) then imply that for all i ∈ I, z k i ĵ = 0 and constraints (33) force y k ĵ = 0 for all ∈ J and the result holds. For all j ∈ J, k ∈ K such that q k j = 0, consider x i = z k ij /q k j and c = y k j /q k j and apply Propostion 4. The result follows.

Consider P roj q,y,s C and P roj q,y D as the feasible regions of the linear relaxations of two MILP formulations-(SDOBSS q,y,s ) and (MIP-p-S q,y )-where we maximize the objective function (34) under the additional requirement that the q variables be binary. Hence, we present (SDOBSS q,y,s ), a security formulation based on (DOBSS q,z,y,s ), (SDOBSS q,y,s )

Max

j∈J k∈K π k (D k (j|p)y k jj + D k (j|u)(q k j -y k jj )) s.t.
(11), ( 12), ( 35), ( 38)

∈J y k j ≤ mq k j ∀j ∈ J, ∀k ∈ K, (40) 0
≤ y k j ≤ q k j ∀j, ∈ J, ∀k ∈ K, (41)
s ∈ R |K| .
And we also present (MIP-p-S q,y ), a security formulation based on (MIP-p-G q,z,y ), (MIP-p-S q,y ) Max

j∈J k∈K π k (D k (j|p)y k jj + D k (j|u)(q k j -y k jj )) s.t. ( 11
), ( 12), ( 37), ( 38), ( 40), (41)

The following corollaries are an immediate consequence of Lemmas 1 and 2.

Corollary 3. P roj q,y,s P(DOBSS q,z,y,s ) ⊆ P(SDOBSS q,y,s ).

Corollary 4. P roj q,y P(MIP-p-G q,z,y ) ⊆ P(MIP-p-S q,y ).

In addition, note that if we restrict (MIP-p-G q,z,y ) to a single type of follower, constraints (14) disappear and one thus obtains the following corollary.

Corollary 5. P roj q,y P(MIP-1-G q,z,y ) =P(MIP-1-S q,y )

The above corollary immediately leads to the following theorem.

Theorem 2. (MIP-1-S q,y ) is a linear description of the convex hull of feasible solutions for the Stackelberg security game with a single type of attacker.

Proof. The result follows from Corollary 5 and from [START_REF] Conitzer | Commitment to correlated strategies[END_REF] showing that (MIP-1-G q,z ) is a linear description for general games.

As in general games, we use Fourier-Motzkin elimination on constraints ( 29) and ( 35) to project out the s variables from formulations (ERASER c,q,s,f ) and (SDOBSS q,y,s ) respectively. This leads to the following two families of inequalities:

(A k (j|p)-A k (j|u))c j +(A k ( |u)-A k ( |p))c +A k (j|u)-A k ( |u) ≤ (1-q k )•M ∀j, ∈ J, ∀k ∈ K, (42)
(A k (j|p) -A k (j|u)) h∈J y k jh + (A k ( |u) -A k ( |p)) h∈J y k h + A k (j|u) -A k ( |u) ≤ (1 -q k ) • M ∀j, ∈ J, ∀k ∈ K, (43)
Replacing constraints ( 29) by ( 42) in (ERASER c,q,s,f ) and ( 35) by ( 43) in (SDOBSS q,y,s ) leads to (ERASER c,q,f ) and (SDOBSS q,y ).

In the same spirit as Proposition 2, we present the following proposition, establishing the tightest values for the big M constants in the formulations seen so far:

Proposition 5. The tightest values for the positive constants M are:

1. In (30), M = max ∈J {D k ( |p), D k ( |u)} -min{D k (j|p), D k (j|u)}, ∀j ∈ J, k ∈ K. 2. In (29), ( 35
), M = max ∈J {A k ( |p), A k ( |u)}-min{A k (j|p), A k (j|u)}, ∀j ∈ J, k ∈ K. 3. In (42), ( 43
), M = max{A k (j|p), A k (j|u)} -min{A k ( |p), A k ( |u)}, ∀j, ∈ J, k ∈ K.

Comparison of the formulations

First, we introduce an additional formulation which we denote by (SDOBSS c,q,y,s,f ). This formulation is equivalent to (SDOBSS q,y,s ), in the sense that the value of their LP relaxations coincide. In this formulation, we introduce variables f k for all k ∈ K to rewrite the objective function so that it matches the objective function of (ERASER c,q,s,f ). We also add variables c for all ∈ J and rewrite constraints (38) as j∈J y k j = c for all ∈ J and all k ∈ K. Using this last condition we can simplify ( 35) to ( 29). The formulation (SDOBSS c,q,y,s,f ) is as follows.

(SDOBSS c,q,y,s,f ) Max 11), ( 12), ( 29), ( 40), (41),

k∈K π k f k s.t. (
f k = j∈J {y k jj (D k (j|p) -D k (j|u))+ q k j D k (j|u)} ∀k ∈ K (44) j∈J y k j = c ∀ ∈ J, ∀k ∈ K, (45) s ∈ R |K| .
Note that P(ERASER c,q,f ) = P roj c,q,f P(ERASER c,q,s,f ) and P(SDOBSS q,y ) = P roj q,y P(SDOBSS q,y,s ).

Proposition 6. P roj c,q,s,f P(SDOBSS c,q,y,s,f ) ⊆ P(ERASER c,q,s,f ). Further, there exist instances for which the inclusion is strict.

Proof. The projection of P(SDOBSS c,q,y,s,f ) onto the (c, q, s, f )-space is obtained by applying Farkas' Lemma. Constraints ( 40)-( 41) and ( 44)-( 45) are the only ones involving variables y k j and are separable by k ∈ K. For a fixed k ∈ K, the projection is given by:

A k = {(c, q, f ) : α(f k - j∈J D k (j|u)q k j ) + ∈J β c + m j∈J γ j q k j + j∈J ∈J δ j q k j ≥ 0 ∀(α, β, γ, δ) : γ, δ ≥ 0, β + γ j + δ ,j ≥ 0 ∀ , j ∈ J : = j, and
α(D k (j|c) -D k (j|u)) + β j + γ j + δ j ≥ 0 ∀j ∈ J} (46)
Consider, for each k ∈ K, the following set B k :

B k = {(c, q, f ) : c ≤ j∈J q k j , ∀ ∈ J, ( 47
) c ≥ 0, ∀ ∈ J, ( 48
) ∈J c ≤ m j∈J q k j , (49)
f k ≤ c j (D k (j|p) -D k (j|u))+ ∈J: =j q k D k ( |p) + q k j D k (j|u) ∀j ∈ J, (50)
q k j ≥ 0 ∀j ∈ J, ∀k ∈ K.}
Let us see that A k ⊆ B k for all k ∈ K. First note that if we set α = 0, the following definitions of the parameters β, γ and δ comply with the conditions in (46):

β = e h , γ = {0} j∈J , δ = {0} ,j∈J , ∀h ∈ J, β = -e , γ = {0} j∈J , δ = {1} j∈J , ∀ ∈ J, β = {-1} ∈J , γ = {1} j∈J , δ = {0} ,j∈J , β = {0} ∈J , γ = {0} j∈J , δ 1 = {e j }, ∀j ∈ J.
Substituting these valid parameters in the generic constraints in A k , produces all of the constraints in B k except (50). Further, for a fixed j ∈ J, consider α = -1, β = 0 and

γ = 1 m (D k ( |p) -D k ( |u)
) for all ∈ J such that = j, β j = D k (j|p) -D k (j|u) and γ j = 0. Finally, set δ j = 0 for all , j ∈ J. This definition of parameters is valid as it satisfies the conditions in (46). Substituting in the generic constraints in A k yields (50).

It remains to show that for all k ∈ K, constraints (50) imply (30) for the tight value of M shown in Proposition 5. The implication holds because

∈J: =j q k D k ( |p) ≤ max ∈J {D k ( |p)} ∈J: =j q k = (1 -q k j ) max ∈J {D k ( |p)} ∀j ∈ J, ∀k ∈ K.
Hence, P roj c,q,s,f P(SDOBSS c,q,y,s,f ) ⊆ P(ERASER c,q,s,f ). To show that the inclusion may . Consider the point defined by q = ( 1 3 , 1 3 , 1 3 ) t , c = (1, 0, 0) t , s = 10 and f = 2/3. Such a point is feasible for (ERASER c,q,s,f ) but violates constraints (50) when j = 2 and is therefore infeasible for P roj c,q,f,s P(SDOBSS c,q,y,s,f ).

Based on Theorem 1 we can present the following theorem comparing the polyhedra P(MIP-p-S q,y ) and P roj q,y P(SDOBSS q,y,s ): Theorem 3. P(MIP-p-S q,y ) ⊆ P(SDOBSS q,y ) = P roj q,y P(SDOBSS q,y,s ).

Proof. The inclusion is a consequence of Theorem 1, the relations between the payoffs described in ( 27) and ( 28) and the relation between the z and y variables described in (33).

To show that the inclusion may be strict, consider the following game. We set m = 2, |J| = 2 and |K| = 1. The reward and penalty payoff matrices for both the defender and the attacker are given by D(

•|p) = [1, 0], D(•|u) = [0, 0], A(•|p) = [0, 0] and A(•|u) = [0, 1].
Additionally, the point with coordinates

c t = (1/2, 1/2), q t = (1/2, 1/2) and y k =   1/4 1/4 1/4 1/4  
has an objective value of 1/4 and is a valid feasible solution of P(SDOBSS q,y ). However, it is not feasible in P(MIP-p-S q,y ) as it does not verify constraints (37) when j = 1 and = 2.

Observe that (MIP-p-S q,y ) can be obtained by applying RLT [START_REF] Sherali | A hierarchy of relaxations and convex hull characterizations for mixed-integer zero-one programming problems[END_REF] to (SDOBSS q,y ). Multiplying both sides of constraints (42) by variable q k and noticing that q k (1 -q k ) = 0, since q k is binary, one obtains constraints that once linearized, by introducing variables y k j , yield (37). Since (ERASER c,q,s,f ) and (f -SDOBSS c,q,s,f ) and (SDOBSS q,y ) and (MIP-p-S q,y ) have the same objective function, the following corollary holds. Corollary 6. v(MIP-p-S q,y ) ≤ v(SDOBSS q,y ) = v(SDOBSS c,q,s,f ) ≤ v(ERASER c,q,s,f ).

Computational experiments for SSGs

Our security experiments are run on randomly generated instances. For each instance, four payoff matrices have to be generated that satisfy

D k (•|p) ≥ D k (•|u) and A k (•|u) ≥ A k (•|p).
We consider two ways of generating these matrices. First, we generate matrices where the values for the penalty matrices (D k (•|u) and A k (•|p)) are randomly generated between 0 and 5 and all values for the reward matrices (D k (•|p) and A k (•|u)) are randomly generated between 5 and 10. We refer to these as matrices with no variability. Second, we consider an alternative where 90% of the values for the penalty matrices are randomly generated between 0 and 5 (between 5 and 10 for the reward matrices) and 10% of the values for the penalty matrices are randomly generated between 0 and 50 (between 50 and 100 for the reward matrices). We refer to these as matrices with variability. We use a solution limit of 3 hours.

A Stackelberg security game instance is defined by |J|, the number of targets, |K| the number of attacker types and m, the number of security resources available to the defender.

Recall from the computational experiments for GSGs that using payoff matrices with variability, amounts to endowing the game with more structure, thus making it somewhat easier to solve. We have encountered the same phenomenon in SSGs. For games whose payoff matrices have variability, we have considered J = {30, 40, 50, 60, 70}, K = {6, 8, 10, 12} and we have allowed m to be either 25%, 50% or 75% of the number of targets. For games whose payoff matrices don't have variability we have had to be less ambitious in order to solve all instances to optimality within the stipulated time limit and have considered J = {10, 20, 30, 40, 50}, K = {2, 4, 6, 8} while still considering m to be either 25%, 50% or 75% of the number of targets. In either case, for each instance size, we generate 5 random instances as described above. In total, we consider 300 randomly generated instances.

We study the behavior of (ERASER c,q,s,f ), (SDOBSS q,y,s ) and (MIP-p-S q,y ). For the sake of clarity, we no longer consider the Fourier-Motzkin formulations (ERASER c,q,f ) and (SDOBSS q,y ). Performance-wise, (ERASER c,q,s,f ) and (SDOBSS q,y,s ) compare to their Fourier-Motzkin formulations in a similar way to how (D2 x,q,s,f ) and (DOBSS q,z,s ) compared to theirs in Section 4 (results not shown). We plot performance profile graphs in Figures 3 and4. Note that for the experiments with variability, (ERASER c,q,s,f ) is the fastest formulation for most of the instances. However, we see that for the more difficult instances, its solution time increases significantly, eventually surpassing the solution time of (MIP-p-S q,y ). This indicates that for these instances (ERASER c,q,s,f ) ceases to be competitive and (MIP-p-S q,y ) is the formulation that solves the fastest. As for the instances whose payoff matrices have no variability, and are thus harder to solve, we observe that (ERASER c,q,s,f ) outperforms the running time of the other two formulations for 80% of the instances. However, for the most difficult instances, (MIP-p-S q,y ) is faster than the other two formulations. For the last 5% of the instances, (ERASER c,q,s,f ) is the worst formulation. In terms of size of the formulations, (ERASER c,q,s,f ) is the formulation with the least number of constraints and variables: O(|J||K|). Observe that (MIP-p-S q,y ) and (SDOBSS q,y,s ) have O(|J| 2 |K|) constraints and variables. Thus, these formulations have larger LP relaxations and thus take longer time to solve than (ERASER c,q,s,f ) does. However, Figures 3 and4 confirm our theoretical findings: (MIP-p-S q,y ) has the tightest LP relaxation and this Based on our results, we observe a trend that indicates that for difficult instances, particularly in the case of payoff matrices with no variability, one could expect (ERASER c,q,s,f ) and (SDOBSS q,y,s ) to perform very poorly compared to (MIP-p-S q,y ). To analyze this, we consider instances where the payoff matrices have no variability and where K = {6, 8, 10, 12}, J = {30, 40, 50, 60, 70} and m is 25%, 50% and 75% of the targets. We generate 5 random instances for each size. In addition, for practical reasons, we consider a time limit of 30 minutes. The computational results for these instances are shown in Figure 5. Note that (MIP-p-S q,y ) is able to solve 95% of the 300 instances within the stipulated time limit, outperforming (SDOBSS q,y,s ) and (ERASER c,q,s,f ), which are only able to solve 56% and 45% of the instances, respectively, within the same time frame. For the 45% of instances which can be solved by the three formulations, we observe that (MIP-p-S q,y ) offers a much tighter percentage optimality gap than the other two formulations. Because of this, the node usage in the B&B tree is significantly smaller in (MIP-p-S q,y ) compared to (ERASER c,q,s,f ) and (SDOBSS q,y,s ). Table 3 records the mean percentage optimality gap at the root node across all the instances for the three formulations under study. Observe that (MIP-p-S q,y ) is significantly tighter than the LP relaxations of the other formulations. We may thus conclude that for the payoff matrices without variability, (MIP-p-S q,y ) is the fastest formulation for the most difficult instances. On the other hand, (ERASER c,q,s,f ) is the fastest formulation when we endow the security game with further structure by allowing matrices to experience variability. Even then, (ERASER c,q,s,f ) looses ground to (MIP-p-S q,y ). This is due to the fact that (MIP-p-S q,y ) has the tightest LP relaxation. The quality of the upper bound obtained from (MIP-p-S q,y ) translates into a smaller B&B tree and this translates into reaching optimality of the integer problem faster in many cases.

Conclusions and future work

In this paper, we consider Stackelberg games in two different settings. We first analyze the general Stackelberg setting, which models a hierarchical competitive game between different agents, and the specific Stackelberg security setting, where an agent must secure subsets of targets from attackers.

In the general setting, we have studied known MILP formulations and have ordered them with respect to the strength of their linear relaxations. We have presented a formal theoretical link between GSG formulations and SSG formulations involving the projection of variables. Exploiting this link has allowed us to i) derive two new SSG MILP formulations (SDOBSS q,y,s ) and (MIP-p-S q,y ); and ii) extend our study of GSG formulations to SSG formulations, leading to a ranking of the security formulations with respect to the strength of their linear relaxations, where (MIP-p-S) has been shown to be the strongest SSG formulation. Further, we have shown its single type of attacker restriction, (MIP-1-S q,y ), to be an ideal formulation.

Our computational studies have shown that (MIP-p-G q,z ) and (MIP-p-S q,y ), the tightest formulations in each setting, are highly competitive with respect to solving time. Further, in the case of (MIP-p-S), we have seen it scales significantly better than competing formulations when tackling instances with no variability in their payoff structure. Formulation (MIP-p-S) represents a significant theoretical and practical improvement over previously existing SSG formulations.

However, the obvious bottleneck, at this time, is solving the tighter but larger LP relaxations for (MIP-p-G q,z ) and (MIP-p-S q,y ). The main challenge is to provide an efficient way of solving these tight formulations. It is our contention that this can be done by exploiting the inherent problem structure in the Stackelberg paradigm to develop either decomposition or cutting plane approaches.

While this paper focuses on the polyhedral analysis of general normal form Stackelberg games and Stacekberg security games, the work of developing efficient algorithms by conducting similar polyhedral analysis of the bilevel interaction could be carried out for Stackelberg games in specific security applications. In particular, extensions to problems that consider multiple attacks by followers, dynamic settings, imperfect information, or non-rational response would be interesting lines of future research.

Figure 1 :

 1 Figure 1: GSGs: |I| ∈ {10, 20, 30}, |J| ∈ {10, 20, 30}, |K| ∈ {2, 4, 6}-without variability.

Figure 2 :

 2 Figure 2: GSGs: |I| ∈ {10, 20, 30}, |J| ∈ {10, 20, 30}, |K| ∈ {2, 4, 6}-with variability.

  be strict, consider the following example where m = 1, |J| = 3 and |K| = 1. Let the reward and penalty matrices for the defender and attacker be D(•|p) = [1, 0, 0], D(•|u) = [0, 0, 0], A(•|p) = [0, 0, 0] and A(•|u) = [0, 0, 0]

Figure 3 :

 3 Figure 3: SSGs: K = {6, 8, 10, 12}, J = {30, 40, 50, 60, 70}-with variability

Figure 4 :

 4 Figure 4: SSGs: K = {2, 4, 6, 8}, J = {10, 20, 30, 40, 50}-without variability

Figure 5 :

 5 Figure 5: SSGs: K = {6, 8, 10, 12}, J = {30, 40, 50, 60, 70}-without variability

Table 2 :

 2 Mean percentage optimality gap at the root node recorded for GSG formulations.

		(D2 x,q,s,f )	(DOBSSq,z,s)	(MIP-p-Gq,z )
	Mean % opt. gap (no variability)	117.68	23.01	9.94
	Mean % opt. gap (with variability)	103.44	40.74	5.17
	Total mean % opt. gap	110.56	31.88	7.56

Table 3 :

 3 Mean percentage optimality gap at the root node recorded for SSG formulations.

Acknowledgements

Casorrán wishes to acknowledge the FNRS for funding his PhD research through a FRIA grant. Fortz and Labbé have been partially supported by the Fonds de la Recherche Scientique -FNRS under Grant(s) no PDR T0098.18. Ordóñez acknowledges the support of CONICYT through grant FONDECYT-1171419 and the Complex Engineering Systems Institute through grant CONICYT-PIA-FB0816. The authors would also like to thank the anonymous reviewers, whose comments have helped to elevate the quality of this paper.


Quentin Breche

Grégory Chagnon

Guilherme Machado

Benjamin Nottelet

Xavier Garric

Edouard Girard

Denis Favier

A non-linear viscoelastic model to describe the mechanical behavior's evolution of biodegradable polymers during hydrolytic degradation

Keywords: PLA-b-PEG-b-PLA, Mechanical behaviour's evolution, Biodegradable, Viscoelasticity, Mechanical model

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Due to its large range of mechanical properties and degradation times, the biodegradable and biocompatible tribloc copolymer PLA-b-PEG-b-PLA is an excellent candidate for tissue engineering applications [START_REF] Tessmar | Customized peg-derived copolymers for tissue-engineering applications[END_REF]. In this case, the polymer is used to be a physical temporary support for the tissular reconstruction called "scaffold". The objective of this scaffold is to be a substitute to injured tissues during healing and to provide a three dimensionnal environment for tissue reconstruction. Moreover, in many cases, the scaffold has to temporarily replace the functions of the tissues, especially the mechanical support. Thus, it is necessary to have a detailed knowledge of the initial mechanical properties and their evolution due to degradation. Therefore, in order to design a polymeric architecture for scaffolds, the development of mechanical constitutive equations that are able to model the initial mechanical properties and that take into account the degradation is essential.

In their initial state, the PLA-based biodegradable materials can present some complex nonlinear behaviors, especially non-linear viscoelasticity [START_REF] Soares | Constitutive framework for biodegradable polymers with applications to biodegradable stents[END_REF]. In order to model these properties, a large class of non-linear viscoelastic models used to describe the mechanical properties of polymeric biomaterials and living tissues is based on single integral models. These models are extensions of the linear viscoelastic model [START_REF] Ward | Mechanical properties of solid polymers[END_REF]. Pipkin and Rogers [START_REF] Pipkin | A non-linear integral representation for viscoelastic behaviour[END_REF] developed a single integral model that can take into account the dependence on strain of stress relaxation modulus. By considering that strain and time variables can be separated, Fung [START_REF] Fung | Biomechanics: mechanical properties of living tissues[END_REF] simplified the Pipkin and Rogers model. In his model, the strain of the linear single integral model is replaced by a non-linear function of strain called "instantaneous elastic response". Nevertheless, the model does not take into account the strain dependence of the time-related response and can be difficult to calibrate [START_REF] Abramowitch | An Improved Method to Analyze the Stress Relaxation of Ligaments Following a Finite Ramp Time Based on the Quasi-Linear Viscoelastic Theory[END_REF]. In order to overcome these problems, Nekouzadeh et al. [START_REF] Nekouzadeh | A simplified approach to quasi-linear viscoelastic modeling[END_REF] developed an alternative approach called "adaptive quasilinear viscoelastic" model that allows to simplify calibration and considers the strain non-linearity of the viscoelastic part.

The properties of PLA-based polymers are changing in use due to degradation. Their main degradation mechanism is the hydrolysis of ester links that breaks the polymeric bonds. It means that, in an aqueous media, a decrease of number average molecular weight, a mass loss and a material erosion can be observed [START_REF] Göpferich | Mechanisms of polymer degradation and erosion[END_REF]. Another important macroscopic effect of polymeric bonds cut is the lost of mechanical properties. Many studies show that Young's modulus and tensile strength decrease after a variably long latency period [START_REF] Tsuji | Autocatalytic hydrolysis of amorphous-made polylactides : effects of l -lactide content , tacticity , and enantiomeric polymer blending[END_REF][START_REF] Karjalainen | Biodegradable lactone copolymers. III. Mechanical properties of c-paprolactone and lactide copolymers after hydrolysis in vitro[END_REF][START_REF] Leroy | Investigation on the properties of linear PLA-poloxamer and star PLA-poloxamine copolymers for temporary biomedical applications[END_REF]. Thus, to predict the evolution of the mechanical properties of scaffolds, the constitutive equation must include hydrolytic damage. Two types of methods can be distinguished to model mechanical properties degradation. The first one consists in modelling the evolution of mechanical parameters such as Young's modulus or tensile strength by a micro-mechanical approach. The entropy spring theory [START_REF] Wang | A phenomenological model for the degradation of biodegradable polymers[END_REF] or the molecular dynamic approach [START_REF] Ding | A molecular dynamics study of Young's modulus change of semicrystalline polymers during degradation by chain scissions[END_REF] are used to model the Young's modulus evolution. By considering scissions as cavities in the polymeric matrix, Samami and Pan [START_REF] Samami | A constitutive law for degrading bioresorbable polymers Journal of the mechanical behavior of biomedical materials[END_REF] modelled the polymer as a continuum solid containing randomly distributed cavities. Then, using existing theories for porous solids, they deduced a constitutive equation that predicts the evolution of Young's modulus and tensile strength. The second method consists in adapting classical hyperelastic [START_REF] Vieira | Mechanical study of PLA-PCL fibers during in vitro degradation[END_REF], viscoelastic models [START_REF] Muliana | Modeling the response of nonlinear viscoelastic biodegradable polymeric stents[END_REF] or viscoplastic model [START_REF] Vieira | Constitutive modeling of biodegradable polymers: Hydrolytic degradation and time-dependent behavior[END_REF] to degradation by assuming that the model constitutive material parameters are dependent of a degradable variable [START_REF] Soares | Constitutive framework for biodegradable polymers with applications to biodegradable stents[END_REF]. These models are only calibrated on simple uniaxial load tensile tests. Breche et al. [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF] performed uniaxial load-relaxation tests and shew that the normalized relaxation function can be considered as an invariant of the degradation. Thus, in the linear viscoelasticity formalism, they expressed the degradation variable as a function of a proportionnality coefficient between the curves at different steps of degradation. Their model permits to describe the mechanical behavior for small deformation i.e. until 2%.

The objective of this paper is to establish a mechanical constitutive equation that is able to fit the PLA-b-PEG-b-PLA uniaxial behavior upon a strain of 2%. To achieve this, a non-linear viscoelastic adaptive quasi-linear viscoelastic (AQLV) model that takes into account the hydrolytic degradation is developed. First, after presenting materials and methods, experimental results of load-relaxation tests on PLA-b-PEG-b-PLA samples are presented. In a second part, the formalism of the AQLV model is introduced. Then, the calibration and the efficiency of the model to fit the experimental data are presented for the non-degraded material. In a last part, the AQLV model is adapted to take into account the hydrolytic degradation by a hydrolytic damage variable. The ability of this model to fit degraded experimental data is finally discussed.

Materials and methods

Material

The material used in this study is the same as the one used in a previous study [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF]. The details about the elaboration of the material are related here. Poly(ethylene glycol) (average Mn 20 000 g/mol), tin(II) 2-ethylhexanoate (Sn(Oct)2, 95%), dichloromethane (DCM), diethyl ether and tetrahydrofuran (THF) were purchased from Sigma-Aldrich (St-Quentin Fallavier, France). D,L-lactide (D,L-LA) was purchased from Purac (Lyon, France).

PLA50-b-PEG-b-PLA50 triblock copolymer was synthesized following a procedure previously described by Leroy et al. [START_REF] Leroy | Investigation on the properties of linear PLA-poloxamer and star PLA-poloxamine copolymers for temporary biomedical applications[END_REF]. Typically, predetermined amounts of D,L-LA and PEG were introduced in a flask. Sn(Oct)2 (0.1 molar % with respect to LA units) was then added. After degassing, the flask was sealed under vacuum and polymerization was allowed to proceed at 110 o C. After 5 days, the copolymer was recovered by dissolution in DCM and precipitation in cold diethyl ether. Finally, the product was dried under reduced pressure to constant mass. The copolymer was obtained with a yield of 90 %. Polymerization degree of each PLA block and molecular weight of the synthesized triblock copolymers were calculated using the following equations:

DP P LA (%) = 1 2 × DP P EG EG LA (1)
M ntribloc (%) = 2 × (DP P LA × 72) + M nP EG (2)
with EG/LA being the ratio of ethylene oxide and lactyl units calculated from 1 H NMR spectra. The polymer films have been manufactured by solvent evaporation. A predefined quantity of polymer was solubilized in acetone, spread in a small dish and placed during 36 hours under an extractor hood for evaporation. Final solvent removal was obtained by further drying in vacuo under 1.10 -3 mBar for 3 days. A 0.5 mm thick film was then obtained.

Methods for mechanical tests

Dogbone tensile specimens were cut in 0.5 mm thick films with a specific punch of 14 mm in length and 2 mm wide gauge length. They were then placed in 20 ml test tubes and submitted to two times of degradation, 1 and 3 weeks, in phosphate buffer solution (PBS) at 37 o C within an agitator/incubator. The mechanical tests were realized in an aqueous media at a controlled temperature of 37 o C via a steel hermetic bath adjustable on the mechanical test machine in order to be in conditions as closed as possible to in vivo ones.

Uniaxial tensile tests on dogbone specimens at three different strain levels (2%, 4% and 6%) with a strain rate during loading of 1% per second followed by a tensile relaxation test are performed. The latter were performed with a Gabo Eplexor mechanical test machine with a load cell of 10 N at 0, 1 and 3 weeks of degradation. The strain field in specimens is measured by a Vic 2D image correlation software. At least two specimens were used for each test. The standard deviations were calculated and indicated on the 6% load curves by error bars as illustrated in figure 1).

Experimental data

In a previous article [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF], the mechanical tests at 2% and 4% have already been presented as well as mass lost, number average molecular weight, and dispersity index. This paper presents the results for larger strains i.e. 6%. The figure 1 presents the stress-strain curves of the material for the tensile tests at a strain of 6% for every steps of degradation. In these curves, the strain is increasing until a strain of about 2% before reaching a plateau. The maximum stress reached in the non-degraded state is of 2.7 MPa. During degradation, the hydrolytic process cleaves the polymeric chains, leading to a decrease of molecular weight. It implies a decrease of the mechanical properties. As can be observed in figure 1, the maximum stress reached decreases even if the the general shape of the curve is conserved. After 3 weeks of degradation, about 20% of the inital non-degraded maximum stress is lost. The curve in 2 and 3 represent the load-relaxation tests for strains at 2% and 4% respectively and have already been presented in in a previous paper [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF]. The figure 4 shows the stress recorded over time during the load-relaxation tests at a strain of 6%. In this last figure, the load part corresponds to figure 1. After the loads, when strain is maintained constant, a stress relaxation is observed. It illustrates that the PLA-b-PEG-b-PLA is a viscoelastic material. A stress steady-state is reached after about 30 seconds, meaning that characteristic relaxation times are quite short. The residual stress was measured at 0.03 MPa at 2% and at 0.08 MPa at 4% in paper [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF] and are measured here at 0.25 MPa at 6% in the non-degraded state. In order to compare the relaxation after the different degradation times, the normalized relaxation is used. It consists in the ratio of the current stress to the static stress before relaxation. If σ m is the static stress reached at the end of the load (corresponding to the dotted vertical lines in figures 2, 3 and 4) and σ relax the stress recorded during relaxation, the equation of the normalized relaxation is:

σ nr = σ relax σ m = σ relax σ relax (t = 0) (3)
In this last equation, t represents the time recorded during the stress relaxation. Thus, the initial time corresponds to the end of the load (represented by the dotted vertical lines in figure 2, 3 and 4). The normalized stress relaxation curves are presented in figure 5 for the different strains. One can observe that curves are not superimposed, meaning that the normalized stress relaxation of the material depends on the strain. This behavior is characteristic of the non-linear viscoelasticity.

Constitutive modelling: an adaptive quasi-linear viscoelastic degradable model

The aim of this section is to develop a viscoelastic degradable model that takes into account the non-linear viscoelasticity of the PLA-b-PEG-b-PLA and that allows to model its mechanical behaviour until a strain of 6%.

Hypothesis

It is considered that the polymer is sufficiently hydrophilic to present a bulk erosion scheme and sufficiently thin to neglect autocatalytic effects [START_REF] Grizzi | Hydrolytic degradation of devices based on poly(DL-lactic acid) sizedependence[END_REF]. Thus, the strain and degradation states are identical in every point of the specimens.

Adaptive quasi-linear viscoelastic model

In a previous study, Breche et al. [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF] developed a simple generalized Maxwell model to predict the behaviour of PLA-b-PEG-b-PLA. The model was efficient for strain lower than 2%. Nevertheless, it was unable to fit the experimental tests upon 2% strain and to take into account the strain dependence of the normalized relaxation function. In the previous paper, a linear viscoelastic Maxwell model was used. In order to improve the modelling, a non-linear viscoelastic model called "adaptive quasi-linear viscoelastic" (AQLV) is considered. In the general case of the AQLV model , the stress is equal to the product of a pure non-linear function k(ε) and an intermediate variable called the "viscoelastic strain" V ε (t) [START_REF] Nekouzadeh | A simplified approach to quasi-linear viscoelastic modeling[END_REF], that is a linear convolution integral.

   σ(t) = k(ε(t))V ε (t) V ε (t) = t 0 g(t -τ ) dε dτ dτ (4)
where g(t) is a "reduced" relaxation function. To take into account the strain dependance of the normalized relaxation functions, degrees of freedom are added to the model by allowing different non-linear behaviours for different "reduced" relaxation functions as suggested by Nekouzadeh et al. [START_REF] Nekouzadeh | A simplified approach to quasi-linear viscoelastic modeling[END_REF].

σ(t) = σ ∞ (ε(t)) + i k i (ε(t))V ε i (t) V ε i (t) = t 0 g i (t -τ ) dε dτ dτ (5)
σ ∞ is a function of strain only representing the residual stress. The "reduced" relaxation functions are chosen to be decreasing exponential functions g i (t) = e -t τ i (ε) [START_REF] Nekouzadeh | A simplified approach to quasi-linear viscoelastic modeling[END_REF] where τ i are characteristic relaxation times. Thus, by analogy with the previous article [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF], the model can be represented in term of hypothetical parallel Maxwell elements. Its rheological scheme is presented in figure 6. In this representation, spring stiffnesses (k i (ε)) and dashpot coefficients (η i (ε)) are functions of the overall strain. For each branch, it exists the following relationships:

   V ε i + V ε i τ i (ε) = ε σ i = k i (ε)V ε i (6)
where :

τ i (ε) = k i (ε) η i (ε) (7)
To simplify, it is considered that the spring stiffnesses and dashpot coefficients are proportionnal to the same nonlinear function of strain f (ε) [START_REF] Nekouzadeh | A simplified approach to quasi-linear viscoelastic modeling[END_REF], meaning that :

τ i (ε) = k i (ε) η i (ε) = k i f (ε) η i f (ε) = k i η i = τ i (8)
Thus, τ i , the characteristic relaxation times, are independent on strain.

k i k 3 k 0 ( ) k 1 ( ) k 2 ( ) ( ) ( ) ηi ( ) η 3 ( ) η 2 ( ) η 1 ( ) Strain: V Strain rate: V / i

Calibration of the model on the non-degraded experimental data

At first, the model is calibrated on the non-degraded experimental data. As seen in previous part (cf. figure 5), the non-linear character of the polymer behavior does not allow to consider the normalized relaxation as an invariant of the strain. Sarver et al. [START_REF] Sarver | Methods for Quasi-Linear Viscoelastic Modeling of Soft Tissue: Application to Incremental Stress-Relaxation Experiments[END_REF] proposed a normalization in the following system of stress coordinates:

σ * nr = σ * relax σ * m = σ relax (ε) -σ ∞ (ε) σ m (ε) -σ ∞ (ε) = Q (9)
Q is called the normalized relaxation factor. In figure 7, the normalized relaxation factors measured at strains of 2%, 4% and 6% are plotted. It can be observed that these curves are nearly similar. Thus, it is considered that it exists only one normalized relaxation factor master curve that is independent on strain. Nekouzadeh et al. [START_REF] Nekouzadeh | A simplified approach to quasi-linear viscoelastic modeling[END_REF], determined the values of k i at every strain points for which a relaxation stress was measured. For this, a complex non-linear least square method that allows to fit simultaneously the relaxation curve at these different strains is used. From the obtained equations, the values of k i are deduced at these particular strain points. Then, the function k i (ε) is fitted from these points to obtain a value at every strain. Sarver et al. [START_REF] Sarver | Methods for Quasi-Linear Viscoelastic Modeling of Soft Tissue: Application to Incremental Stress-Relaxation Experiments[END_REF] shew that if the normalized relaxation factor is an invariant of strain, the calibration of the Fung's QLV model is greatly simplified. It is also the case for the AQLV model. It allows to determine, in a simple way, a value of k i (ε) at every experimental strain. At first, the equation of the master curve must be determined. To this end, the mean of the three normalized relaxation factors curves at strains of 2%, 4% and 6% is calculated. Then, this curve is fitted. Breche et al. [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF] shew that for a PLA-b-PEG-b-PLA, the relaxation curve can be fitted with a three temporal terms Prony's serie. Thus, the master normalized relaxation factors curve is also fitted with three decreasing exponential functions:

Q m = 3 i=1 a i e -t τ i (10)
where a i and τ i are material parameters. The values of the parameters to fit the experimental data curves are presented in table 1. The corresponding calculated relaxation curve is plotted in figure 7. Then, from the values obtained for the tests at 2%, 4% and 6%, and by considering that the residual stress is equal to 0 in the absence of strain, the residual stress σ ∞ (ε) is fitted as a function of the strain in order to obtain a value for every strain. Given the fact that four points need to be fitted, a three order polynomial function is chosen. Thus, σ ∞ (ε) takes the form of a polynomial function of ε, f (ε):

f (ε) = c 0 + c 3 ε 3 + c 2 ε 2 + c 1 ε (11)
c i are the coefficients of the polynomial function and their values are presented in table 2. The fitted curve can be seen in figure 9. In the formalism of the model, the expression of the normalized relaxation factor is:

Q(ε, t) = σ relax (ε) -σ ∞ (ε) σ m (ε) -σ ∞ (ε) = 3 i=1 k i (ε)τ i (e T τ i -1)e -t τ i 1 σ m (ε) -σ ∞ (ε) × dε dτ (12)
σ m (ε) represents the value of the stress in the load measured at a given strain and T represents the duration of the load. Thus, from the values of σ m (ε), and the fitted equations of σ ∞ and Q m , the values of k i can be determined from equations [START_REF] Wang | A phenomenological model for the degradation of biodegradable polymers[END_REF] for every σ m (ε) measured and are presented in figure 8. As for the residual stress, the values of k i must be fitted to obtain an analytical expression of k i (ε) for every strain. In the same way, a three order polynomial function of the form of f (ε) was determined. The coefficients of the function are presented in table 2 for every k i (ε). It perfectly fits the curves except for the first parts of k 2 and k 3 (cf. figure 8 and9). As can be seen in figure 10 that presents the fitted load-relaxation curves and the corresponding experimental curves, the mechanical behavior is perfectly predicted for every strain. The lack of accuracy of the fitting for k 2 and k 3 at the beginning of the curves has thus no consequences on the ability of the model to represent experimental data.

An AQLV degradable model based on the invariance of the normalized relaxation

To take into account the evolution of mechanical properties during degradation, a scalar parameter d h (x, t d , ε) is introduced. It represents the degraded state of the polymer at each material point x and at a given degradation time t d [START_REF] Soares | Constitutive framework for biodegradable polymers with applications to biodegradable stents[END_REF]. As seen in the hypothesis in part 4.1, d h is supposed independent of the location in the specimen. d h has a value between 0 and 1, 0 representing a non-degraded state and 1 a fully degraded state. The figures 11, 12 and 13 represent the normalized relaxation curves for 2%, 4% and 6% respectively at every step of degradation (calculated from the relaxation curves of the figures 2, 3 and 4 and equation 3). For the two first strains, the normalized relaxation curves are perfectly superimposed, meaning that the temporal viscous component remains invariant during the three first weeks of degradation for these levels of strain. Nevertheless, at a strain of 6%, a shift between the normalized relaxation curves is observed. At the third week of degradation, the residual stress of the normalized relaxation curves decreases. However, to simplify the model, it will be considered that the normalized relaxation curves, for every fixed value of strain, are independent on the degradation state. Thus, by admitting, at a given strain ε, that the non-degraded normalized relaxation curve is equal to the normalized relaxation curve at a degraded step, the equation of the degraded relaxation curve can be deduced. For this, equation 3 is used in considering that the relaxation stress depends on the degradation parameter and deformation.

σ relax (d h , ε, t) = σ relax (d h , ε, t = 0) σ relax (d h = 0, ε, t = 0) × σ relax (d h = 0, ε, t) = σ m (d h , ε) σ m (d h = 0, ε) × σ relax (d h = 0, ε, t) (13)
In a load-relaxation test, σ m (d h , ε) is, for a given strain ε and a degradation state d h , the stress reached at the end of the load and corresponds to the time zero of the relaxation. In the AQLV model, the equation of the stress relaxation at a given deformation is :

σ relax (d h , ε, t) = σ ∞ (ε) + 3 i=1 k i (ε(t)) × σ m (d h , ε) σ m (d h = 0, ε) × τ i (e T τ i -1)e -t τ i dε dτ (14)
By combining equation 13 and 14 the general following equation can be obtained:

σ relax (d h , ε, t) = σ ∞ (ε) × σ m (d h , ε) σ m (d h = 0, ε) + 3 i=1 k i (ε(t)) × σ m (d h , ε) σ m (d h = 0, ε) × τ i (e T τ i -1)e -t τ i dε dτ (15)
From the definition of Lemaitre et al. [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF] used by Vieira et al. [START_REF] Vieira | Mechanical study of PLA-PCL fibers during in vitro degradation[END_REF] for hydrolytic damage, the degradation variable is defined as:

d h (ε) = 1 - σ m (d h , ε) σ m (d h = 0, ε) (16)
Consequently, the stress relaxation curve can be expressed as a function of d h :

σ relax (d h , t) = σ ∞ × (1 -d h ) + 3 i=1 k i (ε(t)) × (1 -d h )τ i (e T τ i -1)e -t τ i dε dτ (17)
. Finally, from equation 5, the AQLV degradable model can be written as:

     σ(t) = σ ∞ (ε(t))(1 -d h ) + i k i (ε(t))(1 -d h )V ε i (t) V ε i (t) = t 0 g i (t -τ ) dε dτ dτ (18)
From equation ( 16) the degradation parameters can be calculated as a function of strain for 1 and 3 weeks of degradation. The curves are presented in figure 14 and 15. In the first percent strain of the load, d h suffers important variations. It means that during this transitory part, the non-degraded and the degraded load curves are not perfectly proportionnal. Nevertheless, beyond this part, d h becomes almost constant, meaning that the shape of the curve is conserved in the course of degradation. In order to simplify the model, the hypothesis that d h is constant and independent on strain is used. To determine the value of d h , the mean value between 1% and 6% is taken (black curves in figures 14 and 15). When water diffuses into the polymeric structure, hydrolysis reaction breaks the ester links leading to the chains cleavage. This leads to the creation of oligomers and monomers. The main consequence is the decrease of the number average molecular weight and of the mechanical properties. It will be considered that the evolution of the degradation parameter is directly correlated with the decrease of the number average molecular weight [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF]. In figure 16, the degradation parameter obtained from figures 14 and 15 is plotted as a function of the number average molecular weight for 0, 1 and 3 weeks (black dots). In order to model the evolution of the degradation parameter for every intermediate step, the degradation parameter is fitted as a function of the normalized number average molecular weight. The expression proposed by Samami and Pan [START_REF] Samami | A constitutive law for degrading bioresorbable polymers Journal of the mechanical behavior of biomedical materials[END_REF] for the evolution of tensile strength is used:

d h (ε) = 1 - σ m (d h , ε) σ m (d h = 0, ε) = 1 -(1 -f v ) n (19)
where f v 1 m 3 is the volume fraction of the cavities generated by polymer scissions. f v can be expressed, for early degradation times, as a function of the normalized number average molecular weight [START_REF] Samami | A constitutive law for degrading bioresorbable polymers Journal of the mechanical behavior of biomedical materials[END_REF]. Thus, the expression of the degradation parameter becomes:

d h = 1 -1 -K × M n0 M n -1 n ( 20
)
where M n is the actual number average molecular weight and M n0 is the initial number average molecular weight. K is a fitting parameter that depends on the initial number average molecular weight, initial density, and radius of cavities created by scissions of the polymer. n is also a fitting parameter. The values of these parameters are presented in table 3. They are deduced from the measurement of the normalized number average molecular weight during degradation performed in a previous study [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF]. The fitted function is presented in figure 16. In order to link the evolution's kinetics of normalized number average molecular weight to time, the following expression has been proposed previously [START_REF] Breche | Mechanical behaviour?s evolution of a PLA-b-PEG-b-PLA triblock copolymer during hydrolytic degradation[END_REF]:

M n0 M n = 1 + Kt ( 21
)
where K is a parameter that is proportionnal to the hydrolysis rate and equal to 0.62 in this case. As can be seen in figures 17 and 18 the model perfectly fits both load and relaxation curves at strains of 2% and 4% for every week of degradation. At 6% strain, the model is efficient to predict the mechanical behavior for 0 and 1 week of degradation as can be observed in figure 19. Nevertheless, it starts to be less accurate for the third week of degradation at 6% strain. The beginning of the load part is a bit underestimated and the end overestimated. This result was expected and is due to the approximation that the degradation parameter is constant for every strain, that is not totally the case at three weeks of degradation (cf. figure 15). At this degradation and strain states, the relaxation curve is not perfectly fitted. Indeed, the residual stress is overestimated and its value is twice times larger than the experimental one. It can be explained by the fact that the hypothesis of the invariance of the normalized relaxation function during degradation is not totally verified. It can be seen in figure 13 that at 6%, the normalized relaxation functions are not exactly the same, especially for 3 weeks of degradation.

Conclusions

In order to predict the behavior of a PLA-b-PEG-b-PLA biodegradable polymers during degradation, a non-linear viscoelastic model that takes into account the hydrolytic degradation during the three first weeks of degradation was developed. From an experimental protocol of load-relaxation tests, the non-linear viscoelastic behavior of the PLA-b-PEG-b-PLA was highlighted. An adaptive quasi-linear viscoelastic model was then calibrated on the non-degraded tests. The ability of the model to describe the uniaxial load-relaxation behavior of the non-degraded polymer was demonstrated for every strain lower than 6%. By using the normalized relaxation function invariance during degradation, the adaptive quasi-linear model was adapted to take into account the degradation through a hydrolytic degradation parameter. Then, a physical kinetics evolution of this parameter was used to model every degradation times between zero and three weeks.The model is compared to experimental data for every degradation step and strain. It was shown that the model fits very well the experimental data. Nevertheless, in the most degraded case, the relaxation curve that appears at the load at 6% is not perfectly fitted. This could be improved by choosing the normalized relaxation factor as a degradation invariant rather than the classical normalized relaxation curve but it changes the formalism and the model would become less easy to fit.

Figure 1 :

 1 Figure1: Experimental load curves at 6% strain load at 0, 1, and 3 weeks of degradation.

Figure 2 :

 2 Figure2: Experimental load and relaxation curves at 6% strain load at 0, 1, and 3 weeks of degradation.

Figure 3 :

 3 Figure 3: Experimental load and relaxation curves at 6% strain load at 0, 1, and 3 weeks of degradation.

Figure 4 :Figure 5 :

 45 Figure 4: Experimental load and relaxation curves at 6% strain load at 0, 1, and 3 weeks of degradation.

Figure 6 :

 6 Figure 6: Rheological scheme of the model

Figure 7 :

 7 Figure 7: Normalized relaxation factors for different strain levels and the corresponding fitted master curve

Figure 8 :FitFigure 9 :

 89 Figure 8: Values of the functions ki and their corresponding analytical fit

Figure 10 :

 10 Figure 10: Experimental data versus calibrated model for different strain levels

2 Figure 11 :Figure 12 :

 21112 Figure 11: Normalized relaxations curve at a strain of 2% at different degradation steps

Figure 13 :

 13 Figure 13: Normalized relaxations curve at at a strain of 6% at different degradation steps

Figure 14 :

 14 Figure 14: Degradation parameter as a function of strain after 1 week of degradation

Figure 15 :

 15 Figure 15: Degradation parameter as a function of strain after 3 weeks of degradation

Figure 16 :

 16 Figure 16: Degradation parameter's evolution fitting during degradation

Figure 17 :Figure 18 :

 1718 Figure 17: Experimental load relaxation tests (right) for 0, 1 and 3 weeks of degradation with corresponding load stress-strain curves (left) compared with calibrated AQLV model on the 2% load strain's tests.

Figure 19 :

 19 Figure19: Experimental load relaxation tests (right) for 0, 1 and 3 weeks of degradation with corresponding load stress-strain curves (left) compared with calibrated AQLV model on the 6% load strain's tests.

Table 1 :

 1 Values of the coefficients of the polynomial functions of kε and σ∞

	a 1 (MPa) a 2 (MPa) a 3 (MPa) τ 1 (s) τ 2 (s) τ 3 (s)
	0.4	0.44	0.16	0.22	1.76 12.56

Table 2 :

 2 Values of the coefficients of the polynomial functions of kε and σ∞

		k 1 (MPa)	k 2 (MPa)	k 3 (MPa)	σ ∞ (MPa)
	c 0	2.08 × 10 2	1.50 × 10 2	49.30	0
	2.34 × 10 4 c 2 -5.50 × 10 5 c 1 c 3 4.07 × 10 6	-4.34 × 10 3 -2.10 × 10 3 8.00 × 10 4 4.00 × 10 4 -5.41 × 10 5 -2.72 × 10 5 2.10 × 10 3 2.50 -98.70

Table 3 :

 3 Values of the parameters of the degradation variable kinetic model

	K	n
	0.32 0.29