Alessandro Renzaglia
email: alessandro.renzaglia@inrialpes.fr

Lefteris Doitsidis
email: ldoitsidis@chania.teicrete.gr

Agostino Martinelli
email: agostino.martinelli@inrialpes.fr

Elias B Kosmatopoulos

Multi-Robot 3D Coverage of Unknown Areas

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I Introduction

The use of multi-robot teams has gained a lot of attention in recent years. This is due to the extended capabilities that the teams have to offer comparing to the use of a single robot for the same task. Robot teams can be used in a variety of missions including: surveillance in hostile environments (i.e., areas contaminated with biological, chemical or even nuclear wastes), environmental monitoring (i.e., air quality monitoring, forest monitoring) and law enforcement missions (i.e., border patrol), etc. In all the aforementioned tasks the deployment of limited resources (robots) to optimize the monitoring of the area is the key issue. In order to achieve this purpose, the trajectories of the robots should be designed, in real-time, so that:

(O1) the part of the terrain that is monitored (i.e., is visible) by the robots is maximized;

(O2) for every point in the terrain, the closest robot is as close as possible to that point.

The second objective is necessary for two practical reasons: (a) firstly, the closer is the robot to a point in the terrain the better is, in general, its sensing ability to monitor this point and (b) secondly, in many multi-robot coverage applications there is the necessity of being able to intervene as fast as possible in any of the points of the terrain with at least one robot. The majority of existing approaches for multi-robot surveillance coverage, which concentrate mostly on the 2D case of ground robots, deal only with one of the objectives (O1) or (O2) and in most of them the terrain morphology is considered convex and/or known. In such cases the problem of multi-robot surveillance coverage can be seen to be equivalent to a standard optimization problem where the robot trajectories are generated according to a gradient-descent or gradientdescent-like methodology. However, in the case where it is required that both of the objectives (O1) and (O2) are addressed and the terrain's morphology is nonconvex and unknown, standard optimization tools are not applicable anymore as these tools require full knowledge of an objective function that depends on the unknown terrain's morphology.

To approach this problem, we propose a new solution that is based on the recently introduced Cognitive-based Adaptive Optimization (CAO) algorithm [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF], [START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF]. The main advantage of CAO as compared to standard optimization tools is that it does not require that the objective function to be optimized is explicitly known; CAO instead requires that at each-time instant a value (measurement) of this objective function is available. Then, if it is possible to define an objective function which is available for measurement for every given team configuration, the CAO methodology will be directly applicable to the problem of surveillance coverage treated in this paper. Rigorous arguments establish that, despite the fact that the terrain's morphology is unknown, the CAO methodology shares the same convergence characteristics as those of constrained gradient-descent algorithms (which require perfect knowledge of the terrain's morphology and optimize surveillance coverage subject to the constraints the robot team has to satisfy). As a result, CAO navigates the robots to an arrangement that locally optimizes the surveillance coverage criterion while satisfying physically-imposed constraints such as that the robots should not leave a prespecified area or they should not hit the terrain.

Both objectives (O1) and (O2) can be characterized by two distinct objective functions. In general, one cannot simultaneously optimize both functions, unless the functions share common optima. Hence, the idea is to optimize a combined objective function that strikes a compromise between maximizing visible area and minimizing the distance of the robots to points in the environment. By introducing such objective function, we achieve to render the CAO algorithm applicable to the particular problem of 3D multi-robot surveillance coverage treated in this paper. This objective function depends on the unknown terrain's characteristics and thus its explicit form is not known. However, for any given team configuration the value of this objective function can be directly computed from the robots' sensor measurements, and thus the CAO algorithm can be applied to the problem at hand by using such an objective function.

It has to be emphasized that apart from rendering the optimization problem solvable, the CAO-based approach preserves additional attributes that make it particularly tractable: it can handle in a easily fashion a variety of physical constraints and limitations and it is fast and scalable.

We finally mention that the CAO approach extends the popular Simultaneous Perturbation Stochastic Approximation (SPSA) algorithm [START_REF] Spall | Multivariate stochastic approximation using a simultaneous perturbation gradient approximation[END_REF]. The difference between the SPSA and the CAO approach is that SPSA employs an approximation of the gradient of an appropriate cost function using only the most recent experiments, while the CAO approach employs linear-in-the-parameters approximators that incorporate information of a user specified time window of the past experiments together with the concept of candidate perturbations for efficiently optimizing the unknown function. Comparative evaluations that were performed on complicated optimization problems have shown that CAO exhibits significantly better convergence properties than SPSA [START_REF] Kosmatopoulos | Adaptive fine-tuning of nonlinear control systems with application to the urban traffic control strategy tuc[END_REF][START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF]. Moreover, CAO was shown to exhibit satisfactory (local) convergence characteristics in particular problems where SPSA failed to provide convergent solutions for any choice of its design parameters, [START_REF] Kosmatopoulos | Adaptive fine-tuning of nonlinear control systems with application to the urban traffic control strategy tuc[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF].

It is mentioned that the CAO or the SPSA do not create an approximation or estimation of the obstacles location and geometry; instead, they on-line produce a local approximation of the unknown cost function the robots are called to optimize. For this reason, they require simple -and thus scalable -approximation schemes to be employed.

I.1 Related Work

The majority of approaches for multi-robot surveillance coverage concentrate on objective (O2) described in the previous section. In [START_REF] Cortés | Coverage control for mobile sensing networks[END_REF], the authors present a controller for the coverage with a team of mobile robots of a convex environment, i.e., without obstacles, based on the Voronoi partition. A similar approach, for a convex environment, is proposed in [START_REF] Schwager | Distributed coverage control with sensory feedback for networked robots[END_REF], where additionally the robots estimate a function indicating the relative importance of different areas in the environment, using information from the sensors. A possible approach for non-convex regions is proposed in [START_REF] Pimenta | Sensing and coverage for a network of heterogeneous robots[END_REF]. In this work the Voronoi partition is obtained by using the geodesic distance instead of the Euclidean one taking into account the particular topology of the problem. In (Howard et al., 2002b), the same problem is approached by using the artificial potential field method. Another possible solution for environments which include obstacles is proposed in (Breitenmoser et al., 2010b): the main idea is to combine the classical Voronoi coverage with the Lloyd algorithm and the local path planning algorithm TangentBug. In all the aforementioned works the regions to cover are in 2D. In (Breitenmoser et al., 2010a) the authors approach also the problem of deploying a team of robots on a non-planar surface in 3D space.

As far as it concerns objective (O1) described in the previous section, different solutions have been proposed in the literature. In [START_REF] Ganguli | Maximizing visibility in nonconvex polygons: nonsmooth analysis and gradient algorithm design[END_REF] the authors propose a gradient-based algorithm for the case of a single robot case and they prove that the visible area is almost everywhere a locally Lipschitz function of the observer location. In [START_REF] Ganguli | Visibility-based multi-agent deployment in orthogonal environments[END_REF], an approach for the multi-robot problem is presented based on the assumption that the environment is simply connected. The visibility problem is also related with the Art Gallery Problem where the goal is to find the optimum number of guards in a non-convex environment so that each point of the environment is visible by at least one guard [START_REF] Agarwal | Efficient algorithms for geometric optimization[END_REF], [START_REF] Shermer | Recent results in art galleries[END_REF]. All the aforementioned solutions are based on the hypothesis that a given point can be monitored regardless of its distance from the robot. An incremental algorithm which takes into consideration also a maximum monitoring distance is presented in (Howard et al., 2002a). In [START_REF] Schwager | Optimal coverage for multiple hovering robots with downward facing camera[END_REF], the authors consider the coverage of a 2D region by using a team of hovering robots. In this case, information per pixel is proposed as optimization criterion.

To the best of our knowledge, the problem of considering the two objectives simultaneously to cover a 3D region by using a team of flying robots has never been investigated so far. To do that we propose to use a new stochastic optimization method, the CAO algorithm. The many advantages of using stochastic gradient descent algorithms, like the SPSA algorithm, to approach a sensor-based deployment problem have already been highlighted in (Ny and [START_REF] Ny | Sensor-based robot deployment algorithms[END_REF]. In this work, the authors proposed applications such as: coverage with heterogeneous sensors and source seeking with stochastic wireless connectivity constraints.

The rest of the paper is organized as follows. In section 2 we describe in detail the cognitive based adaptive optimization approach, while in section 3 we formulate the problem for the 3D multi-robot coverage over unknown terrains. In section 4 extensive experimental results are presented and finally in section 5 concluding remarks and future research are discussed.

II The Cognitive-based Adaptive Optimization Approach

The Cognitive-based Adaptive Optimization (CAO) approach [START_REF] Kosmatopoulos | Adaptive fine-tuning of nonlinear control systems with application to the urban traffic control strategy tuc[END_REF][START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] was originally developed and analyzed for the optimization of functions for which an explicit form is unknown but their measurements are available as well as for the adaptive fine-tuning of large-scale nonlinear control systems. In this section, we will describe how the CAO approach can be appropriately adapted and extended so that it is applicable to the problem of multi-robot coverage. More explicitly, let us consider the problem where M robots are involved in a coverage task, attempting to optimize a given coverage criterion. The coverage criterion is a function of the robots' positions or poses (positions and orientations), i.e.,

J k = J x (1) k , . . . , x (M ) k (2.1)
where k = 0, 1, 2, . . . denotes the time-index, J k denotes the value of the coverage criterion at the k-th time-step, x

k , . . . , x (M ) k denote the position/pose vectors of robots 1, . . . , M , respectively, and J is a nonlinear function which dependsapart from the robots' positions/poses -on the particular environment where the robots live; for instance, in the 2D case the function J depends on the location of the various obstacles that are present, while in the 3D case with flying robots monitoring a terrain, the function J depends on the particular terrain morphology.

Due to the dependence of the function J on the particular environment characteristics, the explicit form of the function J is not known in most practical situations; as a result, standard optimization algorithms (e.g., steepest descent) are not applicable to the problem in hand. However, in most practical cases, like the one treated in this paper, the current value of the coverage criterion can be estimated from the robots' sensor measurements. In other words, at each time-step k, an estimate of J k is available through robots' sensor measurements,

J n k = J x (1) k , . . . , x (M ) k + ξ k (2.2)
where J n k denotes the estimate of J k and ξ k denotes the noise introduced in the estimation of J k due to the presence of noise in the robots' sensors. Please note that, although it is natural to assume that the noise sequence ξ k is a stochastic zeromean signal, it is not realistic to assume that it satisfies the typical Additive White Noise Gaussian (AWNG) property even if the robots' sensor noise is AWNG: as J is a nonlinear function of the robots' positions/poses (and thus of the robots' sensor measurements), the AWNG property is typically lost.

Apart from the problem of dealing with a criterion for which an explicit form is not known but only its noisy measurements are available at each time, efficient robot coverage algorithms have additionally to deal with the problem of restricting the robots' positions so that obstacle avoidance as well as robot formation constraints are met. In other words, at each time-instant k, the vectors x (i) k , i = 1, . . . , M should satisfy a set of constraints which, in general, can be represented as follows:

C x (1) k , . . . , x (M ) k ≤ 0 (2.3)
where C is a set of nonlinear functions of the robots' positions/poses. As in the case of J , the function C depends on the particular environment characteristics (e.g., location of obstacles, terrain morphology) and an explicit form of this function may be not known in many practical situations; however, it is natural to assume that the coverage algorithm is provided with information whether a particular selection of robots' positions/poses satisfies or violates the set of constraints (2.3).

Given the mathematical description presented above, the multi-robot coverage problem can be mathematically described as the problem of moving x

(1) k , . . . , x (M ) k
to a set of positions/poses that solves the following constrained optimization problem: minimize (2.1) subject to (2.3) .

(2.4)

As already noticed, the difficulty in solving, in real-time and in real-life situations, the constrained optimization problem (2.4) lies in the fact that explicit forms for the functions J and C are not available. To circumvent this difficulty, the CAO approach, appropriately modified to be applicable to the problem in hand, is adopted. Indeed this algorithm is capable of efficiently dealing with optimization problems for which the explicit forms of the objective function and constraints are not known, but noisy measurements/estimates of these functions are available at each time-step. In the following, we describe the CAO approach as applied to the multi-robot coverage problem described above. It has to be emphasized that the CAO algorithm presented here is an extension of the CAO versions presented and analyzed in [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF]. The main difference is that while [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] address the unconstrained version of the problem (2.4), in the present paper the CAO approach of [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] has to be extended so that it efficiently takes care of the constraint (2.3). In order to do so, the CAO approach of (Kosmatopoulos, 2009, [START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] is augmented by a special -yet simple -projection mechanism. Theorem 1 establishes that the introduction of such a projection mechanism does not destroy the nice properties of the CAO approach of [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF]; as a matter of fact, according to Theorem 1 presented below, the CAO algorithm used in this paper is proven to be approximately a projected gradient-descent algorithm, while the ones of [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] have been established to be approximate unconstrained gradient-descent algorithms.

As a first step, the CAO approach makes use of function approximators for the estimation of the unknown objective function J at each time-instant k according to Ĵk x

(1) k , . . . , x

(M ) k = ϑ τ k φ x (1) k , . . . , x (M ) k .
(2.5)

Here Ĵk x (1) k , . . . , x (M ) k
denotes the approximation/estimation of J generated at the k-th time-step, φ denotes the nonlinear vector of L regressor terms, ϑ k denotes the vector of parameter estimates calculated at the k-th time-instant and L is a positive user-defined integer denoting the size of the function approximator (2.5). The vector φ of regressor terms must be chosen so that it satisfies the socalled Universal Approximation Property [START_REF] Polycarpou | Identification and control of nonlinear systems using neural network models: Design and stability analysis[END_REF], i.e. it must be chosen so that the approximation accuracy of the approximator (2.5) is an increasing function of the approximator's size L. Polynomial approximators, radial basis functions, kernel-based approximators, etc, are known to satisfy such a property, see [START_REF] Polycarpou | Identification and control of nonlinear systems using neural network models: Design and stability analysis[END_REF] and the references therein.

The parameter estimation vector ϑ k is calculated according to

ϑ k = argmin ϑ 1 2 k-1 ℓ=ℓ k J n ℓ -ϑ τ φ x (1) ℓ , . . . , x (M ) ℓ 2 (2.6)
where ℓ k = max{0, k -L -T h } with T h being a user-defined nonnegative integer. Standard least-squares optimization algorithms can be used for the solution of (2.6).

Remark 1 In order for the proposed methodology to guarantee with efficient performance, special attention has to be paid in the selection of the regressor vector φ. Polynomial or polynomial-like regressor vectors as well as sigmoidal regressor vectors can be employed for the construction of φ. The particular choice adopted for the application treated in this paper is described in section III. See [START_REF] Kosmatopoulos | Adaptive fine-tuning of nonlinear control systems with application to the urban traffic control strategy tuc[END_REF][START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] for more details on the design considerations for the regressor vector. ⋄

As soon as the estimator Ĵk is constructed according to (2.5), (2.6), the set of new robots' positions/poses is selected as follows: firstly, a set of N candidate robots' positions/poses is constructed according to 1

x i,j k = x (i) k + α k ζ i,j k , i ∈ {1, . . . , M }, j ∈ {1, . . . , N } , (2.7)
where ζ i,j k is a zero-mean, unity-variance random vector with dimension equal to the dimension of x

(i)
k and α k is a positive real sequence which satisfies the conditions:

lim k→∞ α k = 0, ∞ k=1 α k = ∞, ∞ k=1 α 2 k < ∞ .
(2.8) Among all N candidate new positions x 1,j k , . . . , x M,j k , the ones that correspond to non-feasible positions/poses -i.e., the ones that violate the constraints (2.3) -are neglected and then the new robots' positions/poses are calculated as follows:

x (1) k+1 , . . . , x (M ) k+1 = argmin j ∈ {1, . . . , N } x i,j k not neglected Ĵk x 1,j k , . . . , x M,j k
The idea behind the above logic is simple: at each time-instant a set of many candidate new robots' positions/poses is generated. The candidate, among the ones that provide with a feasible solution, that provides the best estimated value Ĵk of the coverage criterion is selected as the new set of robots' positions/poses.

The random choice for the candidates is essential and crucial for the efficiency of the algorithm, as such a choice guarantees that Ĵk is a reliable and accurate estimate for the unknown function J ; see [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] for more details. On the other hand, the choice of a slowly decaying sequence α k , a typical choice of adaptive gains in stochastic optimization algorithms (see e.g., [START_REF] Bertsekas | Gradient convergence in gradient methods with errors[END_REF]) is essential for filtering out the effects of the noise term ξ k [cf. (2.2)]. The next theorem summarizes the properties of the CAO algorithm described above; as the proof of this theorem is among the same lines as the main results of [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF], only a sketch of the proof is provided.

Theorem 1 Let x (1 * ) , . . . , x (M * ) denote any -local -minimum of the constrained optimization problem (2.4). Let N ≥ 2M × dim x (i) k
and, moreover, the vector φ satisfy the Universal Approximation Property. Assume also that the functions 1 According to [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] J , C are either continuous or discontinuous with a finite number of2 discontinuities. Then, the CAO-based multi-robot coverage algorithm as described above guarantees that the robots' positions/poses x (1) k , . . . , x (M ) k will converge to one of the local minima x (1 * ) , . . . , x (M * ) almost surely, provided that the size L of the regressor vector φ is larger than a lower bound L.

Sketch of the Proof: Let X k denotes the augmented vector of all robots' positions/poses at time-instant k, i.e., the entries of X k are the entries of all vectors x

(1) k , . . . , x (M ) k . Using similar arguments as those of the main results of (Kosmatopoulos, 2009, [START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] it can be seen that at each iteration of the CAO-based algorithm described above, the new vector X k+1 satisfies

X k+1 = Π C {X k -α k (c∇J(X k ) + e k + b k )}
where c is a positive constant; Π C {•} denotes the projection operator onto the set G = {X : C (X) ≤ 0} defined as follows: for any X not satisfying the constraint C (X) ≤ 0, the point X = Π C {X} is the nearest point to X on G, where the norm is defined in the usual Euclidean norm; and e k , b k are two terms that are defined similarly to the respective terms in section III of [START_REF] Spall | Multivariate stochastic approximation using a simultaneous perturbation gradient approximation[END_REF]. By using the same arguments as those in the proof of Proposition 1 of [START_REF] Sadegh | Constrained optimization via stochastic approximation with a simultaneous perturbation gradient approximation[END_REF] [START_REF] Kosmatopoulos | Adaptive fine-tuning of nonlinear control systems with application to the urban traffic control strategy tuc[END_REF][START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF] for more details. ⋄

Remark 4 As an alternative to the CAO approach, the SPSA approach [START_REF] Spall | Multivariate stochastic approximation using a simultaneous perturbation gradient approximation[END_REF] may be employed in multi-robot coverage applications. According to the SPSA approach, the robot positions/poses are updated according to

x (i) k+1 = x (i) k + β k ζ i k , if k is even x (i) k+1 = x (i) k + γ k J n k -J n k-1 ζ (i) k-1 , if k is odd (2.9)
where

ζ (i)
k are zero-mean, unity-variance random vectors and β k , γ k are slowly decaying sequences (similar as the sequence α k ). The SPSA algorithm is computationally simpler than the CAO one, but it does not perform as efficient as the CAO approach as have been demonstrated in a variety of approaches, see [START_REF] Kosmatopoulos | Adaptive fine-tuning of nonlinear control systems with application to the urban traffic control strategy tuc[END_REF][START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF][START_REF] Kosmatopoulos | Large-scale nonlinear control system finetuning through learning[END_REF]. However, extensive simulation experiments have demonstrated that a hybrid scheme which uses SPSA at the first 10-20 time-steps and then switches to the CAO algorithm can have significant improvements over schemes that employ only the CAO algorithm. This is due to the fact that CAO, at its initial steps, may preserve a poor performance because it takes some iterations for the CAO estimator (2.5) in order to come up with a reliable estimate Ĵk of the unknown coverage function J . ⋄

Remark 5

We close this section by mentioning that similarly to the proposed approach, global optimization methods such as simulated annealing and genetic algorithms do not require that the explicit form of the function J is known. However, simulated annealing, genetic algorithms and other similar global optimization methods require that a large amount of different combinations of robots' positions is being evaluated all over the robots' application area. Such a requirement renders these methods practically infeasible as a huge amount of time and energy would have to be spent in order for the robots to visit many different locations all over their application area. Nevertheless, attempting to globally optimize surveillance coverage is practically infeasible as it is an NP-hard problem whose solution requires dense discretization over the space of all possible team configurations and evaluation of all points of the discretized space. ⋄

III CAO for 3D Multi -Robot Coverage over Unknown Terrains

III.1 Problem Definition

In our previous work [START_REF] Renzaglia | Cognitivebased adaptive control for cooperative multi-robot coverage[END_REF][START_REF] Renzaglia | Adaptivebased distributed cooperative multi-robot coverage[END_REF] we have extensively described the case of using the CAO approach for maximizing the monitored area in a given region by using a team of mobile robots in the 2D plane, without any assumption on the topology of the environment. In this section we will extend our approach to the case of robots living in a 3D environment and having fixed orientation. Consider a team of M flying robots that is deployed to monitor an unknown terrain T . Let z denote the unknown height of the terrain at the point (x, y) and assume for simplicity that the terrain T is rectangular along the (x, y)-axes, i.e., x min ≤ x ≤ x max , y min ≤ y ≤ y max . Let P = {x (i) } M i=1 denote the configuration of the robot team, where x (i) denotes the position of the i-th robot. We will say that a point q = (x, y, z) of the terrain is visible if there exists at least one robot so that

• the robot and the point q are connected by a line-of-sight;

• the robot and the point q are at a distance smaller than a given threshold value (defined as the maximum distance the robot's sensor can "see").

Given a particular team configuration P, we let V denote the visible area of the terrain, i.e., V consists of all points q ∈ T that are visible from the robots. We will assume that the robots are equipped with visual sensors together with inertial sensors and/or range sensors; in other words, for each visible point we will assume that the team is able to estimate the terrain's height at this point. A possible way to deploy a robot team satisfying the above, is by using the down-looking-cameraequipped flying robots of [START_REF] Bloesch | Vision based mav navigation in unknown and unstructured environments[END_REF], Weiss et al., 2010) which employ the monocular SLAM algorithm of [START_REF] Klein | Parallel tracking and mapping for small ar workspaces[END_REF].

The main objective for the robot team is to maximize the visible area V. However, this cannot be the only objective for a robot team in a coverage task: trying to maximize the visible area will simply force the robots to climb as high as 3 their visibility threshold allows for. For this reason, while maximizing the visible area is the primary goal of the mission, the team members should be deployed so that for every point in the terrain, the closest robot is as close as possible to that point. In other words, among all possible configurations that maximize the visible area V, the robot team should converge to the one that keeps as small as possible the average distance between each of the robots and the terrain subarea the particular robot is responsible for, where the subarea of the terrain the i-th robot is responsible for is defined as the part of the terrain that (a) is visible by the i-th robot and (b) each point in this subarea is closer to the i-th robot than any other robot of the team. Attempting to keep the robots as close as possible to the subarea they are responsible for is necessary for two practical reasons: (a) firstly, the closer is the robot to a point in the terrain the better is, in general, its sensing ability to monitor this point and (b) secondly, in many multi-robot coverage applications there is the necessity of being able to intervene as fast as possible in any of the points of the terrain with at least one robot.

Having the above reasoning in mind, we define the following combined objective function the robot team has to minimize:

J(P) = q∈V min i∈{1,...,M } x (i) -q 2 dq + K q∈T -V dq (3.10)
where K is a user-defined positive constant and • denote the Euclidean norm.

The first of the terms in above equation is the usual cost function considered in many coverage problem for known 2D environment related to the second objective (minimize the average distance between the robots and the subarea they are responsible for, see [START_REF] Cortés | Coverage control for mobile sensing networks[END_REF]). The second term is related to the invisible area in the terrain ( q∈T -V dq is the total part of the terrain that is not visible by any of the robots).

The positive constant K serves as a weight for giving less or more priority to one of the objectives. For instance, in the case where K = 0, we will have that the robots, in their attempt to minimize their average distance to the subarea they are responsible for, may also seek to minimize the total visible area. It has to be emphasized that the positive constant K should be chosen sufficiently large so that the second term in (3.10) dominates the first term unless no or a negligible part of the terrain remains invisible. In this way, minimization of (3.10) is equivalent to firstly making sure that all, or almost all, of the terrain is visible and then to locate the robots so that their average distance to the subarea they are responsible for is minimized. However, choosing a value for K so that the second term in (3.10) dominates the first term is not straightforward unless the terrain is known. Later in this section we will comment further on how to choose the parameter K for the particular setup considered in this paper. Please also note that whether the CAObased algorithm employing a large K converges to negligible or non-negligible invisible areas depends on the number, the sensing capabilities and maximum height constraints of the robots as well as the terrain's complexity.

The second term q∈T -V dq in (3.10) cannot be, in general, computed in practice; as this term involves the part of the terrain that is not currently visible, its computation requires that the geometry of this part is known or equivalently -as the invisible part changes with the evolution of the team's configuration -that the whole terrain is known. To overcome this problem, instead of minimizing (3.10) the following performance index is actually minimized by the CAO approach:

J(P) = q∈V min i∈{1,...,M } x (i) -q 2 dq -K q∈V dq (3.11)
To see that minimization of (3.11) and (3.10) is equivalent, please note that q∈T -V dq = q∈T dq -q∈V dq and the integral q∈T dq is constant.

Throughout the above analysis, the assumption that K is "sufficiently large" was made. If such an assumption does not hold, the arguments presented above do not hold either. Therefore, the question that naturally arises is what is a value for K that is sufficiently large so that these arguments hold. As a very large choice for K (e.g., K = 10 10 ) can lead to numerical instability problems [switchinglike performance for the algorithm when the numerically-computed invisible area switches from small values to zero], a guide on how to choose K so that such instability problems do not occur should be provided. Extensive simulations with all set-ups considered in the next section (IV. Simulation Results) and with different values for K indicate that it suffices to choose K to be 3 -50 times the parameter max in order to get an efficient performance, where the parameter max can be calculated as follows: let

f (P) = q∈V min i∈{1,...,M } x (i) -q 2 dq, ḡ(P) = q∈T -V I(x, y)dxdy
where I(x, y) denotes the indicator function that is equal to 1 if the point q belongs to the invisible area of the terrain and is zero, otherwise. In other words, the term ḡ(P) would correspond to the total invisible area, if the unknown terrain points (x, y, z) were replaced by (x, y, 1), i.e., if the whole invisible area were flat. Then the parameter max is calculated according to max ≈ sup f (P) sup ḡ(P) i.e., max corresponds to the maximum possible value for f (P) (over all possible feasible team configurations) divided by the maximum value the invisible-areaterm ḡ(P) can take. The sup ḡ(P) is equal to the terrain's area along the (x, y)axes, i.e. sup ḡ(P) = (x max -x min )(y max -y min ). On the other hand an estimate of the term sup f (P) can be produced by running extensive simulations with randomly generated terrains and randomly generated team's configurations. Figure 1 shows the time-histories of the terms f (P) and ḡ(P) for different choices for K and for one of the scenarios described in the simulations section (more precisely, for the scenario described in section IV.2.2). For this particular scenario, we have that sup f (P) ≈ 1000 and sup ḡ(P) = 100 and thus the parameter max can be estimated to be around 10. As exhibited in Figure 1 Figure 1: Time-histories for the terms f (P) = q∈V min i∈{1,...,M } x (i) -q 2 dq (upper plot) and ḡ(P) = q∈T -V I(x, y)dxdy (lower plot) for different values of the parameter K: K = 5 (red), K = 30 (black), K = 100 (blue) and K = 500 (green).

III.2 Algorithm Implementation

An efficient trajectory generation algorithm for optimal coverage must make sure that the physical constraints are also met throughout the whole multi-robot coverage application. Such physical constraints include, but are not limited to, the following ones:

• The robots remain within the terrain's limits, i.e., within [x min , x max ] and [y min , y max ] in the xand y-axes, respectively.

• The robots satisfy a maximum height requirement while they do not hit the terrain, i.e., they remain within [z + d, z max ] along the z-axis, where d denotes the minimum safety distance (along the z-axis) the robots' should be from the terrain and z max denotes the maximum allowable height for the robots.

• The robots do not come closer to the other ones than a minimum allowable safety distance d r .

It is not difficult for someone to see that all the above constraints can be easily cast in the form (2.3) and thus can be handled by the CAO algorithm.

Having defined the optimization problem, a fundamental point for a good behavior of the CAO algorithm is an appropriate choice of the form of the regressor vector φ, introduced in equation (2.5). As mentioned in remark 1, several different choices for its explicit expression are admissible. However, in all different tests for the particular application treated in this paper, it was found that it suffices to choose the regressor vector as follows:

1. choose the size of the function approximator L to be an odd number; 2. select the first term of the regressor vector φ to be the constant term; 3. select randomly the next (L -1)/2 terms of φ to be any 2nd-order terms of the form x

(i) a • x (j)
b [with a, b ∈ {1, . . . , dim(x (i) )}, i, j ∈ {1, . . . , M } randomly-selected positive integers]; 4. select the last (L-1)/2 terms of φ to be any 3rd-order terms of the form x

(i) a • x (k) b • x (j) c [with a, b, c ∈ {1, . . . , dim(x (i) )}, i, k, j ∈ {1, . . . , M } randomly- selected positive integers].
Once the regressor vector φ has been set and once the values of the cost function (3.11) are available for measurement at each time step, it is possible to find at each time step the vector of parameter estimates θ k and thus the approximation of the cost function Ĵk . The other important choice in order to assure the convergence of the algorithm is the expression of the sequence α k , defined in equation (2.7). A typical choice for such a sequence is given by

α k = α (k + 1) η , (3.12)
where α is a positive user-defined constant and η ∈ (0, 0.5).

Remark 6 Please note that the CAO algorithm's computational requirements are dominated by the requirement for solving the least-squares problem (2.6). As the number of free parameters in this optimization problem is L, most popular algorithms for solving least-squares problems have, in the worst case, O(L 3 ) complexity (polynomial complexity with respect to L). Furthermore, as it suffices to choose L around 2M ×dim x (i) k

(see Remark 3) we have that the computational requirements of the proposed algorithm are, at most, O(M 3 ). ⋄

IV Simulation Results

To evaluate the efficiency of the proposed approach, several scenarios were considered with the use of a simulated robot team which was able to move freely at the 3D plane. In all cases studied, the team was homogeneous with exactly the same monitoring capabilities. This assumption has been made for simplification purposes and easier comprehension of the results. The main constraints imposed to the robots are that they remain within the terrain's limits, i.e., within [x min , x max ] and [y min , y max ] in the xand y-axes, respectively. At the same time they have to satisfy a maximum height requirement while they do not hit the terrain, i.e., they remain within [z + d, z max ] along the z-axis. The scenarios considered are terrains with obstacles with same or uneven heights, while for each scenario different values of the parameter α which is responsible for the convergence of the algorithm were tested. Apart from the simulated terrains, a realistic scenario was considered by using a map of a real area [START_REF] Weiss | Intuitive 3d maps for mav terrain exploration and obstacle avoidance[END_REF], extracted with the methodology described in detail in [START_REF] Bloesch | Vision based mav navigation in unknown and unstructured environments[END_REF].

In all experiments reported next, the following choices were made for the algorithm's implementation:

• The CAO parameters N (number of next candidate robots' positions) and L (size of approximator φ) were set equal to 6M and 6M + 1, respectively, where M denotes the robot team's size, while the approximator φ was calculated as described in the previous section. Please note that the above choices for N, L and φ are in accordance to Theorem 1 and Remark 3; moreover, L was set equal to 6M + 1 as it has to be an odd number.

• The parameter K in the cost criterion (3.11) was set equal to 30 which satisfies K ∈ [10 max , 50 max ] for all terrains and team sizes considered in the simulations (see section III.1 for more details on the parameter max ).

• The parameter d (minimum allowable distance from the terrain) was set equal to 0.1, while the robot's were assumed to have unlimited visibility.

• Different choices for the parameters z max (maximum allowable height) and α (magnitude of next candidate robots' positions) were made as these parameters are the most crucial for the algorithm's performance.

• Finally, with the exception of the experiments reported in IV.4 that involve teams of 10 and 20 robots, in all other cases the team comprised 4 robots.

IV.1 Areas with same height obstacles

The first case considered, studies an area sizes 10 by 10 meters, which includes a surface with seven same height randomly placed obstacles. For this area, several scenarios were tested with the robot team starting from different initial positions and heights. In all cases considered the robots had to satisfy a maximum flight height requirement while they did not hit the terrain.

IV.1.1 Scenario 1

In the first scenario studied, all the team members were placed at starting points adjunct to each other, with initial height 0.6 meters. The maximum allowed flight height was 1 meter for all robots. Different values of the expression α were tested for the case of α = 0.3, 0.5, 1 and the respective cost functions are presented in Fig. 2. The initial position for Robot 1 was (0.18, 0.2, 0.6), of Robot 2 was (0.19, 0.2, 0.6), of Robot 3 was (0.2, 0.2, 0.6) and of Robot 4 was (0.21, 0.2, 0.6).

In Fig. 3 successive snapshots of different positions of the robot team for the case of α = 0.3 are presented (different color corresponds to different team member).

The final configuration in all three test cases is presented in Fig. 4. In Table 1 the percentage of the initial and final coverage of the area monitored in all three cases, is presented. It's worth mentioning that the coverage percentage is depended on several factors apart the optimization algorithm i.e. the sensors that might be used in a real implementation. It should be noted that CAO does not converge always to the same swarm configuration, but it converges always to a swarm configuration with similar coverage characteristics which corresponds to similar final J value. In the second scenario studied, the initial positions of the robots were for Robot 1 (9.18, 0.2, 0.4), for Robot 2 (9.19, 0.2, 0.4), for Robot 3 was (0.2, 0.2, 0.4) and for Robot 4 was (0.21, 0.2, 0.4), while the maximum allowed flight height remained the same (1 meter). Different values of the expression α were tested for the case of α = 0.3, 0.5 and the respective cost functions are presented in Fig. 6. In table 2 the percentage of the initial and final coverage of the area monitored in both cases is presented.

IV.2 Areas with uneven obstacle height

The second case considered, studies an area sizes 10 by 10 meters, which includes a surface with seven randomly placed obstacles with uneven height, with maximum value 2 meters. In this test case we have tested several scenarios with the robot team starting from different initial positions and heights. In all cases considered the robots had to satisfy a maximum flight height requirement while they did not hit the terrain.

IV.2.1 Scenario 1

In the first scenario studied for the case of areas with uneven obstacle heights, all the team members were placed at starting points adjunct to each other, with initial height 0.2 meters. The maximum allowed flight height was 1 meter for all robots. The initial positions of Robot 1 was (0.18, 0.2, 0.2), of Robot 2 was (0.19, 0.2, 0.2), of Robot 3 was (0.2, 0.2, 0.2) and of Robot 4 was (0.21, 0.2, 0.2). Different values of the expression α were tested for the case of α = 0.3, 0.5, 1 and the respective cost functions are presented in Fig. 7. The final configuration in all three test cases is presented in Fig. 8, while in Table 3 the percentage of the initial and final coverage of the area monitored in all cases, is presented.

IV.2.2 Scenario 2

In the second scenario studied for uneven surfaces, the maximum allowed flight height was 5 meters for all robots. Different values of the expression α were tested for the case of α = 0.3, 0.5, 1 and the respective cost functions are presented in Fig. 9. In Fig. 10 successive snapshots of different positions of the robot team for the case of α = 0.3 are presented (different color corresponds to different team member). The initial position of Robot 1 was (0.18, 0.2, 0.4), of Robot 2 was (0.19, 0.2, 0.4), of Robot 3 was (0.2, 0.2, 0.4) and of Robot 4 was (0.21, 0.2, 0.4). In table 4 the percentage of the initial and final coverage of the area monitored in all cases, is presented.

IV.3 Cave-like Surface

In the above described set-up and proposed simulations, for simplicity's sake we assumed that the unknown terrain is defined as a set of unique triplets (x, y, z), that is, for each (x, y) the terrain is defined by a unique z-point, i.e. z = f (x, y). In realistic applications there exist cases where there may be more than one z-points (e.g., cases of terrains with buildings, overhangs, ledges, caves, etc). Here we present a similar scenario to show how our method can be applied also for these cases. The simulated environment is a gaussian with a cave (fig. 11). The robots start their mission on the other side with respect to the cave, so at the beginning the it is not visible. that it is minimized until everything is visible.

IV.4 Scalability Issues

To validate the efficiency of our approach in the case of bigger robot teams, we have performed experiments with teams consisting of 10 and 20 members. Our experiments were performed in an area sizes 20 by 20 meters, which includes a surface with fifteen uneven height randomly placed obstacles. The maximum flight height was set to be 2. The basic difference as far as it concerns the computational requirements in the experiments conducted with the teams of 10 and 20 robots, was that the parameters L and N increase linearly according to Theorem 1 and Remark 3; therefore L = 61 and N = 60 in the case of the team with 10 members and L = 121 and N = 120 in the case of the team with 20 members. In the case of 4 robots the best values of J are around 15 which is significantly larger than the values obtained with the bigger teams. In Fig. 14 the cost functions for the case of 10 and 20 robots are presented, while in Fig. 15 we present their final configuration. In table 5 the percentage of the initial and final coverage of the area monitored for a team with 10 & 20 members is presented.

IV.5 Birmensdorf area

To validate our approach in a realistic environment, we used the data which were collected with the use of a miniature quadrotor helicopter specially designed for the needs of the European project sFLY (www.sfly.org). These data correspond to the Birmensdorf area presented in Fig. 16. The area was mapped using a stateof-the-art visual-SLAM algorithm which tracks the pose of the camera while, simultaneously, autonomously, building an incremental map of the surrounding environment. The area mapped using the aforementioned methodology corresponds to the bounded part of Fig. 16. More details about the data and the methodology used, are presented in [START_REF] Bloesch | Vision based mav navigation in unknown and unstructured environments[END_REF] and [START_REF] Weiss | Intuitive 3d maps for mav terrain exploration and obstacle avoidance[END_REF]. The main constraints imposed to the robots are that they remain within the terrain's limits, i.e., within [x min , x max ] and [y min , y max ] in the xand y-axes, respectively. At the same time they have to satisfy a maximum height requirement while they do not hit the terrain, i.e., they remain within [z + d, z max ] along the z-axis. The value of α was equal to 0.3. Several initial configurations for the robot team were tested. In Fig. 17 the cost function of an illustrative scenario is presented, while the final configuration of the team (for the same scenario) is displayed in Fig. 18 (3D view) and in Fig. 19 (2D side view). The initial coverage was 44.49% while the final coverage was 98.55%. In table 6 the final coverage percentage for different initial configuration in the Birmensdorf area, is presented.

V Discussion and Conclusions

A new method for dealing with the problem of performing surveillance coverage in unknown terrain of complex and non-convex morphology has been proposed. The proposed approach has the following key advantages:

• it does not require any a priori knowledge on the environment; • it works in any given environment, without the necessity to make any kind of assumption about its topology;

• it can incorporate any kind of constraints;

• it does not require a knowledge about these constraints since they are learnt during the task execution;

• its complexity is low allowing real time implementations.

The advantages of the proposed methodology make it suitable for real implementations and the results obtained through numerical simulations give us the motivation to adopt the CAO also in other frameworks. We are interested into formulating the same problem in a distributed manner by using different cost functions for the robots in the team. This approach is closer to real world applications since it will not depend into a centralized scheme with all the known disadvantages. Apart from that a decentralized approach will allow us to include communications constraints. We are also interested in incorporating more realistic constraints including sensor limitations. Our aim is to develop a strategy for the surveillance of an unknown urban-like environment with a real MAV swarm. Furthermore, we are interested in enhancing the proposed methodology in order to be able to deal with cases where the team converges to a configuration that fails to cover the entire terrain. It is worth noticing, that although in the 3D case treated here we have never encountered such cases there is always the possibility (due to the local convergence properties of the algorithm) for the methodology to fail to cover the overall area. As a matter of fact, such cases have been encountered in the 2D version of the algorithm reported in [START_REF] Renzaglia | Cognitivebased adaptive control for cooperative multi-robot coverage[END_REF][START_REF] Renzaglia | Adaptivebased distributed cooperative multi-robot coverage[END_REF]) and an extension/enhancement of the proposed methodology is required to efficiently deal with such cases.

Furthermore, we expect that many important tasks in mobile robotics can be approached by CAO-based algorithms: for example coordinated exploration, optimal target tracking, multi-robot localization, and so on. This is basically due to the fact that the CAO approach does not require an a priori knowledge of the environment and it has low complexity. Both these issues are fundamental in mobile robotics.

  it suffices to choose N to be any positive integer larger or equal to 2×[the number of variables being optimized by CAO]. In our case the variables optimized are the robot positions/poses x (1) k , . . . , x (M ) k and thus it suffices for N to satisfy N ≥ 2M × dim x (i) k.

  , the CAO-based algorithm converges to negligible invisible areas for all values of K satisfying K ∈ [3 max , 50 max ] ≡ [30, 500].

Figure 2 :

 2 Figure 2: Cost Functions for α = 0.3, 0.5, 1, in the case described in IV.1.1.

Figure 4 :

 4 Figure 4: Final positions of the robotic teams for α = 0.3 (blue markers), α = 0.5 (red markers), α = 1 (green markers), in the case described in IV.1.1.

Figure 5 :Figure 6 :

 56 Figure 5: Coverage percentage for α = 0.3, 0.5, 1, in the case described in IV.1.1.

Figure 7 :

 7 Figure 7: Cost Functions for α = 0.3, 0.5, 1, in the case described in IV.2.1.

Figure 8 :

 8 Figure 8: Final positions of the robotic teams for α = 0.3 (blue markers), α = 0.5 (red markers), α = 1 (green markers), in the case described in IV.2.1.

Figure 9 :

 9 Figure 9: Cost Functions for α = 0.3, 0.5, 1, in the case described in IV.2.2.

Figure 11 :

 11 Figure 11: A different scenario: the environment is a gaussian with a cave.

Figure 12 :

 12 Figure 12: The team is composed by four robots. Red squares and black circles represent initial and final positions respectively.

Figure 13 :

 13 Figure 13: Fig. (a) shows the behavior of the cost function during the task, fig. (b) the percentage of invisible surface. At the end everything is visible.

Figure 14 :

 14 Figure 14: Comparative cost functions for the case of a robot team with 10 and 20 members.

Figure 15 :

 15 Figure 15: Final configuration of the team with 10 robots (blue circle markers) and the team with 20 robots (red triangle markers).

Figure 16 :

 16 Figure 16: Outdoor flight path through the Birmensdorf area.

Figure 17 :

 17 Figure 17: Cost function for a robotic team performing surveillance coverage in Birmensdorf area.

Figure 18 :

 18 Figure 18: Final configuration of a robot team in a coverage scenario in Birmensdorf area (3D View).

Figure 19 :

 19 Figure 19: Final configuration of a robot team in a coverage scenario in Birmensdorf area (2D View).

Table 1 :

 1 Coverage percentage in the case described in IV.1.1.

		(% of Coverage)
	α	0.3	0.5	1
	Initial Configuration		34.58
	Final Configuration	97.06 97.49 98.59
	IV.1.2 Scenario 2

Table 2 :

 2 Coverage percentage in the case described in IV.1.2.

		(% of Coverage)
	α	0.3	0.5
	Initial Configuration	48.71
	Final Configuration	98.56	97.04

Table 3 :

 3 Coverage percentage in the case described in IV.2.1.

		(% of Coverage)
	α	0.3	0.5	1
	Initial Configuration		29.78
	Final Configuration	98.29 97.76 96.35

Table 4 :

 4 Coverage percentage in the case described in IV.2.2.

		(% of Coverage)
	α	0.3	0.5	1
	Initial Configuration		29.78
	Final Configuration	99.14 98.69 98.36

Table 5 :

 5 Coverage percentage with teams of 10 & 20 members.

		(% of Coverage)
	Team size	10 Members 20 Members
	Initial Configuration	39.62	42.93
	Final Configuration	86.33	90.78

Table 6 :

 6 Coverage percentage for different initial configuration in the Birmensdorf area.

		(% of Coverage)
	Test Case	1	2	3	4	5
	Initial Configuration 44.49 40.49 21.41 57.88 56.81
	Final Configuration	98.55 99.52 98.53 98.94 99.56

Please note that the family of "discontinuous functions with a finite number discontinuities" corresponds to the family of functions that can be approximated with arbitrary accuracy by continuous ones[START_REF] Jin | Neural networks and fuzzy basis functions for functional approximation[END_REF]. For instance, terrains with discontinuities along e.g., a closed or open curve belong to this family of functions and so do the corresponding functions J and C.

Note also that in the ideal case where there are no limits for the robot's maximum height and the robot has unlimited sensing capabilities, it suffices to have a single robot at a very high position to monitor the whole terrain.

Acknowledgments

The research leading to these results has received funding from the European Communities Seventh Framework Programme (FP7/2007-2013) under grant agreement n. 231855 (sFly). The authors would like to acknowledge S.Weiss and D. Scaramuzza from ETHZ/ASL for fruitful discussions and for providing the data set from Birmensdorf area. The authors would like to thank S. A. Chatzichristofis for his help in code optimization for the experimental results. The authors would


P V Milreu
email: milreu@biomserv.univ-lyon1.fr

V Acuña
email: viacuna@biomserv.univ-lyon1.fr

E Birmelé

P Crescenzi

A Marchetti-Spaccamela

M.-F Sagot

L Stougie
email: lstougie@feweb.vu.nl

V Lacroix
email: lacroix@biomserv.univ-lyon1.fr

Enumerating Chemical Organisations in Consistent Metabolic Networks: Complexity and Algorithms

The structural analysis of metabolic networks aims both at understanding the function and the evolution of metabolism. While it is commonly admitted that metabolism is modular, the identification of metabolic modules remains an open topic. Several definitions of what is a module have been proposed. We focus here on the notion of chemical organisations, i.e. sets of molecules which are closed and self-maintaining. We show that finding a reactive organisation is NP-hard even if the network is mass-and flux-consistent and that the hardness comes from blocking cycles. We then propose new algorithms for enumerating chemical organisations that are theoretically more efficient than existing approaches.

Introduction

Until recently, metabolism was analysed via the pathways composing it, which were traditionally established in a non automatic way by experts interested in some specific function (glycolysis for instance, or anaerobic respiration). Metabolic pathways may thus be seen as functional modules. One may wonder, however, about a possible arbitrariness in the definitions of the exact frontiers of each pathway, that may have been introduced by the experts because of some specific objectives they had in mind and of a partial lack of knowledge at the time the limits of the pathways were first drawn. Indeed, the advent of full genome sequences, which has also allowed inferring whole metabolic networks, has revealed the existence and actual use by different organisms of alternative metabolic routes to a same final overall product goal. The question thus rises whether automatic methods for inferring functional modules from whole networks would enable us to preserve the expert knowledge that led to the first pathways, while at the same time providing an insight into alternative functional ones. This is the biological motivation underlying the work presented in this paper. Its computational motivation is to explore one model for metabolic modules, both in terms of the complexity of enumerating such modules and of exact algorithms for performing the enumeration.

Several formal definitions of pathways and modules can be found in the literature on metabolism, the best known of which may be elementary modes [START_REF] Schuster | On elementary flux modes in biochemical reaction systems at steady state[END_REF] or any of its close cousins (see [START_REF] Lacroix | An Introduction to Metabolic Networks and Their Structural Analysis[END_REF][START_REF] Papin | Comparison of network-based pathway analysis methods[END_REF] for a survey). Elementary modes may be informally described as metabolic subnetworks that can function at steady state, meaning that all internal metabolites are produced and consumed in equal rates (that is, nothing accumulates internally). This is a fine definition, but has at least one drawback: it is restricted to the analysis of the system at steady state and does not allow to describe states of the system where metabolites can accumulate. However, such states are relevant as they could correspond to intermediary steps in the evolution of metabolism, or temporary states in the dynamics of metabolism.

As far as we know, two models in the literature enable to study such states. One is Petri nets and the other is a more recent model called chemical organisations. Because it is algebraically easier to manipulate chemical organisations as their formulation follows closely that of (hyper)graphs and matrices, we focus our attention in this paper on chemical organisations. The concept was introduced in 2005 by Peter Dittrich and his group [START_REF] Dittrich | Chemical organisation theory[END_REF], building on earlier work by Fontana and Buss [START_REF] Fontana | The Arrival of the fittest: towards a theory of biological organization[END_REF] and can be used not only for metabolism, but also for any kind of reaction system, including regulatory networks. In this paper, however, we focus exclusively on metabolism.

Chemical organisations are sets of molecules that are self-maintaining and closed (in this paper, we use the terms metabolite and molecule with no distinction). Informally, a self-maintaining set is a set where molecules can accumulate -the feature we were seeking -provided no molecule vanishes. A set is closed if all metabolites produced from reactions for which all the inputs are present in the set will also be present and thus part of the set. By convention, this includes all reactions that take their input from the environment, i.e. are external. All external inputs are therefore considered as being available and used. This introduces a second contrast with elementary modes (EMs). Indeed, EMs may use only part of the externally available inputs. More generally, EMs are not closed.

Finally, as we have already said, the theory of chemical organisations has been proposed for general reaction systems. Its application to metabolic networks raises new specific questions, as the networks have specific properties. They are indeed expected to be mass-consistent (reactions preserve mass) and fluxconsistent (each reaction belongs to at least one elementary mode).

The objectives of this paper are thus twofold. The first is to revisit chemical organisations in the context of mass-and flux-consistent networks. In particular, finding a chemical organisations was shown to be hard in [START_REF] Centler | Computing chemical organizations in biological networks[END_REF]. A legitimate and non trivial question is whether this remains true in biologically more realistic mass-and flux-consistent networks. Section 2 presents the main definitions on chemical organisations and consistency of networks. Section 3 shows that even for consistent networks the enumeration problem is hard. We go however further by identifying the specific structural properties of the network that account for this hardness. Those are discussed in Section 4, while Section 5 fulfills the second objective of this paper. This is to describe a new algorithm that takes advantage of such properties to obtain an exact method that is in all cases theoretically more efficient for consistent networks than the enumeration algorithms presented in [START_REF] Centler | Computing chemical organizations in biological networks[END_REF] because, at best, a smaller part of the solution space needs to be explored.

Preliminaries

A metabolic network, like any reaction system, can be modelled as a weighted directed hypergraph G = (M, R) with M the set of vertices corresponding to the metabolites and R the set of hyperarcs corresponding to the reactions. A directed hyperarc (i.e. a reaction) r ∈ R is an ordered pair of sets of vertices (i.e. metabolites) r = (subs(r), prod(r)) where subs(r) is the set of substrates of r and prod(r) is the set of products of r. For each x in subs(r) (in prod(r)) the weight of x with respect to r denotes the stoichiometric coefficient of x in r, that is, the number of units of x consumed (or produced) when r fires. Note that x can belong to both subs(r) and prod(r); in this case there are two weights associated to x w.r.t. r. Note also that, according to the above definitions, the set of substrates of a reaction r can be empty: in this case, we say that the metabolites in prod(r) are inputs of the network.

Metabolic networks have also been often modelled using matrices [START_REF] Schuster | Detection of elementary flux modes in biochemical networks: a promising tool for pathway analysis and metabolic engineering[END_REF]. The stoichiometric matrix S has |M | rows and |R| columns where S i,j is the stoichiometric coefficients of molecule i in reaction j. S i,j is negative if i is consumed and it is positive if i is produced. We notice here that while the stoichiometric matrix can always be derived from the weighted hypergraph, the reverse is not true. Indeed, metabolites involved as substrates and products of the same reaction cannot be handled in the matrix representation.

For some of the results presented, we also use the concept of the underlying graph of G, which is a directed multigraph with the same set of vertices of G and arcs x → y for every pair of vertices x, y for which there is an hyperarc r such that x ∈ subs(r) and y ∈ prod(r). A reaction is said to be on a path/cycle of the underlying directed graph if any of its (substrate,product)-pairs is an arc of the path/cycle.

In the context of metabolic networks, we say that a flux over the network is the rate at which each reaction occurs. A flux can be represented as a flux vector v ∈ R |R| with v[i] denoting the rate of reaction i. We also define a mass vector m ∈ R |M | with m[j] denoting the mass of metabolite j.

A metabolic network is flux-consistent if there exists a flux vector v > 0 such that Sv = 0 [START_REF] Acuña | Modes and cuts in metabolic networks: Complexity and algorithms[END_REF]. This is the same as saying that every reaction of the network belongs to at least one elementary mode, thus checking for the usefulness of each reaction. For more information on elementary modes, see [START_REF] Schuster | On elementary flux modes in biochemical reaction systems at steady state[END_REF] and [START_REF] Schuster | Detection of elementary flux modes in biochemical networks: a promising tool for pathway analysis and metabolic engineering[END_REF].

A metabolic network is mass-consistent if there exists a mass vector m > 0 such that m T S = 0, where m T denotes the transposed vector. This is the same as saying that there exists some mass distribution for all metabolites such that the whole set of reactions is mass balanced.

An example of a network with two reactions r 1 : a+c → b and r 2 : b+d → d+c that is not mass-consistent is shown in Figure 1. Since reaction r 2 requires that b and c have the same mass, then the mass of a in r 1 has to be 0, which is inconsistent. If we replace r 1 by r ′ 1 : a + c → 2b, then the system becomes mass-consistent (indeed, every positive mass vector with equal components is consistent). We denote by R A ⊆ R the subset of reactions that can be fired when the metabolites in set A ⊆ M are present, i.e., R A = {r ∈ R|subs(r) ⊆ A}. We now introduce the basic definitions that will be used throughout the paper.

Definition 1. A set C ⊆ M is closed if, for all reactions r ∈ R C , prod(r) ⊆ C. Moreover, given a set C ⊆ M , the closure of C, denoted by Cl C , is the smallest closed set H that contains C.
Note that if C is a closed set of molecules, then C must contain all inputs of the network (since the empty set is a subset of C and input reactions therefore belong to R C ). In particular, the closure of the empty set will contain all inputs and whatever can be produced from them.

Definition 2. A set of molecules C ⊆ M is self-maintaining if there is a flux vector v such that: 1. for all reactions r ∈ R C , v[r] > 0; 2. for all reactions r ∈ R C , v[r] = 0; 3. for all molecules i ∈ C, the production rate (Sv)[i] ≥ 0. Definition 3. A set of molecules O ⊆ M is an organisation if it is closed and self-maintaining. O is said to be reactive connected if:
-(reactive) each metabolite in O takes part as substrate or product in at least one reaction inside R O ;

-(connected) for any two molecules x and y in O, there is a path from x to y in the underlying undirected graph.

The motivation to find only reactive connected organisations is illustrated in the example of Figure 2 where the network has 3 connected components and 8 organisations but only 2 of them are reactive connected organisations. The others are just combinations of organisations which cannot directly interact among them. {{s}, {s, p, q, r}, {s, a}, {s, b, c, d}, {s, a, b, c, d}, {s, a, p, q, r}, {s, b, c, d, p, q, r}, {s, a, b, c, d, p, q, r}}, and (c) 2 reactive connected organisations: {{s}, {s, p, q, r}}.

For this reason, from now on we restrict our networks to have only one connected component and some input and output metabolites. For general networks, indeed, we can without loss of generality work on each connected component separately and then combine the results.

Since throughout the paper we need to compute closures of sets, we recall here the forward propagation procedure [START_REF] Romero | Nutrient-related analysis of pathway/genome databases[END_REF] that in an iterative process enables to obtain the closure of a given set C. Informally, this consists in starting from C itself, adding prod(r) for every r ∈ R C , and repeating this procedure until no new metabolites are added to C. The detailed description of the forward propagation algorithm is presented in the Appendix.

As already mentioned, all inputs of the network need to be considered together in order to compute organisations. This is a modelling choice that implies that if one wished to compute organisations for different subsets of the inputs, then it would be necessary to edit the network and recompute the organisations for the subsets of interest.

Chemical Organisations in Consistent Networks

It has been shown that deciding whether a network contains one organisation is NP-complete [START_REF] Centler | Computing chemical organizations in biological networks[END_REF]. However the proof was based on a network that was neither flux-nor mass-consistent. We now characterise organisations in consistent networks. First of all, we observe that it is easy to check whether a set C is an organisation by inspecting the reaction rules to check closure and selfmaintenance using linear programming.

The following theorem shows how to compute two possible organisations.

Theorem 4. If a network is flux-consistent then the whole network and the closure of the empty set are organisations.

Proof. The whole network is always closed. By definition of flux-consistency, we have a flux vector v that covers the whole network, satisfying the condition of selfmaintenance. Therefore the whole network is an organisation. Analogously, if the closure of the empty set produces the whole network then it is an organisation. Otherwise, since every metabolite is produced from the empty set, we can easily obtain a valid flux vector v satisfying the condition of self-maintenance. An algorithm to find a flux for this situation is presented in the Appendix. ⊓ ⊔

In the following, we say that the whole network and the closure of the empty set are trivial organisations.

Observe that the closure of the empty set may not always produce the whole network. An example is given in Figure 1 since the closure of the empty set for that network is {a}.

Theorem 5. If the network is flux-consistent and acyclic, i.e. the underlying directed graph of the hypergraph is acyclic, then the whole network is the only organisation.

Proof. The smallest organisation is given by the closure of the empty set, which can be obtained by applying the forward propagation algorithm to the empty set. As the network is flux-consistent and acyclic, from the inputs any metabolite can be reached, i.e., produced. Hence, the closure of the empty set is the entire network. From the flux-consistency of the network and from Theorem 4, it follows that the smallest organisation is the whole network.

⊓ ⊔

The next result shows that the problem of finding a non trivial organisation in a flux-consistent network is NP-hard. The proof (which is given in the Appendix) is based on a reduction from the 3-SAT problem, which is an appropriate modification of the original reduction given in [START_REF] Centler | Computing chemical organizations in biological networks[END_REF], that showed that finding a reactive organisation in a general reaction system is NP-hard. Theorem 6. Deciding if a flux-consistent network contains a non trivial organisation is NP-hard.

Enumerating Chemical Organisations

Theorem 6 immediately implies that it is not possible to enumerate all organisations in a flux-consistent network in polynomial-time-delay in the size of the network.

We now observe that Theorem 5 indicates that for flux-consistent networks, the difficulty of finding non trivial organisations comes from the presence of cycles in the network. Indeed, as shown in Figure 3 a but needs for this a metabolite b which is not available and therefore blocks the reaction.

In order to find the reactive connected organisations, we need to process cycles every time the forward propagation procedure stops. This simple observation gives an upper bound of 2 k on the number of reactive connected organisations in flux-consistent networks, where k is the number of cycles in the network. In order to proceed we first define cycles formally. By a cycle in the metabolic network we mean a simple directed cycle in the underlying graph. Self-loops are also considered as cycles. Definition 7. A hitting set of a set of cycles is a set of metabolites such that each cycle contains at least one element of the hitting set. Theorem 8. Let H be a hitting set of all the cycles of a directed hypergraph. The set of all reactive connected organisations, denoted as O, is such that

O ⊆ C⊆H {Cl C } Proof. It is sufficient to show that if A is a reactive connected organisation then A = Cl C , where C = A ∩ H.
First observe that, since A is closed and C is a subset of its metabolites, it follows that Cl C ⊆ A.

Let us suppose that A contains vertices which are not in Cl C . We colour these vertices white and the vertices of Cl C black. Consider any white metabolite a 1 . Since A is an organisation, a 1 cannot be vanishing. Moreover, it is not an input of the network as otherwise it would be black. Therefore, there exists a reaction r fired by A that has a 1 as product and has a white substrate a 2 , a 2 = a 1 , otherwise a 1 would again be black by closure.

By iterating the above reasoning, it follows that the subgraph of the underlying directed graph induced by white vertices has minimum in-degree at least 1 and contains a directed cycle. This contradicts the fact that H hits all the cycles. The set of white vertices is therefore empty and A = Cl C .

⊓ ⊔

The bound of the previous theorem is tight. Indeed, an example where the number of organisations reaches 2 |H| is given in Figure 4 where from the input, k metabolites are produced by k independent reactions and all of them are blocked by cycles. Any combination of these independent paths can be de-blocked and produce a new organisation, and therefore we have 2 k organisations. However, some cycles never interrupt the forward propagation procedure as illustrated in Figure 3(a). Other cycles, on the other hand, exhibit structural properties that may lead to a blocking situation. Such cycles are called potentially blocking cycles. A basic solution to find all organisations is to know how to unblock all cycles of the network independently of whether they are potentially blocking cycles or not. Therefore, instead of finding a hitting set for all cycles of the network, it is enough to break all potentially blocking cycles to compute all reactive connected organisations. In order to prove this, we first introduce a more formal definition of potentially blocking cycle. Definition 9. A potentially blocking cycle is a cycle such that there exists a reaction r = ({s 1 , . . . , s h } , {p 1 , . . . , p k }) in the network satisfying the following two conditions: [START_REF] Acuña | Modes and cuts in metabolic networks: Complexity and algorithms[END_REF] there exist i and j such that (s i , p j ) is an edge of the cycle, and (2) there exists ℓ such that s ℓ is not in the cycle.

A potentially blocking cycle may or may not interrupt the forward propagation depending on the metabolites that were produced by the procedure once the cycle is reached. Figure 5(a) shows an example in which the cycle will be traversed, while Figure 5(b) shows an example in which it will block the forward propagation algorithm. Let a i (i = 1, 2, . . . , k) be the set of white vertices such that there exists a reaction r i having a i as product and at least one black substrate. Then r i has also at least one white substrate, which we denote by w i , as otherwise a i would be black by closure. Note that a white vertex does not have to belong to the set of the a i 's, but that this set is not empty as the organisation is connected and the set of white vertices is not empty.

If w i = a i , the selfloop induced by r i is a potentially blocking cycle that contains no vertex of H, leading to a contradiction. Thus we may assume that w i is distinct from a i .

For 1 ≤ i ≤ k, define T i as follows: a vertex w is in T i if either w = w i or there exists a white path starting from w and ending in w i . Up to a reordering, we may assume that

|T 1 | ≤ |T i | for 2 ≤ i ≤ k.
As T 1 is not empty, it has to contain at least one vertex that is a product of a reaction having a black substrate. If that vertex is a 1 , there exists a path from a 1 to w 1 , which yields a white cycle with the edge (w 1 , a 1 ). That cycle is a potentially blocking cycle (because of the reaction r 1 ) which contains no vertex in H, leading to a contradiction.

In other words, T 1 contains a vertex among (a 2 , . . . , a k ), say a 2 . This implies that every path ending in w 2 can be extended to a path ending in w 1 and thus

T 2 ⊆ T 1 . Therefore, by minimality of T 1 , T 1 = T 2 .
This implies that w 1 ∈ T 2 : hence, there exists a white path from w 1 to w 2 . As a 2 ∈ T 1 , there also exists a white path from a 2 to w 1 . Thus, we can construct a white path from a 2 to w 2 . Considering the edge (w 2 , a 2 ), we again obtain a white cycle, which is potentially blocking (because of the reaction r 2 ) and contains no vertex in H, leading to a contradiction.

Thus, the set of white vertices is empty and

A ⊆ Cl C . ⊓ ⊔
Even in the case of the previous theorem, the bound is tight. Indeed, an example where the number of organisations is 2 |H| is, once again, given in Figure 4, since all cycles presented in the example are potentially blocking.

Hitting Set Approach to Enumerate Organisations

Two exact algorithms were proposed in [START_REF] Centler | Computing chemical organizations in biological networks[END_REF] to enumerate organisations. The first one consisted in enumerating all closed sets and then checking for their selfmaintenance, while the second one consisted in enumerating all self-maintaining sets and then checking linear combinations of them in order to obtain closed sets. A third approach was also proposed that was based on the second one but avoided enumerating all self-maintaining sets. This algorithm however was a heuristic not guaranteed to find all self-maintaining sets (and thus organisations). Finally, a variation of the first algorithm was proposed in order to enumerate only reactive connected organisations. This algorithm comes closer to the one we describe later and works as follows. First, the forward propagation of the empty set is computed. Once the procedure is blocked, all possible combinations of metabolites that are connected to the produced set X are considered for addition, in order to obtain further closed sets that include X. At this point, the algorithm recursively continues.

Note that in the above procedures, no concept of blocking cycles has been formally identified and used. Now that we know that the hardness comes from such cycles, two different approaches can be applied in flux-consistent networks. One is to find a global hitting set for all cycles of the network and then, following Theorem 8, to apply the forward propagation procedure on each subset of the hitting set to produce closed sets which together form all candidate organisations and, finally, to check through LP if the candidates are self-maintaining. However, the problem of finding a minimum hitting set for all cycles of a directed graph is NP-hard as indeed it corresponds to the feedback vertex sex (FVS) problem [START_REF] Karp | Reducibility among combinatorial problems[END_REF]. Nevertheless approximation algorithms such as the one described in [START_REF] Seymour | Packing directed circuits fractionally[END_REF] can be used in order to perform this step.

A second possibility is to find a local hitting set. According to Theorem 10 only potentially blocking cycles need to be considered. This is a superset of the blocking cycles that can be identified when the forward propagation procedure stops because it is at this moment that we know we are dealing with actually blocking cycles. A more efficient algorithm to enumerate reactive connected organisations is thus the following one: apply the forward propagation algorithm and once blocked, identify the set B of metabolites that are blocking the closure and find a hitting set that unblocks only the cycles which directly or indirectly involve these blocking metabolites.

In [START_REF] Even | Approximating minimum feedback sets and multicuts in directed graphs[END_REF], the authors presented an approximation algorithm to a generalisation of the FVS problem, called SUBSET-FVS, in which only a subset of the directed cycles in the graph is considered interesting, more specifically the ones that intersects a set of special vertices. In our case, the set of special vertices would be the blocking metabolites locally identified as described in the previous paragraph. The authors in [START_REF] Even | Approximating minimum feedback sets and multicuts in directed graphs[END_REF] gave two approximation algorithms for the SUBSET-FVS problem. The first algorithm achieves an approximation factor of O log 2 |B| . The second achieves an approximation factor of O (min{log T log log T, log n log log n}), where T is the value of the optimum fractional solution of the problem at hand, and n is the number of vertices in the graph.

Before proving that this idea can be correctly used to exactly solve our problem, we need to define the concept of a blocking cycle in relation to a given set C of metabolites. Definition 11. Let C be a set of metabolites. A C-blocking cycle is a cycle of vertices which are not in C such that there exists a reaction r in the cycle whose set of substrates contains at least one metabolite in C.

C-blocking cycles correspond to those which actually stop the forward propagation procedure.

Theorem 12. Let C be a closed set and H a hitting set of the C-blocking cycles of a metabolic network. Let A be a reactive connected organisation whose metabolites contain C. Then either C = A or there exists a non empty subset B of H such that the closure of C ∪ B is still a subset of A.

Proof. Let A be a reactive connected organisation containing C as a subset of its metabolites and let B = A ∩ H. Let us suppose that C = A. To prove the theorem, it is then sufficient to prove that B = A ∩ H is not empty.

We colour the vertices of A \ C in white and those of C in black. Since A is a reactive connected organisation, there exist edges between the white and the black metabolites and some of them go from a black to a white vertex, as otherwise, white vertices would be vanishing. Let (a 1 , . . . , a k ) be the set of white vertices reached by at least one edge coming from a black vertex.

The same argument of the proof of Theorem 10 can now be applied, showing that, as the set of white vertices is not empty, it contains a white C-blocking cycle. Therefore, B is not empty.

⊓ Notice that the size of the hitting set computed by the algorithm is never greater than the number of blocking metabolites. Thus we can guarantee that our algorithm is theoretically better than existing algorithms which consider all blocking metabolites and then test all subsets. In the Appendix, we provide an example in which finding the hitting set gives a better result.

Conclusion

All the results presented in this paper correspond to enumerating closed sets as potential organisations. The problem of enumerating self-maintaining sets is still open. Such an approach could enable us to design a method that enumerates stoichiometrically valid precursor sets, which would be an important follow-up of the work on minimal precursors sets presented in [START_REF] Cottret | Enumerating Precursor Sets of Target Metabolites in a Metabolic Network[END_REF]. Finally, the algorithms introduced in this paper do not take any specific advantage of the fact that the network should be mass-consistent and exploiting this might lead to better algorithms for enumerating self-maintaining sets. from lines 5-8 detects the molecules that are vanishing and increases the fluxes of some reactions that produce them in order to have a self-maintaining set. This is repeated until no molecule is vanishing anymore. The algorithm works because we choose the vanishing molecule to fix based on a previous DFS labelling. This allows us to pick a molecule y knowing that every other molecule that can be reached from y is not vanishing (has a finishing timestamp greater than the one of y). To produce y, we only consume molecules that have finishing timestamp greater than y. At each step, the new vanishing molecule must have a greater finishing time than the previous one and this guarantees that the algorithm will finish and produce a flux vector v in which no molecule vanishes.

⊓ ⊔ Proof of Theorem 6 Theorem 6. Deciding if a flux-consistent network contains a non trivial organisation is NP-hard.

Proof. We reduce our problem from 3-SAT. Given a boolean formula F in 3-CNF with n boolean variables and ℓ clauses, we construct a flux-and mass-consistent reaction network for which the existence of a non-trivial reactive (connected) organisation implies a positive answer to the 3-SAT and vice-versa.

In the reaction network we define the 2n literal metabolites x 1 , x 2 , ..., x n , ¬x 1 , ¬x 2 , ..., ¬x n (¬x is the negation of x), clause metabolites C 1 , C 2 , ..., C ℓ and additional metabolites key, a, input. The hyperarcs (reactions) are defined as:

-Influx: A reaction ∅ → input -Clause reactions: For each clause C h = (A 1 ∨ A 2 ∨ A 3 ), we introduce 3 reactions: key + A j → C h + A j , j = 1, 2, 3. -Key reaction: C 1 + C 2 + ... + C l + input → (l + 1) key.

-Key outflux: A reaction key → ∅ -Variable reactions of two types: For each variable x j , we introduce a type 1 reaction x j + ¬x j → a, and a type 2 reaction a + input → x j + ¬x j + key, j = 1, . . . , n.

An illustration of this transformation is given in Figure 6. The transformation can obviously be done in polynomial time.

Let us first prove that the network is mass-and flux-consistent. To prove the latter, we construct a flux vector v > 0 such that Sv = 0 [START_REF] Acuña | Modes and cuts in metabolic networks: Complexity and algorithms[END_REF]. The balanced flux distribution v constructed has positive flux on every reaction of the network.

Mass-consistency of the network is easily verified by the mass vector that has mass 1 for all metabolites except for a, for which the mass is 2.

Suppose we have a satisfying truth assignment to F . Let K be the set of literals that are true, then we claim that the set containing all metabolites corresponding to K, all clauses, input and key is a non-trivial organisation. Closedness is easily verified. We now show self-maintenance. Let k i be the number of true literals in C i , i = 1, . . . , ℓ. For each true literal in C i we set a flux of 1/k i on the corresponding clause reaction. Together with a flux of at least 1 in the key reaction, the influx and outflux secure self-maintenance.

Reversely, suppose we have a non-trivial reactive (connected) organisation O. First, notice that {input} (the closure of ∅) is a trivial organisation and that it is part of every organisation, hence input ∈ O. As we now explain, this implies that if any organisation contains both x i and ¬x i for any i, then it must contain the whole network, and hence be trivial. The presence of x i and ¬x i triggers a variable type 1 reaction, producing a, which on its turn triggers all variable type 2 reactions producing all literal metabolites and key. All clause reactions are therefore triggered and henceforth the key reaction and the outflux reaction. Hence, O does not contain both literals corresponding to the same variable, and none of the variable reactions will have positive flux.

Since O is reactive, input must be the substrate of some active reaction, hence the key reaction must be active, implying that key ∈ O. Because of self-maintenance, all clause metabolites must be produced. To produce a clause metabolite at least one of its three clause reactions must be fired. Let K be the set of literal metabolites that are part of our reactive (connected) organisation. By closure, they are able to produce all clause metabolites, hence setting their value to true yields a satisfying truth assignment for F .

⊓ ⊔

Local Hitting Set Approach

An illustration of the local hitting set approach is given in Figure 7. In this example, the forward propagation procedure produces the set C and is blocked by two metabolites, w 1 and w 2 . The goal is to find a hitting set that intersects the blocking cycles for these metabolites. In this case, the hitting set has size 1 which is better than checking the combination of C with each of the blocking metabolites.

Fig. 1 .

 1 Fig. 1. Metabolic network whose mass-consistence depends on the stoichiometric coefficients of the first reaction.

Fig. 2 .

 2 Fig.2. A metabolic network with (a) 3 connected components, (b) 8 organisations: {{s}, {s, p, q, r}, {s, a}, {s, b, c, d}, {s, a, b, c, d}, {s, a, p, q, r}, {s, b, c, d, p, q, r}, {s, a, b, c, d, p, q, r}}, and (c) 2 reactive connected organisations: {{s}, {s, p, q, r}}.

  (b), cycles may interrupt the forward propagation if there exists a reaction that can produce a new metabolite

Fig. 3 .

 3 Fig. 3. (a) Non blocking cycle that will be traversed by the forward propagation procedure. (b) Forward propagation blocked by a unreached metabolite b.

Fig. 4 .

 4 Fig. 4. Example with k parallel blocking cycles and 2 k organisations.

Fig. 5 .

 5 Fig. 5. Example of a potentially blocking cycle formed by the vertices a and c and reactions a + b → c and c → a + b. This cycle will be traversed (a) or will block the forward propagation procedure (b) depending on how it is reached by procedure itself.

1 .Fig. 6 .

 16 Fig.6. Finding a non trivial reactive (connected) organisation is NP-hard. The key and input metabolites are presented more than once in the figure for convenience. Edge labels do not stand for stoichiometries but fluxes for a valid flux distribution that covers the whole network.

Fig. 7 .

 7 Fig. 7. Local hitting set approach.

  As in the proof of Theorem 8, it is sufficient to show that, given a reactive connected organisation A, A = Cl C , where C = A∩H (in the following, all paths and cycles are meant directed in the underlying directed graph). Once again, it is easy to see that Cl C ⊆ A. Let us then suppose that A contains vertices which are not in Cl C and let us colour them white and those of Cl C black.

	Theorem 10. Let H be a hitting set of all the potentially blocking cycles of a
	directed hypergraph. The set of all reactive connected organisations, denoted as
	O, is such that
	O ⊆	{Cl C }
	C⊆H
	Proof.

⊔

  Corollary 13. Every reactive connected organisation is included in the set CO returned by the procedure given in Algorithm CCO.Proof. Let A be a reactive connected organisation. It has to contain C 0 as every organisation contains the closure of the empty set. Let C be maximum among the elements of CO which are subsets of A. Then Theorem 12 implies, by maximality of C, that A = C. ⊓ ⊔ Ensure: the set CO of all candidates for being organisations. CO ← {C0} where C0 is the closure of the empty set (Cl {} ) for all elements C in CO which have not been treated before do Compute a hitting set H of the C-blocking cycles for every B ⊂ H do Compute Cl C∪B and add it to CO if it was not present already return CO

	Algorithm CCO(G)
	Require: a metabolic network represented as a hypergraph G = (M, R);

Forward Propagation Algorithm

The forward propagation algorithm computes the closure of a given initial set of metabolites.

Algorithm FP(G, C) Require: a metabolic network represented as a hypergraph G = (M, R), a set C ⊂ M of available metabolites; Ensure: a set X of metabolites that can be reached from C.

Theorem 4. If a network is flux-consistent then the whole network and the closure of the empty set are organisations.

Proof. In order to proceed, we define Sv[x] as the total amount of molecules x that are produced (Sv[x] ≥ 0) or consumed ((Sv[x] < 0) for the flux distribution v. We now present an algorithm to compute the flux vector that fulfills selfmaintenance for the closure of the empty set.

Algorithm FluxForTrivialOrganisation(G) Require: a metabolic network represented as an hypergraph G = (M, R) Ensure: a flux vector v that fulfills the self-maintenance condition for the closure of the empty set.

X ← F P ({})

Apply DFS in the underlying graph and order the vertices by finishing timestamp { ft[v] is the finishing timestamp for a vertex v} v[r] ← 1, for all r ∈ R X while exists some molecule x with Sv[x] < 0 do

x ← molecule with Sv[x] < 0 and minimum finishing timestamp RX ← reactions r in R X with x ∈ prod(r) and all substrates s ∈ subs(r)

In a self-maintaining set nothing can vanish and all reactions that can be fired must be fired. The algorithm finds the closure of the empty set and names it X. We then initialise a flux vector v that fires all reactions in R X . Some of the compounds probably will be vanishing with this flux distribution. The loop