N Malod-Dognin

R Andonov

N Yanev

Maximum Cliques in Protein Structure Comparison

Keywords: protein structure comparison, maximum clique problem, kpartite graphs, combinatorial optimization Compraraison de structures protéiques, problème de clique maximum, graphes k-partis, optimisation combinatoire

Computing the similarity between two protein structures is a crucial task in molecular biology, and has been extensively investigated. Many protein structure comparison methods can be modeled as maximum clique problems in specific k-partite graphs, referred here as alignment graphs.

In this paper, we propose a new protein structure comparison method based on internal distances (DAST) which is posed as a maximum clique problem in an alignment graph. We also design an algorithm (ACF) for solving such maximum clique problems. ACF is first applied in the context of VAST, a software largely used in the National Center for Biotechnology Information, and then in the context of DAST. The obtained results on real protein alignment instances show that our algorithm is more than 37000 times faster than the original VAST clique solver which is based on Bron & Kerbosch algorithm. We furthermore compare ACF with one of the fastest clique finder, recently conceived by Österg ȧrd. On a popular benchmark (the Skolnick set) we observe that ACF is about 20 times faster in average than the Österg ȧrd's algorithm.

Introduction

A fruitful assumption in molecular biology is that proteins of similar threedimensional (3D) structures are likely to share a common function and in most cases derive from a same ancestor. Understanding and computing physical similarity of protein structures is one of the keys for developing protein based medical treatments, and thus it has been extensively investigated [START_REF] Godzik | The structural alignment between two proteins: Is there a unique answer?[END_REF][START_REF] Sierk | Déjà vu all over again: Finding and analyzing protein structure similarities[END_REF]. Evaluating the similarity of two protein structures can be done by finding an optimal (according to some criterions) order-preserving matching (also called alignment) between their components. We show that finding such alignments is equivalent to solving maximum clique problems in specific k-partite graphs referred here as alignment graphs. These graphs could be very large (more than 25000 vertices and 3 × 10 7 edges) when comparing real protein structures. We are not aware of any previous specialized algorithm for solving the maximum clique problem in k-partite graphs. Even very recent general clique finders [START_REF] Konc | An efficient branch-and-bound algorithm for finding a maximum clique[END_REF][START_REF] Tomita | An improved branch and bound algorithm for the maximum clique problem[END_REF] are oriented to notably smaller instances and are not able to solve problems of such size (the available code of [START_REF] Tomita | An improved branch and bound algorithm for the maximum clique problem[END_REF] is limited to graphs with up to 1000 vertices).

For solving the maximum clique problem in this context we conceive an algorithm, denoted by ACF (for Alignment Clique Finder), which profits from the particular structure of the alignment graphs. We furthermore compare ACF to an efficient general clique solver [START_REF] Patric | A fast algorithm for the maximum clique problem[END_REF] and the obtained results clearly demonstrate the usefulness of our dedicated algorithm.

The maximum clique problem

We usually denote an undirected graph by G = (V, E), where V is the set of vertices and E is the set of edges. Two vertices i and j are said to be adjacent if they are connected by an edge of E. A clique of a graph is a subset of its vertex set, such that any two vertices in it are adjacent.

Definition 1

The maximum clique problem (also called maximum cardinality clique problem) is to find a largest, in terms of vertices, clique of an arbitrary undirected graph G, which will be denoted by M CC(G).

The maximum clique problem is one of the first problem shown to be NP-Complete [START_REF] Karp | Reducibility among combinatorial problems[END_REF] and it has been studied extensively in literature. Interested readers can refer to [START_REF] Bomze | The maximum clique problem[END_REF] for a detailed state of the art about the maximum clique problem.

Alignment graphs

In this paper, we focus on grid alike graphs, which we define as follows.

Definition 2 A m × n alignment graph G = (V, E) is a graph in which the vertex set V is depicted by a (m-rows) × (n-columns) array T , where each cell T [i][k] contains at most one vertex i.k from V (note that for both arrays and vertices, the first index stands for the row number, and the second for the column number). Two vertices i.k and j.l can be connected by an edge (i.k, j.l) ∈ E only if i < j and k < l. An example of such alignment graph is given in Fig 2a .   It is easily seen that the m rows form a m-partition of the alignment graph G, and that the n columns also form a n-partition. In the rest of this paper we will use the following notations. A successor of a vertex i.k ∈ V is an element of the set Γ + (i.k) = {j.l ∈ V s.t. (i.k, j.l) ∈ E, i < j and k < l}. V i.k is the subset of V restricted to vertices in rows j, i ≤ j ≤ m, and in columns l, k ≤ l ≤ n. Note that Γ + (i.k) ⊂ V i+1.k+1 . G i.k is the subgraph of G induced by the vertices in V i.k . The cardinality of a vertex set U is |U |.

Relations with protein structure similarity

From a general point of view, two proteins P 1 and P 2 can be represented by their ordered set of components N 1 and N 2 , and estimating their similarity can be done by finding the longest alignment between the elements of N 1 and N 2 . In our approach, such matchings are represented in a

|N 1 | × |N 2 | alignment graph G = (V, E)
, where each row corresponds to an element of N 1 and each column corresponds to an element of N 2 . A vertex i.k is in V (i.e. matching i ↔ k is possible), only if element i ∈ N 1 and k ∈ N 2 are compatible. An edge (i.k, j.l) is in E if and only if (i) i < j and k < l, for order preserving, and (ii) matching i ↔ k is compatible with matching j ↔ l. A feasible matching of P 1 and P 2 is then a clique in G, and the longest alignment corresponds to a maximum clique in G. There is a multitude of alignment methods and they differ mainly by the nature of the elements of N 1 and N 2 and by the compatibility definitions between elements and between pairs of matched elements. At least two protein structure similarity related problems from the literature can be converted into clique problems in alignment graphs : the secondary structure alignment in VAST [START_REF] Gibrat | Surprising similarities in structure comparison[END_REF], and the Contact Map Overlap Maximization problem (CMO) [START_REF] Godzik | Flexible algorithm for direct multiple alignment of protein structures and sequences[END_REF].

VAST, or Vector Alignment Search Tool, is a software for aligning protein 3D structures largely used in the National Center for Biotechnology Information1 . In VAST, N 1 and N 2 contain 3D vectors representing the secondary structure elements (SSE) of P 1 and P 2 . Matching i ↔ k is possible if vectors i and k have similar norms and correspond either both to α-helices or both to β-strands. Finally, matching i ↔ k is compatible with matching j ↔ l only if the couple of vectors (i, j) from P 1 can be well superimposed in 3D-space with the couple of vectors (k, l) from P 2 .

CMO is one of the most reliable and robust measures of protein structure similarity. Comparisons are done by aligning the residues (amino-acids) of two proteins in a way that maximizes the number of common contacts (when two residues that are close in 3D space are matched with two residues that are also close in 3D space). We have already dealt with CMO in [START_REF] Andonov | An efficient lagrangian relaxation for the contact map overlap problem[END_REF], but not by using cliques. Note that a maximum clique formulation in alignment graphs was proposed by Strickland et al. in [START_REF] Strickland | Optimal protein structure alignment using maximum cliques[END_REF], but this formulation differs from ours.

DAST: an improvement of CMO based on internal distances

One of the main drawback of CMO is that in order to maximize the number of common contacts, it also introduces some "errors" like aligning two residues that are close in 3D space with two residues that are remote, as illustrated in Fig 1 . These errors could potentially yield alignments with big root mean square deviations (RMSD) which is not desirable for structures comparison.

To avoid such problems we propose DAST (Distance-based Alignment Search Figure 1: An optimal CMO matching.

Two proteins ( P 1 and P 2 ) are represented by their contact map graphs where the vertices corresponds to the residues and where edges connect residues in contacts (i.e. close). The matching "1 ↔ 1 ′ , 2 ↔ 3 ′ , 4 ↔ 4 ′ ", represented by the arrows, yields two common contacts which is the maximum for the considered case. However, it also matches residues 1 and 4 from P 1 which are in contacts with residues 1 ′ and 4 ′ in P 2 which are remote.

Tool), an alignment method based on internal distances which is modeled in an alignment graph. In DAST, the two proteins P 1 and P 2 are represented by their ordered sets of residues N 1 and N 2 . Two residues i ∈ N 1 and k ∈ N 2 are compatible if they come from the same kind of secondary structure elements (i.e. i and k both come from an α-helix, or from a β-strand) or if both come from a loop. Let us denote by d ij (resp. d k.l ) the euclidean distance between the α-carbons of residues i and j (resp. k and l). Matching i ↔ k is compatible with matching j ↔ l only if |d ijd kl | ≤ τ , where τ is a distance threshold. The longest alignment in terms of residues, in which each couple of residues from P 1 is matched with a couple of residues from P 2 having similar distance relations, corresponds to a maximum clique in G.

Since RM SD = 1 Nm × (|d ij -d kl | 2 )
, where N m is the number of matching pairs "i ↔ k, j ↔ l", the alignments given by DAST have a RMSD of internal distances ≤ τ .

Branch and Bound approach

We have been inspired by [START_REF] Patric | A fast algorithm for the maximum clique problem[END_REF] to propose our own algorithm which is more suitable for solving the maximum clique problem in the previously defined m × n alignment graph G = (V, E). Let Best be the biggest clique found so far (first it is set to ∅), and |M CC(G)| be an over-estimation of |M CC(G)|. By definition, V i+1.k+1 ⊂ V i.k+1 ⊂ V i.k , and similarly

V i+1.k+1 ⊂ V i+1.k ⊂ V i.k .
From these inclusions and from definition2, it is easily seen that for any

G i.k , M CC(G i.k ) is the biggest clique among M CC(G i+1.k ), M CC(G i.k+1 ) and M CC(G i+1.k+1 )
{i.k}, but for the latter only if vertex i.k is adjacent to all vertices in M CC(G i+1.k+1 ). Let C be a (m + 1) × (n + 1) array where

C[i][k] = |M CC(G i.k )| (
values in row m + 1 or column n + 1 are equal to 0). For reasoning purpose, let assume that the upper-bounds in C are exact. If a vertex i.k is adjacent to all vertices in M CC(G i+1.k+1 ), then

C[i][k] = 1 + C[i + 1][k + 1], else C[i][k] = max(C[i][k + 1], C[i + 1][k]). We can deduce that a vertex i.k cannot be in a clique in G i.k which is bigger than Best if C[i + 1][k + 1] < |Best|,
(i.k) ⊂ V i+1.k+1 . Even if C[i + 1][k + 1] ≥ |Best|, if |M CC(Γ + (i.k))| < |Best| then i.k cannot be in a clique in G i.k bigger than Best.
[i][k], C[i][k] is set to max(C[i][k + 1], C[i + 1][k]) if i.k / ∈ V , or to 1 + C[i + 1][k + 1]) if i.k ∈ V .
Note that the order used for visiting the cells in T guaranties that

INRIA when computing the value of C[i][k], the values of C[i + 1][k], C[i][k + 1] and C[i + 1][k + 1] are already computed.
Array C can also be used in function Extend clique() to fasten the maximum clique search. This function is a branch a bound (B&B) search using the following branching rules. Each node of the B&B tree is characterized by a couple (Cli, Cand) where Cli is the clique under construction and Cand is the set of candidate vertices to be added to Cli. Each call to Extend clique({i.k}, Γ + (i.k)) create a new B&B tree which root node is ({i.k}, Γ + (i.k)). The successors of a B&B node (Cli, Cand) are the nodes (Cli {i ′ .k ′ }, Cand Γ + (i ′ .k ′ )), for all vertices i ′ .k ′ ∈ Cand. Branching follows lexicographic increasing order (row first). According to the branching rules, for any given B&B node (Cli, Cand) the following cutting rules holds :

Maximum clique cardinality estimators

Even if the described functions depend on array C, they also use another upperestimator of the cardinality of a maximum clique in an alignment graph. By using the properties of alignment graphs, we developed the following estimators.

Minimum number of rows and columns

Definition 2 implies that there is no edge between vertices from the same row or the same column. This means that in a m × n alignment graph, |M CC(G)| ≤ min(m, n). If the numbers of rows and columns are not computed at the creation of the alignment graph, they can be computed in O(|V |).

Longest increasing subset of vertices

Definition 3 An increasing subset of vertices in an alignment graph

G = {V, E} is an ordered subset {i 1 .k 1 , i 2 .k 2 , . . ., i t .k t } of V , such that ∀j ∈ [1, t-1], i j < i j+1 , k j < k j+1 .

LIS(G) is the longest, in terms of vertices, increasing subset of vertices of G.

Since any two vertices in a clique are adjacent, definition 2 implies that a clique in G is an increasing subset of vertices. However, an increasing subset of vertices is not necessarily a clique (since vertices are not necessarily adjacent), and thus

|M CC(G)| ≤ |LIS(G)|. In a m × n alignment graph G = (V, E), LIS (

Definition 4

The longest increasing subsequence of an arbitrary finite sequence of integers S = "i i , i 2 , . . . , i n " is the longest subsequence S ′ = "i ′ i , i ′ 2 , . . . , i ′ t " of S respecting the original order of S, and such that for all j ∈ [1, t], i ′ j < i ′ j+1 . By example, the longest increasing subsequence of "1,5,2,3" is "1,2,3".

For any given alignment graph G = {V, E}, we can easily reorder the vertex set V , first by increasing order of columns, and second by decreasing order of rows. Let's denote by V ′ this reordered vertex set. Then we can create an integer sequence S corresponding to the row indexes of vertices in V ′ . For example, by using the alignment graph presented in Fig2a, the reordered vertex set V ′ is {4.1, 2.1, 1.1, 3.2, 4.3, 3.3, 2.3, 1.3, 4.4, 3.4, 1.4}, and the corresponding sequence of row indexes S is "4, 2, 1, 3, 4, 3, 2, 1, 4, 3, 1". An increasing subsequence of S will pick at most one number from a column, and thus an increasing subsequence is longest if and only if it covers a maximal number of increasing rows. This proves that solving the longest increasing subsequence in S is equivalent to solving the longest increasing subset of vertices in G. Note that the longest increasing subsequence problem is solvable in time O(l × ln(l)) [START_REF] Fredman | On computing the length of longest increasing subsequences[END_REF], where l denotes the length of the input sequence. In our case, this corresponds to O(|V | × ln(|V |)).

INRIA

Longest increasing path

Definition 5 An increasing path in an alignment G = {V, E} is an in- creasing subset of vertex {i 1 .k 1 , i 2 .k 2 , . . ., i t .k t } such that ∀j ∈ [1, t -1], (i j .k j , i j+1 .k j+1 ) ∈ E. The longest increasing path in G is denoted by LIP (G)
As the increasing path take into account edges between consecutive vertices,

|LIP (G)|, should better estimate M CC(G)|. |LIP (G)| can be computed in O(|V | 2 ) by the following recurrence. Let DP [i][k] be the length of the longest in- creasing path in G i.k containing vertex i.k. DP [i][k] = 1+ max i ′ .k ′ ∈Γ + i.k (DP [i ′ ][k ′ ]).
The sum over all Γ + (i.k)) is done in O(|E|) time complexity, and finding the maximum over all

DP [i][k] is done in O(|V |). This results in a O(|V | + |E|) time complexity for computing |LIP (G)|.
Amongst all of the previously defined estimators, the longest increasing subset of vertices (solved using the longest increasing subsequence) exhibits the best performances and is the one we used for obtaining the results presented in the next section.

Results

All results presented in this section come from real protein structure comparison instances. Our algorithm, denoted by ACF (for Alignment Clique Finder), has been implemented in C and was tested in two different contexts: secondary structure alignments in VAST and residue alignments in DAST. ACF will be compared to Österg ȧrd's algorithm [START_REF] Patric | A fast algorithm for the maximum clique problem[END_REF] (denoted by Östergȧrd ) and to the original VAST clique solver which is based on Bron and Kerbosch's algorithm [START_REF] Bron | Algorithm 457: finding all cliques of an undirected graph[END_REF] (denoted by BK ). Note that BK is not a maximum clique finder but returns all maximal cliques in a graph.

Secondary structures alignments

This section illustrates the behavior of ACF in the context of secondary structure element (SSE) alignments. For this purpose we integrated ACF and Östergȧrd (which code is freely available) in VAST. We afterwards compared them with BK by selecting few large protein chains having between 80 to 90 SSE's (for smaller protein chains the running times of both Östergȧrd and ACF are less than 0.01 sec.). Computations were done on a AMD at 2.4 GHz computer, and the corresponding running times are presented in table 1. We observe that Östergȧrd is 4053 times faster than BK, and that ACF is about 9.3 times faster than Östergȧrd. Although we have chosen large protein chains, the SSE alignment graphs are relatively small (up to 5423 vertices and 551792 edges ). On such graphs the difference between Östergȧrd and ACF performance is not very visible-it will be better illustrated on larger alignment graphs in the next section.

Residues alignment

In this section we compare ACF to Östergȧrd in the context of residue alignments in DAST. Computations were done on a PC with an Intel Core2 processor at 3Ghz, and for both algorithms the computation time was bounded to 5 hours per instance. Secondary structures assignments were done by KAKSI [START_REF] Martin | Protein secondary structure assignment revisited: a detailed analysis of different assignment methods[END_REF], and the threshold distance τ was set to 3 Å. The protein structures come from the well known Skolnick set, described in [START_REF] Lancia | 101 optimal pdb structure alignments: a branch-and-cut algorithm for the maximum contact map overlap problem[END_REF]. It contains 40 protein chains having from 90 to 256 residues, classified in SCOP [START_REF] Andreeva | Data growth and its impact on the scop database: new developments[END_REF] (v1.73) into five families.

INRIA

Amongst the 780 corresponding alignment instances, 164 align protein chains from the same family and will be called "similar". The 616 other instances align protein chains from different families and thus will be called "dissimilar". Characteristics of the corresponding alignment graphs are presented in table 2. Figure 3 compares the running time of ACF to the one of Östergȧrd on the set of 673 instances solved by both algorithms (all instances solved by Östergȧrd were also solved by ACF ). For all instances except one, ACF is significantly faster than Östergȧrd. More precisely, ACF needed 12 hs. 29 min. 56 sec. to solve all these 673 instances, while Östergȧrd needed 260 hs. 10 min. 10 sec. Thus, on the Skolnick set, ACF is about 20 times faster in average than Östergȧrd, (up to 4029 times for some intstances). The ACF time is presented on the x-axis, while the one of Österg ȧrd is on the y-axis. For all instances except one, ACF is faster than Österg ȧrd.

Conclusion and future work

In this paper we introduce a novel protein structure comparison approach DAST, for Distance-based Alignment Search Tool. For any fixed threshold τ , it finds the longest alignment in which each couple of pairs of matched residues shares the same distance relation (+/τ ), and thus the RMSD of the alignment is ≤ τ . This property is not guaranteed by the CMO approach, which inspired initially DAST. From computation standpoint, DAST requires solving the maximum clique problem in a specific k-partite graph. By exploiting the peculiar structure of this graph, we design a new maximum clique solver which significantly outperforms one of the best general maximum clique solver. Our solver was successfully integrated into two protein structure comparison softwares and will be freely available soon. We are currently studying the quality of DAST alignments from practical viewpoint and compare the obtained results with other structure comparison methods.

INRIA

Figure 2 :

 2 Figure 2: A 4 × 4 alignment graph and the visiting order of its array T

  (i) if |Cli| + |Cand| ≤ |Best| then the current branch cannot lead to a clique bigger than |Best| and can be fathomed, (ii) if |M CC(Cand)| ≤ |Best| -|Cli|, then the current branch cannot lead to a clique bigger than |Best|, and (iii) if |M CC(Cand Γ + (i.k))| ≤ |Best| -|Cli| -1, then branching on i.k cannot lead to a clique bigger than |Best|. For any set Cand and any vertex i.k, Cand Γ + (i.k) ⊂ Γ + (i.k) , and Γ + (i.k) ⊂ G i+1.k+1 From these inclusions we can deduce two way of over-estimating |M CC(Cand Γ + (i.k))|. First, by using C[i + 1][k + 1] which over-estimate |M CC(G i+1.k+1 )| and second, by overestimating |M CC(Γ + (i.k))|. All values |M CC(Γ + (i.k))| are computed once for all in Find clique() and thus, only |M CC(Cand)| needs to be computed in each B&B node.

  G) can be computed in O(n×m) times by dynamic programming. However, it is possible by using the longest increasing subsequence to solve LIS(G) in O(|V | × ln(|V |)) times which is more suited in the case of sparse graph like in our protein structure comparison experiments.

Figure 3 :

 3 Figure 3: Running time comparison on the Skolnick set

  and this reasoning still holds if values in C are upper estimations. Another important inclusion is Γ +

Table 1 :

 1 Runing time comparison on secondary structure alignment instancesRuning time comparison of BK, Österg ȧrd and ACF on secondary structure alignment instances for long protein chains (containing from 80 to 90 SSE's). BK is notably slower than the Österg ȧrd's algorithm, which is slightly slower than ACF.

	Instances	BK (sec.)	Österg ȧrd (sec.)	ACF (sec.)
	1k32B	1n6eI	1591.89	1.42	0.09
	1k32B	1n6fB	1546.78	0.01	0.01
	1k32B	1n6fF	1584.25	0.14	0.02
	1n6dD	1k32B	1373.35	0.06	0.01
	1n6dD	1n6eI	1390.27	0.11	0.03
	1n6dD	1n6fB	1328.85	0.65	0.06
	1n6dD	1n6fF	1398.41	0.13	0.05

Table 2 :

 2 DAST alignment graphs characteristics

			array size	-V-	-E-	density	-MCC-
	similar	min	97×97	4018	106373	8.32%	45
	instances	max	256×255	25706	31726150	15.44%	233
	dissimilar	min	97×104	1581	77164	5.76%	12
	instances	max	256×191	21244	16839653	14.13%	48
	All alignment graphs from DAST have small edge density (less than 16%). Similar instances
	are characterized by bigger maximum cliques than the dissimilar instances.

Table 3

 3 compares the number of instances solved by each algorithm on Skolnick set. ACF solved 155 from 164 similar instances, while Östergȧrd solved 128 instances. ACF was able to solve all 616 dissimilar instances, while Östergȧrd solved 545 instances only. Thus, on this popular benchmark set, ACF clearly outperformed Östergȧrd in terms of number of solved instances.

Table 3 :

 3 Number of solved instances comparison

		Österg ȧrd	ACF
	Similar instances (164)	128	155
	Dissimilar instances (616)	545	616
	Total (780)	673	771
	Number of solved instances on Skolnick set: ACF solves 21% more similar instances and 13% more dissimilar instances than Österg ȧrd.

RR n 7053

http://www.ncbi.nlm.nih.gov/Structure/VAST/vast.shtml INRIA

Acknowledgements

This work is a part of ANR project PROTEUS "ANR-06-CIS6-008", Noël Malod-Dognin is supported by the Brittany Region, and Nicola Yanev is supported by the bulgarian project DVU/01/197, South-West University, Blagoevgrad. All computations were done on the Ouest-genopole bioinformatics platform (http://genouest.org). We would like to express our gratitude to J-F Gibrat for numerous helpful discussions and for providing us the source code of VAST.


Ivana Roche
email: ivana.roche@inist.fr

Dominique Besagni
email: dominique.besagni@inist.fr

Claire François
email: claire.francois@inist.fr

Marianne Hörlesberger
email: marianne.hoerlesberger@arcs.ac.at

Edgar L Schiebel
email: edgar.schiebel@arcs.ac.at

Identification and Characterisation of Technological Topics in the Field of Molecular Biology

Keywords: Visualisation and Science Mapping: tools, methods and applications (Theme 5) emerging technologies, bibliometric indicators, terminology evolution, diffusion model, diachronic clustering

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

Background

This paper focuses on methodological approaches for characterising the specific topics within a technological field based on scientific literature data. We introduce a diachronic clustering analysis approach and some bibliometric indicators. The results are visualised with the software-tool Stanalyst® [1]. We are applying our methods to the field "Molecular Biology". This field has grown a great deal in the last decade.

Problem / application

How can we identify and characterise important topics in a set of several thousand articles? Which technological aspects can be detected? Which of them are already established and which of them are new? How are the topics linked to each other? We are trying to answer these questions by applying our bibliometric analysis methods to a set of scientific literature data recorded in the bibliographic database PASCAL.

The data for our study were extracted from PASCAL database in function of their classification categories and keywords.

The diffusion model identifies three categories of terms: the established terms, terms unusual in this topic, and cross section terms. We apply the indicator TFIDF (text frequency inverse document frequency) adapted to our research question and the GINI coefficient, a measure of statistical dispersion most prominently used as a measure of inequality of income distribution or inequality of wealth distribution.

The diachronic cluster analysis is realized with the help of a clustering tool applying first the axial K-means method to produce a non-hierarchical clustering algorithm based on the neuronal formalism of Kohonen's self-organizing maps and then a principal component analysis to map the obtained clusters. This tool is implemented in the information analysis platform Stanalyst®. Considering two successive time periods, our diachronic approach enables us to follow the field time evolution by analysing the two obtained cluster sets and maps.

The applied methods are linked together as the following table shows.

Diffusion Model

Unusual Terms

Established Terms Cross Section Terms

Diachronic

Clusters in the First Period

Cluster Analysis

Clusters in the Second Period new terms terms with roots in the first period terms with roots in the first period

Outcome/findings/results

The diachronic cluster analysis is made for two periods and shows which topics of the second period have roots in the first one and which topics are new in the second period.

On the map of both periods, as we can see in the above figure, it is possible to singularize an interesting dichotomous configuration represented by two very strongly connected cluster networks associating about two thirds of the clusters. However, these networks have quite different characteristics. On the one hand, the themes present in the biggest one are related to the modelling and simulation of biological phenomena. In the other hand, the little network is very homogeneous and deals essentially with instrumentation topics. The remaining clusters are scattered in the map with no significant links with the two above described networks.

Concerning the biggest network, located in the left side of the first period map and in the right side of the second period map, we can observe that structurally it remains the same with two sub-networks brought together by the cluster "Computerized simulation". But its content analysis allows us to find some interesting thematic evolutions. In the first period, its content is focused on both the description of physical characteristics of biological structures and the theoretical aspects related to the physiological process modelling using, for example, neural networks. In the second period, we can detect a refocusing of the couple formed by the clusters "Neural networks" and "Stochastic processes" and an impressive densification of the resulting sub-network with the new associated cluster "Brownian motion".

The smallest network, located in the right side of the first period map, deals with instrumentation techniques applied to measure and therapeutic issues. On the second period map, the network is located in the left side and shows a significant stability despite a reorganization of cluster contents.

1st period cluster map

Conclusion

"Molecular Biology" is a broad field. By applying the complementary methods presented here, it can be characterised and described presenting different views of the features of the field. On the other hand, the two methods point out the relationship of the different topics in our investigated field and their evolution.

  . (2003). Information Retrieval, Suchmodelle und Data-Mining-Verfahren für Textsammlungen und das Web, Heidelberg: dpunkt. POLANCO, X.; FRNAÇOIS, C.; ROYAUTE, J.; BESAGNI, D.; ROCHE, I. (2001). Stanalyst®: An Integrated Environment for Clustering and Mapping Analysis on Science and Technology, In: Proceedings of the 8th ISSI, Sydney, July 16-20, 2001.


Keywords: Electricity, Transmission Capacity Allocation, Market Coupling, ATC-based, flow-based

This paper aims at understanding and assessing the main methods proposed by ETSO and EuroPEX [1] to define and allocate interconnection capacities in Europe. We model these methods and evaluate their technical and economic efficiency on a 7-node network. We find first, that, unsurprisingly, the allocation methods are all the more efficient as the Kirchhoff laws are integrated more precisely in the market clearing. Second, and more surprisingly, we find that the zonal vision of the grid considered by all these methods has important consequences on the technical and economic efficiency as it may cause the capacity limits of powerlines to be exceeded.

I. INTRODUCTION

HE interconnections are fundamental for the European power market as they are the easiest way to stimulate competition between incumbents [START_REF]DG Competition report on energy sector inquiry[END_REF]. That is why the TSOs, following the intention of the European Commission [START_REF]COMMISSION DECISION of 9 November 2006 amending the Annex to Regulation (EC) No 1228[END_REF], have proposed the concept of flow-based allocation method that should allow the market participants to use transmission capacities closer to their physical limits. The market coupling in Europe currently relies on the concept of Available Transfer Capacity (ATC) where the physical Kirchhoff laws are partially ignored. This assumption generally prevents from using all physical interconnection capacity and so from benefiting from more competition and efficiency in the power market. Some TSOs in the Central Western Europe plan to implement in some years the flow-based allocation method. This method should overcome the inefficiencies resulting from the concept of ATC as the market for power with the flowbased allocation method is cleared taking into account explicitly the Kirchhoff laws.

The introduction of this new flow-based concept to allocate transmission capacity in Europe nevertheless raises two main questions. First what is precisely the flow-based allocation method? Indeed, currently, no model of either ATC-or flow-Manuscript received March 8, 2009. This work was supported by the French Regulatory Commission (CRE) under the €nergy Project II.

Vincent Rious and Philippe Dessante are with Supelec, Gif-sur-Yvette, 91190, France (phone: 0033 16985 1531, fax: 0033 16985 1539, vincent.rious@supelec.fr).

based allocation method is publicly available, except their literary descriptions by ETSO and EuroPEX [START_REF] Etso | Development and Implementation of a coordinated model for regional and inter-regional congestion management[END_REF]. Second what are the real gains of the flow-based allocation method compared to the current ATC-based method? The Central Western European TSOs and PXs have already compared the ATC-and flow-based allocation methods [START_REF] Apx | Implementation Study. A report for the MoU signatories on the design of the market coupling solution in the Central West European (CWE) region[END_REF]. This preliminary evaluation results in a quite paradoxical conclusion: the flowbased method might be less efficient than the current ATCbased method.

To answer these two questions, we propose to model, evaluate and compare the ATC-and flow-based methods described by ETSO and EuroPEX [START_REF] Etso | Development and Implementation of a coordinated model for regional and inter-regional congestion management[END_REF]. Modelling these methods prompts us to enter in the intimacy of their design. As a consequence, we can better know the steps needed in the different methods to define the interconnection capacity. Moreover, these models allow us to compare the ATC-and flow-based allocation methods, not only measuring social welfare but also checking that these allocation methods respect the physical limits of the power grid.

The aim of this paper is twofold. Firstly this paper proposes models for the ETSO-EuroPEX methods for defining and allocating interconnection capacities in Europe. Secondly this paper assesses the differences of these methods 1° from a conceptual point of view, from the point of view of its 2° technical and 3° economic efficiency.

This paper is organized as follow. In section II we present the assumptions and data we will use to evaluate the ETSO-EuroPEX methods. In section III, we present and analyse the similarities and differences of these methods. In section IV, we model these methods and evaluate their technical and economic efficiency. Section V concludes.

II. ASSUMPTIONS AND DATA

A. Assumptions

In this paper, we use the five following technical assumptions. 1° The DC approximation is used to calculate the power flows. 2° We don't consider any N-1 or N-k criteria. 3° We focus only on day-ahead exchanges for calculating cross-border capacities. 4° We assume the TSOs are perfectly coordinated mimicking a unique TSO for the operation of the whole interconnected network. Therefore we test only the efficiency of each congestion management scheme. 5° We Real gains from flow-based methods for allocating power transmission capacity in Europe

Vincent Rious, Philippe Dessante T work with a perfect knowledge of the available line capacities and of the generation and load schedule established before day-ahead. We make also the two following economic assumptions. 1° The transmission capacities are allocated jointly with the market clearing through an implicit mechanism. 2° Neither TSO nor generator have market power.

B. Network for simulations

To evaluate the ETSO-EuroPEX congestion management schemes, we use the network already described in [START_REF] Bompard | A comparative analysis of congestion management schemes under a unified framework[END_REF] because the authors used this network too in order to test the efficiency of different congestion management schemes. This network with 7 nodes and 10 lines is cut out in 3 zones as illustrated in Fig. 1. The capacity of each powerline is 100 MW.

The generators are connected only to nodes 1, 3, 4, 5 and 7 Their cost functions are the same ones as those used in [START_REF] Bompard | A comparative analysis of congestion management schemes under a unified framework[END_REF]. A load is connected to each node. We assume that the consumers are not elastic to price.

C. Base case

Generation and load scheduled for every node before dayahead and the resulting line flows describe the base case. A physical margin can then be calculated for each line as the difference between the line capacity and the base case flow. We decide here that the base case corresponds to the situation where the market equilibrium is obtained considering each zone as isolated, that is to say that there is no commercial cross-border exchange (see Fig. 1). Then the market prices are respectively in zone 1, 2 and 3, 31.5 €/MWh, 39.8 €/MWh, and 63.5 €/MWh.

D. Benchmark: nodal pricing and optimal network use

To know the optimal use of the network infrastructures, we consider the equilibrium of a nodal market (see Fig. 2).

Line 3-6 is congested. Zone 1 is then exporting 191 MW, zone 3 importing 186 MW and zone 2 is slightly importing 5 MW. The social cost is then 30403 €. In this section, we present the three main ETSO-EuroPEX congestion management schemes (over the five ones proposed). We identify their similarities and their differences.

A. Presentation of the ETSO-EuroPEX methods

The first main method presented by ETSO and EuroPEX is Available Transfer Capacity (ATC). To ensure security, this method applied on a meshed network takes into account the possible impact of unidentified powerflows resulting by the Kirchhoff laws from transactions between other market zones. It forces the TSOs to display small cross-border capacities compared to their optimal values considering the most constraining situation. This method is widely used for instance on the French borders with Belgium, and Germany [8].

The second main method exposed by ETSO and EuroPEX is "Combined ATC". This method consists in determining some ATC for exchanges related to several borders. This definition allows to offer more important cross-border capacities to the market than with the previous method. The Combined ATC method is used for instance in Germany where there is a joint limit on the exchanges with France, Switzerland, and the Netherlands.

The last main method exposed by ETSO and EuroPEX is exclusively flow-based. The Kirchhoff laws are explicitly considered in the market clearing. It is the Critical Branch flow-based method. In this method, the constraints on the network are not aggregated on the border. Some real powerlines and their physical limits are included in the market clearing because they are critical branches for cross-border exchanges. This method for defining the exchange capacities is currently under study in several regions in Europe because it could offer a maximum of cross-border capacities while ensuring a high level of security for the power system.

B. Similarities and differences in the methods

The methods exposed by ETSO and EuroPEX to define cross-border capacity have differences and similarities. Identifying them is needed to design an analysis framework for evaluating these methods in similar conditions while respecting their particularities.

1) A common point for all the methods from ETSO-EuroPEX, the zonal approach

All the ETSO-EuroPEX methods are zonal congestion management schemes. However, it is the nodal day-ahead modification of the generation pattern that creates, through the Kirchoff laws, additional power flows to the base case ones. To calculate the cross-border capacities in the ETSO-EuroPEX methods, it is necessary to set an assumption on how the day-ahead increase (respectively decrease) of the zonal net exports are distributed and created by the individual increase (resp. decrease) of generators in each zone.

The nodal sharing of the modification of net export of each zone requires to define a matrix called Generation Shift Keys (GSK). This GSK matrix has n raws and z columns, where n is the number of nodes on the network and z the number of market zones cutting out the network. The transposed matrix in Equation 1 below gives approximated GSK values from the base case to the optimum both presented in II. D. We will use the values in the rest of this paper.

The meaning of the GSK matrix is the following one. Let consider the first raw of this transposed matrix. When there is a change in the net export of zone 1 (for instance an increase by 100 MW), 45% of this zonal change comes from generator at node 1 and 55% of the zonal change comes from generator located at node 4. This rationale extends to the other zones.

          = ′ % 50 0 % 50 0 0 0 0 0 0 0 0 % 100 0 0 0 0 0 % 55 0 0 % 45 K GS (1)
ETSO and EuroPEX formalise the concept of a GSK matrix for the Critical Branch method only. At the same time, the GSK matrix generalises two approaches already exposed by ETSO in [START_REF] Etso | Definitions of Transfer Capacities in Liberalised Electricity Markets[END_REF] to calculate ATC. In these two approaches, the ATC were calculated assuming that the variation of zonal net export was proportionally distributed on the generators either related to their remaining generation capacities or related to the power they already generate. The zonal approach is common to all the ETSO-EuroPEX methods. Then, the concept of GSK or a similar one is needed for all the methods, even if ETSO and EuroPEX do not explicitly formulate this need. With the GSK matrix, one can then establish (2) between the real PTDF matrix and a zonal PTDF matrix reflecting the definition of market zones.

PTDF zonal = GSK * PTDF nodal (2)
The GSK matrix may vary between two exports situations. For instance, between two very different export situations, one can expect that more generators will participate in the situation where the export value is the highest. More generators will then be considered in the definition of the GSK matrix in this situation.

A precise evaluation of the GSK matrix is then fundamental. Indeed, if the GSK matrix is wrong, the TSOs cannot ensure that observing contractual congestion in the market clearing corresponds to real physical congestion. An error in the GSK matrix can imply an over-or under-use of the network during contractual congestion. The accuracy of the GSK matrix is also critical for the efficiency of the European congestion management schemes. This point is all the more essential that the regulators may face an asymmetry of information to evaluate the accuracy of the GSK matrix as the TSOs have discretion on setting this matrix.

2) Differences between the ETSO-EuroPEX methods

Besides the characteristics that make each method different from each other, one particular characteristic allows us to classify the ETSO-EuroPEX methods in two categories, 1° the ATC-based methods, and 2° the flow-based method. Obviously, for the two methods namely ATC and Combined ATC, the cross-border capacities are defined relying on the concept of ATC. The other method we study in this paper is the Critical Branch flow-based and as so is classified in the flow-based category. This difference is important because it will force us to model differently the cross-border exchanges for these two categories. For the ATC-based methods, we use a classical equation for any transfer of matter (as water or gas) to link the net (nodal or zonal) production with the interconnection flows. For the flow-based method, we resort on a PTDF (Power Transfer Distribution Factors) matrix that links the net (nodal or zonal) production with the flows on the network lines or interconnections.

The other differences among the ETSO-EuroPEX methods make them different from each other. The two ATC-based methods are different because different assumptions are used to calculate the values of ATC. 1° The ATC are calculated assuming that the available physical capacity on each line is equally shared among the exchanges between adjacent zones. 2° The Combined ATC are calculated giving priority to a zone for export. 3° The Critical Branch flow-based method considers the constraints on each individual line.

IV. MODEL OF THE ETSO-EUROPEX METHOD

From the similarities and differences of the ETSO-EuroPEX methods, we can now develop a standardised model while also integrating the particularities of each method. We first model the ATC-based methods and then the flow-based ones. ( )

A. A uniform model for the ATC-based methods

     ≤ ≤ - ′ ∀ = + ∆ ∀ ′ ′ ′ ′ ′ ∈ ∑ ∑ z z z z z z z z z z i i ATC T ATC z z T P z , zones of couple a , 0 , zone (4)
(5)

Equation ( 3) is the objective function minimising social cost resulting from the change in generation pattern after the day-ahead market clearing. Equation (4) illustrates the equality constraints for each zone. The change in the net export situation of a zone z is equal to the sum of cross-border exchanges with the neighbouring zones z'. And Equation ( 5) stands for the ATC limits.

The different ATC-based methods have in common the three following steps of the calculation of ATC. 1° A nodal base case for generation and load allows to calculate the base case power flow F ij ° flowing through the line between nodes i and j. Only a subset of lines that may be congested is monitored. 2° The TSO must estimate the GSK matrix. 3° The coefficient PTDF ij,zz' defines the additional flow on line ij for an exchange between the two zones z and z'. Ref. [START_REF] Rious | Assessing Available Transfer Capacity on a Realistic European Network: Impact of Assumptions on Wind Power Generation[END_REF] shows that the capacity really available F ij,zz' re. av.

on line ij for an exchange between zones z and z' is given by:

F ij,zz' re. av. = (F ij max -sign(PTDF ij,zz' ) F ij °) (6)
The different ATC-based methods are then characterised by the assumption used to share the available capacity of the monitored lines among the different cross-border exchanges.

1) ATC

The first ATC-based method that we study differs from the other ATC-based methods by the following assumption: the available physical capacity of each powerline ij is equally shared among the k different cross-border exchanges that influence the flow on this line [8]. Ref. [START_REF] Rious | Assessing Available Transfer Capacity on a Realistic European Network: Impact of Assumptions on Wind Power Generation[END_REF] presents this method and we recall it with the three following equations.

F ij,zz' available = (F ij max -sign(PTDF ij,zz' ) F ij °)/k (7) Ex zz',ij max = (F ij max -sign(PTDF ij,zz' ) F ij °)/(k. |PTDF ij,zz' |) (8) ATC zz' = min{(F ij max -sign(PTDF ij,zz' ) F ij °)/(k.|PTDF ij,zz' |), ij} (9) Where F ij,zz'
available is the available physical capacity on line ij for exchanges between the zones z and z' in this method. Equation (7) mathematically expressed the above mentioned assumption. Other things equal, (8) gives the maximal exchange Ex zz',ij max between the zones z and z' that can constrain line ij (while assuming the other lines cannot be congested). Eventually, (9) gives the value of ATC zz' which is the minimal value of Ex zz',ij max as soon as a line is constrained. Because of the assumption about sharing the available capacity of line to saturate the interconnections, it is unlikely that a line be really congested when the cross-border exchanges reach the limit set by these ATC.

2) Combined ATC

The Combined ATC stands for a whole family of definition of ATC-based cross-border capacities. We present here an example where the Combined ATC are defined to maximise in priority the exports of a given zone z. For this zone, we then compute the optimisation described by (10) to (12).

ATC zz' min is a minimal value that can possibly be imposed on some borders because it can be politically difficult to present zero as a maximal value for some cross-borders exchanges.

The ATC corresponding to the exports from the zones other than the z zone are calculated considering that the export exchanges from the priority zone z must always be able to flow through the network. The export exchanges from zone z are then included in a base case specially used to calculate the ATC from the other zones. It may then be needed to give each zone a priority. We then realise the same process as the one described here. These iterative choices of priority for export from the different zones are similar to define approximately a zonal merit order. The zone with the highest priority to export is supposed to have the smallest price, etc. If the realised market equilibrium is to far from the estimation of the zonal merit order established while computing the Combined ATC, it is possible to have a contractual congestion while the network is far to be physically congested. To overcome this difficulty, it is possible to calculate every set of Combined ATC anticipating every order of priority for export from the different market zones. Each set of Combined ATC would then be tested while calculating the market equilibrium and the retained set of Combined ATC would then be the one which maximises the social surplus.

A major difficulty still remains in this method. It is a problem that the TSO as a monopoly has to anticipate the market outcome through the priority order of export from zones. Even if the TSO is regulated. Besides, it is not sure that the TSO himself wants to do this work: it is out of his core business and it places him in a delicate situation vis-à-vis the network users.

B. A uniform model for flow-based methods

Equations ( 13) to (15) model an implicit flow-based congestion management scheme.

( )

∑ ∆ bids the i i i P C min (13) s. t. [ ]      ≤ ∆ + ⋅ ≤ - ∀ = ∆ ∑ max min , border or line 0 l l l i i F P P PTDF F l P (14) (15)
Equation ( 13) is the objective function of our problem. Equation ( 14) illustrates the global equality constraint between generation and load. Equation (15) expresses the constraints on flows. The different flow-based methods differ from each other depending on the form of this last equation where flows can be considered aggregated. However, we consider here only the Critical Branch flow-based.

1) Critical Branch flow-based method

Contrary to the other ETSO-EuroPEX methods, the Critical Branch method distinctively considers the constraints on the individual lines. As a result, in this method, there is no need to define aggregated transmission capacity. However, it is assumed that the influence of generation on flows is zonal, like in the other ETSO-EuroPEX methods. The TSOs then compute a GSK matrix. One can then model this method only changing (15) by ( 16) in the uniform framework for the flowbased method that we presented just above.

[ ] [ ] [

] max

,

l F P P GSK PTDF l line base nodal l ≤ ∆ + ⋅ ⋅ ∀ (16)
The constrained lines being explicitly identified gives three main advantages to this method. First, there is no need of assumptions to anticipate the use of infrastructures after the market clearing, neither through a sharing assumption nor through an assumption about the market outcome itself. Second, this method is more transparent than the other ETSO-EuroPEX methods because the constraints are not aggregated. Lastly, as this method is very close to nodal pricing, it should lead to a more efficient use of the network, although the effect of the zonal view must be evaluated on economic efficiency and on network security.

V. EVALUATION AND NUMERICAL COMPARISON

Now we evaluate the ETSO-EuroPEX methods that we modelled above. First, we verify that the flows resulting from the market outcome respect the real constraints on powerlines. Second, we measure economic efficiency through the social cost on the network presented in II. When needed, we use the GSK matrix presented in III.B.1.

This section is organised as follow. Each subsection studies each of the three main ETSO-EuroPEX methods. We calculate a) the exchange capacities, b) the market equilibrium, and c) the real powerflows resulting from the market outcome.

A. ATC

We evaluate in this subsection the ATC method applying the model of IV.A.1) to the base case in II. To calculate the value of ATC on this system (see table II), the physical available capacity on each line must be shared in 3 because the flow on each line is influenced by three exchanges. For instance, the flow on line 6-7 from 7 to 6 is influenced by the cross-border exchanges from 1 to 2, 2 to 3 and 3 to 1. As a consequence, when the maximal values of ATC are reached, the contractual congestion of a line means that only a third of its available capacity is congested. When the above ATC values are transmitted to the market, the exchanges are constrained by the ATC 1-2, 1-3 and 2-3.

The zonal market prices are then 1° 35 €/MWh, 2° 39 €/MWh and 3° 58 €/MWh. These prices are quite different from the optimal prices that are at maximum 52 €/MWh.

The flows resulting from the market outcome don't congest the power lines. This is because of the sharing assumption. The physical congestion of a line can then happen only if the market simultaneously requires that every exchange saturates this line. This situation is very unlikely as the ATC are limited as soon as the first constraint is reached. In our case, it is possible to congest physically only line 6-7.

Since there is a difference in zone prices whereas the network is not congested, the market outcome is suboptimal with a social overcost of 1794 € compared to the optimum with nodal pricing.

B. Combined ATC

We evaluate in this subsection the combined ATC method applying the model of IV.A.2) to the base case in II. We assume that the TSO sets a minimal value of 1 MW to ATC 2-3 and ATC 3-2 . We then obtain the values of ATC in table III assuming that zone 1 has priority for export. The ATC 2-1 and ATC 3-1 are null because some lines are already congested by the priority exchanges. When the above ATC values are transmitted to the market, the exchanges are constrained by the ATC 1-2, 1-3 and 2-3. The zonal market prices are then 1° 41 €/MWh, 2° 40 €/MWh and 3° 48 €/MWh. These prices are quite close to the prices from nodal pricing. Moreover, line 3-6 is the only one to be congested. The market equilibrium with the Combined ATC with priority to zone 1 is very close to the optimum with nodal pricing because the social overcost is only 9€.

If now, we assume that the TSO makes a mistake in anticipating the market outcome, he may give priority to export from zone 2. Table IV then gives the new values of ATC. The priority given to the export from zone 2, the possible exports from zone 1 are smaller. When the above ATC values are transmitted to the market, the cross-border exchanges are constrained by the ATC 1-2, 1-3 and 2-3. The market price of each zone is then 1° 34 €/MWh, 2° 39 €/MWh and 3° 59 €/MWh. With this market equilibrium, no line is congested. This is because the priority is given to the export from zone 2 while zone 2 is not the cheapest one. We obtain similar results if zone 3 has priority. Since there is difference in zone prices whereas the network is not congested, the market outcome is suboptimal with a social overcost of 2069 € compared to the optimum with nodal pricing. The performances of the Combined ATC method then depend a lot on the assumption giving priority to a zone for export.

C. Critical Branch flow-based

We evaluate in this subsection the Critical Branch flowbased method applying the model of IV.B.1) to the base case in II. With this method, we obtain the following equilibrium prices, 1° 42 €/MWh, 2° 39 €/MWh, 3° 48 €/MWh which are quite close to the prices with nodal pricing. The market equilibrium is close to the optimum with nodal pricing because the social overcost is only 17 €.

Line 3 to 6 constrains the cross-border exchanges. Zone 1 exports 191 MW compared to base case (where there is no cross-border exchange). Zone 2 imports 8 MW. And zone 3 imports 183 MW. The flow on line 3-6 is higher than the maximal capacity by 0.3 MW. The flows can exceed the limit of capacity because of the zonal approximation that consists in estimating the GSK matrix.

The GSK matrix is only an estimation of the participation of generators to the variation of the net export situation of their zone. Once the market equilibrium is established with the Critical Branch flow-based method, it is possible that the realised participation of generators to the net export does not correspond exactly to the anticipated GSK matrix. In this situation, the application of the inaccurate GSK matrix distorts the capacity constraints in the Critical Branch flow-based method. Estimating precisely the GSK matrix corresponding to the reaction of generators to the market price is then important so that the network security can be respected with the Critical Branch flow-based method. Otherwise, the TSO may have to redispatch generation in real time. We haven't seen this problem for the other ETSO-EuroPEX methods because the constraints are not so tight in the other methods.

The Critical Branch flow-based method has also the advantage to reveal precisely where the constrained lines are on the network because the limits on the lines are directly transmitted to the market. The TSOs and the regulators can then directly target their actions toward the lines that are identified as the ones constraining the cross-border exchanges.

Nevertheless, the Critical Branch method has similar drawbacks compared to the other ETSO-EuroPEX methods. First, it is always needed to anticipate to some extent the market equilibrium in order to insert in the Critical Branch method an accurate estimation of the GSK matrix. Second, this method does not always ensure an optimal treatment of constraints on lines that are in a price zone. Indeed, if this internal congestion can be relieved by an internal generator and another one that can be external or internal to the zone, in all the ETSO-EuroPEX methods, this second generator will be systematically external to the zone. Indeed, from the point of view of the algorithm of the Critical Branch method, the joint action of two generators in the same zone (one of them increases its production while the other one decreases it) on any line is null because these two generators have the same zonal PTDF.

VI. CONCLUSIONS AND FUTURE RESEARCHES

In this paper, we have modelled and evaluated the technical and economic efficiency of the main methods proposed by ETSO and EuroPEX to allocate interconnection capacity in Europe. We can draw two main conclusions from the models and analyses of the ATC-and flow-based allocation methods. Firstly, our results about the relative efficiency of the allocation methods are more classical than the preliminary conclusions of the Implementation Study [START_REF] Apx | Implementation Study. A report for the MoU signatories on the design of the market coupling solution in the Central West European (CWE) region[END_REF]. The allocation methods are all the more efficient as the Kirchhoff laws are integrated more precisely in the market clearing. In particular, the gains from a flow-based allocation method are mainly generated when the transmission capacities are precisely related to powerlines rather than related to aggregated border capacities only. Secondly, we find that the paradigm according to which the price zones fit the national boundaries and shared by all the allocation methods of ETSO and Europex has important consequences. The zonal vision of the grid may result in the capacity limits of powerlines being exceeded. Besides, this zonal vision of the grid allows to manage internal congestions with only limited efficiency.

Considering our conclusion, future researches will then focus on two points. First, we will evaluate the influence of the accuracy of the GSK matrix on the technical and economic efficiency of the ETSO-EuroPEX methods. Second, we will study the incentives of the generators and TSOs to make the market converge toward this estimation of the market outcome.

Fig. 1 .

 1 Fig. 1. Base case to test the ETSO-EuroPEX methods

Fig. 2 .

 2 Fig. 2. Optimal use of the network with nodal pricing III. METHODS TO CALCULATE EXCHANGE CAPACITIES

Equations ( 3 )

 3 to (5) model an implicit ATC-based congestion management scheme.

  Table I describes their consumption. Total demand in zone 1 (respectively 2 and 3) is then 330 MW (respectively 170 MW and 530 MW).

TABLE I

 I

				NODAL LOAD
	Node	1	2	3	4	5	6	7
	Load (MW)	125	115	170	90	165	190	255

TABLE II

 II

	ATC VALUES
	ATC 1-2 1-3 2-3 2-1 3-1 3-2
	MW 4	60	1	83	1	39
	Line with contractual congestion 6-7 3-6 6-7 1-3 6-7 3-6

TABLE III

 III

	VALUES OF COMBINED ATC WITH PRIORITY TO ZONE 1
	ATC 1-2 1-3 2-3 2-1 3-1 3-2
	MW 7 178 1	0	0	1
	Line with contractual congestion 3-6 3-6 set 3-6 6-7 set

TABLE IV

 IV

	VALUES OF COMBINED ATC WITH PRIORITY TO ZONE 2
	ATC 1-2 1-3 2-3 2-1 3-1 3-2
	MW 7	48	1 247 0	1
	Line with contractual congestion 6-7 3-6 set 1-3 1-3 set


Ruslan Sadykov
email: ruslan.sadykov@inria.fr

A polynomial algorithm for a simple scheduling problem at cross docking terminals

At cross docking terminals, products from incoming trucks are sorted according to their destinations and transferred to outgoing trucks using a small temporary storage. Such terminals allow companies to reduce storage and transportation costs in supply chain. This paper focuses on the operational activities at cross docking terminals.

We consider the trucks scheduling problem with the objective to minimise the storage usage during the product transfer. We show that a simplification of this NP-hard problem in which the arrival sequences of incoming and outgoing trucks are fixed is polynomially solvable and propose a dynamic programming algorithm for it.

Introduction

Cross docking terminal is a distribution center carrying a considerably reduced amount of stock in contrast to traditional warehouses. Incoming shipments delivered by incoming trucks are unloaded, sorted and loaded onto outgoing trucks waiting at the dock, which forward the shipments to the respective locations within the distribution system. Compared to traditional warehousing, a cost intensive storage and retrieval of goods is eliminated by a synchronization of inbound and outbound flows. An additional advantage of cross docking is efficient usage of truck capacity (i.e. full loads) and the implementation of a good scheduling system [START_REF] Apte | Effective cross dockng for improving distribution efficiencies[END_REF].

In this paper, we consider a simplified cross docking terminal with one receiving, one shipping door and a storage. The incoming trucks fully loaded by products of different types arrive at the receiving door where the products are unloaded. Each truck leaves the platform when it is fully unloaded. Each unloaded product is processed and transferred to the shipping door if the outgoing truck currently staying there demands the products of this type and does not yet have enough of them. Otherwise, the unloaded product is transferred to the temporary storage. Each outgoing truck can leave the door only if it is fully loaded by products moved either directly from incoming trucks or from the storage. The objective is to find the best arrival (or departure) order for the whole set of incoming and outgoing trucks to increase the efficiency of the cross docking terminal.

There are several papers in the literature which dealt with this scheduling problem. Yu and Egbelu [START_REF] Yu | Scheduling of inbound and outbound tracks in cross docking systems with temporary storage[END_REF] developed a model for scheduling incoming and outgoing trucks to minimise the makespan (i.e. the maximum completion time) and proposed several heuristic algorithms for it. A simpler but similar problem was considered by Boysen, Fliedner and Scholl [START_REF] Boysen | Scheduling inbound and outbound trucks at cross docking terminals. OR Spectrum[END_REF]. They have proposed some lower bounds and an exact decomposition approach for it. Maknoon, Baptiste and Kone [START_REF] Maknoon | Optimal loading and unloading policy in cross docking platform[END_REF] concentrated on the objective function which minimizes the storage cost in a simplified setting in which the sequences of incoming and outgoing trucks are fixed and any outgoing truck demands the products of only one type.

We now define the problem formally. The number of incoming and outgoing trucks are, respectively, n and m. The number of different product types is T . We will denote the i-th incoming truck as I i and the o-th outgoing truck as O o . There are several types of products. Each incoming truck I i supplies a it products of type t. Each outgoing truck O o demands b ot products of type t. Let also T o be the set of product types demanded by O o , i.e. T o = {t : b ot > 0}. We suppose that each outgoing truck demands products of at most q different types, i.e. q = max 1≤o≤m | T o |. Note that, for each type t, the total number of supplied products of this type should be not less than the total number of demanded products of this type, otherwise at least one outgoing truck would not be able to depart fully loaded.

A product of type t can be moved directly from an incoming truck at the receiving door to an outgoing truck O o at the shipping door if the number of products of type t which are already in O o is less than b ot . Otherwise the product should be transferred to the temporary storage at a cost c t . One product of type t occupies a volume d t in the storage. The total volume of products in the storage should not exceed its total capacity which is D. The problem consists in finding a policy for unloading and loading products such that the total cost is minimized. For convenience, an equivalent objectif is used which is the maximization of the total cost of the directly transferred products. This problem is NP-hard in the strong sense even for a very restricted case, as shown in Appendix A.

Economically, the storage cost can be also interpreted as the difference between the time needed to transfer a product directly from an incoming to an outgoing truck and via the intermediate storage. So, our objective function is equivalent to the total processing cost. Although this objective is similar to the makespan objective used in [START_REF] Boysen | Scheduling inbound and outbound trucks at cross docking terminals. OR Spectrum[END_REF] and [START_REF] Yu | Scheduling of inbound and outbound tracks in cross docking systems with temporary storage[END_REF], the problems studied there are quite different. The main difference is that we do not allow any concurrent operations. In contrast to this, for example in [START_REF] Yu | Scheduling of inbound and outbound tracks in cross docking systems with temporary storage[END_REF], on can at the same time discharge a product from an incoming truck to the storage and load a product from the storage onto an outgoing truck. As a result, transferring a product via the storage does not always increase the objective function.

In the rest of the paper we concentrate on the problem in which the sequences of incoming and outgoing trucks are fixed. Formally this means that, for each i ∈ {1, . . . , n -1}, the truck I i departs before the arrival of I i+1 , and for each o ∈ {1, . . . , m -1}, the truck O o departs before the arrival of O o+1 .

Note that fixing sequences of incoming and outgoing trucks in problems considered in [START_REF] Boysen | Scheduling inbound and outbound trucks at cross docking terminals. OR Spectrum[END_REF] and [START_REF] Yu | Scheduling of inbound and outbound tracks in cross docking systems with temporary storage[END_REF] make them trivial, and thus polynomially solvable. However, determining the computational complexity of our problem with fixed sequences is not trivial at all. In [START_REF] Maknoon | Optimal loading and unloading policy in cross docking platform[END_REF], Maknoon at al. considered this problem with some restrictions including q = 1, but left the complexity question open.

In this paper, we answer this question by presenting a polynomial dynamic programming algorithm for the problem considered.

Preliminary observations

In the problem, we need to find an optimal policy which is characterized, among others, by the departure order of trucks. Note that the departure order is fixed within the set of incoming or outgoing trucks but not fixed within the set of all trucks. A departure order of all trucks can be represented by a path in the graph depicted in Figure 1. In this graph, each node (i, o) represents the situation in which the trucks I i and O o are at the doors.

The path shown on this picture corresponds to the policy which is characterized by the following departure order of trucks:

(I 1 , I 2 , I 3 , O 1 , O 2 , I 4 , I 5 , O 3 , O 4 , I 6 , . . . , I n-1 , O 6 , . . . , O m-1 , I n , O m ).
Some paths in the departure order graph are not feasible due to the following storage constraints. First, for each t, the number of products of type t in the storage should always be non-negative (all demanded products should be first supplied). Second, the overall volume of products in the storage can never exceed its capacity. We say the a node is infeasible if it does not belong to any feasible path. Formally, a node (i, o) is infeasible if and only if

∃t : i k=1 a kt < o-1 p=1 b pt or T t=1 d t • max 0, i-1 k=1 a kt - o p=1 b pt > D. (1)
From ( 1) it follows that The next fact is very important and serves as the base for the algorithm.

• • • • • • • • • • • • • • • • • • • • • . . .
Observation 1 There exists an optimal policy in which, each time trucks I k and O j are at the doors, for each t, I k transfers directly to O j as many products of type t as possible, i.e. the minimum between the number of products of type t still available in I k and the number of products of type t which are still demanded by O j .

It is easy to see that this observation is correct. Suppose, in an optimal policy, I k does not transfer to O j z products of type t it can and "saves" them for consequent outgoing truck(s). Then O j is obliged to take these z products from the storage. In the modified policy, these z products are transferred directly from I k to O j , and products, transferred directly from I k to consequent outgoing truck(s) in the original policy, are taken from the storage. The cost of the modified policy which comply with the observation do not increase. We will call a policy which complies with Observation 1 direct first. Each path in the departure order graph corresponds to exactly one direct first policy.

The algorithm

We now present a dynamic programming algorithm for the problem which finds an optimal direct first policy.

Each state in our dynamic programming algorithm does not correspond to a node in the departure order graph but to an edge, which a path can follow just after "turning". To clarify the presentation, we first present in Figure 2 the underlying directed graph of the dynamic programming algorithm. For each pair of incoming truck I i and outgoing truck O o , there is the "square node" (i, o) and the "circle node" (i, o). Each node corresponds to a set of states of the algorithm, each edge defines all possible moves between a state in the origin of this edge and a state in the destination of this edge. The paths shown in Figures 1 and2

• • • • • • • • • • • • • • • • • • • • • . . . (1, 1) (4, 1) (4, 3)
(5, 3)

(5, 5) In our dynamic programming algorithm, for each "square node" (i, o), we define the states S out (i, o, {f t } t∈To ), where f t is the number of products of type t transferred directly by truck

(n, 6) (n, m)
I i to O o .
Each "circle node" (i, o) represents the situation in which truck O o has departed, truck O o+1 is going to arrive, truck I i is at the shipping door, and I i had arrived after the departure of the truck O o-1 . This means that I i could ship products directly only to O o .

In our dynamic programming algorithm, for each "circle node", we define the states S inc (i, o, {f t } t∈To ), where f t is the number of products of type t transferred directly by truck

I i to O o .
Now, every policy corresponds to exactly one sequence of states or a states path. Therefore, in the following, we will use the same notation for a states path and the corresponding policy. A path which corresponds to a direct first policy, we will call a direct first path. For each state S out (i, o, {f t } t∈To ) and S inc (i, o, {f t } t∈To ), we keep and update the objective function value of the best path to this state. Let V out (i, o, {f t } t∈To ) and V inc (i, o, {f t } t∈To ) be these values.

To solve the problem, we need to find the best path terminating at a state S out (n, m, {f t } t∈To ) or S out (n, m, {f t } t∈To ).

To simplify the presentation, when there is no ambiguity, we will use the shortened notations S inc (i, o, f ), S out (i, o, f ), V inc (i, o, f ), and V out (i, o, f ), where f is a vector.

Note that, in the algorithm, we handle only the states belonging to at least one direct first path. The total number of such states can be different even for instances of the same size.

Suppose we are in a state S inc (i, o, f ). From this state it is only possible to move to a state S out (i, o , f ), where max{o + 1, f o(i + 1)} ≤ o ≤ lo(i) and f o(n + 1) = m. When we make such a move, truck I i transfers directly to every truck O j , o < j ≤ o , as much products as possible.

The formal procedure for making moves from a state S inc (i, o, f ) is presented in Algorithm 1. The complexity of this procedure is O m(q + ρ) , where ρ is the number of operations needed to check if a state has been already encountered earlier in the algorithm. We will compute ρ later in the paper.

Suppose now we are in a state S out (i, o, f ). From this state it is only possible to move to a state S inc (i , o, f ), where max{i + 1, f i(o + 1)} ≤ i ≤ li(o), and f i(m + 1) = n. When we make such a move, truck O o receive directly from every truck I k , i < k ≤ i , as much products as possible.

The formal procedure for making moves from a state S inc (i, o, f ) is presented in Algorithm 2. The complexity of this procedure is O n(q + ρ) .

In the full algorithm, presented as Algorithm 3, we look through all the created states and make all possible moves from them as described above. To obtain an optimal policy, it suffices to store, for each state, along the value V , the path which gives it. At the end of the algorithm, the best path stored for one of the states S out (n, m, f ) and S inc (n, m, f ) corresponds to an optimal policy.

for t = 1 to T do r[t] ← a it ; 1 for t ∈ T o do r[t] ← a it -f t ; 2 v ← V inc (i, o, f ); 3 for j ← o + 1 to lo(i) do 4 for t ∈ T j do 5 dt[t] ← min{r[t], b jt }; 6 v ← v + c t • dt[t]; 7 if state S out (i, j, {dt[t]} t∈T j ) does not exist then 8 create it: V out (i, j, {dt[t]} t∈T j ) ← -∞; 9 if j ≥ f o(i + 1) and v > V out (i, j, {dt[t]} t∈T j ) then 10 V out (i, j, {dt[t]} t∈T j ) ← v; 11 for t ∈ T j do r[t] ← r[t] -dt[t];
12 Algorithm 1: Algorithm to make moves from a state S inc (i, o, f ). As the base of the reduction for the claim 1, we can take the case i = 2. In this case, we can only have o < o , as, for any fixed o * , there is only one path which terminates at a state S inc (2, o * , f ).

for t ∈ T o do r[t] ← f t ; 1 v ← V out (i, o, f ); 2 for k ← i + 1 to li(o) do 3 for t ∈ T o do 4 dt[t] ← min{b ot -r[t], a kt }; 5 v ← v + c t • dt[t]; 6 if state S inc (k, o, {dt[t]} t∈To ) does not exist then 7 create it: V inc (k, o, {dt[t]} t∈To ) ← -∞; 8 if k ≥ f i(o + 1) and v > V inc (k, o, {dt[t]} t∈To ) then 9 V inc (k, o, {dt[t]} t∈To ) ← v; 10 for t ∈ T o do r[t] ← r[t] + dt[
Analogously, as the base of the reduction of the claim 2, we can take the case o = 2.

Proposition 1 The total number of created states in Algorithm 3 is O(qnm 2 ).

Proof. In the following, when a path contains a state S out (i, o, f ) or a state S inc (i, o, f ), we will say that, for each type t,

• if t ∈ T o ,
ν out i io t * (f ) = max 0, min a it * , b o t * -f - i-1 k=i +1 a kt * .
Therefore, path P contains a triple (i , o , ν out i io t * (f )) inc t * , where i > i . An important observation is that there is at most one value i for which such triple is not t * -canonical. We can say that any triple (i , o , f ) out t * "generates" at most one triple (i , o , f ) inc t * such that i > i which is not t * -canonical. Consider now a direct first path P which contains a triple (i , o , f ) inc t * . In the corresponding policy, truck I i transfers directly to every truck O o , o < o ≤ o , ν inc i o ot * (f ) products of type t * , where

ν inc i o ot * (f ) = max    0, min    b ot * , a i t * -f - o-1 j=o +1 b jt *       .
Therefore, path P contains a triple (i , o , ν inc i o ot * (f )) out t * , where o > o . Again, an important observation is that there is at most one value o for which such triple is not t * -canonical. We can say that any triple (i , o , f ) inc t * "generates" at most one triple (i , o , f ) out t * such that o > o which is not t * -canonical. We now fix a value o * . From the above analysis, it follows that every t *canonical triple (i, o, f

) inc t * or (i, o, f ) out t * such that o ≤ o * "generates" at most one non-t * -canonical triple (i , o * , f ) inc t * or (i , o * , f ) out t * such that i > i. Therefore, the total number of different triples (i, o * , f ) inc t * or (i, o * , f ) out t * contained in paths in P inc (i, o * , f ) and P out (i, o * , f ) is O(nm).
A straightforward analysis now indicates that the total number of states in the algorithm does not exceed O((nm 2 ) q ). However, it is possible to get rid of the exponent q by the following reasoning.

We fix values i * and o * , and define a complete lexicographic order for states S out (i We now return to the question of estimation of ρ, which is the number of operations needed to check whether a state S out (i, o, f ) or S inc (i, o, f ) has been already created. From the proof of Proposition 1, remember that, given fixed values i * and o * , the number of different states S out (i * , o * , f ) is not more than O(qnm). As these states can be lexicographically ordered, the storage and the search of these states can be done using a binary tree. Therefore, to check whether a state S out (i, o, f ) or, analogously, a state S inc (i, o, f ), has been already created and retrieve the value V out (i, o, f ) or V inc (i, o, f ) we need O log(qnm) operations. This concludes the estimation of the complexity of the algorithm.

Proposition 2

The complexity of the dynamic programming algorithm is O qnm 2 (n + m)(q + log m + log n) .

Conclusion

In this paper, we presented a polynomial dynamic programming algorithm for a scheduling problem with fixed sequences of incoming and outgoing trucks to minimise the storage cost in cross docking terminal. This algorithm allows to determine the computational complexity of the problem for the first time.

Although polynomial, the complexity of the algorithm is high. However, a theoretical estimation of the number of created states can be much higher than what we have in reality. Therefore, we hope that the algorithm may be used in practice especially when coupled with some bounding procedure and when q is small. As an indirect evidence of this, we can mention good experimental results reported in [START_REF] Maknoon | Optimal loading and unloading policy in cross docking platform[END_REF] for an algorithm proposed for the case q = 1 of our problem.

One interesting direction for a future research is to try to find a linear programming formulation for the problem, since it was shown to be polynomially solvable. Such a formulation would help a lot in developing methods for solving more practical generalisation of the problem in which the truck sequences are not fixed.

Given an instances of the 3-partition problem, we now define the corresponding instance of our cross docking problem. There are 3n incoming, 4n outgoing trucks (3n of the first type, n of the second type) and two types of products. The supplies and demands are the following:

a i1 = 1, i = 1, . . . , 3n, a i2 = 2n + r i , i = 1, . . . , 3n, b i1 = 1, i = 1, . . . , 3n, b i2 = 0, i = 1, . . . , 3n, b i1 = 0, i = 3n + 1, . . . , 4n, b i2 = 6n + B, i = 3n + 1, . . . , 4n.
We claim that there exists a 3-partition if and only if at most n products are transferred via the storage.

Suppose that there exists a 3-partition {A 1 , A 2 , . . . , A n }, where A j = {i j1 , i j2 , i j3 }. Then, the trucks are sequenced in n groups. Group j, 1 ≤ j ≤ n, includes set A j of incoming trucks and outgoing trucks O 2j-1 , O 2j and O 3n+j . The departure order of group j is the following

O 2j-1 , I ij1 , I ij2 , O 3n+j , I ij3 , O 2j .
As r ij1 + r ij2 + r ij3 = B, the incoming trucks transfer 6n + B products of type 2 directly to truck O 3n+j , I ij1 transfers 1 product of type 1 directly to O 2j-1 , and I ij3 transfers 1 product of type 1 directly to O 2j . Only 1 product of type 1 is transferred to the storage from I ij2 . These transfers are depicted in Figure 3. As there are n groups, n products in total are transferred to the storage and then put to the outgoing trucks O 2n+1 , . . . , O 3n , which are sequenced at the end.

I ij2 I ij3 I ij1 storage O 3n+j O 2j-1 O 2j 1 1 1 r i j 1 +
2 n ri j2 + 2n Suppose now there is a sequence of trucks such that at most n products are transferred through the storage. Then, every outgoing truck O i , i = 3n + 1, . . . , 4n must receive products directly from at least 3 incoming trucks.

r i j 3 + 2 n
Otherwise it would receive from the storage at least 2n products of type 2. Therefore, at least 1 product of type 1 goes to the storage while each such truck O i is supplied, and the total number of products of type 1 transferred through the storage is at least n, meaning that all products of type 2 should be transferred directly. Then, each incoming truck can transfer directly products of type 2 to exactly one outgoing truck. Otherwise, between two outgoing trucks of type 2, only one outgoing truck can receive a product of type 1 directly, and the total number of products of type 1 transferred through the storage would exceed n. We conclude that there should exist a partition of incoming trucks into triples {A 1 , A 2 , . . . , A n } such that i∈Aj r i = B. Otherwise there would exist a triple A j such that r ij1 + r ij2 + r ij3 < B, and the outgoing truck which is supplied by the incoming trucks in A j would need to take at least one product of type 2 from the storage.

Figure 1 :

 1 Figure 1: The departure order graph and an example of the path

Figure 2 :

 2 Figure 2: The underlying directed graph for the dynamic programming algorithm

Figure 3 :

 3 Figure 3: Product transfers within a group j of trucks

  correspond to the same policy. Again, the infeasible nodes are shown in black.

					outgoing trucks
		1	2	3	4	5	6	• • •	m
		1
		2
		3
	incoming trucks	4 5
		6
		. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .
		n

  t]; Algorithm to make moves from a state S out (i, o, f ). all trucks O o , o ≤ o < o * , in policy P . As P and P are direct first policies, there are two possible cases: • o = o and f t > f t . Then, as P -∈ P inc (i * , o , f ) and P -∈ P inc (i * , o , f ), by claim 2, f t ≥ f t for all types t ∈ T o . • o < o . In both cases, as P and P are direct first policies, for every t, truck I i * does not transfer more products of type t to all trucks O o , o ≤ o < o * in policy P than to all trucks O o , o ≤ o < o * , in policy P , and f t ≤ f t for all types t ∈ T o * . Analogously, if we suppose that claim 1 is true for i < i * and o ≤ o * , we can prove the claim 2 for i = i * and o = o * .

	11
	Algorithm 2:

  the path contains triple (i, o, f ) out Consider a direct first path P which contains a triple (i , o , f ) out t * . In the corresponding policy, truck O o receives directly from every truck I i , i < i ≤ i ,

	t	or (i, o, f ) inc t ;
	• if t ∈ T o , the path contains triple (i, o, 0) out t	or (i, o, 0) inc t .
	Let now fix type t ν out i io t

* . We will call triples (i, o, 0) out t * , (i, o, b ot * ) out t * , (i, o, a it * ) inc t * , and (i, o, 0) inc t * t * -canonical. There are O(nm) such triples in total. * (f ) products of type t * , where

  , o * , f ). A state S out (i * , o * , f ) is lexicographically smaller than a state S out (i * , o * , f ), f = f , if and only if f t ≤ f t for all t ∈ T o . It is always possible to compare in this way two such states, as we cannot have f t < f t and f t > f t for any two types t , t ∈ T o * by Lemma 1. Now, the total number of states S out (i * , o * , f ) is equal to the sum, for every type t ∈ T o * , of different triples (i * , o * , f ) out t , as, in order to pass from a state S out (i * , o * , f ) to the lexicographically next state, at least one of the values f t , t ∈ T o * should be increased. Analogously, the total number of states S inc (i * , o * , f ) is equal to the sum, for every type t ∈ T o * , of different triples (i * , o * , f ) inc t . Consequently, for every fixed value o * , the total number of states S out (i, o * , f ) and S inc (i, o * , f ) is equal to O(qnm). From this, we conclude that the total number of states is O(qnm 2 ).

*

run Algorithm 1 for all created states S inc (i, o, f );

run Algorithm 2 for all created states S inc (i, o, f );

We will now determine the complexity of the algorithm. The critical point here is estimating the total number of created states in the algorithm.

To do it, we will need the following lemma. But first we introduce some additional notations. Let P out (i, o, f ) and P inc (i, o, f ) be the sets of paths which teminate at states S out (i, o, f ) and S inc (i, o, f ). Let also P out (i, o, f, o , f ) be the set of paths which terminate at state S out (i, o, f ) and contain state S inc (i, o , f ) such that o < o. Similarly, let P inc (i, o, f, i , f ) be the set of paths with terminate at state S inc (i, o, f ) and contains state S out (i , o, f ) such that i < i. Finally, for a path P , we will denote as P -its sub-path in which the terminating state of P is excluded.

Lemma 1

1. For any two direct first paths P ∈ P out (i, o, f ) and

2. Analogously, for any two direct first paths P ∈ P inc (i, o, f ) and

Proof. We will prove this lemma by induction. Suppose that the claim 2 is true for i ≤ i * and o < o * . We will prove the claim 1 for i = i * and o = o * . Without loss of generality, let P ∈ P out (i * , o * , f , o , f ) and P ∈ P out (i * , o * , f , o , f ), f = f . As f t < f t , truck O o * receives less products of type t from I i * in policy P than in policy P . Therefore, I i * transfers to trucks O o , o ≤ o < o * more products of type t in policy P than to

Appendix A NP-hardness proof

Here we consider the general problem in which the sequences of trucks are not fixed. We show the this problem is NP-hard in the strong sense even for the case in which each incoming truck supplies products of at most two types, each outgoing truck demands products of at most one type, all the storage costs and product volumes are unitary, and the storage capacity is unlimited. We will perform a reduction from the 3-partition problem.

Remember that, in the 3-partition problem, we are given an integer B and a set of 3n integers r 1 , r 2 , . . . , r 3n such that 3n i=1 r i = Bn and B/4 < r i < B/2 for each i. We need to decide whether there exists a partition of the set of indexes {1, 2, . . . , 3n} into n sets {A 1 , A 2 , . . . , A n } such that i∈Aj r i = B, ∀j = 1, . . . , n. Note that, if such a partition exists, each subset A j contains exactly 3 indexes.


Mfoihaya Bedja

Isabelle Magne
email: isabelle.magne@edf.fr

Martine Souques

Jacques Lambrozo

Laurent Le Brusquet

Gilles Fleury

Alain Azoulay

Stéphane Ruszczynski

French Population Exposure to 50 Hz Magnetic Fields: Intermediate Results

INTRODUCTION

In 1979, a study indicated that childhood cancers might be linked with extremely low frequency (ELF) magnetic field [START_REF] Wertheimer | Electrical wiring configurations and childhood cancer[END_REF]. Leaning on numerous epidemiological studies, the International Agency for Research on Cancer classified in 2001 ELF MF as "possibly carcinogenic to human". These conclusions were based on a statistical association found in some epidemiological studies, unconfirmed by experimental results, between childhood leukaemia risk and a mean exposure over 24h higher than 0.4µT.

One remaining -and critical-question about these epidemiological results is the exposure assessment. Another relevant question is to know, in the overall exposure of people, the relative contribution of all the possible field sources in the daily environment.

Starting from a previous study (50 Hz magnetic field levels in houses, half of them close to power lines [START_REF] Clinard | Residential magnetic field measurements in France: comparison of indoor and outdoor measurements[END_REF]), the French Ministry of Health initiated in 2006 a large study of the individual exposure of a randomly selected sample of the French population.

The total exposure database will contain 1000 measurements for children and 1000 measurements for adults. The analysis of results is still ongoing, and should be published in 2010. This paper presents intermediate results, based on the data already analysed.

MATERIALS AND METHODS

The exposure data were collected during 3 measurement campaign between the beginning of 2007 and the end of 2008.

The volunteers were randomly selected from phone lists. Each of them has worn an EMDEX II (Enertech, USA) MF recorder during 24h, and has simultaneously filled in a timetable describing his activities. In addition, all volunteers answered to a questionnaire about his socioprofessional data and his house (one of the parents answered for their child).

The present paper gives the intermediate results on 990 persons : 437 children (0-14 years) and 553 adults.

RESULTS

The observed arithmetic and geometric means are respectively 0.11µT and 0.03µT for children, 0.15µT and 0.04µT for adults.

The statistical analysis has led to the following results :

-The search for the most discriminating MF descriptors from a hierarchical clustering classification followed by CART method (Classification And Regression Tree)led to a distribution of the studied sample into three groups for each type of population (adults and children). The most discriminating descriptors are o the arithmetic and geometric means and the Rate Change of metric standardized (RCMS) for both the two populations, o the maximum and the median values for the adults, o the third quartile and the standard deviation for the children.

-The factors that lead a person to be more exposed has been identified by linear and logistic regression methods or non-parametric regression models: o for both the two populations: the time spent in railway transport, the time spent using a computer, to sleep near a clock radio, to live in a collective building, o the age for the children, o to live in a city of more then 2.000 inhabitants and to have an individual electric water heating for the adults.

CONCLUSIONS

The most important factor influencing the personal exposure is the clock radio, which was identified in 120 series. The question whether these measurements are representative of personal exposure will be discussed.

The influence of proximity to electric transmission and distribution facilities will be analysed.

ACKNOWLEDGMENTS

This study was funded by the Ministry of Health and Solidarities and realized by Supélec, with the collaboration of EDF and RTE.