Yoann Fouquet
email: yoann.fouquet@hds.utc.fr

Dritan Nace
email: dritan.nace@hds.utc.fr

Michal Piόro

Michael Poss
email: michael.poss@hds.utc.fr

Flow Adjustmenta Flexible Routing Strategy for Demand Protection Against Multiple Partial Link Failures

Keywords: survivable networks, multiple partial link failures, multicommodity flow, linear and mixed integer programming, path generation. I

In this paper, we study a flexible routing strategy for demand protection and a corresponding optimization problem for networks that permanently experience fluctuations of the capacity available on their links. This is an important and novel topic as limited link availability is a fundamental feature of wireless networks; yet majority of work in survivable network design is restricted to total failures of single links. Hence, protection against partial failures of multiple links is considered as congestion avoidance. We assume a given finite set of network states. Each state is characterized by a vector of link availability coefficients specifying, for each link, the fraction of its nominal (maximum) capacity available in this state, and by a traffic coefficients vector specifying, for each demand, the proportion of its nominal traffic to be realized in the considered state. Our routing strategy allows for adjustment (thinning or thickening) of the reference pathflows. For a given nominal value x of a path-flow, its thickening is limited to x where  is a given constant greater than or equal to 1. Thus, in each state, the value of every pathflow can range from 0 to times its reference value. It turns out that the corresponding link cost minimization problem (where link capacities and state-dependant path-flows are decision variables) is NP-hard. We present a non-compact linear programming formulation of the problem together with a solution algorithm based on path generation. We illustrate the effectiveness of the introduced routing strategy by presenting numerical results for a set of representative network examples.

INTRODUCTION

In this paper, we study a routing strategy for demand protection called Flow-Adjustment Routing (FAR) together with the corresponding network optimization problem. FAR aims at protecting traffic demands against Multiple Partial Link Failures (MPLF) by appropriate adjustment of the reference (nominal) path-flows. Strictly speaking, FAR is not a rerouting strategy because no paths can be created to restore the affected flows. Instead, the nominal path-flows are thinned or thickened according to the current state of the links. To account for MPLF, we consider a nominal state and a set of link availability states (containing the nominal state) and a finite set A fractional numbers. For each link e ∈ E and each availability state ∈ , the corresponding link availability coefficient ∈ A specifies the fraction of the nominal capacity of link e that is available in state s ( . FAR assumes that the nominal path-flows (i.e., path-flows defined for the nominal state with all links fully available) can be thinned or thickened to adapt to capacity fluctuations but not restored on paths not used in the nominal state. On top of this, we assume that the demand volumes to be realized in the availability states different from the nominal state are possibly reduced as compared to the nominal traffic. In the following, we use the terms availability states and failure states interchangeably.

For the so-specified setting, we study the corresponding optimization problem referred to as Flow-Adjustment Problem (FAP). FAP aims at minimizing the cost of links under the following constraints: (i) the nominal path-flows and the state-dependent path-flows satisfy the nominal and the state-dependent traffic demand requirement, respectively, (ii) path-flows obey the maximum thickening assumption of the protection strategy, and (iii) for each link and each availability state, the link load does not exceed the currently available link capacity. As FAP is NP-hard (because a special case of FAP shown to be NP-Hard is studied by Tomaszewski et al. [START_REF] Tomaszewski | On the complexity of resilient network design[END_REF]), its Linear Programming (LP) formulations (see Section III) are unavoidably non-compact and require path generation to be solved to optimality. For generating paths we need to apply a mixed-integer pricing problem-here the difficulty of FAP is manifested.

The introduced model is original. In fact, not much work has been done in survivable network optimization under the MPLF assumption. To the best of our knowledge, only the so called global rerouting (restoring flows from scratch in surviving capacity) has been studied in this context [START_REF] Piόro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF] a strategy quite impractical in the considered wireless framework. It happens that multiple failures have so far been considered in survivable network design almost always (i.e., besides global rerouting) assuming total link failures ( ∈ { }). The work done here is substantial (again, see [START_REF] Piόro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF][5] and the references therein) but this case is, unfortunately, not relevant for our framework.

The rest of this paper is organized as follows. In the next section, we give motivation behind the considered flowadjustment strategy. In Section III, we introduce the notation and formulate the basic FAP model. Section IV describes a path generation algorithm necessary to find optimal solutions for FAP. In Section V, we present results of a numerical study that illustrates a potential value of the considered protection strategy in terms of the cost of the resulting network in comparison with global rerouting. In Section VI, we present the performance of the restoration process using simulations. Finally, in Section VII, we present concluding remarks.

II. FLOW ADJUSTMENT-A SIMPLE PROTECTION STRATEGY FOR MULTIPLE PARTIAL FAILURES

The flow-adjustment routing strategy for MPLF proposed in this paper extends the Flow Thinning Routing strategy (FTR) presented in [START_REF] Piόro | An optimization model for communication networks with partial multiple link failures[END_REF] [START_REF] Piόro | On Protected Traffic Routing in Wireless Networks with Partial Multiple Link Failures[END_REF]. Both FTR and FAR are inspired by the idea of Elastic Rerouting (ER) presented in [START_REF] Shinko | A study on a distributed rerouting scheme[END_REF] and Path Diversity Protection (PDP) [START_REF] Piόro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF], [START_REF] Orlowski | On the complexity of column generation in network design with path-based survivability mechanisms[END_REF]. We recall that both ER and PDP assume total link failures, i.e., failure scenarios which admit only binary vectors characterizing the availability states (referred to as the failure states in this case) in (i.e., A = {0,1}), that is, when links fail, then their entire capacity is lost. PDP is a protection strategy where in a failure state ∈ path-flows through the failing links are simply disconnected, and the surviving path-flows must be sufficient to realize the demand volumes, possibly decreased with respect to nominal demand volumes. ER allows for decreasing the flows of unaffected demands, as well as for increasing (to a certain threshold) path-flows of affected demands.

In this paper, we follow the concept of ER, but no distinction will be made between unaffected and affected demands. Hence, our strategy is allowed to thin or thick the nominal path-flow on any routing path. Thus, FAR assumes that each demand is in general routed over several paths, not necessarily disjoint, with over-dimensioned nominal pathflows to ensure an assumed level of demand survivability. Contrary to conventional end-to-end restoration strategies, no flow is rerouted to recover from a failure. In fact, it seems that allowing only for thinning makes the resulting protection strategy less complicated for network implementations. Our strategy is an extension of FTR in which the path-flows are allowed only to be decreased. Instead, FAR admits both thinning and thickening allowing for less costly networks as compared with FTR. Yet, FAR and FTR are similar from the mathematical modeling point of view, especially the mixed-integer pricing problem.

To clarify the idea of how FAR works, we wish to emphasize a few points. First, note that a common pool of link capacities is used for the nominal path-flows and in the remaining availability states. Hence, the selected routing paths carrying the flows are dimensioned so that the total traffic realized by the demand's path-flows could in general be greater than the nominal traffic. The most important feature of FAR is handling partial failures without any flow rerouting at all. In other words, no paths besides the nominal paths are used for handling other availability states. Therefore, the proposed approach results in using a sort of limited dynamic routing, adapted to the network states. To summarize, for each demand there is a fixed set of nominal routing paths carrying nominal flows. In an availability state in general only a part of the total nominal demand flow will be realized on these paths, depending on the available capacity and the required demand restoration ratio. Consequently, the restoration will be practically done by thickening some of the path-flows, and no new (re)routing paths are allowed.

III. FAP: FLOW-ADJUSTMENT PROBLEM

The basic problem considered in this paper is referred to as FAP and is as follows. We minimize the cost of link capacity assuming that in the nominal state of network operation all demand volumes are realized by means of (nominal) path-flows. When the network is subject to a failure from a given set of failure states (we assume that a failure state consists of partial failures of multiple links) then the demand volumes, possibly reduced, are realized for the duration of the failure state by appropriate thinning or thickening of the nominal flows. The detailed formulation of FAP will be given in Subsection III-B.

A. Notation

The considered network is modeled using a graph ( , undirected or directed, composed of a set of nodes and a set of links . In the sequel, we will always consider directed graphs unless stated explicitly otherwise. Thus, each link ∈ represents a directed pair ( of some nodes ∈ , and is assigned a non-negative unit capacity cost which is a parameter, and a capacity reservation which is an optimization variable. The cost of the network is given by the quantity ∑ ∈ . The demands are represented by the set . Each demand ∈ is associated with a directed pair of nodes (end-nodes of ) ( ( ; ( ) for some ( ( ∈ ; a volume (a parameter) has to be sent from ( to ( (demand volumes and link capacities are expressed in the same units). Also, each demand is assigned a set of admissible paths composed of selected elementary paths from ( to ( in graph . (Recall that an elementary path does not traverse any node more than once.) Paths from are used to realize the demand volume by means of path-flows ∈ , which are optimization variables. The given sets of admissible paths define the link-path incidence coefficients ∈ ∈ ∈ , where if path ∈ traverses link ∈ (i.e., if ∈ , treating the paths as subsets of links:

), and if path ∈ does not traverse link . It is important to note that the sets of admissible paths ∈ , are parameters in the FAP problem formulation considered in the sequel, although in general it assumes that all possible elementary paths can potentially be used if this is required to reach the optimum.

Network links are subject to (partial) failures. The failure-less state (with all links fully available) is called the nominal state and is denoted by 0, and the flows ∈ are referred to as nominal flows. The set of failure states is denoted by (we call set the failure scenario). Each failure state ∈ is specified by a vector of link availability ratios ( ∈ , where for ∈ . The link capacity available in state is assumed to be equal to . Hence, is the fraction of capacity that survives in state . Certainly, in state 0 all the ratios are equal to 1. In general, in a failure state more than one link can have its availability ratio less that 1 so in fact we are considering multiple partial failures of the network links.

When a failure state ∈ affects the network, in general not all the nominal flows ∈ ∈ can be realized anymore as the available link capacity is decreased. It is assumed that the demand volumes to be realized in state ∈ can be decreased, to values ∈ . The demand volumes assumed for a failure state must be realized by means of the nominal flows that are appropriately thinned (decreased) or thickened (increased) so to fit to the reduced link capacity.

∑ ∈ ( [ ] ∑ ∈ ∈ ( [ ] ∑ ∑ ∈ ∈ ∈ ( [ ] ∑ ∈ ∈ ∈ ( [ ] ∑ ∑ ∈ ∈ ∈ ∈ ( [ ] ∈ ∈ ∈ ( (
In the formulation, the quantities in brackets on the lefthand sides are dual variables associated with the constraints (see Section I.2.2 in [START_REF] Orlowski | On the complexity of column generation in network design with path-based survivability mechanisms[END_REF]). All these variables are, by assumption, nonnegative. Objective (1a) minimizes the total cost of links. Constraint (1b) makes sure that for each demand its paths have jointly sufficient capacity to satisfy the demand volume assumed for the nominal state. Constraint (1c) does not allow the nominal link loads to exceed the nominal link capacities. Next, constraint (1d) assures that in each failure state ∈ , the adjusted flows are sufficient to realize the (possibly reduced) volume of each demand ∈ . Then, constraint (1e) makes sure that in each state ∈ , the surviving capacity of each link ∈ is not exceeded. Finally, constraint (1f) relates the state-dependent path-flows to the reference flows.

IV. PATH GENERATION

Formulation ( 1) is a link-path LP formulation. It is noncompact because of exponentially many path-flow variables as we potentially consider all possible paths in the admissible path-sets ∈ . Thus, in order to consider all possible routing paths in graph ( in FAP, we need to apply the technique of linear programming known as column generation [START_REF] Lasdon | Optimization Theory for Large Systems[END_REF], [START_REF] Minoux | Mathematical Programming: Theory and Algorithms[END_REF], called Path Generation (PG) in our context. With PG, starting from some initial path-sets ∈ , we generate new paths (corresponding to variables/columns ∈ ∈ and ∈ ∈ ∈ ), and iteratively add them to the path-sets. As discussed below, this is done by solving an appropriate pricing problem using, as parameters, optimal dual variables associated with constraints (1b)-(1f), i.e., an optimal solution of the problem dual to FAP formulated in the next subsection.

A. Dual LP formulation of FAP

The problem dual to LP (1) is as follows (see for example [START_REF] Piόro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF]):

∑ ( ∑ ∈ ) ∈ ( ∑ ∈ ∈ ( ∑ ∈ ∑ ∈ ∈ ∈ ( ∑ ∈ ∈ ∈ ∈ ( ( Let (
denote the problem defined by (2) for a given set of admissible paths ⋃ ∈ and let ( denote the (dual) polyhedron defined by constraints (2b)-(2e). Essentially, path generation is related to the so called dual separation problem. The problem consists in trying to find, for each demand ∈ , a path ( between ( and ( such that adding the constraints (2c)-(2d) for ( to formulation (2) makes the current optimal solution of ( infeasible for ( { ( } (see for example [START_REF] Piόro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF]). If there is no demand for which a path exists, then the current set is sufficient to find the optimum of FAP admitting all possible admissible paths. The problem of finding paths ( ∈ is called the pricing problem.

B. Pricing problem

Suppose form an optimal solution of the dual problem ( defined by [START_REF] Piόro | An optimization model for communication networks with partial multiple link failures[END_REF]. Then, to be sure that and form an optimal solution also for the dual problem with all possible paths, we need to check whether for each demand ∈ , and for each path in graph ( between the nodes ( and ( , there exist nonnegative numbers ∈ such that:

∑ ∈ | | ( | | ∈ ( ∈ ( Above, | | ∑ ∈
and | | ∑ ∈ ∈ denote the state-dependent dual length of path (note that to simplify the notation we do not use the link-path incidence coefficients here). Certainly, for all paths ∈ the above inequalities fulfilled by are a part of the considered optimal solution of (2). So the question now is how to find a path q outside the assumed set , if any, for which inequalities (3) are infeasible with respect to variables ∈ . In fact, in path generation it is advantageous (to speed up the PG algorithm, see Subsection III-D) to find, for each ∈ , not only a path (if any) in ̂ that just separates the current dual solution , but rather a path ∈ ̂ for which the dual constraints (2c) and (2d) are maximally violated by the considered dual solution (see [START_REF] Lasdon | Optimization Theory for Large Systems[END_REF], [START_REF] Minoux | Mathematical Programming: Theory and Algorithms[END_REF], [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF]). This is especially true when finding such a path is not substantially more complex than finding an arbitrary path that violates the dual constraints. In our setting, the (negative) measure of violation of the dual constraints corresponding to path is equal to

‖ ‖ ( | | ∑ | | ∈ ) ( ∑ ∈ ) (
provided it is negative. The quantity ‖ ‖, is referred to as the generalized dual length of path q. Thus, the pricing problem for demand ∈ can be specified as:

∈ ̂

‖ ‖ (

The so-formulated pricing problem ( 6) is difficult because of the particular form of the dual length ‖ ‖. Nevertheless, the problem can be stated as a binary program using formulation [START_REF] Lasdon | Optimization Theory for Large Systems[END_REF]. In the formulation, binary variables ∈ , specify the path we are looking for: { ∈ }. Binary variables ∈ in turn, identify the set corresponding to the so defined path { ∈ }. Besides, ( and ( denote the sets of all links outgoing from, and all links incoming to, respectively, node ∈ , and is the originating node of the considered demand and is its terminating node. where is a "big M" constant. Constraints (7b) and (7c) assure that variables that are equal to 1 form a path from to . Constraints (7d) force each variable ∈ to be equal to 1 when the length, with respect to , of the path defined by variables is smaller than . If we for a while assume that also = 0 when the length, with respect to , of path is greater than or equal to , then it is clear that the objective function computes the value of the quantity

‖ ‖ (| | ∑ (| | ∈
(see definition ( 5)). The optimal solution of ( 7) defines the optimal path { ∈ } we are looking for, with L equal to its reduced cost.

C. Path Generation algorithm

The path generation algorithm for FAP is as follows: Algorithm 1. Flow-adjustment problem Step 1: Initialization Define initial feasible path-sets ∈ .

Step 2: Solving

Solve the dual problem ( given by (2) to obtain dual variables and .

Step 3: Update

For each ∈ solve the pricing problem [START_REF] Lasdon | Optimization Theory for Large Systems[END_REF]. If the optimal objective L is negative, then add the resulting path q to the path-sets . If for any demand no new path has been added then stop: the resulting path-sets are sufficient to solve FAP to optimality. Otherwise, go to Step 2. End In this section, we present the results of a computational study illustrating the performance of FAR. The undirected network instances used in our tests and listed in Table 1 are taken from SNDlib [START_REF] Orlowski | SNDlib 1.0survivable network design library[END_REF]. In the following subsections, we will compare the Global Rerouting strategy (GR), the FTR strategy, and the FAR strategy in terms of cost-effectiveness. Recall that in a non-nominal link availability state, GR is allowed to restore flows for all demands in the available capacity from scratch, and that FTR is equivalent to FAR with . For the experiments reported in this section, we assume a uniform availability ratio for all the affected links in a given availability state, i.e., for all ∈ and ∈ (where denotes the subset of links affected in state s), and ∈ ∈ . Also, we assume the demand satisfaction ratio , i.e., for all ∈ and ∈ . In the comparisons, (respectively and will denote the optimal value of the network cost for FAR (respectively, the network cost for FTR and GR) with the availability ratio fixed to . GapFA (respectively, GapGR) gives the relative gain of the cost indicated by (respectively,

) with respect to for a given availability ratio : and .

A. Minimum network cost for single partial link failures

We first consider the single partial link failures, as a reference to compare FAR with other rerouting strategies. We assume that each link can fail but one at a time, and when it fails its availability ratio is equal to . Hence, the availability states can be identified with the links. In Figures 2345, we present GapFA and GapGR as a function of for 19 selected cases, varying from 95% to 5%. As illustrated by Figures 2 to 5, the overall cost of the network obtained when allowing thickening is substantially lower than that obtained using flow-thinning only, and this remains true for all considered availability scenarios and all networks.

It turns out that network topology and demand density characteristics have a significant impact on the efficiency of the flow-adjustment routing strategy. Net_10 and Net_14 have a high number of demands compared to their number of nodes (Net_10 has a full traffic matrix and Net_14 has a full undirected traffic matrix). For these networks the link cost optimized for FAR is close to the cost optimized for GR already for relatively low values of parameter (1.75 to 2). In contrast, Net_11 and Net_17 have a limited number of demands compared to their number of nodes and the optimal network cost for FAR requires a high value of parameter to meet the results of GR. We finally notice that for a sparse network with a high number of demands (as Net_17), allowing path thickening has a limited impact on the network cost for .

B. Minimum network cost for single partial node failures

We now consider single partial node failures (each node can fail but one at a time). This kind of failures represents a local perturbation in a network affecting several links at the same time. In this case, the availability states correspond to the nodes, and the link capacities in a failure state are given by if e is incident to the failing node, and to otherwise. In Figures 6789, we present GapFA and GapGR as We notice that the graphs have similar shapes for all the considered networks. Indeed, the gain in the network cost using the flow-adjustment strategy increases until a certain threshold (which seems not to be dependent on the value of ), and next decreases to 0% for full node failures ( =5%.) This shows that the flow-thinning strategy remains efficient for node failures for the extreme cases as sufficiently low or high perturbation levels. However, the gap between these two strategies and the global rerouting strategy is significant and reaches 80% to 90%, depending on the network, for a low availability ratio =5%. (Recall that for , the network cost for FAR converges to that of GR.) Figure 6. Relative gap between FAR, FTR and GR strategies for Net_10. Observe that whereas sparse networks (as Net_14 and Net_17) display a decrease of the network cost from 8% to 10% with (compared to flow-thinning), dense networks (as Net_11) have a limited decrease of the cost (maximum 5%). In fully meshed networks (as Net_10), it is virtually impossible to improve the network cost-the gain is less than 1%. However, further investigations are required to analyze this effect.

C. Range of the perturbation

Finally, we study the range of the perturbation, i.e., the number of paths which have their flow adjusted when a failure occurs.

We have two routing paths update processes: thinning and thickening. show the average number of paths that will have their flow adjusted, for all failure states with a fixed capacity availability ratio set to 50% for single link failures (Figure 10) and for single node failures (Figure 11). In these figures, the percentage of paths which has to be thinned (resp. thickened) is represented in solid (resp. dashed) bars.

The figures show that FAR requires modifying flows on almost every path when a failure occurs. The percentage of "thickening operations" depends on the value of and exhibits a similar behavior for all networks. Figure 12 displays the number of distinct paths used by the other strategies: FTR and GR. We first notice that FTR requires less paths than FAR, and the gap increases with . Now we compare the number of paths used by FAR and GR. When is low, FAR requires less paths than GR. The difference decreases for higher values of , and for , for all networks, GR uses less paths than FAR to achieve an optimal solution. This is mainly due to the fact that FAR uses paths with very low value in the nominal state to be able to increase them in a failure state. In regard to scalability and management cost, the above observations suggest that FAR should be used with small values of . Another way to limit the number of paths per demand used by this scheme is to modify the formulation (1) by adding specific constraints restricting the number of routing paths per demand, (which would surely have an extra computational cost).

VI. FAILURE RECOVERY PROCESS

The main interest in the flow-thinning strategy stems from the simplicity of the path-flow handling process. Indeed, the reaction to a particular availability state basically consists in decreasing the flow on (some of) the perturbed paths. In the following, we assume that a signaling protocol sends a message from the end nodes of a link to the source nodes of the paths passing through this link when its capacity decreases and reaches a certain threshold. Hence, the time required to recover from a failure state is equal to the time required for the signaling messages to traverse the longest path from the end nodes of the perturbed links to the source node of the disturbed demands.

In the following, we study the FAR restoration process.

A. Simulation of the FAR flow adjustment process

The flow adjustment process in FAR is not as simple as in FTR. This process is composed of two simultaneous stages managed by the source and the destination of the traffic demands [START_REF]Resource ReSerVation Protocol (RSVP) --Version 1 Functional Specification[END_REF]. The source nodes of the demands will first decrease the flow of some concerned paths, in order to make room for enlarging the flows on some other paths. Next, the destination nodes can increase the flow of the latter paths. Finally, we do not ensure full synchronization as this can be time consuming; hence, some routing paths may need a certain time to become fully operational. This will depend on the number of paths to restore, and on the quantity of flow required to be added to these paths.

B. Simulation results

We study the evolution of the ratio of demands that is perturbed when a failure occurs, and the time required for complete restoration by doing a basic simulation. This simulation is called "basic" for 2 reasons: we first assume that the buffer size is not modified by the perturbation (which is not correct for TCP where a failure is managed using congestion-avoidance algorithm) and we also consider that the time to pass through any link is the same, equal to 1 unit.

Figure 13 reports the ratio of perturbed demands for a fixed link failure as a function of the time, for network Net_14 with

. With FAR there are about 20% of demands perturbed for this link failure situation. Then, the bandwidth of some paths is decreased to allow the bandwidth of some other paths to be increased. This leads to further demands to be temporarily perturbed with about 55% of perturbed demands after 3 units of time. Nevertheless, the process converges very fast comparing to GR which reroutes all demands, perturbed or not. Finally, we notice that the ratio of perturbed demands remains too high compared to that with the flow-thinning strategy.

VII. CONCLUSION AND FUTURE WORK

In this paper, we have studied flow adjustment routinga routing strategy for demand protection which allows for limited path-flow thickening on top of a commonly understood flow thinning mechanism. FAR is capable of handling multiple partial link failures that are observed in wireless networks or in upper layers of wired communication networks. FAR requires less capacity in the network as compared with a pure Flow-Thinning Strategy, and is cost-wise very close to Global Rerouting for some networks. Our computational results show that FAR achieves encouraging results in terms of the investment cost, at least for the sets of availability states consisting of partial single link failures or single node failures.

We have noticed that the FAR restoration process affects a large number of path flows in the network and may lead to additional perturbation at the very beginning of the process if thinning and thickening operations are not coordinated. Hence, some non-perturbed demands may be affected for a short time when some bandwidth needs to be released from its paths. Making a distinction between paths of perturbed and non-perturbed demands could help to avoid this. On one hand, paths of non-perturbed demands should not be allowed to decrease their total flow more than their traffic requirement. On the other hand, paths of perturbed demands should be able to thin or thick their bandwidth. This leads to a restricted flow-adjustment routing strategy which is nothing else than an extension of the elastic robust rerouting strategy [START_REF] Fouquet | Elastic routing: a distributed variant, implementation issues, and numerical results[END_REF], with some constraints relaxed, applied to multiple partial link failures. Moreover, in its current form, FAR and FTR to a less degree, could require significant implementation effort as a large volume of paths are concerned and an important quantity of information is required to maintain all the connections at the routing nodes; this may impact the scalability of the strategy. These issues need further investigation and will be the subject of future work.

Figure 1 .

 1 Figure 1. Path generation algorithm for FAP V. NUMERICAL RESULTS

Figure 2 .

 2 Figure 2. Relative gap between FAR, FTR and GR strategies for Net_10.

Figure 3 .

 3 Figure 3. Relative gap between FAR, FTR and GR strategies for Net_11.

Figure 4 .

 4 Figure 4. Relative gap between FAR, FTR and GR strategies for Net_14.

Figure 5 .

 5 Figure 5. Relative gap between FAR, FTR and GR strategies for Net_17.

  for 19 selected cases, varying from 95% to 5%.

Figure 7 .

 7 Figure 7. Relative gap between FAR, FTR and GR strategies for Net_11.

Figure 8 .

 8 Figure 8. Relative gap between FAR, FTR and GR strategies for Net_14.

Figure 9 .

 9 Figure 9. Relative gap between FAR, FTR and GR strategies for Net_17.

Figure 10 .

 10 Figure 10. Fraction of paths to be thinned (solid bar) and thickened (dashed bar) with respect to the total number of paths per link failures.

Figure 11 .

 11 Figure 11. Fraction of paths to be thinned (solid bar) and thickened (dashed bar) with respect to the total number of paths per node failures.

Figure 12 .

 12 Figure 12. Average number of distincts paths used by FAR, FTR and GR per single link failures.

Figure 13 .

 13 Figure 13. Evolution of the ratio of perturbed demands for Net_14 (α=0,5).

  τ =1.1) FA (τ =1.5) FA (τ =2) FA (τ =6)

  The thinned nominal flows for state ∈ are denoted by ∈ ∈ . These flows are allocated to the admissible paths for the duration of the failure state.

							The
	state-dependent	flows	∈	∈	∈	are
	optimization variables. The relation, for each demand ∈
	, is defined by a given demand-dependent flow control
	factor , where	. Note that when	= 1, then, the
	reference path-flows of demand can only be thinned. Still,
	when	, they can be as well thickened, but only up to
	the factor so that,		∈	∈	∈ .
	B. Formulation of FAP
	FAP assumes that in each failure state ∈ , only a part
	of the nominal volume	has to be realized for each
	demand ∈ . This is achieved by thinning or thickening
	the nominal flows	to values	( ∈	so that the
	links capacities	∈ available in state are not
	exceeded. For given sets of admissible paths	∈ ,
	problem FAP can be represented by the following path-flow
	linear programming (LP) formulation involving variables
	( ∈	∈	,	( ∈	∈	∈	and
	( ∈	:

TABLE I

 I

	.	NETWORK DESCRIPTION
	Network	Nodes	Links Demands
	Net_10 (dfn-bwin)	10	45	90
	Net_11 (di-yuan)	11	42	22
	Net_14 (nobel-us)	14	21	91
	Net_17 (nobel-germany)	17	26	121

ACKNOWLEDGMENT

The work of French authors was carried out in the framework of the Labex MS2T, which was funded by the French Government, through the program "Investments for the future" managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02). M. Pioro was supported by National Science Centre (Poland) under grant 2011/01/B/ST7/02967 "Integer programming models for joint optimization of link capacity assignment, transmission scheduling, and routing in fair multicommodity flow networks."


Reliability analysis with ill-known probabilities and dependencies

In this work, we are interested in estimating the uncertainty that a system will work, given the probabilities that its components will be either failing or working. Usually, this estimation is done by propagating the component uncertainties through a structure function. This problem is well-studied when component probabilities are known and when the component are assumed to be independent. In this paper, we first recall some results obtained when considering that probabilities are unknown and when one either assume component independence or no knowledge about the components dependence. We also explore some structures not studied up to now, such as consecutive k-out-of-n structures.

INTRODUCTION

In many cases, a system reliability, here assumed to be binary (working or failing) can be estimated from the reliabilities of its components and from a structure function linking the binary states (working or failing) of those components to the system state. Usually, this estimation is done by assuming:

• that probabilities of component failures are precisely known or given;

• that the components behavior are independent.

The first assumption (preciseness) is quite strong, especially when few or no data concerning the components are available, or when one has to model some expert opinion about the component. In these cases, using a precise probability usually means one has to add some assumption not supported by available evidence (e.g., using maximum entropy principle [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF]). An alternative is to include the imprecision in the uncertainty model, by considering probability bounds.

The second assumption is in general more likely to hold, however the analyst may wish to relax this assumption, considering the case where the possible dependencies between components are unknown or only partially known.

Both issues have been investigated, in general settings, by imprecise probability theories [START_REF] Walley | Statistical reasoning with imprecise probabilities[END_REF][START_REF] Couso | A survey of concepts of independence for imprecise probabilities[END_REF]. However, the specific problem of assessing a system reliability under such conservative assumptions has only been explored in a very few works [START_REF] Utkin | Reliability models of m-out-of-n systems under incomplete information[END_REF][START_REF] Berleant | Bounding the times to failure of 2-component systems[END_REF][START_REF] Pedroni | Uncertainty analysis in fault tree models with dependent basic events[END_REF][START_REF] Fetz | Probability bounds for series systems with variables constrained by sets of probability measures[END_REF], and the case of partially specified independence in even less [START_REF] Hill | Inequality-based reliability estimates for complex systems[END_REF][START_REF] Troffaes | Applying the imprecise dirichlet model in cases with partial observations and dependencies in failure data[END_REF]. In this paper, we recall some of the main results of these previous papers, setting them in a general framework. We also provide some preliminary results about consecutive k-out-of-n systems [START_REF] Lambiris | Exact reliability formulas for linear & circular consecutive-k-out-of-n: F systems[END_REF], that have not been studied yet within an imprecise probabilistic framework.

PRELIMINARIES AND NOTATIONS

We consider that we have a set of components X 1 , . . . , X N , whose values are binary and described by domain X = {1, 0} where 1 usually stands for working and 0 for not working. We will denote by X N = × N i=1 X the Cartesian product describing the set of all possible system states. We will denote by x = (x 1 , . . . , x N ) ∈ X N a state of the system.

The uncertainty about whether a given component X i is working is given by an imprecisely specified probability, that is by two bounds

p i = p(X i = 1 i ) and p i = p(X i = 1 i )
Note that precise probabilities are retrieved when p i = p i .

Example 1. The assessment "component X i has a probability of working that is between 0.8 and 0.9" corresponds to p i = 0.8, p i = 0.9 For a given system, the structure function φ : X N → {0, 1} maps each system state x ∈ X N to 1 if the system works in this state, and 0 if the system fails in this state. φ -1 (0) and φ -1 (1) ⊆ X N will respectively denote the set of states for which the system fails and the set of states for which it works (note that they form a partition of X N ). We also assume that the system is coherent, that is if x ≥ x ′ (element-wise), then φ(x) ≥ φ(x), meaning that if more components work, then the system can only improve.

GENERIC PROBLEM

The generic problem we aim to tackle is to estimate the uncertainty bounds of φ -1 (1), that is p(φ -1 (1)) = 1-p(φ -1 (0)) and p(φ -1 (1)) = 1-p(φ -1 (0)), given our knowledge about the component uncertainties. Provided the bounds p i , p i are all we know, the problem of estimating p(φ -1 (1)) can be expressed as:

min p x∈X N ,φ(x)=1 p(x) (1)
under the constraints

p i ≤ x∈X N ,x i =1 p(x) ≤ p i , ∀i ∈ [1, N ] (2) x∈X 1:N p(x) = 1, p(x) ≥ 0 ∀x ∈ X N .
The upper bound is obtained by taking the max instead of the min. The above linear program expresses that we search for the joint probability that minimizes p(φ -1 (1)) that respects the marginals p i , p i , without making any assumption about the dependence. It is interesting to note that under such assumptions, the obtained inference on p(φ -1 (1)) will be imprecise, even if one starts from precise marginal probabilities p i = p i .

Although the above problem is a linear program, the number of variables increases exponentially with the number of components, making it NP-hard to solve (as is the fundamental reliability problem (Myers 2010, Ch 3.)) for generic systems. It is therefore interesting to seek specific cases displaying a lower computational complexity.

In the sequel, we will detail some results for such specific system structures, both for the case of independent components and unknown independence. We will stick to the easier but often considered case of components with identical uncertainty (i.e., p i = p w and p i = p w for all i). Also note, that, as the system is coherent, the minimum in ( 1) is obtained by considering p i = p i for every i, that is we can replace (2) by

p i = x∈X 1:N ,x i =1 p(x)
for every i. This makes the problem slightly simpler to solve, and basically says that we can consider it as a precise case problem where component precise probabilities of working are equal to their lower bounds. The same is true for the upper bounds, this time taking the upper probabilities p i as precise values.

THE CASE OF INDEPENDENT COMPONENTS

In the case of independence, we have that

p(x) = i,x i =1 i p i i,x i =0 i (1 -p i ), (3)
therefore in this case obtaining upper and lower bounds p(φ -1 (1)), p(φ -1 (1)) simply consists in replacing the probability that a component will be working by the appropriate bound in Eq. ( 3), for instance take p i = p i to compute p(φ -1 (1)).

Even in this case, computing p(φ -1 (1)), p(φ -1 (1)) in general remains hard. However, there are special configurations of systems for which easy-to-use formulas exist. In this paper, we will be interested in the following types of structures:

• k-out-of-n:F (k/n : F ) systems where a system of n components fails if k of them fails.; • Consecutive linear k-out-of-n:F (L : k/n : F ) systems, where components are linear and the system fails if k consecutive elements in the line fail. Such

k/n : F p(φ -1 (1)) = k-1 i=0 n i (p w ) n-i (1 -p w ) i p(φ -1 (1)) = k-1 i=0 n i (p w ) n-i (1 -p w ) i L : k/n : F p(φ -1 (1)) = n i=0 n-i•k i (-1) i (p w (1 -p w ) k ) i -(1 -p w ) k k-1 i=0 n-k(i+1) i (-1) i (p w (1 -p w ) k ) i p(φ -1 (1)) = n i=0 n-i•k i (-1) i (p w (1 -p w ) k ) i -(1 -p w ) k k-1 i=0 n-k(i+1) i (-1) i (p w (1 -p w ) k ) i C : k/n : F p(φ -1 (1)) = n i=0 n-i•k i (-1) i (p w (1 -p w ) k ) i +k k-1 i=0 n-k(i+1)-1 i (-1) i (p w (1 -p w ) k ) i+1 -(1 -p w ) n p(φ -1 (1)) = n i=0 n-i•k i (-1) i (p w (1 -p w ) k ) i +k k-1 i=0 n-k(i+1)-1 i (-1) i (p w (1 -p w ) k ) i+1 -(1 -p w ) n
Table 1. Reliability bound formulas in the independent case systems can model, for instance, communication relays [START_REF] Chiang | Relayed communication via consecutive-kout-of-n: F system[END_REF] or pumping stations; • Consecutive circular k-out-of-n:F (C : k/n : F ) systems, where components are placed around a circle and the system fails if k consecutive elements in the circle fail. Such systems can model, for instance, particle accelerators [START_REF] Kao | Computing reliability from warranty[END_REF].

All these systems generalize well-known series and parallel systems, as these latter systems correspond respectively to 1-out-of-n:F and n-out-of-n:F systems (whether they are consecutive or not). It should also be noted that, once n, k and the component uncertainty are fixed, a L : k/n : F system is more reliable than a C : k/n : F system, itself more reliable than a k/n : F system (in the sense that p(φ -1 (1)) increases and is the highest for L : k/n : F system). It means that formulas for k/n : F systems can be used as pessimistic bounds for the other systems. k/n : F , L : k/n : F and C : k/n : F systems under an assumption of component independence have been well studied in the literature [START_REF] Lambiris | Exact reliability formulas for linear & circular consecutive-k-out-of-n: F systems[END_REF][START_REF] Utkin | Reliability models of m-out-of-n systems under incomplete information[END_REF][START_REF] Chao | Survey of reliability studies of consecutivek-out-of-n: F and related systems[END_REF][START_REF] Eryilmaz | Review of recent advances in reliability of consecutive k-out-of-n and related systems[END_REF], hence obtaining lower and upper bounds case simply consists in replacing precise probabilities by appropriate bounds. Table 1 summarizes the formulas for the different systems. Note that, although trivial, the formulas for L : k/n : F and C : k/n : F systems considering imprecise probabilities have not been investigated before.

Example 2. Consider components such that [p w , p w ] = [0.95, 0.99] and 2/4 systems. Then, using the formulas of Table 1, we obtain the following bounds:

2/4 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.9859, 0.9993]; C : 2/4 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.9905, 0.9996]. L : 2/4 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.9927, 0.9997];

THE CASE OF UNKNOWN INDEPENDENCE

In this section, we interest ourselves to the case of unknown independence, that is the case where we have to solve the optimization program given by ( 1) and (2). For the k/n : F systems, we can rely on previous results [START_REF] Utkin | Reliability models of m-out-of-n systems under incomplete information[END_REF] indicating that

p(φ -1 (1)) = max(0, (n -k + 1)p w + k -n); p(φ -1 (1)) = min(1, kp w ).
In the case of series and parallel systems, we retrieve the well-known Frechet bounds [START_REF] Pedroni | Uncertainty analysis in fault tree models with dependent basic events[END_REF][START_REF] Fréchet | Généralisation du théoreme des probabilités totales[END_REF]. The cases of L : k/n : F and C : k/n : F systems have not been investigated up to now, but obtaining bounds under an assumption of unknown independence for such systems is harder than for the assumption of independence. However, the two following propositions show that, for L : 2/3 : F systems and C : 2/4 : F systems (C : 2/3 : F systems are equivalent to 2/3 : F systems), closed formulas can be obtained.

Proposition 1. Given component uncertainty p w , p w and unknown independence, the lower and upper bounds p(φ -1 (1)), p(φ -1 (1)) for a L : 2/3 : F system are

p(φ -1 (1)) = p w p(φ -1 (1)) = min(1, 2p w )
Example 3. Consider components such that [p w , p w ] = [0.95, 0.99] and 2/3 : F , L : 2/3 : F systems with an assumption of unknown independence. Then, using the corresponding formulas for these systems we obtain : 2/3 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.9, 1];

L : 2/3 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.95, 1].
The lower bound of the L : 2/3 : F being indeed slightly higher than the bound of the 2/3 : F system. Proposition 2. Given component uncertainty p w , p w and unknown independence, the lower and upper bounds p(φ -1 (1)), p(φ -1 (1)) for a C : 2/4 : F system are p(φ -1 (1)) = max(0, 2p w -1)

p(φ -1 (1)) = min(1, 2p w )
Example 4. Consider components such that [p w , p w ] = [0.95, 0.99] and 2/4 : F , C : 2/4 : F systems with an assumption of unknown independence. Then, using the corresponding formulas for these systems we obtain :

2/4 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.85, 1]; C : 2/4 : F [p(φ -1 (1)), p(φ -1 (1))] = [0.9, 1].
The lower bound of the C : 2/4 : F being indeed slightly higher than the bound of the 2/4 : F system.

In comparison with the formulas given in Table 1 for the independent case, the above formulas look extremely simple and are straightforward to use. This suggests that the case of unknown independence for consecutive systems may actually be simpler to handle. We also think that Proposition 2 will be more easy to extend to generic k/n systems then Proposition 1, due to the symmetry of circular systems that is not present in linear systems, where components at the beginning and end of the line play particular roles. It would also be interesting to relax the assumption that our knowledge about each component is identical (i.e., p i = p w for all i ∈ [1, N ]).

Related to this last remark, it should also be noted that bounds for k/n : F and for C : k/n : F could be used to provide approximations of the L : k/n : F systems.

CONCLUSION

We have recalled some results regarding the evaluation of lower and upper probabilities that a given system will work. In particular, we have settled them as a generic optimization problem, and have mentioned some structures for which this evaluation can be done through a closed-form formula.

Additionally, we have started to investigate the case of consecutive k-out-of-n:F systems, which to our knowledge have never been studied in the framework of imprecise probabilities. Our results are still preliminary for the case of unknown independence, but suggest that relatively easy formulas can be obtained for such systems as well. As such systems are used in a lot of modern applications, finding such formulas would be interesting, at least to have an idea of the extreme bounds one can obtain for such systems.

In future works, we intend to study how to integrate some known dependency information [START_REF] Hunter | An upper bound for the probability of a union[END_REF]) in the constrained problem (2) so that it can still be easily solved, and in particular for which case computations can still be achieved through closed form formulas.

Another interesting future line of work is to study other aspects of consecutive kout-of-n systems when probabilities or dependencies are ill-known, such as for instance importance measures [START_REF] Zhu | Patterns of the birnbaum importance in linear consecutive-k-out-of-n systems[END_REF], multi-state versions of the systems [START_REF] Radwan | Bounds for increasing multi-state consecutive k-out-of-r-from-n: F system with equal components probabilities[END_REF] or design optimization [START_REF] Shingyoch | Proposal of simulated annealing algorithms for optimal arrangement in a circular consecutivek-out-of-n: F system[END_REF]).

ACKNOWLEDGMENT

This work was carried out and funded in the framework of the Labex MS2T. It was supported by the French Government, through the program "Investments for the future" managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02).


Hugo Bazille

Jacques Nicolas

Computational Protein Design: trying an Answer Set Programming approach to solve the problem

Proteins are macromolecules made of a chain of amino-acids. The combinatorial nature of the space of possible protein conformations makes computer-aided protein study a major research field in bioinformatics. The problem of computational protein design aims at finding the best protein conformation to perform a given task. This problem can be reduced to an optimization problem, looking for the minimum of an energy function depending on the amino-acid interactions in the protein.

We have designed a model based on Answer Set Programming. The CPD problem may be easily modeled as an ASP program but a practical implementation able to work on real-sized instances has never been published. We have raised the main source of difficulty for current ASP solvers and ran a series of benchmarks highlighting the importance of finding a good upper bound estimation of the target minimum energy to reduce the amount of combinatorial search. Our solution clearly outperforms a direct ASP implementation without this estimation and has comparable performances with respect to SAT-based approaches. It remains less efficient that the recent approach by cost function networks of D. Allouche & al., showing there exists still some place for improving the optimization component in ASP with more dynamical strategies.

Introduction

Proteins are essential compounds of living organisms, implied in almost all structural, catalytic, sensory, and regulatory functions. Proteins are amino-acid sequences with many different functions, mostly determined by their three dimensional structure. The study of these structures is thus an important field in biology with applications in various fields such as medicine, biotechnology, synthetic biology... [START_REF] Pleiss | Protein design in metabolic engineering and synthetic biology[END_REF]. Computer-assisted study of proteins offers opportunities to mimic the evolution and create new mutations or new structures in proteins [START_REF] Peisajovich | Protein engineers turned evolutionists[END_REF].

The goal of computational protein design (CPD) is precisely to find among a collection of proteins those that most likely target a function. It is sometimes referred as the inverse folding problem: whereas protein folding is looking for the 3D structure of a given sequence, protein design searches the amino-acid sequences that would fold into a given 3D structure. As they are 20 possible amino-acids for each position in a protein sequence, each of them accepting several structural variants, the number of combinations to be tested is out of range of any experimental approach, even for short sequences. Consequently, it became a strategic simulation challenge and numerous works have demonstrated the power of computer-aided protein design [START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Gainza | Protein design using continuous rotamers[END_REF][START_REF] Hallen | Dead-end elimination with perturbations (deeper): A provable protein design algorithm with continuous sidechain and backbone flexibility[END_REF][START_REF] Hu | A gradient-directed monte carlo approach for protein design[END_REF][START_REF] Gg | Improved prediction of protein side-chain conformations with scwrl4[END_REF][START_REF] Ollikainen | Sat-based protein design[END_REF]. Among the most striking results of this approach, one can mention the production of antimalarial drugs from the engineered bacteria E. coli [START_REF] Vincent Jj Martin | Engineering a mevalonate pathway in escherichia coli for production of terpenoids[END_REF] or the development of the most efficient computationally designed enzyme for the Kemp elimination to date [START_REF] Hk | Iterative approach to computational enzyme design[END_REF].

In CPD, choosing the best amino-acid sequences to perform some function is formulated as an optimization problem. Even with the most recent advances in CPD, techniques still need to be improved: many approximations are made in order to make this problem solvable and more realistic models are needed. It is also necessary to generate a limited collection of solutions close to the optimum since, due to model approximations, the optimum of the combinatorial problem is not always the best protein in practice.

This paper presents a specific approach to solve the CPD problem, Answer Set Programming (ASP) [START_REF] Lifschitz | What is answer set programming[END_REF]. This work is based on the previous results of [START_REF] Filipe | Protein design using answer set programming[END_REF], with a number of enrichments. In order to evaluate this work, two comparison studies were carried out: the first one quantifies the progress with respect to [START_REF] Filipe | Protein design using answer set programming[END_REF], and the second one is a more demanding task addressing the most advanced results in the domain [START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Allouche | Computational protein design as a cost function network optimization problem[END_REF][START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF]. The ASP model performances do not reach the level of the best current approaches and a first discussion on the advantages and drawbacks of the ASP approach for solving the CPD problem is provided.

We start by giving the biological background of this work and defining the CPD problem with a fixed rigid backbone and discrete sets of possible rotamers in proteins. Then, we give an overview of state of the art, with various techniques and paradigms. Next, we recall basic ASP notions and give our encodings of the CPD problem. Results are described in a last section, together with the analysis of the different encodings. We conclude by different perspectives that seem interesting to be studied in a future work.

The computational protein design problem (CPD)

The protein design problem is an optimization problem on the conformation (geometrical structure) of a protein whose structure is partially known.

Elements of protein structure

Proteins are biological macromolecules made of several amino acids linked together by peptide bonds. There are mainly 20 different amino acids with a common organization illustrated in figure 1. The 'R' represents the side-chain that makes the amino acid unique. An amino acid in a polypeptide is also called a residue. The sequence of carbon, oxygen and nitrogen atoms without the side chain defines a 3D structure called backbone. The geometry of the protein is relatively rigid and generally defined by three dihedral angles (ϕ, ψ and ω) for each amino acid along the backbone. A dihedral angle is the angle between two planes defined by a sequence of four atoms. Furthermore, each side-chain has up to four degrees of freedom: the dihedral angles χ (see an example figure 2).

The backbone structure is highly constrained (the 6 atoms implied in the peptide bond C α C(O)N HC α form a plane). This paper assumes that its structure has no degree of flexibility. In contrast, side chains have different possible conformations called rotamers. Whenever there is a degree of freedom, there is an infinite number of possible rotamers, as the molecules can continuously rotate around the axis. However, it is often sufficient to consider a representative finite set of rotamers, which can be determined by a statistical analysis of the actual conformations. For each amino acid, there are from 1 to 30 of most common rotamers that are listed in public libraries. We have used the 2010 version [START_REF] Shapovalov | A smoothed backbonedependent rotamer library for proteins derived from adaptive kernel density estimates and regressions[END_REF] of the Dunbrack backbone-dependent rotamer library [START_REF] Dunbrack | Backbone-dependent rotamer library for proteins application to side-chain prediction[END_REF].

Protein design through energy minimization

The aim of CPD is to find a protein that will perform a desired function. This function depends on the backbone structure and on the rotamer configurations. For a given backbone structure, the protein tends to adopt a global stable configuration of minimal energy. Therefore, the practical goal of CPD is to find a sequence of residues and a conformation that folds into a defined backbone and minimizes the protein energy, which can be expressed as an optimization problem. The computation of realistic energy functions is a whole topic by itself and is not the subject of this paper. Energy values are part of the benchmark data sets we have used and they are considered here as input facts. The energy minimization depends only on the energy of interaction between rotamers or between the backbone and rotamers. Large scale protein design uses pairwise approximation: the energy function is supposed to be pairwise decomposable on residues. Moreover, some design positions are chosen to be mutated. With i r the rotamer at position i, E(i r ) the energy of interaction between rotamer i and the backbone, and E(i r , j r ) the energy of interaction between rotamers r at position i and r at position j, the formula to minimize is:

E = i E(i r ) + i j,j<i E(i r , j r ) 1 2 3
The energy functions E(i r ) and E(i r , j r ) are based on empirical measures and depend on many parameters such as van der Waals potentials, electrostatics,. . . Protein design is a NP-hard problem [START_REF] Pierce | Protein design is np-hard[END_REF] associated to a huge hypothesis space. Indeed, a protein with m residues and a mean number of n rotamers per amino acid gives rise to (20 × n) m possible conformations. Using reasonable parameters such as m = 100 and n = 10 leads to ≈ 10 234 possible solutions. Variants of this problem are also studied, such as protein design with continuous rotamers or with backbone flexibility [START_REF] Hallen | Dead-end elimination with perturbations (deeper): A provable protein design algorithm with continuous sidechain and backbone flexibility[END_REF].

State of the art for CPD

In this section, we introduce a few state of the art algorithms to solve this problem. A first part introduces algorithms that are allowed to find approximations of the best solution and are supposed to scale to larger proteins. A second part introduces algorithms that are designed to find exactly the best solution.

CPD with approximations

Genetic algorithm based approaches have been tried on the CPD problem [START_REF] Luis | Using genetic algorithm to design protein sequence[END_REF]. Each element of the population is a sequence of amino-acids with their rotamer conformations. A cross-over operator is defined by exchanging amino-acids at a position in two elements of the population. A mutation operator is defined by introducing a new amino-acid at some position. The algorithm discards the conformations of higher energy and a new population is generated by applying operators on the remaining rotamers. This knowledge-poor approach was applied early on CPD [START_REF] David | De novo protein design using pairwise potentials and a genetic algorithm[END_REF] but other stochastic methods are now preferred for this task.

Rosetta [START_REF] Leaver-Fay | Chapter 19 -rosetta3: An object-oriented software suite for the simulation and design of macromolecules[END_REF] is a large open source package that is the most representative of this tendency. It uses the Dunbrack backbone-dependant library to reduce the set of conformations and a Monte-Carlo based method to sample this space [START_REF] Hu | A gradient-directed monte carlo approach for protein design[END_REF][START_REF] Kaufmann | Practically useful: What the rosetta protein modeling suite can do for you[END_REF]. It takes into account several factors to build a more precise energy function: van der Waals potentials, electrostatics... Even if is not used as a solver, Rosetta offers many functionalities and was useful in our first series of experiments to compute the needed interaction energy values.

SCWRL4 [START_REF] Gg | Improved prediction of protein side-chain conformations with scwrl4[END_REF] is a program used to solve a restriction of the protein design problem: side-chain positioning. In this problem, amino acids are fixed and only the side-chain conformation is varying. This assumption greatly reduces the search space. An interaction graph that represents the side-chain placement problem is first created. The graph is decomposed into trees on which a branchand-bound search can be performed. A rotamer collision detection step reduces the risk to produce physically impossible solutions. Most of the time, tree search algorithms are looking for exact solutions.

Exact search algorithms

As the search space of the CPD problem can be seen as a decision tree, some of the early methods to explore it have been based on graph traversal algo-rithms such as A* [START_REF] Andrew R Leach | Exploring the conformational space of protein side chains using dead-end elimination and the A* algorithm[END_REF]. The A* algorithm needs a heuristic function to perform the search. This function is calculated by adding in a first part the energy for all nodes already assigned and in a second part a suboptimal function over unassigned nodes. For each of node i r , the suboptimal function is E(i r ) + js assigned E(i r , j s )+ k unassigned min r E(i r , k r ). The interested reader may consult [START_REF] Hallen | Dead-end elimination with perturbations (deeper): A provable protein design algorithm with continuous sidechain and backbone flexibility[END_REF] to look at the current results achieved by such methods.

The problem of minimizing the total free energy of a protein conformation can also be represented as an Integer Linear Program (ILP). Indeed, binary variables q i (r i ) can represent the presence of a rotamer at place i. The presence of exactly one rotamer at each position is expressed by the formula ri q i (r i ) = 1. The energy minimization is then represented by a minimization term

min i ri q i (r i ) E(r i ) + j =i rj E i,j (r i , r j ) q i (r i ) q j (r j ).
This approach is extended in [START_REF] Zhu | Mixed-integer linear programming algorithm for a computational protein design problem[END_REF] where only some variables have integer domains (Mixed ILP Programming). CPD can also be formulated as a satisfiability problem with minimization of a formula. For each rotamer and each position, a boolean variable indicates if it is chosen or not. This approach enables to use SAT solvers, and effective encodings with improved branch and bounds algorithms are presented in [START_REF] Ollikainen | Sat-based protein design[END_REF]. Basically, given a position i and a set of possible rotamers i 1 . . . i r the problem of choosing one rotamer can be encoded by formula

Φ i = (i 1 ∨ . . . ∨ i r ) ∧ s =s ¬(i s ∧ i s ).
Then, a global formula can be established: Φ = i Φ i . Clauses of the form ¬(i s ∧ j s ) derive from the computation of incompatible pairs of rotamers. In practice, hierarchical decisions are made: first the amino acid is assigned, then the dihedral angles (side chain). These programs are quite efficient, provided they take into account only three possibilities by dihedral angle.

The most recent advance in solving the CPD problem proposes to model it as a weighted constraint satisfaction problem [START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Allouche | Computational protein design as a cost function network optimization problem[END_REF][START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF], using cost function networks. A Cost Function Network (CFN) is a pair (X, W ) where X is a set of n variables and W is a set of cost functions. Each variable i ∈ X has a finite domain D i of values. A value r ∈ D i is denoted i r . For a set of variables S ⊆ X, D S denotes the cartesian product of the domains of the variables in S. A cost function w S ∈ W , with scope S ⊆ X, is a function w S :

D S → [0, k]
where k is a maximum integer cost used for forbidden assignments. There is a natural correspondence between the formalism of CFN and the elements of CPD: each variable i corresponds to a position to be redesigned and each value r in domain D i corresponds to a possible rotamer r at position i. Furthermore, each unary cost function w(i r ) represents the energy of interaction rotamer-backbone, and each binary cost function w(i r , j r ) represents the energy of interaction between two rotamers. This approach has led to major improvements and it dramatically outperforms previous works for exact results. It allows to solve big instances with tens of positions and thousands of rotamers.

We have studied yet another framework, Answer Set Programming (ASP), which has been tried only once to our knowledge, in a Master thesis [START_REF] Filipe | Protein design using answer set programming[END_REF], where a direct translation of the problem specification has been tested, together with several variants such a as a separated treatment of positive and negative energy score values. The ASP solvers appeared to be able to solve only a very low number of residues. Since ASP has exhibited interesting performances in a number of combinatorial optimization problems, we have tried to understand the reasons of this relative failure. An earlier paper points to the difficulty of solving the related problem of protein folding on a lattice [START_REF] Dovier | A comparison of clp (fd) and asp solutions to np-complete problems[END_REF] but the ASP solvers have made progress since then and the interest of CPD is to offer a combinatorial problem where optimization plays a major role. We have emphasized the specific features of the problem that make it difficult and looked at better encodings, achieving clear progress with respect to the way its resolution could be improved.

4 Modeling CPD using Answer Set Programming (ASP)

A bit of syntax and the solver clingo4

Answer Set Programming (ASP) is a form of declarative logic programming using the the stable models semantics. It is a Boolean constraint solving framework designed for Knowledge Representation and reasoning and for complex optimization problems: solvers for ASP programs implement heuristics in order to solve efficiently minimization queries over weighted atoms.

This document follows the current standard of ASP language1 . Given a set of atoms A, a normal logic program over A is a finite set of normal rules (a 0 is called the head and the rest is the body):

a 0 : -a 1 ; . . . ; a m ; not a m+1 ; . . . ; not a n , where a i ∈ A for all i, . It is possible to use predicates of any arity and first order variables in ASP, provided that these variables are defined over a finite domain. Each variable will be replaced during a grounding phase by its possible values in the Herbrand universe of the program, leading to a fully instantiated program.

It is possible to write integrity constraints, which are rules with no head: : -a 1 ; . . . a m ; not a m+1 ; . . . not a n .

A cardinality rule is of the form a 0 : -l{a 1 ; . . . ; a m ; not a m+1 ; . . . ; not a n }u It allows to control the cardinality of sets of atoms, with lower bound l and upper bound u. Similarly, a choice rule is of the form: l{a 1 ; . . . ; a m }u : -a m+1 ; . . . ; a n ; not a n+1 ; . . . ; not a p , meaning that if the body holds, at least l and at most u atoms of the head have to hold. Such rules are expanded in a number of normal rules quadratic in m.

ASP systems include statements to express cost functions [START_REF] Gebser | On the implementation of weight constraint rules in conflict-driven asp solvers[END_REF] and multicriteria optimization. With w i the weight of literal l i and p i its priority level, the minimize statement writes #minimize {w 1 @p i : l 1 . . . ; w n @p n : l n }. The solver minimizes the sum of the weighted literals, starting by the highest priority level.

The first ASP solvers were quite "monolithic": first the grounding generated a propositional program, then a solver computed its stable models. The solvers use enhancements of the Davis-Putman-Logemann-Loveland (DPLL) algorithm, analyzing and learning conflicts in case of failure to prune the search space. In the recent versions of the ASP system clingo4, two scripting languages have been integrated with ASP in a common environment in order to achieve complex reasoning processes: Python and Lua [START_REF] Gebser | Clingo= asp+ control: Preliminary report[END_REF]. On the declarative side, it becomes possible to define procedures allowing to instantiate different logic programs with different parameters (directive #program) and they may be grounded and solved at any time using built in procedures ground() and solve(). External rules ("i. e. "volatile rules ") may also be added/removed in a program (#external).

We have chosen clingo4 because it is one of the most efficient ASP solver to date and because it offers large possibilities to adjust the solving strategy. The user can easily design its program by incremental refinements or produce high level interactive solvers. A client-server ASP solving process encoded in Python is presented in [START_REF] Gebser | Clingo= asp+ control: Preliminary report[END_REF]. Scripting has been used in our work to combine the search for approximate solutions and local exploration. Since clingo4 is still evolving, this kind of study can help its designers to improve the optimization component.

Algorithms and ASP models to solve the CPD

DEE pruning algorithms Like in every combinatorial problem, the main issue in CPD is the very fast increase of the size of the search space with respect to the size of the protein. The goal of Dead End Elimination (DEE ) algorithms is to remove the choice of some residue at some place because it always gives worse results than other residues at this place. They can also be extended to eliminate pairs of residues, triplets... Numerous DEE algorithms are available that have different complexities and efficiencies We have used two of them, simple split and double Golstein, offering a good tradeoff with respect to these parameters.

Simple split eliminates a rotamer i r if ∃k, ∀v, ∃i t ,

E(i r ) -E(i t ) + (E(i r , k v ) -E(i t , k v )) + j =k =i min u [E(i r , j u ) -E(i t , j u )] > 0 (1)
If p is the number of positions and n the number of rotamers at some position, the complexity of the algorithm is O(p 2 n 3 ). The Double Goldstein DEE eliminates a pair of rotamers (i 1r 1 , i 2r 2 ) if there exists a pair (i 1t 1 , i 2t 2 ) such that

(E(i 1r 1 ) + E(i 2r 2 ) + E(i 1r 1 , i 2r 2 )) -(E(i 1t 1 ) + E(i 2t 2 ) + E(i 1t 1 , i 2t 2 )) + j =i1j =i2 min u (E(i r , j u ) -E(i t , j u )) > 0 (2)
The complexity of the algorithm is O(p 3 n 5 ), a noticeable increase of complexity that allows to eliminate a high number of pairs. In practice, this algorithm is launched after the first one to decrease the value of n.

Search for exact results

The basic ASP code for the CPD problem is introduced as in [START_REF] Filipe | Protein design using answer set programming[END_REF] in order to serve as a reference point for the improvements that are tested thereafter. The set of positions, of possible residues and possible rotamers at each position have been represented by the following facts : position(Position in the protein). possibleResidue(Position,Aminoacid). possibleRotamer(Position,Aminoacid,Rotamer,EnergyB).

EnergyB is the interaction energy with the backbone. Then the interaction energies between rotamers are tabulated from Rosetta or Osprey computations: interEnergy(Pos1,AA1,Rot1,Pos2,AA2,Rot2,EnergyR).

Given these facts as input data, the CPD problem can be expressed in 3 rules. For each position, one residue and one rotamer are chosen to be in the solution. It is recorded in predicate residue(Pos,AA) and rotamer(Pos,AA,Rotamer).

Algorithm 1 Exact encoding

% Exactly one residue must be assigned to each position (a) 1{residue(Pos,Amin) : possibleResidue(Pos,Amin)}1 :-position(Pos). % Exactly one rotamer must be assigned to each residue (b) 1{rotamer(Pos,Amin,Id) : possibleRotamer(Pos,Amin,Id,Energy)}1

:-residue(Pos,Amin), position(Pos).

% Minimize the sum of the energies (c) #minimize{ EnergyB@1,Pos,Id: rotamer(Pos,Amin,Id),possibleRotamer(Pos,Id,EnergyB) ; EnergyR@1,Pos1,Pos2: rotamer(Pos1,Amin1,Id1),rotamer(Pos2,Amin2,Id2), interEnergy(Pos1,Amin1,Id1,Pos2,Amin2,Id2,EnergyR)}.

-Rule (a) chooses a residue among the possible residues at each position; -Rule (b) chooses a rotamer among possible rotamers for each chosen residue.

-Rule (c) minimizes the sum of the recorded energies predicates. The weight of each predicate is the value of its Energy parameter.

As for the SAT specification of CPD, a file containing the list of pruned rotamer pairs computed by DEE algorithms before the search is added to this program. An eliminated pair is taken into account by an integrity constraint: :-rotamer(Pos1,Id1), rotamer(Pos2,Id2).

This short program computes exact results but is limited to very small instances as runtime may quickly become too long. The optimization step is crucial in CPD and we have studied the way it is managed by clingo to understand the limitations and improve the code.

Transformation in an equivalent problem with better properties

Different minimization algorithms are implemented in the ASP solver (branch and bound, unsatisfiable core...) but they all require positive values. This means that programs are automatically rewritten to fulfill this constraint in order to transform the minimization problem into an equivalent problem with only positive values and lower bound 0. This rewriting algorithm is quite simple: each time there is a minimization statement of the form #minimize { V:p }., V < 0, the solver transforms it in #minimize { -V: not p}. It appears that the transformation is far from being optimal in our problem. Let us consider a short example:

a;-4 b;0 c;-1 3 0 -1 1
In this schema, it is possible to choose rotamer a or b at position 1 and rotamer c or d at position 2. The minimal solution of this problem is the pair (a, d), a configuration with total energy -6 (-4+-2+0). Solving this problem with (an adaptation of) the naive encoding given in section 4.2 leads to the good answer, but the optimization value reached in the transformed system is 2. In details, weights are added to the energy of the solution depending on the following conditions: d) is chosen. On this example, the final optimization value is not too far from the lower bound 0. On real instances with many negative values it is no more the case in general and since the choice of a rotamer entails that all others are not chosen and it is difficult for the solver to find good lower bounds.

4 if a is not chosen, 0 if b is chosen, 1 if c is not chosen, 2 if d is not chosen 3 if the pair (a, c) is chosen, 0 if (a, d) is chosen, 1 if (b, c) is not chosen, 1 if (b,
To bypass this issue, the problem is transformed such that minimization is carried over positive values only and the solver can estimate a much better lower bound. This issue may be stated as a constraint satisfaction problem and we apply a technique imported from this field, constraint propagation, to reduce the set of possible values. This is achieved by looking for two local consistency properties, soft arc consistency, which is propagating energy constraints on pairs of rotamers, and soft node consistency, which is propagating energy constraints on a single rotamer with the backbone. The principle is used in the cost function network approach [START_REF] Cooper | Arc consistency for soft constraints[END_REF][START_REF] Larrosa | Solving weighted {CSP} by maintaining arc consistency[END_REF] where it appeared to be a major source of efficiency in solving the CPD problem. The idea is the following one: for each position i, let N i,1 , . . . , N i,r be the interaction energy of rotamers i 1 , . . . , i r with the backbone, and let N i = min j N i,j . Then, the energy N i may be considered as a fixed cost E ∅ associated to any choice in the position and it is only necessary to search for the minimization on the energies N i,1 -N i , . . . , N i,r -N i that are all positive. A similar process can be done for the energies of interaction between rotamers (see also the hierarchical approach below).

In the previous example, the minimum energy of interaction with the backbone is -4 at position 1 and -2 at position 2. Thus any choice over positions 1 and 2 includes a fixed cost of interaction with the backbone equal to -4+-2= -6, which is the best possible lower bound in this case, and the cost of a, b, c and d may be respectively replaced by 0, 4, 1 and 0 for the minimization.

In terms of the ASP program, this means that all predicates possibleRotamer(P os, AA, Rotamer, Energy). will be changed in possibleRotamer(P os, AA, Rotamer, Energy -M inEnergy at P os).

It is a bit more complex since the energies of interaction between rotamers have also to be transformed. For instance in the previous problem, the minimum energy of interaction between rotamers once b has been chosen is -1 and thus the interaction energy -1 may be added as a fixed cost for the choice of b.

Add a hierarchy to the encoding

In a branch and bound search, the sooner a cut can be made, the better it is. We propose to cluster rotamers with similar properties at some position in order to be able to do more cuts sooner. Fortunately, there is a natural way to achieve this clustering in constant time: regroup the rotamers by amino-acid since they share in general similar interaction energies.

Let i be a position of an amino-acid A with rotamers i A1 . . . i Ar , having energies of interaction with the backbone E(A 1 ) . . . E(A r ). If N = min E(A 1 ) . . . E(A r ).

Then energy E (A) = N is associated to the choice of A and energies E (A 1 ) = E(A 1 ) -N . . . E (A r ) = E(A r ) -N to the choice of i A1 . . . i Ar . It allows to replace r cuts at the rotamer level by a single cut at the amino-acid level. The same principle may be applied for pairs of rotamers. For amino-acid A at position i with conformations r and B at j with conformations r ,one defines (i A , j B ) = argmin r,r E(i r , j r ). Then energies E(i r , j r ) are replaced by energies

E (i A , j B ) and E (i r , j r ) = E(i r , j r ) -E(i A , j B ).
Algorithms 2 and 3 correspond to a sample of transformation steps described previously to get an equivalent program, adapted for the hierarchical encoding. Algorithm 4 gives an overview of the whole process. The lower bound is iteratively improved as much as possible for every position. Its complexity is O(p 2 n 3 ) (with p the number of positions and n the number of rotamers by position), the same complexity that most used DEE algorithms. More on arc consistency enforcing algorithms may be found in [START_REF] Bessière | An optimal coarse-grained arc consistency algorithm[END_REF].

Algorithm 2 Energy transfer from rotamers to a residue. E ← min E(Ar) for all rotamers Ar at position i do

E(iA r ) ← E(iA r ) -E end for E(iA) ← E

Enumeration of ε-solutions

As we work on an abstraction of a biological problem, some approximations are made and solutions of the optimization problem may not be solutions of the CPD problem. Thus solutions close to the optimum are interesting too. In order to produce a set of good solutions with a fixed divergence from the global minimum solution, we assume the minimum value E best has been produced and restart the algorithm with a modified DEE pruning including a parameter ε in the right Algorithm 3 Energy transfer from pairs to rotamers.

for all position j do E ← min E(iA r , jB r ) for all rotamers B r at position j do E(iA r , jB r ) ← E(iA r , jB r ) -E end for E(iA r ) ← E(iA r ) + E end for members of the equations in 4.2. It is also possible to adapt DEE algorithms to discard rotamer i r if a lower bound of the energy E lower (i r ) when i r is chosen is greater than E best + ε. The cost of obtaining lower bounds increases when the required precision increases. A possible lower bound is

E lower (i r ) = E(i r )+ j =i min r [E(j r )+E(i r , j r )]+ j =i k>j,k =i min r ,r" [E(k r" , j r )]
Computing this lower bound for all rotamers at all possible place has complexity O(n 3 p 3 ) with n the number of rotamers and p the number of positions. Similar criteria to eliminate pairs of rotamers have complexity O(n 5 p 4 ). After pruning the search space, the search itself can be launched. An ASP program doing this task may be constituted of the program Exact encoding with an additional rule stating that the sum of energies must not be greater than E best +ε. A 2% variation of E best may increase the number of solutions to several thousands.

Results

Benchmarks have been run on two different datasets: the first one is extracted from [START_REF] Filipe | Protein design using answer set programming[END_REF], which contained initially 120 instances. They are based on 12 instances of 10 proteins , considering 10, 12, 15 and 17 positions and three different aminoacid sets (fixed, only hydrophobic or all aminoacids). We retained the 40 most difficult problems. It serves to compare the progress made between [START_REF] Filipe | Protein design using answer set programming[END_REF] and our own work. The second dataset is extracted from [START_REF] Allouche | Computational protein design as an optimization problem[END_REF] and contains 47 instances. It allows us to compare our results to the last advances in the domain. The solver clingo4 has been run with default parameters. Trying dfferent solver options for clasp (particularly for the choice of heuristics) has shown no significant improvement. All runs have been performed on a Intel Xeon W3520 quad-core, 2.66 GHz. All reported times are in seconds.

DEE

The efficiency of DEE algorithms has been evaluated on the first dataset. Results are presented in table 1 (columns 4 to 7). For each instance, we perform 3 successive runs of Goldstein DEE (a fast algorithm) that allows to run then 5 times Simple Split (a complex algorithm) on a highly reduced dat set. Note that more runs can be necessary to reach a fix point but it would not further improve the total search time. Several facts are interesting here:

-First, except for very small instances, the percentage of eliminations decreases with the size of the instance, as shown in figure 3. The same phenomena applies for the elimination of pairs . This may be explained by the fact that there is a "forall" in the elimination criterion increasing the stringency of the criterion with the number of possibilities.

13

-The second fact is that elimination takes a long time on larger instances: due to the cubic complexity it is not possible to launch the double Goldstein DEE on largest instances in a reasonable time. However, even if double Goldstein cannot be applied, simple split DEE discards many pairs (see column 6).

Exact searches

The first series of results illustrates the progress made over [START_REF] Filipe | Protein design using answer set programming[END_REF] (columns 8 and 9 of table 1). A time limit has been set to 3 hours. A "-" denotes a problem that could not been solved due to time limit and a "!" denotes a problem that could not been solved due to insufficient memory. Many more instances have be solved, and none crashes because of memory. The difficult problems occur around 1800 rotamers instead of 100 in the previous approach. With the preprocessing we have implemented, the ASP solver is able to infer a good lower bound and it dramatically reduces the search space. Note that the simple idea consisting of reducing the energy computation to a spatial neighborhood of residues does not work in practice to break the complexity.

The second series of instances is taken from [START_REF] Allouche | Computational protein design as an optimization problem[END_REF] and has been built using Osprey2.0 [START_REF] Gainza | Osprey: protein design with ensembles, flexibility, and provable algorithms[END_REF]. These instances are far more complex, as they deal with several tens of positions to be redesigned and use a different library of rotamers that cannot be easily eliminated by DEE. Then, the threshold for instances that can not be solved is lower. The solver Toulbar2 presented in [START_REF] Allouche | Computational protein design as an optimization problem[END_REF] has by far the best performances. Instances (not in this table) that could not be solved with it could not be solved with other solvers. We wanted to compare ASP with Protsat2 [START_REF] Ollikainen | Sat-based protein design[END_REF], but the repository is empty and we could not reach the authors. Then, we compared to maxhs, another core-based SAT-solver, also used in benchmarks in [START_REF] Allouche | Computational protein design as an optimization problem[END_REF]. Maxhs is not very effective: many instances could not be solved. The most interesting comparison is with Cplex, a sophisticated solver for ILP. On some instances, our ASP-based solver outperforms Cplex, however on larger instances, Cplex takes the lead. It means that something in our approach makes it less scalable.

Enumeration of neighboring solutions

We ran solvable instances from the first series characterized by 17 positions to be redesigned and between 500 and 700 rotamers. For each instance, we have enumerated the solutions for ε ranging from 500 to 4500 with step 500. It represents about 0.5 to 8% of the free energy for these instances.

Figure 4 represents the number of solutions and time needed to produce them. They grow exponentially with ε. The main fact is that it is easy to get a big set of candidate proteins in small time if the optimal solution is easily found. On larger instances, the search for a set of proteins becomes hard to solve.

Explaining the differences

The first problem is the search for a lower bound. Our transformation algorithms are not the best possible ones: on big instances, the lower bound is far from the The ASP solver clingo4 chosen for this study had the advantage of proposing primitives to access the grounding and solving phases through a dedicated API. It is no more a black box but it still remains a "grey" box in the version we used (4.3.0): even if the user can parameterize a lot of features (optimization strategy, policy of restart. . . ), enforcing some procedure such as node and arc consistencies each time a choice is made was not possible. Systems such as Toulbar are relatively slow in exploring the search space, but they offer a better preprocessing and the possibility to interact more accurately with the solver.

Finally, it would require more work to study the influence of built-in ASP configurations on the search (frumpy, handy, crafty...).

Conclusion

The goal of this study was to evaluate the potential of ASP in solving the CPD problem. Our work lies in the scope of exact searches and can be compared to many different techniques: ILP, 01QP, MMRF, A* searches, CFN . . . We have chosen for the benchmarks some of the most recent and best approaches, which inspired our own work. We have established a first milestone on the subject, and are confident that it can still be largely improved.

Trying to reach performances of CFN is an interesting challenge, and may mobilize different competences. We also studied different approaches and proposed different techniques to quickly get approximated results that could be redirected as new constraints in the solver. Refining the use of ASP solver heuristics and improving the search for lower and upper bounds are the main perspectives of this research. We can also pinpoint that in practice, due to the approximations and the choice of energy function, the most stable structure may slightly differ from the computed optimum.The extension to neighboring solutions can be easily implemented in ASP by launching another solve process with an extended search after the search of the optimal solution, but this strategy suffers from the disadvantages of exact search: the search for a lower bound is worse than in many other approaches and then the search time grows much faster.

Among other variants, it may also be interesting to study the modeling of the CPD problem with a flexible backbone in ASP. On this problem, the best results are given by software such as Osprey. The most difficult challenge for ASP would be to solve the CPD problem with a continuous space of rotamers.

A great interest in trying different approaches to solve problems such as CPD may lead to new general ideas to solve difficult problems and to ameliorate solvers by merging the best ideas of each approach. As an interesting by-product of this work, we have identified some current limitations of the clingo's minimization process and this should help improving the optimization part in the solver.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Structure of an amino acid

Fig. 3 .

 3 Fig. 3. Bigger the instances are, smaller the ratio of eliminated rotamers or pairs is

Fig. 4 .

 4 Fig. 4. Number of solutions and solving time depending on ε

Table 1 .

 1 Results of algorithms performed on the first series of instances

	Protein Size AA # Rotamers # Rotamers Discarded pairs DEE time Time in Time in
			clusters	before	after			[10] our work
	1BE9	17 15	All Hydro All Hydro	2547 529 2063 435	1720 130 1023 56	270811 1903 122977 270	2919 3.2 1063 1.3	! -! 1.719	-0.66 155 0.05
	1I92	17 15	All Hydro All Hydro	2648 593 2055 455	1728 186 816 93	288332 4318 86194 1042	3262 6.0 806 2.1	! -! -	-1.26 1174 0.10
	1MFG	17 15	All Hydro All Hydro	2608 560 1972 410	1798 188 854 74	299418 3717 106835 485	3145 5.3 759 1.3	! -! 679	-0.51 336 0.07
	1N7F	17 15	All Hydro All Hydro	2548 544 1934 402	1525 136 648 57	211868 1811 59369 202	2473 4.2 477 1.2	! -! 66	-0.21 15.6 0.047
	1QAU	17 15	All Hydro All Hydro	2643 562 2023 421	1776 196 895 97	288165 4177 101761 1021	1598 6.2 1228 1.8	! -! 2295	-0.73 98 0.10
	1RZX	17 15	All Hydro All Hydro	2722 554 2121 404	1748 167 941 56	257158 3044 94310 205	3398 4.5 883 1.1	! -! 36	-0.37 906 0.05
	1TP3	17 15	All Hydro All Hydro	2582 536 2142 447	1732 142 1037 63	282007 2253 135604 409	3145 3.8 1057 1.7	! -! 463	-0.19 403 0.054
	2EGN	17 15	All Hydro All Hydro	2548 546 2005 413	1570 150 804 75	240733 2516 78377 655	2478 3.7 573 1.4	! -! 280	-0.32 1038 0.039
	2FNE	17 15	All Hydro All Hydro	2736 577 2144 438	1909 197 1060 98	325100 4639 132397 1093	3927 6.1 1158 1.7	! -! 4976	-0.49 154 0.10
	2GZV	17 15	All Hydro All Hydro	2934 613 2247 456	1991 189 1122 87	348785 3823 132910 696	5243 7.5 1340 2.1	! -! -	-0.62 1242 0.09

Table 2 .

 2 Comparison with the most recent tools in CPD

	Protein Positions # Rotamers ASP(s) Toulbar2(s) Cplex(s) maxhs
	2TRX	11	410	0.4	0.1	2.6	4086
	1HZ5	12	427	0.3	0.1	7.6	5695
	1PGB	11	438	2.7	0.1	3.6	5209
	1MJC	28	440	151	0.1	4.1	3698
	1UBI	13	498	209	0.2	139	-
	1CSK	30	508	2030	0.1	9.6	-
	1SHF	30	527	-	0.1	8.6	-
	2PCY	18	598	2589	0.2	26.9	-
	1SHG	28	613	366	0.2	39.4	-
	1NXB	24	625	-	0.2	17	-
	1FNA	38	887	-	0.5	121	-
	1CSP	30	1026	146	0.84	1264	-
	1BK2	24	1089	6.8	0.65	125	-
	1LZ1	59	1202	-	1.5	1084	-
	1FYN	23	2110	-	2.8	3136	-
	1CM1	17	2242	-	3.3	473	-

https://www.mat.unical.it/aspcomp2013/ASPStandardization

http://sourceforge.net/projects/protsat/

Algorithm 4 Enforcing consistency

for all position i do for all amino-acid A at position i do for all position j = i do for all amino-acid B at position j do for all rotamers Ar at position i do from residue to rotamers(j,B) E(iA r ) ← E(iA r ) + E(iA) end for E(iA) ← 0 end for for all rotamer r at position j do for all rotamers B r at position j do from rotamer to pairs(j,r ,i) E(iA r , jB r ) ← E(iA r , jB r ) + E(iA r ) end for E(iA r ) ← 0 end for end for for all rotamer Ar at position i do Energy transfer from pairs to rotamer(i,Ar) end for Energy transfer from rotamers to residue(i,A) for all position j = i do for all rotamer B r at position j do Energy transfer from pairs to rotamer(j,Br) end for for all amino-acid B at position j do Energy transfer from rotamers to residue(j,B) end for end for E ← min E(iA) for all amino-acids A at position i do from residues to the lower bound

Acknowledgments

Authors would like to thank I. Lynce and F. Gouveia at INESC-ID Lisboa and T. Schiex, S. Barbe and D. Allouche at INRA Toulouse for having made this study possible with their own work and having kindly provided their benchmarks. We also thank the reviewers for their constructive suggestions.