Hadrien Cambazard

Nicolas Catusse

ROADEF CHALLENGE 2014: A MODELING APPROACH Rolling stock unit management on railway sites

Introduction

We report here our analysis of the problem presented in [START_REF] Franc | Roadef/euro 2014 challenge, trains don't vanish !final phase, rolling stock unit management on railway sites[END_REF] and propose a methodology strongly based on modeling with MIP and CP technologies. Due to the complexity of the problem formulation, we believe that a robust engineering is easier to achieve by relying on models than dedicated code. Two core modeling ideas are presented for relating the daily maintenances limit and the linked departures to an assignment based MIP model. Additionnaly, a variant of the maximum matching problem lying at the heart of the problem is shown to be NP-Complete. Intermediate experimental results are given along the way to support the ideas reported. Section 2 gives an overview of the proposed methodology. The two main components of this approach relying on a MIP model and a CP model are discussed in section 3 and 4 respectively.

Overview of the proposed method

We propose to decompose the problem into two natural stages of decisions.

Firstly, We suggest to handle the assignment of the arrival's trains to the departures as well a the choice of platforms for all events with a MIP solver. The rationale is twofolds: this first problem has a strong matching component for which MIP technologies are efficient and the objective function can be nearly entirely expressed on this first part. Two key modeling ideas are presented to integrate daily maintenances limits and linked-arrivals naturally in the matching based MIP model.

Secondly, we believe that the routing of trains through the station can be computed without the need of really questionning the assignment first step. This was the case on the first problem formulation. From an engineering point of view as well as efficiency, we believe that CP is the best technology to implement this part and adress the routing under time-windows constraints.

We were not able to fully implement the proposed methodology but will give preliminary experimental results that support our analysis of the problem.

Assignment of trains, platforms and maintenances

We chose to give the most important lines of the MIP model used and skip the details that we believe do not prevent to understand the approach (and reproduce it) but would significantly burden the readibility due the complexity of the problem.

Notations

We use the notations introduced in [START_REF] Franc | Roadef/euro 2014 challenge, trains don't vanish !final phase, rolling stock unit management on railway sites[END_REF] regarding the problem definition. Additionally, a bipartite graph is denoted G(A, B, E) and the set of edges incident to a vertex x ∈ A ∪ B is denoted Γ (x). If G is oriented the notation Γ is extended so that Γ + (x) and Γ -(x) are the set of vertices that are incident to arcs respectively leaving and entering vertex x: Γ + (x) = {y | (x, y) ∈ E} and Γ -(x) = {y | (y, x) ∈ E}. This notation is generalized to sets and for X ⊆ A, we denote by

Γ + (X) = ∪ x∈X Γ + (x). Finally for Y ⊆ B we denote Γ -1 (Y ) = {x ∈ A | Γ + (x) ⊆ Y } the set
of vertices of X that have all their neighours in Y .

Model restricted to unlinked arrivals

In this section we design a model under the assumption that all arrivals are known and do not depend on decisions made for earlier departures i.e we assume that A does not contain any linked-arrivals. Under this assumption, the category and all associated features of an arrival are fully known. We then show how the model presented for unlinked arrivals can be generalized to linked arrivals.

Assignment of arrivals to departures.

We first present the variables encoding the assignment of arrivals to departures:

y j ∈ {0, 1} : indicates if departure j ∈ D is covered (assigned to an arrival). w i ∈ {0, 1} : indicates if arrival i ∈ A is covered (if yes it must go on a platform). x ij ∈ {0, 1} : indicates whether arrival i ∈ A ∪ T is assigned to departure j ∈ D. u gj ∈ {0, 1} : indicates if departure j ∈ LD is of group category g.
The key assignment variables are therefore the x ij variables. The possible assignment of a train (an initial or an arriving train) to a departure is a-priori limited by three elements (chronology, maximum maintenance capacities and categories). For a given pair of train i ∈ A ∪ T and departure j ∈ D, we directly state x ij = 0 (infeasible assignment) if one of the three conditions hold:

-Chronology: The time between the arrival and departure is not large enough to accomodate for any maintenance and junction/disjunction needed as well as routing through the station. The routing is represented by a unique constant b 1 representing a time buffer that we estimate initially.

arrT ime i +

1 remDBMt i <reqDBMj maintT imeD catt i + 1 remDBMt i <reqDBMj maintT imeT catt i + 1 i∈Jarr disjT ime × (|jaList i | -1) + 1 j∈J dep (junT ime × (|jdList j | -1) + minAsbT ime) + b 1 > depT ime j
Note that the full disjunction/junction time is included whenever i/j are respectively part of a joint arrival/departure. This can be refined (or removed in particular to use the MIP for a lower bound) when the assignmengt i to j can potentially be done without these operations (their neighbhoring rolling units are potentially compatible in the corresponding convoys).

-Maximum maintenance capacity: The trains don't have enough maintenance capacity even if maintenance is done.

maxDBM catt i < reqDBM j or maxTBM catt i < reqTBM j
-Categories: Categories are not compatible.

cat ti ∈ compCatDep j
Constraints (1)-( 2) channels x, w and y variables and also enforce an arrival to be assigned to at most one departure; (2) channels u and y variables enforcing exactly one group category for a covered departure belonging to a joint; (3) channels x and u making sure that the arrivals assigned to a joint departure have the same corresponding category.

(1)

j∈D x ij ≤ w i ∀i ∈ A ∪ T (2) i∈A∪T x ij = y j ∀j ∈ D (3) g∈G u gj = y j ∀j ∈ J dep (4) i∈A|catGroupt i =g x ij ≤ |jdList j |u gj ∀j ∈ J dep , ∀g ∈ G 2.
Assignment of arrivals/departures to platforms.

We turn our attention to the assignment of platforms. We use another time buffer b 2 corresponding to a minimum time that the rolling unit should be able to stay on the platform without causing a crash. The following decision variables are added:

pl ik ∈ {0, 1}: indicates if platform k ∈ P is assigned to arrival i ∈ A.
pl jk ∈ {0, 1}: indicates if platform k ∈ P is assigned to departure j ∈ D. l j ≥ 0 : total length of the convoy made for departure j ∈ D.

ct jc ∈ {0, 1}: indicates if category c ∈ C is assigned to departure j ∈ D.
The assignment of a platform k ∈ P to an arrival i ∈ A is a-priori restricted by four elements. We state pl ik = 0 if:

-Length:

length k < length ti -Categories: cat ti / ∈ compCatRes k -Imposed consumption: l ∈ I and [beg l , end l -1] ∩ [arrT ime i , arrT ime i + minResT ime + b 2 ] = ∅ -Arrival sequence: the sequence of i (arrSeq i ) is not connected to platform k.
Similarly Imposed consumption and Departure sequence can be used to prefix pl jk to 0. The model is now shown below. Constraints [START_REF] Franc | Roadef/euro 2014 challenge, trains don't vanish !final phase, rolling stock unit management on railway sites[END_REF] states that a single category is assigned to a covered departure; (6) channel departures category with the arrival assigned; (7) and (8) enforces a single platform to be chosen for the a covered arrival/departure; (9) channels the length of the train of a departure and to the one of the chosen arrival; (10) ensures the departure's platform is long enough; (11) ensures that a departure's platform is compatible with its category; (12) make sure that all arrivals in the same joint arrival are assigned to the same platform.

(5

) c∈C ct jc = y j ∀j ∈ D (6) ct jc ≥ x ij ∀j ∈ D, ∀c ∈ C, ∀i ∈ A|cat ti = c, (7) k∈P pl ik = w i ∀i ∈ A (8) k∈P pl jk = y j ∀j ∈ D (9) l j ≥ i∈A∪T length ti x ij ∀j ∈ D (10) length k pl jk + L(1 -pl jk ) ≥ l j ∀j ∈ D, ∀k ∈ P (11) pl jk ≤ 1 -ct jc ∀j ∈ D, ∀k ∈ P, ∀c ∈ C|c / ∈ compCatRes k (12) pl j1k = pl j2k ∀j ∈ J arr , ∀(j 1 , j 2 ) ∈ jaList j , ∀k ∈ P
Finally two platforms can be shared by two events (an event is either an arrival or a departure) only if the time gap inbetween is large enough to free the platform. A platform k can not be shared by two events (a,b) (ordered chronologically) if :

-a ∈ D, b ∈ D : depT ime b -depT ime a < idealDwell b + b 2 -a ∈ A, b ∈ A : arrivT ime b -arrivT ime a < idealDwell a + b 2 -a ∈ A, b ∈ D : depT ime b -arrivT ime a < idealDwell a +idealDwell b +b 2 -a ∈ D, b ∈ A : arrivT ime b -depT ime a < b 2 Fig. 2.
An exemple a maximum flow problem encoding the choice of days for maintenances. On this example arrival i is performing two maintenances (time and distance).

So finally ∀k ∈ P, ∀(a, b) ∈ D ∪ A × D ∪ A, such that event a can not share a platform with event b, we simply state:

pl ak + pl bk ≤ 1
The LP formulation can be strenghtened by stating these constraints over the maximum cliques of the incompatibility graph (but cplex does that automatically).

3.

Modeling daily maintenances limit. We know if an arrival will require maintenance operations when it is assigned to a given departure, and we don't need to know the exact day it is performed to state the daily limits. This idea is also proposed by [START_REF] Buljubašić | Team s18. In Challenge ROADEF/EURO 2014[END_REF] and stated similarly with constraints (13). We enforce the daily limits on all possible slidingtimewindows of all sizes (in {1, . . . , nbDays}) over the horizon:

(13) ∀d a ∈ {1, . . . , nbDays}, ∀d b ∈ {d a , . . . , nbDays}, (i,j)∈M (da,d b ) o ij x ij ≤ (d b -d a + 1) × maxM aint
We claim that the days of maintenances can be found afterwards by solving a maximum flow problem in a bipartite graph where arrivals are on the left side (an arrival i has an incoming flow o ij ) and possible days are on the right, each with a capacity of maxMaint. Figure 2 shows an example. We now give a proof of this claim.

Let's consider a bipartite flow graph G(s, t, A, B, E) (figure 2 shows an example of such a graph) defined by a bipartite graph G(A, B, E ′ ) (edges are oriented from A to B) as well as two additionnal vertices, a source s and a sink t, the additionnal edges (s, i), ∀ i ∈ A and (j, t), ∀ j ∈ B. Moreover a single capacity Capa is added on all edges leading to the sink i.e all (j, t), ∀ j ∈ B. We recall from 3.1 that for

Y ⊆ B, Γ -1 (Y ) = {x ∈ A | Γ + (x) ⊆ Y }.
Lemma 1: Given a bipartite flow graph and an entering flow f i for each edge (s, i)

for all i ∈ A, the flow is feasible if and only if i∈Γ -1 (Y ) f i ≤ Capa × |Y | for all Y ⊆ B.
Proof : The Hall's theorem states that a perfect matching (matching saturating A) exists if and only if for all X ⊆ A we have |X| ≤ |Γ + (X)|. This can be generalized to flows. Typically the flow (defined by the f i ) is feasible in the bipartite flow graph if and only if ∀ X ⊆ A, we have i∈X f i ≤ Capa × |Γ + (X)|. To prove it, we apply Hall's theorem in the bipartite graph where each vertex of i ∈ A is replicated f i times and each vertex j ∈ B is replicated Capa times. A feasible flow gives a maximum matching in the graph with replicated vertices by dividing the flow in units. Reversely by combining units, we can obtain the flows. Now for any X ⊆ A we can state Y = Γ + (X) (and so Γ -1 (Y ) = X) and by assumption we have

i∈X f i = i∈Γ -1 (Y ) f i ≤ Capa × |Y | = Capa × |Γ + (X)|.
Therefore there is a feasible flow since Hall's theorem is verified. Reversely, if we assume a feasible flow we obtain the condition from Hall's theorem.

A convex bipartite flow graph is a bipartite flow graph such that the underlying bipartite graph is convex over the set of vertices A. This means that there is an ordering of the vertices of A such that all neigbhors Γ + (i) of a vertex i ∈ A are consecutives (figure 2 shows an example of such a convex bipartite flow graph).

Lemma 2: Given a convex bipartite flow graph and an entering flow f i for each edge (s, i), ∀i ∈ A, the flow is feasible if and only if i∈Γ -1 (W ) f i ≤ Capa×|W | for all consecutive subset W of vertices of B.

Proof : Consider any subset Y of B, Y is made of a number k of subsets of con- secutives vertices so that Y = Y 1 ∪ Y 2 ∪ . . . ∪ Y k . For each Y l (l ≤ k) we have i∈Γ -1 (Y l ) f i ≤ Capa × |Y l | by assumption. Since the graph is convex, the neigb- hors of any two parts of Y are disjoint so that Γ -1 (Y a ) ∩ Γ -1 (Y b ) = ∅ for all pairs (Y a , Y b ) of distinct parts of Y (a ≤ k, b ≤ k and a = b). As a result we have i∈Γ -1 (Y ) f i = k l=1
i∈Γ -1 (Y l ) f i and we can state:

i∈Γ -1 (Y ) f i = k l=1 i∈Γ -1 (Y l ) f i ≤ k l=1 Capa × |Y l | = Capa × |Y |
So by Lemma 1, the flow is feasible. Reversely, if the flow is feasible, the condition (on all subsets Y ) of Lemma 1 holds and consequently the condition (on all con-secutive subsets W ) of Lemma 2.

Back to our problem, the bipartite graph of the arrivals and days is a convex bipartite flow graph (the set of possible maintenance days for an arrival is made of consecutive days without any holes and we set Capa to maxMaint). The conditions stated in the MIP is the one of Lemma 2 and is enough to make sure feasible maintenance days are always found afterwards by solving the corresponding flow. The non-satisfied preferred platform assignment cost as well as the non-satisfied train reuse cost are directly encoded based on the pl and x variables. Similarly the over-maintenance cost of an arrival/departure pair (i, j) that require a maintenance is known and counted when the corresponding x ij is set to 1. Modeling the train junction and disjunction operation costs requires more effort and intermediate variables must be introduced to represent each pair of consecutive train in a convoy. We do not give the full details here. We mention however that our model remains a relaxation of this cost as it does not take into account the need for junctions/disjunctions when maintenances have to be performed.

5.

The MIP formulation presented here has a very strong relaxation in practice. The reason is that the core problem is an assignment problem and that the side constraints seem relatively easy to take into account for this benchmark. Despite the size of the formulation, the results 1 shown in Table 1 prove that the MIP model without linked arrivals can be solved to optimality (which is not really needed in practice) for the B instances. Half of the instances are solved without branching.

Table 1. Performances for solving to optimality (precision 10 -3 ) the B instances when ignoring the linked arrivals. Buffers (in seconds) were set to b1 = 2 × revT ime + 400 and b2 = 30.

instance B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 Total time (s) 36 37 308 90 386 329 6 3 156 4 306 7 Search space (nodes) 0 0 2249 0 0 1016 0 0 0 2 2481 0

Model generalized to linked arrivals

Generalizing the matching.

Consider the matching problem modelled with the x ij variables and lying at the heart of the previous MIP formulation. We denote the underlying bipartite graph by G = (A ∪ I, D, E). An edge (i, j) belongs to E if the aforementionned criteria (chronology, maximum maintenance capacities and categories) are satisfied for the pair (i, j).

The presence of linked-arrivals (LA ⊆ A) alter the problem by introducing dependencies. Typically, the possible departures Γ (i) that can be assigned to an arrival i can now depend on the assignment made on an earlier departure j referred here as a linked departure. We formalize the problem as follow: For each vertex x i of D, we define Γ (x i ) = {Γ (a,yi) (x i ) | a ∈ Γ (y i )} so that each Γ (a,yi) (x i ) ⊆ Γ (x i ) represents the possible neigbhors of x i if a is matched to y i . THE DEPENDENT MATCHING PROBLEM (DMP) is to determine if there is a complete matching M in G so that for all i = 1, . . . , k:

-If x i is matched to b of B and y i is matched to a of A, we have b ∈ Γ (a,yi) (x i ).

-If y i is unmatched and x i is matched to b, we have b ∈ Γ (x i ). Theorem: DMP is NP-Complete Proof : It is easy to see that DMP belongs to the class NP. We show a reduction from the SAT [START_REF] Stephen | The complexity of theorem-proving procedures[END_REF] problem to DMP. Let ψ = C 1 ...C p be a boolean expression in conjunctive normal form and (x 1 , ..., x q ) the variables of ψ. We define the bipartite graph as follows:

-For each literal x i and each clause C j wherein x i appears, we create a gadget, named x i C j , consisting of a K 2,2 graph. Notice that this gadget has exactly two possible complete matchings: two parallel or two crossing edges. In the following, we will code the value x i true with two parallel edges and x i true with two crossing edges (see Fig. 4).

-For each clause C j with n literals, we create a gadget, named C j , consisting of a K n,1 graph, such as vertices on the left represent every literals in C j . In order to diffuse the value of a literal in every clauses, we will build a chain with conditions between every gadgets on the same literal. We add the following condition on the lower right vertex (vertex b in Fig. 5) of the literal gadget: if the edge connected with this vertex is the crossing one, we add a dependance with the lower left vertex (vertex 4 in Fig. 5) of the next literal gadget (corresponding to the next occurence of this literal in a clause) in the chain to remove the parallel edge from its neighborhood ( Γ (4) = {Γ (1,b) (4) = {c}, Γ (2,b) (4) = {d}} in Fig. 5). We also create another dependance in the same way for parallel edges. Notice that with these conditions, if one literal gadget is set to true or false, every other literal gadget on the same literal must be set to the same value to obtain a complete matching. Then in order to avoid the selection of literal without the proper value in the clause, we add a condition on the higher right vertex (vertex a in Fig. 5) of the corresponding literal gadget: if the clause C j has a positive literal x i (resp. a negative literal x i ) and the edge connected with this vertex is the crossing one (resp. the parallel one), we remove the edge between x i (resp. x i ) and C j in the clause gadget ( Γ (x 1 ) = {Γ (1,a) (x 1 ) = {C j }, Γ (2,a) (x 1 ) = ∅} in Fig. 5).

We now show that the DMP problem B with the conditions C has a solution if and only if ψ is satisfiable.

(i) Assume that ψ is satisfiable. For each clause C j choose a single literal y (y = x i or y = x i ) such that y is true. If y = x i (resp. y = x i ) then include the edge (C j , (x i , j)) (resp. (C j , (x i , j))) in the matching M . For each literal gadgets where y appears, set two parallel edges if x i is true, and two crossing edges if x i is true. Finally, if some literal gadgets aren't already set, we can choose to set them to the value true with two parallel edges. Thus, M is a complete matching. For each literal y, every gadgets with y is either set with two parallel edges or two crossing edges, but not both, depending of the value of y in ψ. And since we set one value for each clause gadget and this value stays the same on literal gadget, all conditions are respected.

(ii) Let M be a solution to the DMP problem. The truth value of x i is assigned as follows: if literal gadgets of x i are set to parallel edges (resp. crossing edges), x i is given the value true (resp. false). The conditions imply that the value of x i is well defined, x i and x i cannot both be true since all gadgets for x i must be set to the same type of complete matching. Since M is a complete matching, every vertex C j from clause gadget is matched with some vertex y. From the construction, the literal y is true and the clause C j is satisfied. Thus ψ is satisfiable.

Solution 1: An extended formulation.

A solution of the original problem can involve chains of events alternating arrivals and departures: i 1 , j 1 , i 2 , j 2 , . . . i c , j c using the same rolling stock unit. We denote by Ω the set of all such possible chains and suggest to replace the x ij variables of the previous MIP model by variables c k ∈ {0, 1} indicating if the k-th chain is used in the matching. A similar idea was also investigated by [START_REF] Thorlund | Team j10[END_REF][START_REF] Létocart | Team s20[END_REF]. The key idea is therefore that one can see DMP as a matching problem in a hyper-graph where all possible chains are represented (as edges). This would allow to preserve the matching structure of the previous MIP and give the rationale for simply replacing x ij by c k . Now in our problem, maintenances can be performed at different positions along the chain and this information must be represented. The k-th chain is denoted chain k = (i 1 , md 1 , mt 1 , j 1 , i 2 , md 2 , mt 2 , j 2 , . . . , i c , md c , mt c , j c ) where i x are arrival nodes, j x are departures nodes and md x (resp. mt x ) are 0/1 values denoting if a distance (resp. time) maintenance is performed between the two events i x and j x . All features (such as category, group category, length, remaining DBM/TBM at any stage) of a given chain are known. Note that the previous model was simply limited to chains (i, j) of size 2 represented by the x ij variables. The MIP model is modified by replacing the x ij variables with the appropriate c k variables in all constraints. Typically constraints (1), ( 2), ( 4), ( 6), (9) now becomes:

(1 ′ ) k∈Ω|i∈chain k c k ≤ w i ∀i ∈ A ∪ T (2 ′ ) k∈Ω|j∈chain k c k = y j ∀j ∈ D (4 ′ ) k∈Ω|catGroup chain k =g c k ≤ |jdList j |u gj ∀j ∈ J D, ∀g ∈ G (6 ′ ) ct jc ≥ c k ∀c ∈ C, ∀k ∈ Ω|cat chain k = c, (9 ′ ) l j ≥ k∈Ω|j∈chain k length chain k c k ∀j ∈ D
Since the distance/time maintenances operations are known on a chain, the daily maintenances limit can be stated on the c k variables. Typically for a pair (i r , j r ) in chain k , we simply set o irjr = md r + mt r . So again, we outline the fact that the exact days of maintenances are not represented in the pricing problem:

(13 ′ ) ∀d a ∈ {1, . . . ,
o irjr c k ≤ (d b -d a + 1) × maxM aint
As Ω is of exponential size, this model is meant to be solved using column generation techniques and Branch and Price. Each category leads to an independent path related problem. At most four choices are possible regarding the maintenance between an arrival i r and a departure j r : none, time, distance, both time and distance. Each choice can be modeled with a different arc between i r and j r (the four arcs are not necessarily present depending of the time available between i r and j r ). The objective function is expressed on the use of vertices and edges:

the cost for using an arrival vertex is related to the dual variable of (1 ′ )

the cost for using a departure vertex is related to the dual variable of (2 ′ )

the cost for using a maintenance arc between a given pair of arrival and departure (i r , j r ) is related to the dual variables of (13 ′ )

a constant cost related to the category itself is due to the dual variables of (4 ′ ), (6 ′ ), ( 9′ )

The problem is to find a shortest path in this graph satisfying two ressource constraints on the level of distance and maintenance which are reset depending of the type of arc used (indicating if some maintenance is performed). We believe it can be modeled as a MIP and do not go here into more details. Solution 2 : A greedy approach.

The currently implemented solution proceeds greedily by forbidding linked-arrivals as long as their corresponding linked-departure is not known. We simply solve the previous MIP model step by step so that chains of linked events have a length of at most 2i at the i-th iteration. After each step, potentially knew linked-departures are known (their features are known) and the linked-arrivals can now be included in the model. This is repeated as long as new arrivals are assigned.

The results presented in Table 2 shows that this approach is able to actually provide solutions that involve dependencies but that many departures remain uncovered. We strongly believe that the column generation model given previously has the potential to accurately deal with this issue. The size of the formulation of the column generation model might be actually smaller and require less memory as not all chains of size 2 might be needed to identify near-optimal solutions. The station can be modelled by an oriented graph G(V, E) where vertices are associated to gates, arcs to the transition of a ressource (with a potential transition time trT ime).

Track groups are typically represented by complete bipartite graphs. Some ressources such as platforms allow the trains to change direction. If such a move is possible on a ressource between two gates A and B, it is modelled in the graph with two additionnal vertices A ′ and B ′ as shown on Figure 6. Routing in the station is subject to forbidden time-windows expressing the usage of resources by other actors. This is expressed as forbidden time-windows on vertices i ∈ V . We denote by I i the set of such forbidden time-windows on vertex i : The following variables are used:

I i =
t i ∈ {0, . . . , H} -∪ |Ii| k=1 I k i : the time at which the train is entering vertex (gate) i.

-s i ∈ {k | (i, k) ∈ E} ∪ {i} : successor of i in the path (s i = i if i is not in the path).
The model2 is as follows:

(1)

t si = t i + d isi ∀i ∈ V -{b} (2) ALLDIFFERENT(s 1 , . . . , s n ) ∀i ∈ V (3) NOCYCLE(s 1 , . . . , s n ) ∀i ∈ V (4) t i ∈ {0, . . . , H} -∪ |Ii| k=1 I k i ∀i ∈ V
The constraint t si = t i + d isi can be stated with two intermediate variables y i and z i representing respectively the time to reach the successor of i (y i = t si ) and the distance to the successor of i (z i = d isi ). t si = t i + d isi is stated using two ELEMENT constraints (on a table of variables and on a table of constants i.e the distances): ELEMENT(y i , [t 1 , . . . t n ], s i ) ELEMENT(z i , [d i1 , . . . d in ], s i ) y i = t i + z i Strengthening Propagation. Since distances can be null (d ii = 0), the shortest path from the source a to all nodes is not actually propagated leading to poor filtering regarding the time-windows. Alternative models using the predecessor variables, min/max constraints to inform about the shortest path from a to all vertices and to any vertex to p can be added and help filtering vertices that can not be part of a feasible path.

Branching. The route is built in a constructive manner selecting vertices from a along the shortest (in number of edges) route leading to b, propagation at each step helps getting rid of invalid successors.

Overall routing strategy.

Our routing component simply computes routes one by one from the earliest to latest arrival and cancels the arrival if the CP solver is unable to identify a feasible route for the train. A train might require several routes if intermediate resources need to be reached (typically for maintenances or simply for waiting the proper time to be sent to the departure platform). Many heuristics can be added here for the selection of trains such as the ones proposed by [START_REF] Geiger Martin | Team s1. In Challenge ROADEF/EURO 2014[END_REF][START_REF] Buljubašić | Team s18. In Challenge ROADEF/EURO 2014[END_REF].

Conclusion and future work

A natural problem-decomposition was presented and require limited dedicated algorithms. The router part (see Figure 1) is the only part where a specific strategy has to be designed. We showed that a large part of the problem (larger than expected) related to assignment of arrivals/departures, platforms and in particular the daily maintenances limits can be efficiently handled with MIP. A key idea for integrating efficiently linkedarrivals in such an approach is suggested using a column generation framework. This remains to be evaluated experimentaly but the column generation model can turn out to be smaller in size and should handle linked-arrivals very accurately.

The proposed methodology only holds if the routing can be considered separatly from the assignment. This was the case on the first problem formulation but remains an open question now that track group conflicts are considered as hard constraints. In any case a stronger relaxation of the routing should be included in the assignment part of the model (the two buffers b 1 and b 2 included are a very crude relaxation).

Fig. 1 .

 1 Fig. 1. Overall view of the proposed solution approach.

  Let's denote by M (d a , d b ) the set of pairs (arrival, departure) (i, j) that require o ij ∈ {1, 2} maintenance operations (distance and/or time) and that takes place in a time-window of days [d a , d b ]. A pair (i, j) belongs to M (d a , d b ) if arrT ime i and depT ime j belong to the interval of days [d a , d b ].

4 .

 4 Modeling the objective function. We can model all elements of the objective function to the exception of the platform usage costs. The uncovered arrivals/departures and unused initial trains costs are easily modeled based on the y and w variables.

The

  Dependent Matching Problem (DMP): Consider a bipartite graph G(A, B, E) and let D = {x 1 , . . . , x k } ⊆ A and C = {y 1 , . . . , y k } ⊆ B be two subsets of respectively A and B so that |C| = |D|. D stands for Dependencies and C for Conditions.

Fig. 3 .

 3 Fig. 3. An instance of DMP (left), a complete matching M1 solution (middle), a complete matching M2 which do not satisfy the dependencies.

Fig. 4 .

 4 Fig. 4. Gadget for the literal xi (left), set to true (middle) and false (right).

Fig. 5 .

 5 Fig. 5. Literal gadget set to false (upper left), diffusion to its clause gadget Cj (upper right), diffusion to the next literal gadget (lower left).

Fig. 6 .

 6 Fig. 6. An exemple structure for performing a reverse move.

Fig. 7 .

 7 Fig. 7. Additionnal vertices expressing crossing issues in track groups.

  nbDays}, ∀d b ∈ {d a , . . . , nbDays}, (ir,jr)∈c k |k∈Ω ∧ (ir,jr)∈M (da,d b )

Table 2 .

 2 Results obtained on the first part when handling greedily (step by step) the linked arrivals.

	instance	B1	B2	B3	B4	B5	B6	B7	B8	B9 B10	B11	B12
	Cost (s) 169365 103547 42608 150372 182576 164165 73768 73768 253704 79641 322471 261639
	Uncovered Arr	0	0	0	1	7	2	0	0	0	14	44	66
	Uncovered Dep	86	81	14	52	80	62	55	55	121	51	195	149
	Total time (s)	47	54 487	250 1311 1229	26	22	453	10	101	27
	4 Routing
	4.1 Modeling the station and hardness of routing

Experiments were performed on a MacBook, running MacOSX 10.5.8, running on a

2GHz Intel Core 2 Duo CPU, with 2 GB 667 MHz DDR2 SDRAM.

We assume that NOCYCLE allows single loops of type si = i.


Jingxian Lu

Christophe Dousson

Benoit Radier
email: benoit.radier@orange-ftgroup.comfrancine.krief@labri.fr

Francine Krief

Towards an autonomic network architecture for self-healing in telecommunications networks

Keywords: Self-healing, Model-based diagnosis, Causal graph

We propose a solution to achieve a global self-healing process based on local knowledge of each network element using a model-based diagnosis approach. Then we illustrate this solution in an IMS platform and detail the corresponding causal graph and diagnosis process.

Introduction

With the increased complexity of system installation, maintenance, configuration and recovering, an increase in intelligence in the system is required and the concept of autonomic system will come in.

An autonomic system manages itself and becomes autonomic without human interventions. It has four distinct objectives: Self-configuring, Self-healing, Selfoptimizing and Self-protecting. The closed control loop proposed by IBM is implemented by autonomic managers, which control managed resources [START_REF]IBM: An architectural blueprint for autonomic computing[END_REF]. With the development of autonomic system, the identifying, tracing, and determining of the cause of failures and the repairing action became more and more important to the operators in order to satisfy the needs of customers. Thus, the purpose of our studies is to introduce self-healing in operator networks.

Model-based diagnosis for self-healing

The causal graph is a general, intuitive and effective representation to model the consistent causalities governing the function of the system to monitor [START_REF] Poole | Normality and faults in logic-based diagnosis[END_REF]. It is supposed to be an acyclic graph that limits the expressiveness of the representation, however, this restriction is important to reduce the complexity of reasoning [START_REF] Console | A theory of diagnosis for incomplete causal models[END_REF]. The causal graph structure in this paper is based on five kinds of nodes and two kinds of arcs which are shown in Figure 1.

The diagnosis process consists in collecting all the current observations, and then, determines all the primary causes which could explain them. To achieve this process, each node has a unique state as following: Guilty (a test result or an active alarm or a proved cause), Innocent (a test result or an absent alarm or an exonerated cause), Suspect (a test to run or an unknown cause to investigate), Unknown (the default state).

Then, applying the causal propagation rules that follows until stability of the node states over the whole graph gives the diagnosis and repairing action: R1 : if a node is guilty, then all its "should cause" successors are guilty. R1' : if a node is guilty, at least one of its predecessors is guilty. R2 : if a node is innocent, then all its "should cause" predecessors are innocent. R2' : if all the predecessors are innocent, the node is innocent. R3 : if a node is guilty, then all its predecessors become suspect. R4 : if a node is suspect, then all its "should cause" predecessors are suspect. R4' : if a node is suspect, then all its "should cause" successors are suspect. R5 : if a test node is suspect, then the test is performed. Depending of the result, the state becomes guilty or innocent. R6 : if a repairing action is guilty, then it should be processed.

Application to IMS service

The autonomic network architecture used a distributed Knowledge Plane (KP) to interconnect different network elements [START_REF]ANA: Autonomic Network Architecture[END_REF][START_REF] Mbaye | Autonomic Networking : The Knowledge Plane[END_REF]. We define the ARM (Autonomic Resource Manager) to manage network element types and ANM (Autonomic Network Manager) to manage the ARM plane and to control the network [START_REF] Strassner | FOCALE -A novel autonomic networking architecture[END_REF][START_REF] Strassner | The design of an autonomic communication element to manage future internet services[END_REF].

We illustrate the previous diagnosis process by applying our approach to the case of a SIP-based VoIP service delivered over an autonomic management architecture (Fig. 2). Let's suppose that there is the following problem: I-CSCF sends the user's IP Multimedia Public Identity (IMPU) to SLF which can not find the corresponding user information.

According to the above problem, the global diagnosis relies on the causal graph shown in the Figure 3. We have divided it into four regions. Each region is a subgraph which corresponds the ARM or the ANM. The operators validate each subgraph in each ARM to decide if it could be involved in a global diagnosis when the equipment are integrated in the network. ARM allows the administrator to understand in a smooth way how the equipment could interact with contextual information. If the ARM is compatible with the operators' expectation, the subgraph is integrated in the KP of ARM and ANM by connecting the different subgraphs together. As a consequence, a primary cause in a subgraph could become an intermediate cause in the implicit global graph. It's why we need to standardize the interface (showed in the Fig. 3 Running the diagnosis process (Section 2) on the example of IMPU alarm is procedure in the following steps by applying the specified rules:

as arrows between

1. (R3): As "Alarm unknown IMPU in SLF" is guilty (alarm present), then "IMPU unknown in SLF" and "IMPU barring in SLF" are suspect. 2. (R4): "IMPU unknown in SLF" is suspect, then its predecessor "Unknown client in SLF (No account)" is set to suspect. 3. (R4): "Unknown client in SLF (No account)" is suspect, then "Unknown client in HSS (No account)" is set to suspect. 4. (R4'): "Unknown client in HSS (No account)" is suspect, then "Unknown client in AS (No account)" and "IMPU unknown in AS" are set to suspect. 5. (R2): "Alarm unknown IMPU in AS" is innocent (no alarm presents), then "IMPU unknown in AS", "Unknown client in AS (No account)" and "Unknown client in HSS (No account)" are set to innocent. 6. (R2'): Then "Unknown client in SLF (No account)" is set to innocent. 7. (R2,R2',R4): "IMPU barring in SLF" is suspect, and the same process will be launched as steps 2 to 6. Then "Client barring in SLF" is set to innocent. 8. (R1'): As "Unknown client in SLF (No account)" and "Client barring in SLF" are innocent, then "SLF no synchronization" is guilty. 9. (R6): The repairing action "Resynchronization SLF" is performed.

Conclusion

We have proposed a solution to achieve a global network diagnosis with local knowledge of each network element using causal graphs validated by the operators. It permits to connect different necessary processes (subgraphs) to achieve an end-to-end diagnosis and to perform the repairing actions.

The further research should aim to specify the interface between elements or domains which corresponds to the connection between different local diagnosers and consider how the performance of the global diagnosis process could be adapted to our autonomous architecture by defining when and how to send information. Then, we will implement our self-healing approach on an operational IMS platform.

Fig. 1 .

 1 Fig. 1. Nodes and arcs of causal graph

Fig. 2 .Fig. 3 .

 23 Fig. 2. The case of a SIP-based VoIP service


Ahlem Bouchahda-Ben

Nhan Tekaya

Adel Thanh

Faten Bouhoula

Labbene-Ayachi

Ahlem Bouchahda-Ben Tekaya

Nhan Le Thanh

Adel Bouhoula

Faten Labbene-Ayachi

An Access Control Model for Web Databases

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

Introduction

Nowadays, Web applications depend more and more on the back-end database to provide much more functionalities. Containing valuable business information, these systems are highly interesting to attackers and special care needs to be taken to prevent any malicious access to this database layer. Access control is the primary means of attack prevention for databases. But as long as web databases cannot identify their real users, they cannot supply them with proper authorization. Database Views are another means of unauthorized access restriction as they can define the only part of a database relevant to a user. But for web databases they are useless. In fact, web applications are run from the user's browser windows. The browser does not directly connect to the database, but instead transfer a request to a web server who processes the request and if an access to the underlying database is needed, transfers it to the application server which performs a transaction to the database. It implies that the database does not identify the real user who accesses it and so, traditional identity-based mechanisms for performing access control are useless for web databases. Further, a DBMS can not handle users who access it indirectly via the application server, no user-based access control can be applied since the only recognized user is the user of the application server and for most of the web applications this is the user with very high privileges. Databases can no longer differentiate between transactions of different application users. The principle of minimal privilege is violated. It is impossible to authorize the web application user with appropriate privileges at the database level: all application users have access to the same data. Restrictions on what authenticated users are allowed to do are not properly enforced. Attackers can exploit these flaws to view sensitive data, or use unauthorized functions. So, no more fine-grained access control to the database exists and authorization can be provided only at the application level. To protect web databases from attacks, access control policies should be based on a strong model that is implemented by the DBMS. Checking for authorization should be done on every attempt to access secure information. For that end, we extend the NIST RBAC (Role-Based Access Control) standard [START_REF]American national standard for information technology, role based access control[END_REF]. The central idea of RBAC + is including the concepts of application, application profile and sub-application session when controlling the access to web databases. The application profile is necessary to track the user behavior throughout a whole session and mainly to prevent business logic violation attacks by enforcing access control. The RBAC + monitors transactions issued by users and malicious transactions are viewed as intrusion behaviors. If a malicious transaction is identified, the RBAC + cancels the transaction before it succeeds, thus minimize damage caused by malicious transactions. The rest of the paper is organized as follows. We present the related work and discuss it in Section 2. In section 3, we present an overview of our model. We define formally and detail our model in Section 4. Access control policies are presented in section 5. We conclude our work and present future work in Section 6.

Related Work

The problem of access control to databases accessible over the web is very important one. This problem is well known to the web application developers and security consultants, but little existing work has addressed it. Gertz et al. in [START_REF] Gertz | Security Re-engineering for Databases: Concepts and Techniques[END_REF] pose this problem and presented some fundamental concepts and techniques that help administrators and security personnel to gradually evaluate and improve the security of a database. Also, Roichman in [START_REF] Roichman | Intrusion prevention and detection for web databases[END_REF] proposed a method that uses the databases' built-in access control mechanisms enhanced with Parameterized Views and adapts them to work with web applications in order to prevent intrusions. He defines also the concept of session vector to represent the application fingerprints used in an application session in order to detect intrusions. This concept is similar to our application profile concept with the difference that the session vector represents the session profile while the application profile represents an execution way of the application. Beyond these two approaches, to protect web databases from attacks of malicious users, two main approaches exist. The first consist on using ad hoc tools specifically oriented to the detection of specific kinds of attacks like SQL injection [START_REF] Halfond | A classification of sql-injection attacks and countermeasures[END_REF]. The second consists on using Intrusion Detection Systems (IDSs). [START_REF] Vieira | Detection of malicious transactions in dbms[END_REF] presents a database IDS that uses the profile of the transactions implemented by database applications (authorized transactions) to identify user attempts to execute unauthorized transactions. Although we believe that database IDSs should play an important role in database security, we have to point out that the web application's access to databases remains untraceable. Further, an IDS can not overcome the absence of web database internal access control and the uselessness of views as a means of access restriction. Moreover, with the assumption that the attack does not go unnoticed, IDSs focus on detecting attacks after the malicious user has accessed the DB with all the damage it could cause. What we propose is strengthening access control and continuously monitor users. Consequently, the majority of attacks can be stopped from the access control stage and the IDS will be used to detect attacks that have escaped the access control stage. Intrusion detection without enforcing access control is not as efficient and effective. IDS is a complement but can not, alone, protect DB from attacks.

Overview of our approach

As accesses to the data occur through several layers, starting with a person then the application, which, in turn, performs operations on the database, correlating anomalous behaviors with a person is not a trivial task. To address these problems, we propose RBAC + , an extension of RBAC able to detect malicious users and stopping the attack before it succeeds. Assuming that the database management system (DBMS) has an RBAC model in place, the key idea of our approach is as follows. We create application profiles that represent all the possible execution paths of the application. Given the permissions necessary to the execution of an application and the set of roles that the underlying database user (DBU) is authorized for, we calculate for each pair (application, DBU) the subset of roles to activate in a web user session, called sub-application session. It is called so because, in the context of a web application, a web user session is included in a database session. Hence, a sub-application session contains only the permissions really needed to fulfill exactly the tasks it was created for, and so we take advantage of all RBAC assets such as least privilege and separation of duty. A sub-application session allows to the DBMS distinguishing between web users working with the database. It will also allow distinguishing between the requests of different web users that belong to the same database session. When the web user logs in, the SQL queries that he submits are associated with a database session, an application and the underlying database user that issued them. All queries belonging to a sub-application session must match an application execution path else the access is denied because the action to be executed is illegitimate. We assume that the user is identified at each tier. The recent tendency in the architecture of web applications allows preserving the identity of the real user through the middle tier As example, Oracle9i introduced n-tier authentication [START_REF]Oracle corporation: Oracle9i database concepts release 2[END_REF], i.e., that is "lightweight session". Now, when an employee wants to attack enterprise resources, he, for example, can submit an SQL injection attack. But because his database privileges are limited only to legitimate actions, an SQL injection will be entirely mitigated or at least, its effect is strongly limited. The importance of our solution is that it enforces access control based on business application logic rather than primitive reads and writes. A user's ability to access and manipulate data is typically dependent of the application function the user executes thus reducing drastically attacks against databases and in particular, business logic violation attacks because an action may be legitimate on its own but illegitimate in the context of a whole session. Take the example of an online retailer application. If the intruder can submit the insert statement into the Orders table without submitting an insert into the Credit Card table, then he can buy goods without paying. This business rule violation can be detected only at the session level since each statement by itself is a legitimate one. Databases cannot prevent them because the existing database access control can grant or revoke access to resources only according to the accessor identity/role. It cannot rely on the business logic of an organization. [START_REF] Halfond | A classification of sql-injection attacks and countermeasures[END_REF] The RBAC + Core Model

We, now, introduce a rigorous definition of the model. The purpose is to provide a comprehensive definition of the components, thus including all the aspects of the model. The general structure of the model is illustrated in Figure 1. We use the graphical representation adopted in RBAC. In particular, APPS, AP and SASES represent the sets of applications, application profiles and sub-application session respectively.

Application Profile

An application has many possible ways of execution each of which is called an application profile representing one valid execution sequence of the application (sequences of selects, inserts, updates, deletes). It consists of a sequence of nodes such that from each of its nodes there is an edge to the next node in the sequence. Each node in the path represent a SQL statement. It has one start node and one end node where the application execution starts and finishes, respectively. The other nodes in the path are called internal nodes. An application profile may include cycles representing the repetitive execution of sets of commands. Obviously, the set of application profiles must cover the different database application functionalities. To that end, we choose to analyze the application code which implicitly contains a policy that allows for distinguishing legitimate and malicious queries. Also, the source code contains enough information to infer models of the expected, legitimate sequences of SQL queries generated by the application.

Definition 1. (Application Profile). An application profile AP is a binary vector with the length equal to the number of permissions (PRMS) in the DBMS

, where the i th bit is 1 if the application profile needs the permission p i , else bit i is 0.

p i ∈ P RM S.
We also define :

-AAP : AP P → AP , the mapping of an application onto its corresponding application profiles. Formally, AP P prof iles(app) = {ap ∈ AP |(ap, app) ∈ AAP }. -RAP ⊆ ROLES×AP , a many-to-many mapping Role-to-Application profile Assignment relation. -AP roles : AP → ROLES, the mapping of an application profile to a set of roles. Formally, AP roles(ap) = {r ∈ ROLES|(r, ap) ∈ RAP }. -RAA ⊆ ROLES × AP P , a many-to-many mapping Role-to-Application Assignment relation. -AP P roles : AP P → ROLES, the mapping of an application to a set of roles. Formally, AP P roles(app) = {r ∈ ROLES|(r, app) ∈ RAA}.

Note that, AP roles(ap) ⊆ AP P roles(app) with {ap ∈ AP |(ap, app) ∈ AAP }

Sub-application session

An application session is composed of all the transactions that an application runs on behalf of all its users. A sub-application session (SASES) is the subset of transactions related to one user. Formally, we define :

Definition 2. (Sub-application session) We define:

app sas : AP P → 2 SASES . The mapping of an application onto a set of sub-application sessions. -session sas : SESSION S → 2 SASES . The mapping of a session onto a set of sub-application sessions.

Permissions

In our model, permissions are associated with roles and with application profiles. Applications are then associated with the appropriate roles based on the set of permissions assigned to application profiles.

Definition 3. (Permissions) The set of permissions PRMS is defined as P RM S = 2 (OP S×OBJ) . We also define:

-P AA ⊆ P RM S × AP , a many-to-many mapping Permission-to-Application profile Assignment relation. -AP perms : AP → 2 P RM S , the mapping of an application profile onto a set of permissions. Formally, AP perms(ap) = {p ∈ P RM S|(p, ap) ∈ P AA}.

Users

Each user is associated with a set of applications he/she is authorized to execute.

Definition 4. (Users) We define:

-AA ⊆ AP P S × U SERS, a many-to-many mapping application-to-user assignment relation. -U SER AssignedApps : U SERS → 2 AP P S , the mapping of a user to a set of applications. Formally, U SER AssignedApps(u) = {u ∈ U SERS|(app, u) ∈ AA}.

Sessions

When a user logs in, a new session is activated and a number of roles are selected to be included in the session role set. Formally: Definition 5. (Sessions): We define:

session user: SESSION S → U SERS, the mapping from a session s to the user of s. -session roles : SESSION S → 2 ROLES , the mapping of session s onto a set of roles. Formally: session roles(s) ⊆ {r ∈ ROLES|(session U ser(s), r) ∈ U A}. -session applications : SESSION S → 2 AP P S , the mapping of session s onto a set of applications. -avail app roles : (SESSION S, AP P S) → 2 ROLES , the mapping of a session and an application onto a set of roles. Formally, avail app roles(s, app) ⊆ {r ∈ ROLES|r = session roles(s) ∩ app roles(app)} -avail app prms : (SESSION S, AP P S) → 2 P RM S , the permissions available to an application in a session. Formally, avail app prms(s, app) = ∪ r∈avail app roles(s,app) assigned permissions(r). An sql query is a set of permissions. the above function is repeated as many permissions as the sql query requires permissions to be executed.

To support a differentiated and adequate protection of the information stored in a database made available over the web, the access control policies must be flexible enough to support a spectrum of web-based applications, such as e-commerce web applications or health care web-applications. To that end, we specify two access control policies. The first consists on monitoring all the SQL statements submitted by an user and the second consists on monitoring only special statements that we call "critical points". A critical point may be an SQL statement manipulating sensitive data or an SQL statement modifying data (insert, delete, update). When a malicious transaction is detected, the transaction is rolled back.

Policy 1

Under this policy, every command executed must match a profile. When the first command of the transaction is executed the tool searches for all the application profiles starting with that same command, which are marked as candidate profiles for the current transaction. The next command executed is then compared with the second command of these candidate profiles. If it corresponds, the access is granted, either, access is denied. Only those who match remain candidate profiles. This process is executed over and over until the transaction reaches its end or there are no more candidate profiles for that transaction. In this latter case the transaction is identified as malicious. In practice, to detect malicious transactions, we implement the following generic algorithm over the application profiles:

1: while (1) do 2:

for each new command submitted do obtain list of authorized application profiles starting with this command for each valid (authorized) trans. for the user do 7:

if the current command represents a valid successor node in the application profile then 8:

the command is valid

Fig. 1 .

 1 Fig. 1. Core RBAC +

3 :

 3 if user does not have any active transaction then {command is the 1st command in a new transaction } then 4:

  ⟨U, is, app, p, o⟩ ∈ U SERS × SASES × AP P S × OP S × OBJ. ar can be satisfied if (p, o) ∈ avail app prms(s, a) and is ∈ session sas(s).

	5 Access control policies
	Application profiles representing valid application execution paths are used to
	detect unauthorized SQL statements, which are seen as invalid sequences of SQL
	commands.
	Definition 6. (Authorization control function): An access request ar is a tuple
	ar =

end for 18: end while

Policy 2

Under this policy, the SQL statements submitted by the user called user context are stored and the access is granted until he submits a critical point. In this case, the tool searches an application profile that corresponds to the user context. If one or more corresponding application profiles are found, the access is granted. Otherwise, the transaction is rolled back.

Conclusion and Future Work

Database security problems seriously persecute web-based applications that rely on web databases storing invaluable data. In this paper we have presented RBAC + , an extension of the RBAC model addressing access control requirements for RBAC-administered web databases. We do not only monitor DB users to detect potential attacks, but timely stop the attacks when they are detected to minimize losses caused by the attacks. As future work, we will focus on the implementation of a prototype of the proposed system and its experimentation against a variety of simulated intrusions to prove its effectiveness and efficiency.


Anh Thu

Phan Ho

Bao An

Hoang Mai

Wadih Sawaya

Patrick Bas

Patrick Bas Authentication

A-T Phan Ho

B-A Mai Hoang

AUTHENTICATION USING GRAPHICAL CODES: OPTIMISATION OF THE PRINT AND SCAN CHANNELS

Keywords: Authentication, Hypothesis testing, minimax game, print and scan models

published or not. The documents may come

INTRODUCTION

The authentication of printed materials, which consists in bringing a forensic evidence that materials are genuine, is a huge challenge nowadays and can be used to prevent forgeries of valuable items, such as identity documents, or products, such as drugs by securing their associated packages. Authentication can be done either 1) by characterizing the "fingerprint" of the package, for example by recording the random patterns of the fiber of the package or the paper [START_REF] Haist | Optical detection of random features for high security applications[END_REF], but such a system is practically heavy to deploy since each product needs to be linked to its high definition capture stored in a database, or 2) by relying on the degradation induced by the interaction between the materials and a physical process such as printing, marking, embossing or carving. Since both the physical process and the matter are stochastic, the interaction between the two entities can be considered as a Physically Unclonable Function (PUF) [START_REF] Suh | Physical unclonable functions for device authentication and secret key generation[END_REF] that cannot be reproduced by the forger and can consequently be used to perform authentication.

We study here the authentication system first proposed by [START_REF] Picard | Digital authentication with copy-detection patterns[END_REF], and based on the use of a printed binary code to perform authentication. The whole system is depicted in Figure 1: the legitimate source prints an original secret random code x N (a binary matrix of N elements) on a document or a package and the receiver scans it to perform authentication as the opponent may have manufactured his forgery and generates his copied code. In order to generate this copy, the opponent observes the printed and scanned version y N of x N (step 1 in Figure 1) and extracts a binary code xN (step 2). This step is unavoidable due to the fact that an industrial offset printer can only use binary input to generate dots of inks. With this second print and scan process, the copied code z N generated by the opponent (step 3) has a different distribution than the initial printed and scanned code y N . This distinction drives the authentication system (step 4).

PRINT AND SCAN MODELS AND AUTHENTICATION

Principle of the authentication system

The legitimate receiver observes a gray level code o N , and we assume that the random observed sequence (O N | x N ) (conditioned to the secretly shared binary code x N ) is independent and identically distributed (i.i.d. ). The Neyman Pearson test is expressed as:

L = log P (o N | x N , H 1 ) P (o N | x N , H 0 ) H1 ≷ H0 λ. (1)
where H 0 is the hypothesis that the observed sequence comes from an original source with distribution P (O | x , H 0 ), and

Asymptotic expressions of α and β

Before considering our optimization game in section 3, we recall [START_REF] Ho | Document Authentication Using Graphical Codes: Impacts of the Channel Model[END_REF], where a method to compute reliably type I error probability α (the probability to consider a genuine code as a copy) and type II error probability β (the probability to not detect a copy) has been presented. Contrary to the Gaussian approximation of L which provides inaccurate error probabilities when the threshold λ in (1) is far from its the mean value, this solution uses the Chernoff bound [START_REF] Dembo | Large deviations techniques and applications[END_REF] as very small error probabilities of type I and II may be desired [START_REF] Gallager | Information theory and reliable communication[END_REF]. For i.i.d. random sequences, the distribution of the random variable L = i ℓ(O i /x i ) depends on the origin (H 0 or H 1 ) of the observed code o N , and for any real number s, the semi-invariant moment generating function of each ℓ(

O i /x i ) is µ ℓ (s; H j ) = x=0,1 µ ℓ/x (s; H j ) = x=0,1
log E O|x,Hj e sℓ(O/x) . Type I and II errors may then be tightly expressed for sufficiently large N (with approximately N/2 white and N/2 black dots), as:

α = Pr(L ≥ λ | H 0 ), → N →∞ 1 s0 √ N πµ ′′ ℓ (s0; H0) e N 2 [µℓ(s0; H0)-s0µ ′ ℓ (s0 ; H0)] .
(2) and

β = Pr(L ≤ λ | H 1 ), → N →∞ 1 |s1| √ N πµ ′′ ℓ ((s1; H1) e N 2 [µℓ(s1; H1)-s1µ ′ ℓ (s1 ; H1)] .
(3) where µ ′ ℓ (s j ; H j ) and µ ′′ ℓ (s j ; H j ) are respectively the first and second derivatives of µ ℓ (s; H j ) at value sj such that

N 2 µ ′ ℓ (s; H j ) = λ.

Models for the print and scan channel

In this paper we use two different families of distributions to model the print and scan channel, but the general methodology of this paper does not depend on the model and can still be applied. The first one is the Generalized Gaussian distribution which has been chosen because it can model both sparse and dense distributions. The second one is the Lognormal distribution since it has been previously shown by Baras and Cayre [7] that this distribution is an accurate model of the print and scan channel. Note that other print and scan models based on the gamma transfer function or additive noise with input dependent variance can be found in [START_REF] Lin | Distortion modeling and invariant extraction for digital image print-and-scan process[END_REF].

The distribution T V |x (v | x) modeling the physical device, i.e. the association of a printer with a scanner, can be written as (for x taking binary values):

-For the Generalized Gaussian distribution:

p(v | x) = b 2aΓ(1/b) e -(|v-m(x)|/a) b , (4)
where Γ(•) is the gamma function, m(x) the mean, and parameters a can be computed for a given variance

σ 2 (x) = var[V | x]: a = σ(x)Γ(1/b)/Γ(3/b). (5)
The parameter b is used to control the sparsity of the the distribution, for example when b = 1 the distribution is Laplacian, b = 2 the distribution is Gaussian, and b → +∞ the distribution is uniform.

-For the Lognormal distribution:

p(v | x) = 1 vs(x) √ 2π e - (log v-m(x)) 2 2s 2 (x) , (6)
where log(V | x) has a Gaussian distribution with mean m(x) and variance s 2 (x). The mode of the distribution is M (x) = e m(x)-s 2 (x) , and the variance is given by σ 2 (x ) = (e s 2 (x) -1)e 2m(x)+s(x) 2 .

To provide values within [0, . . . , 255] to model a scanning process, we quantize and truncate distributions ( 4) and [START_REF] Gallager | Information theory and reliable communication[END_REF] . Each channel is parametrized by 4 parameters, 2 per each type of dots. For the Generalized Gaussian distribution the parameters are m(0) and σ(0) for black dots and m(1) and σ(1) for white dots. The Lognormal distribution can be parametrized by the standard deviations σ(0) and σ(1) and the modes M (0) and M (1) respectively for black and white dots.

Figure 2

OPTIMAL CONFIGURATIONS FOR AUTHENTICATION

This authentication problem can be seen as a game where the main goal of the designer of the authentication system is, for a given false alarm probability α, to find a channel that minimizes the probability of misdetection β.

Practically this means that the channel can be chosen by using a given quality of paper, an ink of appropriate density and/or by adopting a given resolution. For example if the le- gitimate source wants to decrease the noise variance, he can choose to use oversampling to replicate the dots, on the contrary if the legitimate source wants to increase the noise variance, he can use a paper of lesser quality. It is important to recall that because the opponent will have to print a binary version of its observation, and because a printing device at this very high resolution can only print binary images, the opponent will in any case have to print with decoding errors after estimation X.

We analyze two scenarios described below: • The legitimate source and the opponent have identical printing devices (by devices we mean printer, ink, paper, scanner), practically this means that they use exactly the same printing setup. In this case the legitimate source will try to look for the channel C such that for a given α, the legitimate party will have a probability of misdetection β * such that:

β * = min C β(α). (7)
In this case, the opponent is defined to be passive. • The opponent can modify its printing channel C o , practically it means that he can modify one or several parameters of his printing setup. Actually, we assume that he changes the variance of its noise given that it will be the most efficient way for him to confuse the receiver. The opponent thus tries to maximize the probability of misdetection by choosing his adequate printing channel, whereas the legitimate source will adopt a printing channel C l which minimizes the probability of misdetection. We end up with so-called a min-max game in game theory, where the optimal β * is the solution of:

β * = min C l max Co β(α). (8)
In this case the opponent is active since he tries to adapt his strategy in order to degrade the authentication performance.

For the Generalized Gaussian model and the Lognormal model, we assume that, respectively, the means m(0) and m(1), and the modes M (0) and M (1) are constant for all the players in the different channels (which implies that the scanning process has the same calibration for the two types of images). We assume also that variances of black and whites dots are equal at each channel and denote them σ 2 m and σ 2 o for main and opponent respectively..

Passive opponent

Here the opponent undergoes a channel identical to the main channel. The only parameter of the optimization problem ( 7) is consequently σ m .

Active opponent

In this setup the opponent can tune his variance σ 2 o to confuse the receiver with the higher β. Figure 5.a shows the evolutions of β w.r.t σ o for different σ m when a Generalized Gaussian channel is assumed. We can see that in each case it's in the opponent interest to optimize his channel.

Analysis

When facing a passive opponent, it is not surprising to notice that in each case β is important whenever σ m is very small, i.e. when the print and scan noise is negligible hence the estimation of the original code by the opponent is easy; or very large; i.e. when the print and scan noise is so important that the original and forgery become equally noisy. The legitimate source will consequently avoid a channel that generates noise of very small or very large variance.

For an active opponent, the active scenario offers a saddle point satisfying (8) either for Generalized Gaussian or Lognormal distribution. This means that even if the adversary owns ideally perfect print and scan devices (σ o → 0, o N = x N ), it is not to his advantage to use it since the authentication is still efficient due to the decoding errors he will create by generating the binary code XN .

Another general remark is to notice that the optimal opponent parameters are very close to the optimal parameters of the passive scenario, which means that the adversary has little room to maneuver when choosing his best attack (see Figures 4 (a) and5 (b,c)) and nearly no room when the noise is close to uniform (b = 6).

For Generalized Gaussian distribution, it is important to notice that for distributions of same variance, dense distributions yields to better authentication performance than sparse distributions for both scenarios (see Figures 4 (a) and 5 (b)). This is due to the fact that a distribution close to uniform tend to create a bigger overlap between the two decision regions than a sparse distribution that will generate codes mainly lying in the original one.

For the Lognormal distribution we can notice that the authentication performances are similar for different values of modes, both for a passive and an active opponent (see Figures 4 (b) and 5 (c)). However, the larger the difference, the larger the optimal standard deviation, which means that it is in designer strategy to force the opponent to generate decoding errors in this case.

CONCLUSIONS AND PERSPECTIVES

In this paper we have proposed to cast the problem of authentication of printed documents using binary codes into an optimization game between the legitimate source and the opponent, each player potentially tries to select the best print and scan channel to minimize/maximize his authentication performance. This game was possible by considering accurate computations of the type I and type II probability errors and by using additive stochastic processes to model the print and scan channel.

We have shown that for both the Generalized and Lognormal distributions the game can be tractable, and that it is in the interest of the legitimate source to adopt a channel which is close to the uniform distribution.

Our future work will consist in finding, using information theoretic arguments, what can be the best additive channel for this setup.

Fig. 1 .

 1 Fig. 1. Principle of authentication using graphical codes.

Fig. 2 .

 2 Fig. 2. Example of a 20x20 code which is printed and scanned by an opponent following a Generalized Gaussian distribution for b = 1 (first row), b = 2 (second row) and b = 6 (third row). Columns (a), (b), (c), (d) represent respectively X N , XN , Y N and Z N . Main and opponent channels are identical with m(0) = 50, m(1) = 150, σ(0) = 40, σ(1) = 40 .

Fig. 3 .

 3 Fig. 3. Representation of the print and scan model for the black dots (on the left) and the white parts of the paper (on the right) for different standard deviations σ(0) = σ(1) = σ with M (0) = 70 and M (1) = 150 for the Lognormal distribution.

Figure 4 .MFig. 4 .

 44 Fig. 4. Evolution of the probability of non detection w.r.t the standard deviation of the channel (α = 10 -6 ) for the Generalized Gaussian distribution (a), and Lognormal distribution (b).

Figures 5 .

 5 b and 5.c shows the evolution of the best opponent strategy max σo β w.r.t σ m . By comparing it with Figure4, we can see that the opponent's probability of non detection can be multiplied by one or several orders of magnitude for the Generalized Gaussian distribution (×10 6 for b = 1, ×10 5 for b = 2) or for the Lognormal distribution (×105 for each mode separation) but stays the same when the distribution is close to uniform (b = 6).

MFig. 5 .

 5 Fig. 5. Evolution of opponent strategy β for the Generalized Gaussian distribution for b = 2 (a), and the best opponent strategy max(β) w.r.t the standard deviation of the channel (α = 10 -6 ) for the Generalized Gaussian distribution (b) and the for the Lognormal distribution (c).

ACKNOWLEDGEMENTS

This work was partly supported by the National French project ANR-10-CORD-019 "Estampille".


Bao An

Mai Hoang

Wadih Sawaya

Patrick Bas

IMAGE MODEL AND PRINTED DOCUMENT AUTHENTICATION: A THEORETICAL ANALYSIS

This paper combines the principles of statistical estimation and hypothesis testing to analyze the impact of parameter estimation on an authentication system based on graphical codes. The studied authentication system uses the fact that a code, once printed, undergoes a stochastic and non invertible alteration. A statistical test applies a likelihood ratio between the model of the authentic printed and scanned image and the model of the reproduced one, with the particularity here that the later model is unknown. The proposed solution consists in using an optimal estimation of the image model coming from observed fake codes in order to perform the likelihood test. Using a second order expansion, we derive a linear relation between the quadratic error of the estimated parameters and the probability of type II error. We are then able to formulate analytically and practically the error spread region of the Receiver Operating Characteristic (ROC) curves, and to compute the average authentication performance when the receiver has to estimate the opponent print and scan channel.

Index Terms-Authentication, Statistical estimation, Hypothesis testing, Printed documents

STUDIED PROBLEM

In the world of global exchanges people can use graphical code as a way to perform authentication of physical products such as documents, goods, drugs,.... Our authentication system is based on printed graphical 2D codes using very high resolution printers (2400dpi). Each printed and scanned set of dots (a dot being a binary element) suffers from a stochastic non-invertible noise which makes the reproduction of the original graphical code impossible [START_REF] Picard | Towards fraud-proof id documents using multiple data hiding technologies and biometrics[END_REF][START_REF] Villan | Multilevel 2 d bar codes: toward high-capacity storage modules for multimedia security and management[END_REF][START_REF] Nguyen | Printer technology authentication from micrometric scan of a single printed dot[END_REF][START_REF] Nguyen | Printer identification from micro-metric scale printing[END_REF] (see in Fig. 1). The opponent's goal is then to reproduce a printed and scanned code similar to the original printed one, using a printer that will also generate a non-invertible noise.

The goal of the receiver is to reject copied codes. In previous works [START_REF] Phan Ho | Document authentication using graphical codes: impacts of the channel model[END_REF][START_REF] Phan Ho | Document authentication using graphical codes: Reliable performance analysis and channel optimization[END_REF] we modeled the mentioned authentication system as a hypothesis testing problem and we derived tight bounds on its performance. The authentication system works as follow: a binary authentication image is constructed from a Fig. 1. Left: an original printed and scanned graphical code. Right: a re-printed and scanned (forged) graphical code. randomly chosen binary sequence x N in {0, 1} N by the legitimate source. It is shared secretly with the legitimate receiver and published as a noisy version y N modeling the original printed and scanned code (see Fig. 1 on the left). The opponent observes y N and tries to estimate the original image obtaining x N (see [START_REF] Baras | Towards a realistic channel model for security analysis of authentication using graphical codes[END_REF] in Fig. 2). He then prints it to create a forged observable noisy image z N hopping that it will be accepted by the receiver as coming from the legitimate source (see Fig. 1 on the right). The observed images y N and z N are 8 bits grey level images. In practice, this attack will be used to create false documents or fake packages that could be considered as authentic.

The whole physical process, precisely printing and scanning devices used by the legitimate parts (see [START_REF] Andersen | Asymptotic properties of conditional maximum-likelihood estimators[END_REF] in Fig. 2) and by the counterfeiter (see [START_REF] Bhattacharya | Asymptotic statistics[END_REF] in Fig. 2), are respectively modeled by probability distributions conditioned to the original data P Y /X,Θ and P Z/X, Θ, Θ and Θ are sets of parameters, taken in Ξ, specifying the devices in each case. As pointed in [START_REF] Phan Ho | Document authentication using graphical codes: impacts of the channel model[END_REF], the receiver observes one of the two possible images y N or z N and have to decide whether it comes from the legitimate source or not (see [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] in Fig. 2), supposed that the models P Y /X,Θ and P Z/X,Θ are known. The print and scan process in this particular setup has been modeled by an AWGN channel or an additive i.i.d. lognormal noise in [START_REF] Baras | Towards a realistic channel model for security analysis of authentication using graphical codes[END_REF].

Authentication here is based on classical Neyman-Pearson test (NP-test) (see ( [START_REF] Lehmann | The Fisher, Neyman-Pearson theories of testing hypotheses: One theory or two[END_REF][START_REF] Christensen | Testing Fisher, Neyman Pearson, and Bayes[END_REF])) in which the receiver considers two hypothesis H 0 and H 1 . The former hypothesis attests authenticity, i.e. that the received sequence is generated by P Y /X,Θ and the latter one unveils a fake code, i.e. that the observed sequence is driven from P Z/X, Θ. Performances are evaluated computing the probability of type I error (rejecting hypothesis H 0 while actually the observed sequence comes from the legitimate source) and the probability of type II error (accepting hypothesis H 0 while it is actually a fake).

In this paper we extend our analysis to the case where the receiver doesn't know the true parameter Θ related to the opponent print and scan process, but establishes a test statistic using estimated ones obtained by a maximum likelihood based algorithm. The estimated parameters are computed from several codes identified previously as fake codes which represent a set of printed and scanned dots driven from P Z/X, Θ (see [START_REF] Christensen | Testing Fisher, Neyman Pearson, and Bayes[END_REF] in Fig. 2).

We derive a linear approximation relating the probability of type II error to the quadratic error on the estimated parameters. This approximation helps us to express analytically the error spread of the authentication performance, and hence to evaluate its average when the receiver doesn't know perfectly the opponent print and scan process.

MODELS AND ESTIMATION

Printing and scanning processes

In our global authentication model, P Y /X,Θ represents a grey level distribution of the authentic image conditioned to the knowledge of both the authentication dots and the parameters governing the legitimate print and scan process.

When performing a detection to obtain an estimated sequence XN , the opponent undergoes errors coming from the facts that he is not able to infer the original code. It is important to note that the opponent will have to print a binary version of its observation because a printing device at this very high resolution can only print binary images.

These errors include probabilities P e,W when confusing an original white dot with a black and P e,B when confusing an original black dot (generated by bit 0) with a white dot (generated by bit 1).

Given P Z/X, Θ a grey level distribution of the forged image conditioned to the knowledge of both the authentication dots and the parameters governing opponent print and scan process T Z/X, Θ, we have (equivalently for X = 0 or 1):

P Z/X, Θ(Z = v/X = 0(1), Θ) = (1 -P e,B(W ) )T Z/ X, Θ(v/ X = 0(1), Θ) + P e,B(W ) T Z/ X, Θ(v/ X = 1(0), Θ) (1)

Error spread region for the estimated parameters Θ

Although our analysis in this paper does not depend on the estimation method, we consider maximum likelihood estimation (MLE) in order to achieve optimal estimation. In MLE, we have E[ Θ] = Θ and the estimated set of parameters are jointly normally distributed N ( Θ; V Θ) where V Θ is their covariance matrix which can be computed either practically or analytically using the Fisher information ( [START_REF] Andersen | Asymptotic properties of conditional maximum-likelihood estimators[END_REF][START_REF] Bhattacharya | Asymptotic statistics[END_REF][START_REF] Van Der Vaart | Asymptotic statistics[END_REF]). Under these assumptions the covariance matrix will help us to provide a measure of how the estimated parameters spread w.r.t the true value. The quadratic form of the error (the variation of the estimation) is chi-squared distributed and:

ρ( Θ) = ( Θ -Θ) T V -1 Θ ( Θ -Θ) ∼ χ 2 κ . (2)
where χ 2 κ is the chi-squared distribution with κ degree of freedom ( [START_REF] Bhattacharya | Asymptotic statistics[END_REF][START_REF] Van Der Vaart | Asymptotic statistics[END_REF][START_REF] Wonnacott | Introductory statistics[END_REF]). Herein, κ is the number of parameters that govern the print and scan model. One may observe that ρ( Θ) = cte is an ellipsoid in the κ-dimensional space. Using the property of ρ( Θ) in (2), we can compute the probability that the κ-dimensional error vector lies between two ellipsoids. The error spread region for the estimated parameters is given by: R = Θ :

χ 2 κ,γ1 ≤ ρ( Θ) ≤ χ 2 κ,γ2 , (3)
where χ 2 κ,γ1 and χ 2 κ,γ2 are critical levels w.r.t γ 1 and γ 2 , i.e, P r[ρ( Θ) ≤ χ 2 κ,γ1 ] = γ 1 and P r[ρ( Θ) ≤ χ 2 κ,γ2 ] = γ 2 .

AUTHENTICATION PERFORMANCE

In this section we express reliably the probabilities of type I error α and probability of type II error β used to assess the performance of the authentication system in the NP-test setup.

In [START_REF] Phan Ho | Document authentication using graphical codes: impacts of the channel model[END_REF][START_REF] Phan Ho | Document authentication using graphical codes: Reliable performance analysis and channel optimization[END_REF], we indicate several methods to approximate β while keeping α fixed and show that an Asymptotic Expression (AE) provides a reliable approximation for β. We present a brief summary for the AE method.

Given a real number s the Chernoff bounds on α and β may be expressed as:

α = Pr(L ≥ λ/H 0 ) ≤ e -sλ g L/H 0 (s) for any s > 0, (4)
β = Pr(L ≤ λ/H 1 ) ≤ e -sλ g L/H 1 (s) for any s < 0,

where L is the log-likelihood ratio of P Z N /X N over P Y N /X N and λ is a threshold. The function g L/Hj (s) (j = 0, 1) is the moment generating function of L given the hypothesis H j :

g L/H j (s) = E L/H j e sL . (6)
Without loss of generality, we can use the original code with the same number of bits 0 (N b ) and bits 1 (N w ), i.e., N b = N w = N 2 with N is the length of the codeword. For large N , ( 4) and ( 5) become very tight and specifically with an additional correcting factor we obtain for i.i.d. samples:

α → N →∞ 1 |s 0 | N πµ ′′ ℓ (s 0 ) exp N 2 [µ ℓ (s 0 ) -s 0 µ ′ ℓ (s 0 )] . (7)
β → N →∞ 1 |s 1 | N πµ ′′ ℓ (s 1 ) exp N 2 [µ ℓ (s 1 ) -s 1 µ ′ ℓ (s 1 )] . (8)
for s 0 > 0 and s 1 < 0. Here, given hypothesis H j , µ ℓ (s j ) = log g ℓ/H j (s j ) is the cumulant generating function of P Z/X P Y /X at the point s j which is the solution of N 2 µ ′ ℓ (s) = λ, and note that one can easily prove that s 1 = s 0 -1.

INTERPLAY BETWEEN CHANNEL ESTIMATION AND AUTHENTICATION PERFORMANCE

In this section, we analyze how the set of estimated parameters impacts the performance of the probability of type II error β(α, Θ) for a fixed value of α. Precisely, we relate the error spread region of log β(α, Θ) to the error spread region of Θ defined in [START_REF] Bhattacharya | Asymptotic statistics[END_REF]. For large enough N , the changes of correcting factors

1 |sj | √ N πµ ′′ ℓ (sj )
in ( 7) and ( 8) are imperceptible and we drop their analysis. Due to the limitation of the paper, we will present our analysis using only one estimated parameter. The extension of this analysis for vectors of estimated parameters will be addressed in future works.

Analytical Analysis

We want to show the linear tendency of the scatter of log β(α, Θ) w.r.t ρ( Θ) by using Taylor expansion. For the case of one parameter, let θ be the estimated version of the true parameter θ t ∈ Θ, and θ m ∈ Θ be the parameter of the legitimate model, and define:

β * (θ) = 2 N log β(α, θ), (9)
A Taylor expansion gives:

β * (θ) ∼ = β * (θ t ) + △θ ∂β * (θ) ∂θ θt + (△θ) 2 2 ∂ 2 β * (θ) ∂θ 2 θt + ... (10
) where △θ = (θθ t ).

We aim now to compute the first and second derivatives of β * (θ) expressed for θ = θ t . Because α is fixed and we choose the same threshold for detection, we have:

E P0 l(θ)e l(θ)s0(θ) E P0 e l(θ)s0(θ) = E P1 l(θ)e l(θ)s1(θ) E P1 e l(θ)s1(θ) , (11)
with P 0 ≡ P Y /X,θm , P 1 ≡ P Z/X,θt , l ≡ l(θ) = log

P Z/X,θ
P Y /X,θm . This leads to:

∂β * (θ) ∂θ = s 1 (θ) E P1 ∂l ∂θ e l(θ)s1(θ) E P1 e l(θ)s1(θ) - E P0 ∂l ∂θ e l(θ)s0(θ) E P0 e l(θ)s0(θ) (12)
Notably, when θ = θ t , we have s 1 (θ t ) = s 0 (θ t ) -1 and so it can be proved that for every function f (θ),

E P0 f (θ t )e l(θt)s0(θt) = E P1 f (θ t )e l(θt)s1(θt) (13)
Therefore:

∂β * (θ) ∂θ θt = 0 (14)
The equality ( 14) is not surprising since the NP-test is known to reach the optimum when applied on the true parameter. Now if we denote:

E 1 = E P1 e l(θt)s1(θt) E 2 = E P1 l(θ t )e l(θt)s1(θt) E 3 = E P1 l 2 (θ t )e l(θt)s1(θt) E 4 = E P1 ∂l(θ)
∂θ θt e l(θt)s1(θt)

E 5 = E P1 ∂l(θ) ∂θ θt 2 e l(θt)s1(θt) E 6 = E P1 l(θ t ) ∂l(θ) ∂θ θt e l(θt)s1(θt) (15)
and given R i1 = E i /E 1 for all i = 2, 3, ..., 6 we obtain the second derivative of β * at θ = θ t as

∂ 2 β * (θ) ∂θ 2 θt = [R 61 -R 41 R 21 ] 2 R 31 -R 2 21 -1 +R 2 41 -R 51 , (16)
based on the fact that

2 N ∂ 2 log β(α, θ) ∂θ 2 θt = ∂ 2 β * (θ) ∂θ 2 θt . (17)
If now we call γ(α, θ t ) = ∂ 2 log β(α,θ) ∂θ 2 θt × Var(θ) 2 the slope of analytical linear expression, then from ( 10), ( 14) and ( 17), it yields

log β(α, θ) ∼ = log β(α, θ t ) + γ(α, θ t )ρ(θ) (18)
where Var(θ) is the variance of the estimated parameter θ and ρ(θ) = (△θ) 2 Var(θ) is the variation of the estimation. Moreover, from the property of NP-test, β(α, θ) ≥ β(α, θ t ) for all θ so γ(α, θ t ) is always nonnegative. Using ( 2) and (18), we show that log β(α, θ) follows a shifted and scaled χ 2 distribution and we are now able to derive a spread error region (see 4.2).

Numerical results

In order to perform our analysis, we have to construct a MLE scheme for parameter estimation. It is known that the Expectation Maximization (EM) algorithm is an iterative method for finding maximum likelihood. Without loss of generality we assume that T Z/ X=0, Θ and T Z/ X=1, Θ are modeled by truncated discrete normal distributions with Θ = (μ b , σ2 b , μw , σ2 w ) such that P Z/X, Θ is a mixture of two truncated Gaussians (see ( 1)). We then develop an EM algorithm ( [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF]) for this particular mixture to estimate the set of unknown parameters. In Fig. 3, we suppose that only μw is unknown and we run EM algorithm N iter times using each time N obs observations and obtain a set of μw . The scatter plot of Fig. 3 represents the computed values of log β coming from AE method. It is compared with the analytical expression (18) and the statistical linear regression.

In Fig. 4, we analyze the impact of the estimation error on the ROC curves. We select a 95% confidence error region for μw , i.e., ρ( μw ) is bounded by two critical levels χ 2 1,0.025 and χ 2 1,0.975 such that P r[ρ( μw ) ≤ χ 2 1,0.025 ] = 0.025 and P r[ρ( μw ) ≤ χ 2 1,0.975 ] = 0.975, and we thus obtain a corresponding 95% confidence error region for log β(α, μw ). We then derive two critical ROC curves C 0.025 min and C 0.975 max computed analytically from χ 2 1,0.025 and χ 2 1,0.975 and we choose the mean value for ρ( μw ) to find the mean ROC curve C mean . We then compare C 0.025 min , C 0.975 max and C mean with the three ones computed from the dataset of ρ( μw ) and we observe that our approximation is accurate.

CONCLUSIONS

This paper analyzes the impact of parameters estimation on the performance of document authentication.

-We show experimentally and theoretically an analytical linear relation between the variation of the estimated parameters and the logarithm of the corresponding probabilities of type II error.

-We are able to predict the average authentication loss when performing the NP-test of the estimated distribution for the opponent's channel.

-The proposed analysis is not impacted by the nature of the noise, and can be applied for different memoryless channels that are more realistic to model the printing process, such as the lognormal distribution [START_REF] Baras | Towards a realistic channel model for security analysis of authentication using graphical codes[END_REF].

-Although we particularly focus on document authentication, our analysis in this paper generally can be applied for a larger class of forensics problems.

Fig. 2 .

 2 Fig. 2. Principle of authentication using graphical codes.

Fig. 3 .

 3 Fig. 3. Comparison between log β(ρ) for the true parameter (horizontal line) and the one (dots) from estimated opponent's parameters. Linear regression (stars); analytical expression (straight line). Critical value χ 2 1,0.025 and χ 2 1,0.975 are represented by vertical lines. N obs = 3.10 3 , α = 10 -16 , N iter = 5.10 3 , N = 2.10 3 .

Fig. 4 .

 4 Fig. 4. Comparison between three analytical ROC curves C 0.025 min (squares), C 0.975 max (stars) and C mean (circles) with min ROC curve (dash line), max ROC curve (straight line) and mean ROC curve (straight line with triangles) computed from N iter = 5000 data of μw .