Clémence Frioux

Torsten Schaub

Sebastian Schellhorn

Anne Siegel

Under consideration for publication in Theory and Practice of Logic Programming 1 Hybrid Metabolic Network Completion

Metabolic networks play a crucial role in biology since they capture all chemical reactions in an organism. While there are networks of high quality for many model organisms, networks for less studied organisms are often of poor quality and suffer from incompleteness. To this end, we introduced in previous work an ASP-based approach to metabolic network completion. Although this qualitative approach allows for restoring moderately degraded networks, it fails to restore highly degraded ones. This is because it ignores quantitative constraints capturing reaction rates. To address this problem, we propose a hybrid approach to metabolic network completion that integrates our qualitative ASP approach with quantitative means for capturing reaction rates. We begin by formally reconciling existing stoichiometric and topological approaches to network completion in a unified formalism. With it, we develop a hybrid ASP encoding and rely upon the theory reasoning capacities of the ASP system clingo for solving the resulting logic program with linear constraints over reals. We empirically evaluate our approach by means of the metabolic network of Escherichia coli. Our analysis shows that our novel approach yields greatly superior results than obtainable from purely qualitative or quantitative approaches. Under consideration in Theory and Practice of Logic Programming (TPLP).

1 Introduction Among all biological processes occurring in a cell, metabolic networks are in charge of transforming input nutrients into both energy and output nutrients necessary for the functioning of other cells. In other words, they capture all chemical reactions occurring in an organism. In biology, such networks are crucial from a fundamental and technological point of view to estimate and control the capability of organisms to produce certain products. Metabolic networks of high quality exist for many model organisms. In addition, recent technological advances enable their semi-automatic generation for many less studied organisms, also described as non-model organisms. However, the resulting metabolic networks are usually of poor quality, due to error-prone, genome-based construction processes and a lack of (human) resources. As a consequence, they usually suffer from substantial incompleteness. The common fix is to fill the gaps by completing a draft network by borrowing chemical pathways from reference networks of well studied organisms until the augmented network provides the measured functionality.

In previous work [START_REF] Schaub | Metabolic network expansion with ASP[END_REF], we introduced a logical approach to metabolic network completion by drawing on the work in [START_REF] Handorf | Expanding metabolic networks: Scopes of compounds, robustness, and evolution[END_REF]. We formulated the problem as a qualitative combinatorial (optimization) problem and solved it with Answer Set Programming (ASP [START_REF] Baral | Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF]. The basic idea is that reactions apply only if all their reactants are available, either as nutrients or provided by other metabolic reactions. Starting from given nutrients, referred to as seeds, this allows for extending a metabolic network by successively adding operable reactions and their products. The set of compounds in the resulting network is called the scope of the seeds and represents all compounds that can principally be synthesized from the seeds. In metabolic network completion, we query a database of metabolic reactions looking for (minimal) sets of reactions that can restore an observed bio-synthetic behavior. This is usually expressed by requiring that certain target compounds are in the scope of some given seeds. For instance, in the follow-up work in [START_REF] Collet | Extending the metabolic network of Ectocarpus siliculosus using answer set programming[END_REF][START_REF] Prigent | The genome-scale metabolic network of ectocarpus siliculosus (ectogem): a resource to study brown algal physiology and beyond[END_REF], we successfully applied our ASP-based approach to the reconstruction of the metabolic network of the macro-algae Ectocarpus siliculosus, using the collection of reference networks Metacyc [START_REF] Caspi | The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases[END_REF].

We evidenced in [START_REF] Prigent | Meneco, a Topology-Based Gap-Filling Tool Applicable to Degraded Genome-Wide Metabolic Networks[END_REF] that our ASP-based method partly restores the bio-synthetic capabilities of a large proportion of moderately degraded networks: it fails to restore the ones of both some moderately degraded and most of highly degraded metabolic networks. The main reason for this is that our purely qualitative approach misses quantitative constraints accounting for the law of mass conservation, a major hypothesis about metabolic networks. This law stipulates that each internal metabolite of a network must balance its production rate with its consumption rate at the steady state of the system. Such rates are given by the weighted sums of all reaction rates consuming or producing a metabolite, respectively. This calculation is captured by the stoichiometry 1 of the involved reactions. Hence, the qualitative ASP-based approach fails to tell apart solution candidates with correct and incorrect stoichiometry and therefore reports inaccurate results for some degraded networks.

We address this by proposing a hybrid approach to metabolic network completion that integrates our qualitative ASP approach with quantitative techniques from Flux Balance Analysis (FBA 2 [START_REF] Maranas | Optimization methods in metabolic networks[END_REF], the state-of-the-art quantitative approach for capturing reaction rates in metabolic networks. We accomplish this by taking advantage of recently developed theory reasoning capacities for the ASP system clingo [START_REF] Gebser | Theory solving made easy with clingo 5[END_REF]. More precisely, we use an extension of clingo with linear constraints over reals, as dealt with in Linear Programming (LP [START_REF] Dantzig | Linear Programming and Extensions[END_REF]. This extension provides us with an extended ASP modeling language as well as a generic interface to alternative LP solvers, viz. cplex and lpsolve, for dealing with linear constraints. We empirically evaluate our approach by means of the metabolic network of Escherichia coli. Our analysis shows that our novel approach yields superior results than obtainable from purely qualitative or quantitative approaches. Moreover, our hybrid application provides a first evaluation of the theory extensions of the ASP system clingo with linear constraints over reals in a non-trivial setting.

Metabolic Network Completion

Metabolism is the sum of all chemical reactions occurring within an organism. As the products of a reaction may be reused as reactants, reactions can be chained to complex chemical pathways. Such complex pathways are described by a metabolic network.

We represent a metabolic network as a labeled directed bipartite graph G = (R ∪ M, E, s), where R and M are sets of nodes standing for reactions and compounds (also called metabolites), respectively. When (m, r) ∈ E or (r, m) ∈ E for m ∈ M and r ∈ R, the metabolite m is called a reactant or product of reaction r, respectively. Metabolites and reactions nodes can both have multiple ingoing and outgoing edges. More formally, for any r ∈ R, define rcts(r

) = {m ∈ M | (m, r) ∈ E} and prds(r) = {m ∈ M | (r, m) ∈ E}.
The edge labeling s : E → R gives the stoichiometric coefficients of a reaction's reactants and products, respectively, i.e., their relative quantities involved in the reaction. Finally, the activity rate of reactions is bound by lower and upper bounds, denoted by lb r ∈ R + 0 and ub r ∈ R + 0 for r ∈ R, respectively. Whenever clear from the context, we refer to metabolic networks with G (or G ′ , etc) and denote the associated reactions and compounds with M and R (or M ′ , R ′ etc), respectively.

We distinguish a set S ⊆ M of compounds as initiation seeds, that is, compounds initially present due to experimental evidence. Another set of compounds is assumed to be activated by default. These boundary compounds are defined as:

S b (G) = {m ∈ M | r ∈ R, m ∈ prds(r), rcts(r) = ∅}.
For simplicity, we assume that all boundary compounds are seeds: S b (G) ⊆ S. Note that follow-up concepts like reachability and activity in network completion are independent of this assumption.

For illustration, consider the metabolic network in Fig. 1. The network consists of 9 reactions, r s1 , r s2 , r e and r 0 to r 5 , and 8 compounds, A, . . . , F , S 1 , S 2 and S 3 . Here, S = {S 1 , S 2 , S 3 }, S 1 and S 2 being the two boundary compounds of the network. Dashed rectangle describes the boundary of the system, outside of which is the environment of the organism. Consider reaction r 4 : E → 2C transforming one unit of E into two units of C (stoichiometric coefficients of 1 are omitted in the graphical representation; cf. Fig. 1). We have rcts(r 4 ) = {E}, prds(r 4 ) = {C}, along with s(E, r 4 ) = 1 and s(r 4 , C) = 2.

In biology, several concepts have been introduced to model the activation of reaction fluxes in metabolic networks, or to synthesize metabolic compounds. To model this, we introduce a function active that given a metabolic network G takes a set of seeds S ⊆ M and returns a set of activated reactions active G (S) ⊆ R. With it, metabolic network completion is about ensuring that a set of target reactions (reaction r 5 in Fig. 1) is activated from seed compounds in S by possibly extending the metabolic network with reactions from a reference network (cf. shaded part in Fig. 2).

Formally, given a metabolic network G = (R∪M, E, s), a set S ⊆ M of seed compounds Fig. 2: Metabolic network completion problem. The purpose of its solving is to select the minimal number of reactions from a database (dashed shaded reactions) such that activation of target reaction r 5 is restored from boundary and/or initiation seeds. There are three formalisms for activation of target reaction: stoichiometric, topological and hybrid.

such that S b (G) ⊆ S, a set R T ⊆ R of target reactions, and a reference network (R ′ ∪ M ′ , E ′ , s ′ ), the metabolic network completion problem is to find a set R ′′ ⊆ R ′ \ R of reactions of minimal size such that R T ⊆ active G ′′ (S) where 3 G ′′ = ((R ∪ R ′′ ) ∪ (M ∪ M ′′ ), E ∪ E ′′ , s ′′ ) , (1)
M ′′ = {m ∈ M ′ | r ∈ R ′′ , m ∈ rcts(r) ∪ prds(r)} , (2)
E ′′ = E ′ ∩ ((M ′′ × R ′′ ) ∪ (R ′′ × M ′′ )), and (3)
s ′′ = s ∪ s ′ . (4)
We call R ′′ a completion of (R∪M, E, s) from (R ′ ∪M ′ , E ′ , s ′ ) wrt S and R T . Our concept of activation allows different biological paradigms to be captured. Accordingly, different formulations of metabolic network completion can be characterized: the stoichiometric, the relaxed stoichiometric, the topological, and the hybrid one. We elaborate upon their formal characterizations in the following sections.

Stoichiometric Metabolic Network Completion

The first activation semantics has been introduced in the context of Flux Balance Analysis capturing reaction flux distributions of metabolic networks at steady state. In this paradigm, each reaction r is associated with a metabolic flux value, expressed as a real variable v r confined by the minimum and maximum rates:

lb r ≤ v r ≤ ub r for r ∈ R. (5)
Flux distributions are formalized in terms of a system of equations relying on the stoichiometric coefficients of reactions. Reaction stoichiometries are governed by the law of mass conservation under a steady state assumption; in other words, the mass of the system remains constant over the reaction. The input and output fluxes of reactions consuming and producing a metabolite are balanced.

r∈R s(r, m) • v r + r∈R -s(m, r) • v r = 0 for m ∈ M. (6)
Given a target reaction r T ∈ R T , a metabolic network G = (R ∪ M, E, s) and a set of seeds S, stoichiometric activation is defined as follows:

r T ∈ active s G (S) iff v rT > 0 and ( 5) and ( 6) hold for M and R.

Note that the condition v rT > 0 strengthens the flux condition for r T ∈ R in the second part. More generally, observe that activated target reactions are not directly related to the network's seeds S. However, the activation of targets highly depends on the boundary compounds in S b (G) for which ( 6) is always satisfied and thus initiates the fluxes. Since boundary compounds are produced by at least one reaction without prerequisite, an arbitrary amount might be produced. Therefore, the incoming flux value always balances the sum of the flux values associated to outgoing edges. Intuitively, boundary compounds are nutrients that are expected to be available in the system for the consumption by the metabolic network, thus initiating the reactions within. In our draft network G, consisting of all non-dashed nodes and edges depicted in Fig. 2 (viz. reactions r s1 , r s2 , r e and r 0 to r 5 and compounds A, . . . , F , S 1 , S 2 , and S 3 and r 5 the single target reaction) and the reference network G ′ , consisting of the shaded part of Fig 2, (viz. reactions r 6 to r 9 and metabolite G) a strict stoichiometry-based completion aims to obtain a solution with r 5 ∈ active s G ′′ ({S 1 , S 2 , S 3 }) where v r5 is maximal. This can be achieved by adding the completion R ′′ 1 = {r 6 , r 9 } (Fig. 3). The cycle made of compounds E, C, D and the boundary seed S 2 is already balanced and notably self-activated. Indeed, initiation of D and E producibility requires the producibility of C (in addition to the presence of the boundary seed S 2 ) that itself depends on D and E. Yet, according the flux conditions, that models steady state conditions, the cycle is activated. Such self-activation of cyclic pathways is an inherent problem of purely stoichiometric approaches to network completion. This is a drawback of the semantics because the effective activation of the cycle requires the additional (and unchecked) condition that at least one of the compounds was present as the initial state of the system. This could be the case provided there exist another way to enable the production of one or several components of the cycle (here an activable Fig. 3: Solution to metabolic network completion under stoichiometric activation hypothesis in order to satisfy Equations ( 5), ( 6) and ( 7). Within this network, there exists at least one flux distribution which activates r 5 . Fig. 4: Solution to metabolic network completion under relaxed stoichiometric activation hypothesis in order to satisfy Equations ( 5), ( 8) and ( 9). Notice that within this completed network, there exist no flux distribution allowing the reaction r 5 to be activated.

reaction producing E for instance) [START_REF] Prigent | Meneco, a Topology-Based Gap-Filling Tool Applicable to Degraded Genome-Wide Metabolic Networks[END_REF]. The instance of Equation ( 6) controlling the reaction rates related to metabolite

C is 2 • v r4 -v r2 -v r5 = 0.
To solve metabolic network completion with flux-balance activated reactions, Linear Programming can be used to maximize the flux rate v rT provided that the linear constraints are satisfied. Nonetheless, this problem turns out to be hard to solve in practice and existing approaches scale poorly to real-life applications (cf. [START_REF] Orth | Systematizing the generation of missing metabolic knowledge[END_REF].

This motivated the use of approximate methods. The relaxed problem is obtained by weakening the mass-balance equation ( 6) as follows:

r∈R s(r, m) • v r + r∈R -s(m, r) • v r ≥ 0 for m ∈ M. (8)
This lets us define the concept of relaxed stoichiometric activation:

r T ∈ active r G (S) iff v rT > 0 and ( 5) and ( 8) hold for M and R.

The resulting problem can now be efficiently solved with Linear Programming [START_REF] Satish Kumar | Optimization based automated curation of metabolic reconstructions[END_REF]. Existing systems addressing strict stoichiometric network completion either cannot guarantee optimal solutions [START_REF] Latendresse | Efficiently gap-filling reaction networks[END_REF] or do not support a focus on specific target reactions [START_REF] Thiele | fastGapFill: efficient gap filling in metabolic networks[END_REF]. Other approaches either partially relax the problem [START_REF] Vitkin | MIRAGE: a functional genomics-based approach for metabolic network model reconstruction and its application to cyanobacteria networks[END_REF] or solve the relaxed problem based on Equation ( 8), like the popular system gapfill (Satish [START_REF] Satish Kumar | Optimization based automated curation of metabolic reconstructions[END_REF]. Applied to the network of Fig. 2, the minimal completion under the relaxed stoichiometric activation is R ′′ 1 = {r 6 } (Fig. 4) but does not carry flux because of the accumulation of metabolite G, allowed by Equation (8). Note however that for strict steady-state modeling an a posteriori verification of solutions is needed to warrant the exact mass-balance equation ( 6).

Topological Metabolic Network Completion

A qualitative approach to metabolic network completion relies on the topology of networks for capturing the activation of reactions. Given a metabolic network G, a reaction r ∈ R is activated from a set of seeds S if all reactants in rcts(r) are reachable from S. Moreover, a metabolite m ∈ M is reachable from S if m ∈ S or if m ∈ prds(r) for some reaction r ∈ R where all m ′ ∈ rcts(r) are reachable from S. The scope of S, written Σ G (S), is the closure of compounds reachable from S. In this setting, topological activation of reactions from a set of seeds S is defined as follows:

r T ∈ active t G (S) iff rcts(r T ) ⊆ Σ G (S). ( 10
)
Note that this semantics avoids self-activated cycles by imposing an external entry sufficient to initiate all cycles (S 3 is not enough to activate the cycle as it does not activate one of its reaction on its own). The resulting network completion problem can be expressed as a combinatorial optimization problem and effectively solved with ASP [START_REF] Schaub | Metabolic network expansion with ASP[END_REF].

For illustration, consider again the draft and reference networks G and G ′ in Fig. 1 and Fig. 2. We get Σ G ({S 1 , S 2 , S 3 }) = {S 1 , S 2 , S 3 , B}, indicating that target reaction r 5 is not activated from the seeds with the draft network because A and C, its reactants, are not reachable. This changes once the network is completed. Valid minimal completions are R ′′ 2 = {r 6 , r 7 } (Fig. 5) and

R ′′ 3 = {r 6 , r 8 } (Fig. 6) because r 5 ∈ active t G ′′ i ({S 1 , S 2 }) since {A, C} ⊆ Σ G ′′ i ({S 1 , S 2 }) for all extended networks G ′′ i obtained from completions R ′′ i of G for i ∈ {2, 3}.
Relevant elements from the reference network are given in dashed gray.

Hybrid Metabolic Network Completion

The idea of hybrid metabolic network completion is to combine the two previous activation semantics: the topological one accounts for a well-founded initiation of the system from the seeds and the stoichiometric one warrants its mass-balance. We thus aim at network completions that are both topologically functional and flux balanced (without suffering from self-activated cycles). More precisely, a reaction r T ∈ R T is hybridly activated from a set S of seeds in a network G, if both criteria apply:

r T ∈ active h G (S) iff r T ∈ active s G (S) and r T ∈ active t G (S). (11)
Applying this to our example in Fig. 2, we get the (minimal) hybrid solutions 11) (that is Equations ( 5), ( 6), ( 7) and ( 10)). 11) (that is Equations ( 5), ( 6), ( 7) and ( 10)).

R ′′ 4 = 2 S 3 S 2 S 1 F A B C D E G r 0 r 1 r 2 r 3
{r 6 , r 7 , r 9 } (Fig. 7) and R ′′ 5 = {r 6 , r 8 , r 9 } (Fig. 8). Both (topologically) initiate paths of reactions from the seeds to the target, ie.

r 5 ∈ active t G ′′ i ({S 1 , S 2 , S 3 }) since {A, C} ⊆ Σ G ′′ i ({S 1 , S 2 , S 3 }) for both extended networks G ′′ i obtained from completions R ′′ i of G for i ∈ {4, 5}.
Both solutions are as well stoichiometrically valid and balance the amount of every metabolite, hence we also have r 5 ∈ active s

G ′′ i ({S 1 , S 2 , S 3 }).

Union of Metabolic Network Completions

As depicted in the toy examples for the topological (Fig. 5 and Fig. 6) and hybrid (Fig. 7 and Fig. 8) activation, several minimal solutions to one metabolic network completion problem may exist. There might be dozens of minimal completions, depending on the degradation of the original draft network, hence leading to difficulties for biologists and bioinformaticians to discriminate the individual results. One solution to facilitate this curation task is to provide, in addition to the enumeration of solutions, their union. This has been done previously for the topological completion [START_REF] Prigent | Meneco, a Topology-Based Gap-Filling Tool Applicable to Degraded Genome-Wide Metabolic Networks[END_REF].

Notably, the concept of "union of solutions" is particularly relevant from the biological perspective since it provides in a single view all possible reactions that could be inserted in a solution to the network completion problem. Additionally, verifying the union according to the desired (stoichiometric and hybrid) activation semantics, offers a way to analyze the quality of approximation methods (topological and relaxed-stoichiometric ones). If individual solutions contradict a definition of activation that the union satisfies, it suggests that the family of reactions contained in the union, although possibly nonminimal, may be of interest. Thus providing merit to the approximation method and their results.

Importantly, we notice that the operation of performing the union of solutions is stable with the concept of activation, although it can contradict the minimality of the size of completion. Indeed, the union of solutions to the topological network completion problem is itself a (non-minimal) solution to the topological completion problem. Similarly, the union of minimal stoichiometric solutions always displays the stoichiometric activation of the target reaction(s). In fact, adding an arbitrary set of reactions to a metabolic network still maintains stoichiometric activation, since flux distribution for the newly added reactions may be set to zero. Consequently, the union of minimal hybrid solutions always displays the hybrid activation in the target reaction(s).

The following theorems (Theorems 1, 2 and 3) are a formalization of the stability of the union of solutions with respect to the three concepts of activation.

The union

G = G 1 ∪ G 2 of two metabolic networks G 1 = (R 1 ∪ M 1 , E 1 , s 1 ) and G 2 = (R 2 ∪ M 2 , E 2 , s 2 ) is defined by G = (R ∪ M, E, s), (12)
R = R 1 ∪ R 2 , (13)
M = M 1 ∪ M 2 , (14)
E = E 1 ∪ E 2 , (15)
s = s 1 ∪ s 2 . ( 16
) Theorem 1. Let G 1 and G 2 be metabolic networks. If R T ⊆ active t G1 (S), then R T ⊆ active t G1∪G2 (S).
Proof. The proof is given by monotonicity of the union and the monotonicity of the closure. Thus it can never be case that having more reactions disables reachability. More formal, R T ⊆ active t G1 (S) holds iff rcts(r T ) ⊆ Σ G1 (S). Furthermore, we have Σ G1 (S) ⊆ Σ G1∪G2 (S) by the definition of the closure. This implies rcts(r T ) ⊆ Σ G1∪G2 (S). Finally, we have R T ⊆ active t G1∪G2 (S).

Theorem 2. Let G 1 and G 2 be metabolic networks. If R T ⊆ active s G1 (S), then R T ⊆ active s G1∪G2 (S).

Proof. First, we define following bijective functions

f :R 1 → {1, . . . , l} ⊆ N, r → f (r) = i g :M 1 → {1, . . . , k} ⊆ N, m → g(m) = j f ′ :R 1 ∪ R 2 → {1, . . . , l ′ } ⊆ N, r → f ′ (r) = f (r) , if f (r) is defined i , otherwise g ′ : M 1 ∪ M 2 → {1, . . . , k ′ } ⊆ N m → g ′ (m) = g(m) , if g(m) is defined j , otherwise for k = |M 1 |, l = |R 1 |, k ′ = |M 1 ∪ M 2 | and l ′ = |R 1 ∪ R 2 | regarding G 1 and G 1 ∪ G 2 ,
respectively. Now, we rewrite the system of ( 6) regarding G 1 as a matrix equation Av = 0 of form

   a 11 . . . a 1l . . . . . . . . . a k1 . . . a kl       v 1 . . . v l    =    0 . . . 0   
where A is a k × l matrix with coefficients

a g(m)f (r) =        s 1 (r, m) , (r, m) ∈ E 1 -s 1 (m, r) , (m, r) ∈ E 1 0 , otherwise and v consists of variables v f (r) for r ∈ R 1 . By L = {v | Av = 0}
we denote the set of solutions induced by Av = 0. Furthermore, we represent the system of linear equations of ( 6) regarding

G 1 ∪ G 2 as a matrix equation A ′ v ′ = 0 of form            a 11 . . .
k ′ l+1 . . . a k ′ l ′                       v 1 . . . v l v l+1 . . . v l ′            =    0 . . . 0    where A ′ is a k ′ × l ′ matrix with coefficients a g ′ (m)f ′ (r) =        s(r, m) , (r, m) ∈ E 1 ∪ E 2 -s(m, r) , (m, r) ∈ E 1 ∪ E 2 0 , otherwise
where s = s 1 ∪ s 2 and v ′ consists of variables v f ′ (r) of ( 6) for r ∈ R 1 ∪ R 2 . Note that A ′ can always be written in this form, since switching columns and rows will not change solutions. By L ′ = {v ′ | A ′ v ′ = 0} we denote the set of solutions induced by ), we know that the homogeneous system of linear equations from ( 6) regarding G 1 is non-trivial satisfiable, which finally implies that R T ⊆ active s G1∪G2 (S). Theorem 3. Let G 1 and G 2 be metabolic networks. If R T ⊆ active h G1 (S), then R T ⊆ active h G1∪G2 (S). Proof. Follows directly by the definition of hybrid activation together with Theorem 1 and Theorem 2. More formal, R T ⊆ active h G1 (S) holds iff R T ⊆ active t G1 (S) and R T ⊆ active s G1 (S). From Theorem 1 and R T ⊆ active t G1 (S) follows R T ⊆ active t G1∪G2 (S). Analogously, from Theorem 2 and R T ⊆ active s G1 (S) follows R T ⊆ active s G1∪G2 (S). Finally, this implies R T ⊆ active h G1∪G2 (S). In particular, studying the union in case of topological modeling can pinpoint interesting cases. Individual solutions satisfying the topological activation can additionally satisfy the stoichiometric and thus the hybrid activation semantics. A union including such a solution will also adhere to the hybrid standard. In some cases, the union of solutions will display the stoichiometric activation whereas the individual solutions only satisfy the topological activation. Fig. 9 to Fig. 11 display an example of topological metabolic network completions that do not satisfy stoichiometric (and hybrid) activation whereas ), but contradicts minimality and carries no flux r 4 ∈ active s G ({S}), due to accumulation of compound E that contradicts Eq. 6. their union does. Fig. 12 to Fig. 14 provide an example of minimal topological completions that do not satisfy stoichiometric (and hybrid) activation and for which the union does not satisfy it either.

A ′ v ′ = 0. Since A ′ v ′ = 0 is homogeneous, L ⊆ L ′ holds by extending L with zeros for v f ′ (r) with r ∈ R 2 \ R 1 . Thus {v | v ∈ L, ∀r T ∈ R T , v f (rT ) > 0} ⊆ {v | v ∈ L ′ , ∀r T ∈ R T , v f ′ (rT ) > 0} by extending the first set with zeros for v f ′ (r) with r ∈ R 2 \ R 1 . From R T ⊆ active s G1 (S
∪ R 2 = {r 2 , r 3 }. G = G 1 ∪ G 2 sat- isfies r 4 ∈ active t G ({S}
Both observations induce that in general we cannot derive anything about activation of reactions in a graph resulting from the union of two or more graphs. And similarly, we cannot infer about the activation of reactions in subgraphs arbitrarily derived from a graph in which these reactions are activated.

Answer Set Programming with Linear Constraints

For encoding our hybrid problem, we rely upon the theory reasoning capacities of the ASP system clingo that allows us to extend ASP with linear constraints over reals (as addressed in Linear Programming). We confine ourselves below to features relevant to our application and refer the interested reader for details to [START_REF] Gebser | Theory solving made easy with clingo 5[END_REF].

As usual, a logic program consists of rules of the form a 0 : -a 1 ,... , a m , not a m+1 ,... , not a n where each a i is either a (regular) atom of form p(t 1 ,...,t k ) where all t i are terms or a linear constraint atom of form4 '&sum{w 1 *x 1 ;. . .;w l *x l } <= k' that stands for the linear constraint

w 1 • x 1 + • • • + w l • x l ≤ k.
All w i and k are finite sequences of digits with at most one dot5 and represent real-valued coefficients w i and k. Similarly all x i stand for the real-valued variables x i . As usual, not denotes (default) negation. A rule is called a fact if n = 0.

Semantically, a logic program induces a set of stable models, being distinguished models of the program determined by stable models semantics [START_REF] Gelfond | Classical negation in logic programs and disjunctive databases[END_REF]. Such a stable model X is an LC-stable model of a logic program P ,6 if there is an assignment of reals to all real-valued variables occurring in P that (i) satisfies all linear constraints associated with linear constraint atoms in P being in X and (ii) falsifies all linear constraints associated with linear constraint atoms in P being not in X. For instance, the (non-ground) logic program containing the fact 'a("1.5").' along with the rule '&sum{R*x} <= 7 :-a(R).' has the stable model {a("1.5"), &sum{"1.5"*x}<=7}. This model is LC-stable since there is an assignment, e.g. {x → 4.2}, that satisfies the associated linear constraint '1.5 * x ≤ 7'. We regard the stable model along with a satisfying real-valued assignment as a solution to a logic program containing linear constraint atoms. For a more detailed introduction of ASP extended with linear constraints, illustrated with more complex examples, we refer the interested reader to [START_REF] Janhunen | Clingo goes linear constraints over reals and integers[END_REF].

To ease the use of ASP in practice, several extensions have been developed. First of all, rules with variables are viewed as shorthands for the set of their ground instances. Further language constructs include conditional literals and cardinality constraints [START_REF] Simons | Extending and implementing the stable model semantics[END_REF]. The former are of the form a:b 1 ,...,b m , the latter can be written as s{d 1 ;...;d n }t, where a and b i are possibly default-negated (regular) literals and each d j is a conditional literal; s and t provide optional lower and upper bounds on the number of satisfied literals in the cardinality constraint. We refer to b 1 ,...,b m as a condition. The practical value of both constructs becomes apparent when used with variables. For instance, a conditional literal like a(X):b(X) in a rule's antecedent expands to the conjunction of all instances of a(X) for which the corresponding instance of b(X) holds. Similarly, 2{a(X):b(X)}4 true whenever at least two and at most four instances of a(X) (subject to b(X)) are true. Finally, objective functions minimizing the sum of weights w i subject to condition c i are expressed as #minimize{w 1 :c 1 ;. . .;w n :c n }.

In the same way, the syntax of linear constraints offers several convenience features. As above, elements in linear constraint atoms can be conditioned, viz.

'&sum{w 1 *x 1 :c 1 ;...;w l *x l :c n } <= k' where each c i is a condition. Moreover, the theory language for linear constraints offers a domain declaration for real variables, '&dom{lb..ub} = x' expressing that all values of x must lie between lb and ub. And finally the maximization (or minimization) of an objective function can be expressed with &maximize{w 1 *x 1 :c 1 ;...;w l *x l :c n } (by minimize). The full theory grammar for linear constraints over reals is available at https://potassco.org.

Solving Hybrid Metabolic Network Completion

In this section, we present our hybrid approach to metabolic network completion. We start with a factual representation of problem instances. A metabolic network G with a typing function t : M ∪ R → {d,r,s,t}, indicating the origin of the respective entities, is represented as follows:

F (G, t) = {metabolite(m,t(m)) | m ∈ M } ∪ {reaction(r,t(r)) | r ∈ R} ∪ {bounds(r,lb r ,ub r )| r ∈ R} ∪ {objective(r,t(r))| r ∈ R} ∪ {reversible(r) | r ∈ R, rcts(r) ∩ prds(r) = ∅} ∪ {rct(m,s(m, r),r,t(r))| r ∈ R, m ∈ rcts(r)} ∪ {prd(m,s(r, m),r,t(r))| r ∈ R, m ∈ prds(r)}
While most predicates should be self-explanatory, we mention that reversible identifies bidirectional reactions. Only one direction is explicitly represented in our fact format. The four types d, r, s, and t tell us whether an entity stems from the draft or reference network, or belongs to the seeds or targets.

In a metabolic network completion problem, we consider a draft network G = (R ∪ M, E, s), a set S of seed compounds, a set R T of target reactions, and a reference network

G ′ = (R ′ ∪ M ′ , E ′ , s ′
). An instance of this problem is represented by the set of facts F (G, t) ∪ F (G ′ , t ′ ). In it, a key role is played by the typing functions that differentiate the various components:

t(n) =    d, if n ∈ (M \ (T ∪ S)) ∪ (R \ (R S b ∪ R T )) s, if n ∈ S ∪ R S b t, if n ∈ T ∪ R T and t ′ (n) = r, where T = {m ∈ rcts(r) | r ∈ R T } is the set of target compounds and R S b = {r ∈ R | m ∈ S b (G), m ∈ prds(r)}
is the set of reactions related to boundary seeds.

Our encoding of hybrid metabolic network completion is given in Listing 1. Roughly, the first 10 lines lead to a set of candidate reactions for completing the draft network. Their topological validity is checked in lines 12-16 with regular ASP, the stoichiometric one in lines 18-24 in terms of linear constraints. (Lines 1-16 constitute a revision of the encoding in [START_REF] Schaub | Metabolic network expansion with ASP[END_REF].) The last two lines pose a hybrid optimization problem, first minimizing the size of the completion and then maximizing the flux of the target reactions.

In more detail, we begin by defining the auxiliary predicate edge/4 representing directed edges between compounds connected by a reaction. With it, we calculate in Line 4 and 5 the scope Σ G (S) of the draft network G from the seed compounds in S; it is captured by all instances of scope(M,d). This scope is then extended in Line 7/8 via the reference network G ′ to delineate all possibly producible compounds. We draw on this in 24 & sum { R } > "0" : -reaction (R , t ).

26 & maximize { R : o b j e c t i v e(R , t ) }. 27 # minimize { 1 , R : c o m p l e t i o n( R ) }.
Listing 1: Encoding of hybrid metabolic network completion Line 10 when choosing the reactions R ′′ of the completion (cf. Section 2) by restricting their choice to reactions from the reference network whose reactants are producible. This amounts to a topological search space reduction.

The reactions in R ′′ are then used in lines 12-14 to compute the scope Σ G ′′ (S) of the completed network. And R ′′ constitutes a topologically valid completion if all targets in T are producible by the expanded draft network G ′′ : Line 16 checks whether T ⊆ Σ G ′′ (S) holds, which is equivalent to R T ⊆ active t G ′′ (S). Similarly, R ′′ is checked for stoichiometric validity in lines 18-24. For simplicity, we associate reactions with their rate and let their identifiers take real values. Accordingly, Line 18 accounts for ( 5) by imposing lower and upper bounds on each reaction rate. The mass-balance equation ( 6) is enforced for each metabolite M in lines 20-22; it checks whether the sum of products of stoichiometric coefficients and reaction rates equals zero, viz. IS*IR, -OS*OR, IS'*IR', and -OS'*OR'. Reactions IR, OR and IR', OR' belong to the draft and reference network, respectively, and correspond to R ∪ R ′′ . Finally, by enforcing r T > 0 for r T ∈ R T in Line 24, we make sure that R T ⊆ active s G ′′ (S). In all, our encoding ensures that the set R ′′ of reactions chosen in Line 10 induces an augmented network G ′′ in which all targets are activated both topologically as well as stoichiometrically, and is optimal wrt the hybrid optimization criteria.

System and Experiments

In this section, we introduce fluto, our new system for hybrid metabolic network completion, and empirically evaluate its performance. The system relies on the hybrid encoding described in Section 4 along with the hybrid solving capacities of clingo [START_REF] Gebser | Theory solving made easy with clingo 5[END_REF] for implementing the combination of ASP and LP. We use clingo 5.2.0 incorporating as LP solvers either cplex 12.7.0.0 or lpsolve 5.5.2.5 via their respective Python P P P P P P P interfaces. We describe the details of the underlying solving techniques in a separate paper and focus below on application-specific aspects.

The output of fluto consists of two parts. First, the completion R ′′ , given by instances of predicate completion, and second, an assignment of floats to (metabolic flux variables v r for) all r ∈ R ∪ R ′′ . In our example, we get R ′′ = {completion(r 6 ), completion(r 8 ), completion(r 9 )} and {r s1 = 49999.5, r 9 = 49999.5, r 3 = 49999.5, r 2 = 49999.5, r e = 99999.0, r 6 = 49999.5, r 5 = 49999.5, r 4 = 49999.5}.

Variables assigned 0 are omitted. Note the flux value r 8 = 0 even though r 8 ∈ R ′′ . This is to avoid the self-activation of cycle C, D and E. By choosing r 8 , we ensure that the cycle has been externally initiated at some point but activation of r 8 is not necessary at the current steady state.

We analyze (i) the impact of different system configurations (ii) the quality of fluto's approach to metabolic network completion, and (iii) compare the quality of fluto's solutions with other approaches. To have a realistic setting, we use degradations of a functioning metabolic network of Escherichia coli [START_REF] Reed | An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR)[END_REF] comprising 1075 reactions. The network was randomly degraded by 10, 20, 30 and 40 percent, creating 10 networks for each degradation by removing reactions until the target reactions were inactive according to Flux Variability Analysis [START_REF] Becker | Quantitative Prediction of Cellular Metabolism with Constraint-based Models: The COBRA Toolbox[END_REF]. 90 target reactions with varied reactants were randomly chosen for each network, yielding 3600 problem instances in total [START_REF] Prigent | Meneco, a Topology-Based Gap-Filling Tool Applicable to Degraded Genome-Wide Metabolic Networks[END_REF]. The reference network consists of reactions of the original metabolic network.

We ran each benchmark on a Xeon E5520 2.4 GHz processor under Linux limiting RAM to 20 GB. At first, we investigate two alternative optimization strategies for computing completions of minimum size. The first one, branch-and-bound (bb), iteratively produces solutions of better quality until the optimum is found and the other, unsatisfiable core (usc), relies on successively identifying and relaxing unsatisfiable cores until an optimal solution is obtained. Note that we are not only interested in optimal solutions but if unavailable also solutions activating target reactions without trivially restoring the whole reference network. In clingo, bb naturally produces these solutions in contrast to usc. Therefore, we use usc with stratification [START_REF] Ansótegui | SAT-based MaxSAT algorithms[END_REF], which provides at least some suboptimal solutions.

System configurations

The configuration space of fluto is huge. In addition to its own parameters, the ones of clingo and the respective LP solver amplify the number of options. We thus concentrate on distinguished features revealing an impact in our experiments.

The first focus are two system options controlling the hybrid solving nature of fluto. First, prop-n controls the frequency of LP propagation: the consistency of linear constraints is only checked if n% of atoms are decided. Second, the fluto option core-n invokes the irreducible inconsistent set algorithm [START_REF] Ostrowski | ASP modulo CSP: The clingcon system[END_REF] whenever n% of atoms are decided. This algorithm extracts a minimal set of conflicting linear constraints for a given conflict. Note that the second parameter depends on the first one, since conflict analysis may only be invoked if the LP solver found an inconsistency.

The default is to use core-100, prop-0, and use LP solver cplex 7 . This allows us to detect conflicts among the linear constraints as soon as possible and only perform expensive conflict analysis on the full assignment.

To get an overview, we conducted a preliminary experiment using bb and usc with fluto's default configuration on the 10, 20, and 30 percent degraded networks, 2700 instances in total, limiting execution time to 20 minutes. For our performance experiments, we selected at random three networks with at least one instance for which bb and usc could find the optimum in 100 to 600 seconds. With the resulting 270 medium to hard instances, we examined the cross product of values n ∈ {0, 25, 50, 75, 100} for core-n and prop-n, respectively, limiting time to 600 seconds.

Table 1 and Table 2 display the results using bb and usc respectively. The columns increase the value for prop-n and the rows for core-n in steps of 25, i.e., LP propagation becomes less frequent from left to right, and conflict minimization from top to bottom. The first value in each cell is the average runtime in seconds and the value in brackets shows the number of timeouts. The shade of the cells depends on the average runtime, i.e., the darker the cell, the less performant the combination of propagation and conflict minimization heuristics.

Table 1 shows that propagation and conflict minimization heuristics have an overall small impact on the performance of bb optimization. Since bb relies on iterating solutions and learns weaker constraints, only pertaining to the best known bound, while optimizing, the improvement step is less constraint compared to usc. Due to this, conflicts are more likely to appear later on in the optimization process allowing for less impact of frequent LP propagation and conflict minimization Nevertheless, we see a slight perfor- mance improvement of propagating and conflict minimizing for every partial assignment (prop-0, core-0) compared to only on full assignments (prop-100, core-100). To prove the optimum, the solver is still required to cover the whole search space. For this purpose, early pruning and conflict minimization may be effective. Furthermore, we see the best average runtime in the area prop-0-50 at core-75. That indicates a good tradeoff between the better quality conflicts which prune the search effectively and the overhead of the costly conflict minimization. There is no clear best configuration, but prop-25 and core-75 shows the best tradeoff between average runtime and number of timeouts.

usc on the other hand (Table 2), clearly benefits from early propagation and conflict minimization. The area prop-0-75 and core-0-50 has the lowest average runtime and number of timeouts, best among them prop-25 and core-25 with the lowest timeouts and average runtime that is not significantly different from the best value. usc aims at quickly identifying unsatisfiable partial assignments and learning structural constraints building upon each other, which is enhanced by frequent conflict detection and minimization. Disabling LP propagation on partial assignments with usc leads to the overall worst performance and we also see deterioration with core-75 and core-100 in the interval prop-0-50. Overall, usc is more effective than bb for the instances and we see a benefit in early LP propagation and conflict minimization as well as in fine-tuning the heuristics at which point both are applied. Now, we focus on the portfolio configurations of clingo. Those configurations were crafted by experts to enhance the solving performance of problems with certain attributes. To examine their impact, we take the best result for bb (prop-25 and core-75) and usc (prop-25 and core-25), and employ the following clingo options: FR Refers to clingo's configuration frumpy that uses more conservative defaults. JP Refers to clingo's configuration jumpy that uses more aggressive defaults. TW Refers to clingo's configuration tweety that is geared toward typical ASP problems. TR Refers to clingo's configuration trendy that is geared toward industrial problems. CR Refers to clingo's configuration crafty that is geared towards crafted problems. HD Refers to clingo's configuration handy that is geared towards larger problems.

For more information on clingo's configurations, see [START_REF] Gebser | Progress in clasp series 3[END_REF].

Table 3 shows the average runtime in seconds (t) and number of timeouts (to) for all six configurations using bb and usc on the same 270 instances. Even though CR has slightly higher average runtime for bb compared to TW, it is the overall best configuration. This configuration is geared towards problems with an inherent structure compared to randomly generated benchmarks which fits with the metabolic network completion problem at hand since the data is taken from an existing bacteria. Interestingly, bb performs worse under more specific configurations and favors moderate once like TW and CR. This might be due to the changing nature of improvement steps as the optimization process goes on from finding any random solutions to an unsatisfiability proof in the Table 5: Results using best system options.

end. usc on the other hand, benefits from a more structural heuristics in CR and more conservative defaults in FR which allow the solver to explore and collect conflicts instead of frequently restarting and forgetting.

Solution quality

Now, we examine the quality of the solutions provided by fluto. Table 4 gives the number of solutions (#sols) and optima (#opts) obtained by fluto (f) in its default setting within 20 minutes for bb, usc and the best of both (bb+usc), individually for each degradation and overall. The default setting for fluto includes the default configurations for clingo and cplex. The data was obtained in our preliminary experiment using networks with 10, 20, and 30 percent degradation. For 94.3% of the instances fluto(bb+usc) found a solution within the time limit and 82.3% of them were optimal. We observe that bb provides overall more useful solutions but usc acquires more optima, which was to be expected by the nature of the optimization techniques. Additionally, each technique finds solutions to problem instances where the other exceeds the time limit, underlining the merit of using both in tandem. Column verified shows the quality of solutions provided by fluto. Each obtained best solution was checked with cobrapy 0.3.2 [START_REF] Ebrahim | COBRApy: COnstraints-Based Reconstruction and Analysis for Python[END_REF], a renowned system implementing an FBA-based gold standard (for verification only). All solutions found by fluto could be verified by cobrapy. In detail, fluto found a smallest set of reactions completing the draft network for 77.6%, a suboptimal solution for 16.7%, and no solution for 5.6% of the problem instances. Finally, we change the system configuration and examine how fluto scales on harder instances. To this end, we use the best configurations from Section 5.1, prop-25, core-75 and CR for bb, and prop-25, core-25 and CR for usc, and rerun the experiment on all 3600 instances. The results are shown in (Table 5). f(vbs) denotes the virtual best results, meaning for each problem instance the best known solution among the two configurations was verified. For 20% and 30% degradation, we obtain additional 29 and 68 solutions and 41 and 193 optima, respectively. Overall, we find solutions for 92.4% out of the 3600 instances and 72.1% of them are optimal. The number of solutions decreases slightly and the number of optima more drastically with higher degradation. The results show that fluto is capable of finding correct completions for even highly degraded networks for most of the instances in reasonable time.

Comparison to other approaches

We compare the quality of fluto with meneco 1.4.3 [START_REF] Prigent | Meneco, a Topology-Based Gap-Filling Tool Applicable to Degraded Genome-Wide Metabolic Networks[END_REF] and gapfill8 (Satish [START_REF] Satish Kumar | Optimization based automated curation of metabolic reconstructions[END_REF].9 Both meneco and gapfill are systems for metabolic network completion. While meneco pursues the topological approach, gapfill applies the relaxed stoichiometric variant using Inequation (8). We performed an enumeration of all minimal solutions to the completion problem under the topological (meneco), the relaxed stoichiometric (gapfill ), and hybrid (fluto) activation semantics for the 10 percent degraded networks of the benchmark set (900 instances to be completed). First, we compare the quality of individual solutions of fluto and meneco.10 Results are displayed in Table 6. The first two rows give the minimum, average and maximum number of solutions per instance, and reactions per solution, respectively, for fluto and meneco. While fluto finds 19% more solutions on average and twice as many maximum solutions per instance compared to meneco, the numbers of reactions in minimal solutions of both tools are similar. The next four rows pertain to the solution quality as established by cobrapy. First, what percent solutions over all instances could be verified, second, what percent of instances had verified solutions exclusively, third, how many instances had no verified solutions at all, and finally, percent of instances where only a portion of solutions could be verified. All of fluto's solutions could be verified, compared to the 72.04% of meneco across all solutions and 72.94% of instances that were correctly solved. Interestingly, meneco achieves hybrid activation in some but not all solutions for 0.45% (4) of the instances. fluto does not only improve upon the quality of meneco, but also provides more solution per instances without increasing the number of relevant reactions significantly.

To empirically evaluate the properties established in Section 2.4, and be able to compare to gapfill, for which only the union of reactions was available, we examine the union of minimal solutions provided by all three systems and present the results in Table 7. The four rows show, first, for what percent of instances the union of solutions could be verified, second, how many instances had only verified solutions and their union was also verified, third, the percentage of instances where the union of solutions displayed activation of the target reactions even though all individual solutions did not provide that, and forth, instances where the solutions were partly verifiable and their union could also be verified. While again 100% of fluto's solutions could be verified, only 73.3% and 6.2% are obtained for meneco and gapfill, respectively, for 10 percent degraded networks. As reflected by the results, the ignorance of meneco regarding stoichiometry leads to possibly unbalanced networks. Still, the union of solutions provided a useful set of reactions in almost three quarters of the instances, showing merit in the topological approximation of the metabolic network completion problem. On the other hand, the simplified view of gapfill in terms of stoichiometry misguides the search for possible completions and eventually leads to unbalanced networks even in the union. Moreover, gapfill 's ignorance of network topology results in self-activated cycles. By exploiting both topology and stoichiometry, fluto avoids such cycles while still satisfying the stoichiometric activation criteria. The results support the observations made in Section 2.4. For both fluto and meneco all instances, for which the complete solution set could be verified, the union is also verifiable, as well as all unions for instances where meneco established hybrid activation for a fraction of solutions.

Discussion

We presented the first hybrid approach to metabolic network completion by combining topological and stoichiometric constraints in a uniform setting. To this end, we elaborated a formal framework capturing different semantics for the activation of reactions. Based upon these formal foundations, we developed a hybrid ASP encoding reconciling disparate approaches to network completion. The resulting system, fluto, thus combines the advantages of both approaches and yields greatly superior results compared to purely quantitative or qualitative existing systems. Our experiments show that fluto scales to more highly degraded networks and produces useful solutions in reasonable time. In fact, all of fluto's solutions passed the biological gold standard. The exploitation of the network's topology guides the solver to more likely completion candidates, and furthermore avoids self-activated cycles, as obtained in FBA-based approaches. Also, unlike other systems, fluto allows for establishing optimality and address the strict stoichiometric completion problem without approximation.

fluto takes advantage of the hybrid reasoning capacities of the ASP system clingo for extending logic programs with linear constraints over reals. This provides us with a practically relevant application scenario for evaluating this hybrid form of ASP. To us, the most surprising empirical result was the observation that domain-specific heuristic allow for boosting unsatisfiable core based optimization. So far, such heuristics have only been known to improve satisfiability-oriented reasoning modes, and usually hampered unsatisfiability-oriented ones (cf. [START_REF] Gebser | Progress in clasp series 3[END_REF]).

Fig. 1 :

 1 Fig.1: Example of a metabolic network. Compounds and reactions are depicted by circles and rectangles respectively. Dashed reactions are reactions involving the boundary between the organism's metabolism and its environment. r 5 is the target reaction. S 1 and S 2 are boundary (and initiation) seeds. S 3 is assumed to be an initiation seed. Numbers on arrows describe the stoichiometry of reaction (default value is 1).

Fig. 5 :Fig. 6 :

 56 Fig.5: First solution to metabolic network completion under topological activation hypothesis satisfying Equation (10). The production of C cannot be explained by a self-activated cycle and requires an external source of compounds via S 3 and reaction r 7 .

Fig. 7 :

 7 Fig.7: First solution to metabolic network completion under hybrid activation hypothesis satisfying Equation[START_REF]R0 " , d ). reaction (" R1[END_REF] (that is Equations (5), (6), (7) and (10)).

Fig. 8 :

 8 Fig.8: Second solution to metabolic network completion under hybrid activation hypothesis satisfying Equation[START_REF]R0 " , d ). reaction (" R1[END_REF] (that is Equations (5), (6), (7) and (10)).

Fig. 9 :Fig. 10 :Fig. 11 :Fig. 12 :Fig. 13 :

 910111213 Fig.9: Topological completion R 1 = {r 2 } satisfies r 4 ∈ active t G1 ({S}), but carries no flux, due to accumulation of compound B that contradicts Eq. 6.S

1

  edge (R ,M ,N , T ) : -reaction (R , T ) , rct (M ,_ ,R , T ) , prd (N ,_ ,R , T ). 2 edge (R ,M ,N , T ) : -reaction (R , T ) , rct (N ,_ ,R , T ) , prd (M ,_ ,R , T ) , r e v e r s i b l e( R ).

4

  scope (M , d ) : -m e t a b o l i t e(M , s ). 5 scope (M , d ) : -edge (R ,_ ,M , T ) , T != r , scope (N , d ): edge (R ,N ,_ ,T ') , N != M , T '!= r . 7 scope (M , x ) : -scope (M , d ). 8 scope (M , x ) : -edge (R ,_ ,M , _ ) , scope (N , x ): edge (R ,N ,_ , _ ) , N != M . 10 { c o m p l e t i o n( R ) : edge (R ,M ,N , r ) , scope (N , x ) , scope (M , x ) }. 12 scope (M , c ) : -scope (M , d ). 13 scope (M , c ) : -edge (R ,_ ,M , T ) , T != r , scope (N , c ): edge (R ,N ,_ ,T ') , T '!= r , N != M . 14 scope (M , c ) : -c o m p l e t i o n( R ) , edge (R ,_ ,M , r ) , scope (N , c ): edge (R ,N ,_ , r ) , N != M . 16 : -m e t a b o l i t e(M , t ) , not scope (M , c ). 18 & dom { L .. U } = R : -bounds (R ,L , U ). 20 & sum { IS * IR : prd (M , IS , IR , T ) , T != r ; IS '* IR ' : prd (M , IS ' , IR ' , r ) , c o m p l e t i o n( IR '); 21 -OS * OR : rct (M , OS , OR , T ) , T != r ; -OS '* OR ' : rct (M , OS ' , OR ' , r ) , c o m p l e t i o n( OR ') 22 } = "0" : -m e t a b o l i t e(M , _ ).

Table 1 :

 1 Comparison of propagation and core minimization heuristics for bb.

Table 2 :

 2 Comparison of propagation and core minimization heuristics for usc.

Table 3 :

 3 Comparison of clingo's portfolio configurations for bb and usc.

Table 4 :

 4 Comparison of qualitative results.

		f(vbs)	verified
	degradation #sols #opts f(vbs)
	10% (900)	900	900	900
	20% (900)	896	855	896
	30% (900)	848	575	848
	40% (900)	681	68	681
	all (3600)	3325	2398	3325

Table 6 :

 6 Comparison of fluto and meneco solutions for 10 percent degraded networks.

		fluto meneco gapfill
	verified union 100% 73.39% 6.20%
	verified union of verified solutions 100% 72.94%	NA
	verified union of unverified solutions	0% 0.00%	NA
	verified union of partially verified solutions	0% 0.45%	NA

Table 7 :

 7 Comparison of fluto, meneco and gapfill unions for 10 percent degraded networks.

In clingo, theory atoms are preceded by '&'.

In the input language of clingo, such sequences must be quoted to avoid clashes.

This corresponds to the definition of T -stable models using a strict interpretation of theory atoms[START_REF] Gebser | Theory solving made easy with clingo 5[END_REF], and letting T be the theory of linear constraints over reals.

Update of 2011-09-23 see http://www.maranasgroup.com/software.htm

The results for meneco and gapfill are taken from previous work[START_REF] Prigent | Meneco, a Topology-Based Gap-Filling Tool Applicable to Degraded Genome-Wide Metabolic Networks[END_REF], where they were run to completion with no time limit.

There was no data available for the individual solutions of gapfill.

Acknowledgments This work was partially funded by DFG grant SCHA 550/9 and 11 and benefited from the support of the French Government via the National Research Agency investment expenditure program IDEALG ANR-10-BTBR-04.

Appendix A Factual representation of example metabolic network

The factual representation of the metabolic network in Fig. 2 is given in Listing 2.

1 metabolite (" S1 " , s ). metabolite (" S2 " , s ). metabolite (" S3 " , s ). 2 metabolite (" a " , t ). metabolite (" b " , d ). metabolite (" c " , t ). 3 metabolite (" d " , d ). metabolite (" e " , d ). metabolite (" f " , d ). 4 metabolite (" g " , r ).

6 reaction (" R_importS 1 " , s ).

reaction (" R_importS2 " , s ). 7 reversible (" R_importS1 ").

reversibl e (" R_importS 2 "). 8 prd (" S1 " ,"1" ," R_importS1 " , s ). prd (" S2 " ,"1" ," R_importS2 " , s ). 9 reaction (" R_exportF " , d ).

rct (" f " ,"1" ," R_exportF " , d ).


Xutong Liu

Yu-Zhen Janice Chen

John C S Lui

Konstantin Avrachenkov

Graphlet Count Estimation via Convolutional Neural Networks

1 Introduction. Graphlets are defined as k-node connected induced subgraph patterns. For an undirected graph, 3-node graphlets include close triangle ( ) and open triangle ( ). When k = 4, there are six different types of graphlets, e.g., tailed-triangle ( ) and clique ( ) are two possible 4-node graphlets. The number of each graphlet, called graphlet count, is a signature which characterizes the local network structure of a given graph. Graphlet count plays a prominent role in network analysis of many fields, most notably bioinformatics [START_REF] Prulj | Biological network comparison using graphlet degree distribution[END_REF] and social science [START_REF] Chen | A general framework for estimating graphlet statistics via random walk[END_REF].

However, enumerating exact graphlet count is inherently difficult and computational expensive because the number of graphlets grows exponentially large as the graph size and/or graphlet size k grow [START_REF] Chen | A general framework for estimating graphlet statistics via random walk[END_REF]. To deal with this difficulty, many sampling methods were proposed for estimating graphlet count with bounded error [START_REF] Bressan | Counting graphlets: Space vs time[END_REF][START_REF] Chen | A general framework for estimating graphlet statistics via random walk[END_REF][START_REF] Rahman | Graft: An efficient graphlet counting method for large graph analysis[END_REF]. Nevertheless, these methods require large number of samples to be statistically reliable, which is still computationally demanding. Moreover, they have to repeat laborious counting procedure even if a new graph is similar or exactly the same as previous studied graphs.

Intuitively, learning from historic graphs can make estimation more accurate and avoid many repetitive counting to reduce computational cost. Based on this idea, we propose a convolutional neural network (CNN) framework and two preprocessing techniques to estimate graphlet count. 1 Extensive experiments on two types of random graphs and real world biochemistry graphs show that our framework can offer substantial speedup on estimating graphlet count of new graphs with high accuracy.

2 Method. Given a set of undirected graphs and a particular type of k-node graphlet, our objective is to develop a CNN which will be trained using part of dataset with known graphlet counts. After training, the CNN can quickly and accurately predict graphlet counts of other unseen graph samples in the set. Our framework takes the graph adjacency matrix as input and outputs the graphlet count of the input graph. Let us define some notations for our CNN. Let O (l) ∈ R N (l) ×N (l) ×C (l) be the output tensor at layer l, where l = 0, 1, 2, 3, N (l) denotes the width (and height) along each channel and C (l) denotes the channel size. Let O (l) i, j,t be the (i, j) th element along the t th channel. We assign O (0) as the graph adjacency matrix. Mathematically, our CNN structure can be described as follows:

O (l) i, j,t = ReLU(W (l) t • O (l-1) [i : i + H (l) -1, j : j + H (l) -1, : ] + b (l) t ), l = 1, 2, (1) O (3) = ReLU(Flatten(O (2) ) T W (3) + b (3) ). (2)
Equation ( 1) corresponds to two convolution layers. Each layer applies C (l+1) filters over the input feature map O (l-1) , and the t th filter is parameterized by a trainable 3D weight tensor

W (l) t ∈ R H (l) ×H (l) ×C (l)
, where H (l) denotes the width (and height) of the filter. [a : b, c : d, :] is a slicing function which extracts subset of elements indexing from a to b in width, c to d in height and all in channel to form a new tensor. • is the sum of element wise product of two tensors. After adding bias term b (l) t , we apply ReLU (max(0, x)) as the activation function to obtain the output feature map O (l) . Equation ( 2) is associated with the fully connected layer. It flattens the output O (2) into a column vector, applies W (3) , b (3) and ReLU to obtain the estimated graphlet count. Finally, our CNN is trained with back propagation and mean squared error as the loss function.

The above CNN structure inherits the learning power for local structural information of graphs. However, we still need to address the following challenges: (1) The input adjacency matrix is not consistent because graphs in the training set may have different sizes. [START_REF] Bressan | Counting graphlets: Space vs time[END_REF] In practice, real world network dataset may not contain sufficient amount of graph samples for training, which will cause overfitting problem. To address these challenges, we introduce two preprocessing techniques: Adjacency Matrix Zero Padding. To preserve edge connectivity information of all training graphs, we consider the largest graph in the training set, and use its dimension (say N) as the dimension of the input adjacency matrix (N × N). For other graphs in the training set, we take each adjacency matrix and pad it with zero till we have an input matrix of dimension N × N. This solves the varying input size problem. Swapping Augmentation. To acquire sufficient data for training, we take advantage of the graph isomorphism property, where a graph can be expressed by different input adjacency matrices having the same underlying network structure. Our approach is to randomly pick indices i and j, then swap the i th row with j th row and i th column with j th column of the adjacency matrix. We can repeat the swapping operation for each graph m times to create m more training data. Analogous to flipping or rotation of images, we improve CNN's generalization ability and thus improve the accuracy of our model.

3 Data and Metric. Here, we introduce our testing datasets, benchmarking works, and evaluation metrics. Random Graph. We synthesize datasets with two random graph models: random geometric graph (RGG) and Erdos-Renyi (ER) graph. A RGG is constructed by placing nodes uniformly at random in a unit cube and connecting two nodes by an edge if and only if their distance is within a given radius r. In a ER graph, the edge between every two nodes exists with probability p. In each synthetic dataset, we have 3000 training graphs, 300 validation graphs, and 300 testing graphs. Empirical Network. We test on three real world biochemistry datasets: MUTAG [START_REF] Vishwanathan | Graph kernels[END_REF], NCI1 and NCI109 [START_REF] Wale | Comparison of descriptor spaces for chemical compound retrieval and classification[END_REF]. MUTAG dataset contains 188 mutagenic compound graphs. NCI1 and NCI109 each has 4110 and 4127 chemical compound graphs tested on lung and ovarian cancer cells respectively. For MUTAG, we use swapping augmentation to increase the number of training samples. We also apply adjacency matrix zero padding to make all graphs in each dataset have the same size. Benchmark. We compare CNN with three existing frameworks: GRAFT [START_REF] Rahman | Graft: An efficient graphlet counting method for large graph analysis[END_REF], CC2 [START_REF] Bressan | Counting graphlets: Space vs time[END_REF], GUISE [START_REF] Bhuiyan | Guise: Uniform sampling of graphlets for large graph analysis[END_REF], which are based on edge sampling, color coding, and Markov Chain Monte Carlo method respectively. Relative Error. Let c i be the ground truth graphlet count of sample graph i, c i be its estimated count, and there are S samples in the dataset. We compute the mean absolute error of the estimations, mae = Σ S i=1 |c ic i |/S, and mean of ground truth counts, µ = Σ S i=1 c i /S. We take relative error as e = mae/µ . Speed. To ensure a fair comparison, we do not choose running time as the performance metric since it highly depends on hardware and implementation (e.g. running on GPU/CPU). Instead, we measure the number of arithmetic operations they use. For CNN model, we compute the number of floating-point operations (FLOPs). For benchmarking works, we calculate the number of comparison operations in the algorithms.

4 Result. We test our framework on random graph datasets. For approximating 4clique counts, our CNN model achieves less than 8% relative error on 50-node RGGs with radius 0.45 and less than 5% relative error on 50-node ER graphs with edge existing probability 0.5. We also train our CNN models for estimating 4-path ( ), 3-star ( ), 5-path ( ) on the empirical biochemistry datasets. The relative errors on all three datasets are less than 20% of the ground truth counts. For estimating 4-path on MUTAG dataset, our model performs especially well making only 6% relative error. To compare the speed of our CNN with existing methods, the number of arithmetic operations used are calculated. For a fair comparison, we tune the number of iterations for all benchmarking sampling methods, so that they obtain as close relative errors to that of CNN as possible. Figure 1 (a,c) shows that the numbers of arithmetic operations used by GRAFT, CC2, or GUISE are significantly more than that used by CNN. This result demonstrates that our CNN based graphlet count estimation approach offers remarkable speedup on predicting graphlet counts while still maintaining high accuracy.

Fig. 1 .

 1 Fig. 1. Comparison of the number of arithmetic operations used for estimating 4-clique counts, tailed-triangle counts on 50-node ER graphs with edge existing probability 0.5. (a, c) The number of operations used by each framework. (b, d) The relative error each framework makes.


Aurélien Cabarbaye

Titouan Verdu

Miguel Olszanecki Jacson

Gautier Barth

Hattenberger

Accurate control law for low-cost UAV

Present article proposes a method to design a control law and the associated observer to stabilize a UAV. The design is only based on the UAV geometry, mass and propulsion system characteristics which do not require any expensive facilities or software to be obtained. The control inputs are the longitudinal airspeed, the roll angle and the slope angle to ease the guidance control whether manual or automatic. The resulting control only relies on the information provided by a 6 Degrees of Freedom (DoF) Inertial measurement unit (IMU) that makes it suitable for implementation on very basic autopilot board (cf. Paparazzi Chimera, Pixhawk XRacer, ArduPilot APM...). The propeller acts indeed as an airspeed probe which makes additional sensor unnecessary. This low-cost implementation makes it of particular interest for large UAVs fleet.

INTRODUCTION

The ENAC UAV team has developed, since 2003, the Paparazzi UAV (Unmanned Aerial Vehicle) open-source drone project. It enables to convert quickly and easily an electric powered fixed-wing airframe in a semi-automatic drone. However, just like alternative autopilot, the default control law consists of various PID which gains must be tweaked manually during the first flight tests [START_REF]Fixedwing pid tuning guide[END_REF] [START_REF]Subsystem/control[END_REF]. Such method turns out to present some issues:

• It is time-consuming

• It requires good piloting capabilities

• It is hazardous since it does not offer any protection against unstable modes

• It requires some practices to identify which gain must be modified and how much it may be tweaked.

• the thrust control does not consider the complex behaviour of the propulsion system.

• The resulting gain adjustment does not usually offer the shortest response time or oscillation attenuation. * Email address: contact aurelien.cabarbaye@enac.fr

To solve those issues, the present article proposes a more accurate control. This latter is designed from the aeroplane model detailed in section 2, page 1 whose procurement only relies on free data and calculation tools. The control law is exposed in section 3, page 4. This control depends on the state vector which is estimated by the observer, presented in section 4, page 5. The observer is based on the already known aeroplane model of section 2. It takes thus into account the aeroplane dynamics which makes it more precise than other methods that only rely on IMU data (e.g. complementary, Kalman filters [START_REF] Jossué | Quaternion kalman filter for inertial measurement units[END_REF]). However, as these alternatives, it is designed to only require information from a low-grade IMU to be implementable on a wide range of low cost existing autopilot boards (e.g. Paparazzi Apogee [4], Pixhawk XRacer [5]) which is of great interest to make up a drone swarm [START_REF] Sanders | Drone swarms[END_REF]. The resulting control law and observer are tested by simulation in section 5, page 6

The whole process is implemented on the airplane shown on figure 1 Figure 1: Experiment airplane

AIRCRAFT MODEL

The model of the UAV is constructed from its dynamics [START_REF] Etkin | Dynamics of Atmospheric Flight[END_REF]:

10th International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia.

F X -mg sin θ = m uE + qw E -rv E F Y + mg cos θ sin φ = m vE + ru E -pw E F Z + mg cos θ cos φ = m ẇE + pv E -qu E Q L = I x ṗ -I zx ṙ + qr (I z -I y ) -I zx pq Q M = I y q + rp (I x -I z ) + I zx p 2 -r 2 Q N = I z ṙ -I zx ṗ + pq (I y -I x ) + I zx qr p = φ -ψ sin θ q = θ cos φ + ψ cos θ sin φ r = ψ cos θ cos φ -θ sin φ φ = p + (q sin φ + r cos φ) tan θ θ = q cos φ -r sin φ ψ = (q sin φ + r cos φ) sec θ ẋE = u E cos θ cos ψ + v E (sin φ sin θ cos ψ -cos φ sin ψ) + w E (cos φ sin θ cos ψ + sin φ sin ψ) ẏE = u E cos θ sin ψ + v E (sin φ sin θ sin ψ + cos φ cos ψ) + w E (cos φ sin θ sin ψ -sin φ cos ψ) żE = -u E sin θ + v E sin φ cos θ + w E cos φ cos θ u E = u + W x v E = v + W y w E = w + W z (1)
where F X , F Y and F Z are the forces, u, v and w are the velocities, p (roll), q (pitch)and r (yaw) are the rotation speed respectively and W x , W y and W z are the wind speed on the airplane axis X (upstream), Y (right wing direction), Assumption 4 Roll, pitch and yaw rate angles are likely to be small However, the yaw angle (ψ) is not bounded. Therefore, stabilization and guidance must be addressed separately as it is currently done in paparazzi [START_REF]Subsystem/control[END_REF].

To solve the issues exposed in section 1, only the stabilization must be modified. Therefore, the state vector is reduced to:

X =             u v w p q r ϕ θ             (2)
Assumption 5 v << 1 w << u which imply that the sideslip angle and the angle of attack are small: β = tan -1 v u << 1 (symmetric flight) and α = tan -1 w u << 1 (unstall flight). Therefore, the aerodynamics should present an almost linear behaviour.

Thanks to assumptions 3 4 and 5, the differences between a linearized model and the real dynamics should be small. The stabilization of the UAV around a flight configuration is therefore based on its linearized model. This latter is obtained in two steps. The main body is first studied. Then the propulsion system is modelled.

Main body

Athena Vortex Lattice (AVL) is an open source aerodynamics simulation software based on the Vortex Lattice Methods [START_REF] Drela | Avl 3.36 user primer[END_REF]. It assumes an inviscid flow which does not allow good zero-lift drag estimation. However, it seems very accurate to assess aerodynamics evolutions [START_REF] Lopez | Validation of software for the calculation ofaerodynamic coefficients: with a focus on the software package tornado[END_REF]. It is therefore used to generate the linear model of the UAV glider base.

AVL requires information about the shape and the moment of inertia of the aeroplane.

The aeroplane shape is defined with its main aerodynamics characteristic dimensions as shown on figure 2 and with the airfoil shapes.

If the airfoil shapes are unknown, they can be reverse engineered with a profile gauge with sufficient accuracy.

Thanks to assumption 2, the moment of inertia can be easily measured following the Bifilar Pendulum methodology [START_REF] Jardin | Optimized measurements of uav mass moment of inertia with a bifilar pendulum[END_REF] as shown on figure 3 The model provided by AVL is of the form:

Ẋ = AX + BU (3)
Since the experiment airplane only have two aerodynamic controls (elevons), the control input is defined as: It must be noticed that the control surface deflections are managed in degree by AVL rather than in radian for the attitude angles and rates. This particularity is kept unchanged in the following of the article. The resulting elevons deflection δ lef t elevon and δ right elevon is computed as follows:

U = δ elevator δ aileron (
δ lef t elevon = δ elevator -δ aileron (5)
δ right elevon = δ elevator + δ aileron (6)

Propulsion system

The propulsion system is composed of a propeller and an electric motor.

Propeller Propeller behaviour can be described by the force F and torque Q produced. F and Q can be very well assessed from the theory resulting in the mixt between Momentum theory and blade element theory [START_REF] Johnson | Helicopter Theory[END_REF]. Thanks to assumption 5, the impact of the radial airspeed is neglected. The equations can be thus simplified as follow:

F = Kf ω ω 2 + Kf χ ωU (7)
Q = Kq ω ω 2 + Kq u U 2 + Kq χ ωU ( 8
)
where ω is the propeller rotation speed and Kf ω , Kf χ , Kq ω , Kq u , Kqχ are constants that must be estimated from wind tunnel test. Library of such test results are available for a wide range of propellers [START_REF] Brandt | UIUC Propeller Database[END_REF]. The simplified modelled is adjusted to the test results obtained for the APC 8x6E mounted on the experiments aeroplane. The force estimation, as well as the test data, are shown on figure 4. The model estimates very well the experimental data except for those performed in hover (on the top left-hand side). This is because the induced air velocity generated by the propeller is no more negligible compared to the general airspeed. However, this situation is very unlikely to happen in normal flight, so it does not call into question equation [START_REF] Etkin | Dynamics of Atmospheric Flight[END_REF]. The torque estimation and its corresponding test data are shown on figure 5. The very little difference between experimental and estimated data validates equation ( 8). The propeller proposed model can be ultimately validated with the propeller efficiency: η = F •u Q•ω shown on figure 6. The results confirm once more the suitability of equations ( 7) and [START_REF] Drela | Avl 3.36 user primer[END_REF].

Since these equations are very similar, propeller generated force and torque are shown together on figure 7. It can be seen that torque varies in a linear fashion with the force which gives the approximation:

∆Q = a Q/F ∆F (9)
where a Q/F is a constant.

Electric motor

The electric motor behaviour is assessed thanks to two equations which gives the relation between the electrical and the mechanical parameters [START_REF] Drela | First-Order DC Electric Motor Model[END_REF]:

ω = K V (V -RI) (10)
I = K V Q + I 0 ( 11
)
where V and I are respectively the voltage and the current applied to the motor and K V , I 0 and R are respectively the motor velocity constant, zero load current and electric resistance.

The model can be complexified to assess even better the motor behaviour but the adopted one seems to be sufficient [START_REF] Drela | Second-Order DC Electric Motor Model[END_REF].

Motor constants K V , I 0 and R can be easily retrieved if not provided by the manufacturer [START_REF] Drela | Measurement of Brushed DC Electric Motor Constants[END_REF].

Figure 7: torque function of force

CONTROL

The aim of the control is to stabilise the aeroplane to make it as easy to steer as possible whether in manual or automatic guidance mode (cf. AU T O1, AU T O2 [START_REF]Subsystem/control[END_REF]). To do so, the control must integrate a stall protection, fly the aeroplane symmetric and be easily integrated into existing autopilot.

Control idea

The control states are designed first. Usually, the thrust or the power are directly controlled by the pilot in closed loop fashion. I order to ease the control, a closed loop control is built to stabilise the longitudinal velocity u.

Most aeroplane lateral controls are built to set the roll angle ϕ to a reference value. This enables an effective path following since there is a direct relation between ϕ and the turning radius. ϕ is therefore chosen as the second control state.

For the longitudinal control, there is no consensus. Aircraft are usually controlled in pitch rate, in vertical acceleration or in pitch angle [START_REF] Hamel | In-Flight Simulators and Fly-by-Wire/Light Demonstrators: A Historical Account of International Aeronautical Research[END_REF]. A control on a derivative (cf. pitch rate, vertical acceleration) imposes more work to the guidance control. A control on pitch angle is not very precise because of the angle of attack variations during the flight. Therefore, the climb slope γ is chosen as third control state. Thanks to assumption 5, γ follows:

γ = θ -α ≈ θ - w u ( 12
)
Applying the small disturbance linearization method [START_REF] Etkin | Dynamics of Atmospheric Flight[END_REF]:

∆γ = ∆θ - ∆w u 0 + w 0 u 2 0 ∆u = ∆θ - 1 u 0 ∆w + α 0 u 0 ∆u (13)
Therefore, the output vector becomes:

Y c =   u ϕ γ   = C c X ( 14
)
where

C c =   1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 α u0 0 -1 u0 0 0 0 0 1  
Three controls input are required to control the three outputs. F is chosen for the third one to removes the nonlinearity of equations 7. Therefore, from equation (4) U becomes:

U =   δ elevator δ aileron F   ( 15
)
And equation ( 3) gives

Ẋ = A c X + B c U ( 16
)
10th International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia.

where

A c = A and B c =             B 1 M 0 0 - a Q/F Ixx 0 0 0 0            
Notice that the relation between force and torque of equation [START_REF] Etkin | Dynamics of Atmospheric Flight[END_REF] has been considered to compensate the torque as soon as it is applied, without the delay that would have been otherwise induced by the dynamics.

The controllability matrix

C m = B c A c B c A 2 c B c ... A 7 c B c
has a rank equal to the length of A c , so the system is controllable.

Control law

The control law is based on the method of the steady state tracking which consists in stabilising first all the states of X and then to alter some of them to converge the output to the desired value [START_REF] Williams | Linear State-Space Control Systems[END_REF].

The control law is of the form:

U = -KX + GR (17)
where K and G are two gain matrices ad R is the desired output.

Combining equations ( 16) and ( 17) leads to:

Ẋ = (A c -B c K) X + B c GR (18)
K is first chosen to make all the states stable. Therefore, the poles p c of the matrix (A c -B c K) are set negatives thanks to the pole placement method. In order to offer good stability properties, the poles are chosen from A c poles. Imaginary parts are kept unchanged, positive real parts are substituted by their opposite and non-sufficiently negative real parts are changed by better ones. Then the convergence of the output is satisfied as follows:

lim t→∞ (A c -B c K) X (t) + B c GR = 0 lim t→∞ Y c (t) = lim t→∞ C c X (t) = R (19)
The system [START_REF] Chen | Analog and Digital Control System Design: Transfer-Function, State-Space, and Algebraic Methods[END_REF] leads to:

G = -C c (A c -B c K) -1 B c -1
(20)

Control saturation and trim

Saturations can be added to consider the servo travel limit and the propulsion maximum power. A trim can also be added to balance the model errors. It is simply implemented modifying the values of the control required to fly in the condition of the linearization. The actual control input U trim is:

U trim =   δ elevator + δ 0elevator δ aileron + δ 0aileron F + F 0   (21)
The deflections of the ailerons are obtained from equations ( 5) and ( 6). Resolving quadratic equation ( 7) to get U trim [START_REF] Jossué | Quaternion kalman filter for inertial measurement units[END_REF] leads to ω. Then combining equations ( 8), ( 10) and ( 11) leads to:

V = ω K V + RK V Kq ω ω 2 + RK V Kq u U 2 + RK V Kq χ ωU + RI 0 ( 22
)
where U = u + U 0 and U 0 is the airspeed of the linearization. Noting that 0 ≤ V ≤ V batt and δ min ≤ δ elevon ≤ δ max , the control applied becomes:

  V δ right elevon δ lef t elevon   applied = min     V batt δ max δ max   , max     V δ right elevon δ lef t elevon   ,   0 δ min δ min       (23)

Stall protection

A stall protection can be added to the control input γ. To balance the weight, the lift must be [START_REF] Raymer | Aircraft design: A conceptual approach[END_REF]:

L cos (ϕ) = mg cos (γ) (24)
which gives

γ max = cos -1 1 2 ρU 2 SCl cos (ϕ) mg ( 25
)
where γ max is the limit of the desired output γ when 1 2 ρU 2 SCl cos (ϕ) ≤ mg, verifying γ max ∈ -π 2 , 0 . 1 2 ρU 2 is the dynamic pressure, S is the wing surface and Cl is the lift coefficient. Cl can be for example fixed to fly, as for general aviation, at least at 1.3 • V s in normal flight and at least at 1.1 • V s in final approach (V s is the stall speed). That is to say: Cl = Clmax 1.14 in normal flight and Cl = Clmax 1.05 in final approach.

OBSERVER

A good estimation of the state vector X is required by the control law (cf. equation [START_REF] Williams | Linear State-Space Control Systems[END_REF]). An observer is therefore built to provide it.

Baseline idea

To be as low cost as possible to implement, the observer only relies on the information provided by a 6 DoF IMU. The output vector Y o of the observer is therefore:

Y o =         u vẇ p q r         = C o X + D o U (
O m =       C o C o A c C o A 2 c ... C o A 7 c      
has a rank equal to the length of A c , so the system is observable.

The dynamics of the simulated observer is [START_REF] Chen | Analog and Digital Control System Design: Transfer-Function, State-Space, and Algebraic Methods[END_REF]:

Ẋ = A c X + B c U + L Y o -Ŷo (27)
where Ŷo = C o X + D o U From equations ( 16), ( 26) and ( 27), the estimation error X = X -X dynamics is:

Ẋ = Ẋ -Ẋ = (A c -LC o ) X (28)
To urge the X to zero, matrix L is chosen following the pole placement method with all the poles p o of (A c -LC o ) negatives.

Observer improvements

As it will be shown in section 5, page 6, the approximation introduced in equation ( 7) produces an error of estimation. To improve the observer, the non-linearity of the propulsion system is considered. Equation ( 27) is thus modified as follows:

Ẋ = A c X + B o U + EQ + L Y o -Ŷo (29)
where

B o =             B 1 M 0 0 0 0 0 0 0             and E =             0 0 0 -1 Ixx 0 0 0 0            
and matrix L obtained in section 4.1, page 5 is kept unchanged. Q and F are computed thanks to equations ( 7) and [START_REF] Drela | Avl 3.36 user primer[END_REF] where ω is computed solving quadratic equation ( 22), where V is the result of equation (23).

SIMULATION

A simulation is performed for the experiments aeroplane. The linearization has been made for u = 10 m • s -1 in trimmed steady flight. The resulting data considered for the simulation are exposed in Appendix A:, page 8.

Control results

The control input of equation ( 17) is fixed to:

R =   5 -0.5 0.5   (30)
The initial state is fixed arbitrarily to:

X 0 =             1 0 -1 0 0 0 0.5 0             (31)
States evolution is shown on figure 8. All the states are stable

0 5 10 15 20 time (s) -2 0 2 4 Speed (m/s

Observer results

The observer is based on the information given by an IMU. This latter is supposed to be very noisy. To simulate it properly, a white noise has been added to each state of equation (26). The resulting collected information is shown on figure 11 The integration of the observer starts with the esti-

X0 =             0 0 0 0 0 0 0 0             . ( 32
)
Figure 12 shows what would have been the estimation error of the observer if it has been kept in its baseline form defined in equation ( 27). It can be noticed that if it seems to 2) are precisely estimated. u estimation is very little bias even though no specific sensor has been used (e.g. pitot tube). That can be explained by the propeller which depends a lot on airspeed u (cf. equation ( 7)) and therefore acts as an alternative airspeed probe. Moreover, the estimation errors converge very rapidly to zero (less than a second) which makes it able to compensate for most of wind gusts. If a faster convergence is needed, the poles p o can be increased.

CONCLUSION

The aircraft model has been built in section 2 in two steps. First, the linearized main aeroplane body model is obtained from AVL software. Then the non-linear propulsion system 10th International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia. model is adjusted from freely available wind tunnel tests. A simplification of the propeller model enables to build a linear overall model that is used as a starting point for the control law design in section 3. This control law focusses on airspeed u and angles φ and γ to ease the guidance whether manual or automatic. The states knowledge required by the control is provided by an observer designed in section 4. This observer is based on the non-linear model resulting from the mix between the linear main body model and the non-linear propulsion model, rather than on the linear model used for the control law design. This enables to track much better the states and in particular, the longitudinal speed since the propeller plays the part of an airspeed probe. In addition, the observer high convergence speed makes it suitable to deal with wind gusts. Tests of control law and observer are performed by simulation in section 5. Results seem to prove the suitability of the proposed method.

ZAssumption 1 Assumption 2 Assumption 3

 123 (downward), m and I =   I x I xy I xz I xy I y I yz I xz I yz I z   are the airplane mass and moments of inertia. g is the gravity acceleration. x, y, and z are the earth reference frame axis. φ (roll), θ (pitch) and ψ (yaw) are the euler angles between the earth and the airplane reference frames.The model can be simplified with the following assumptions: Thanks to the relative symmetries around the airplane axis, I = The gusts are supposed to be slow so the wind parameters W x , W y and W z are constant. Roll and slope angles are likely to be small, in the range :(ϕ, φ) ∈ [-0.5; 0.5] 2

  4) 10th International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia.

Figure 2 :

 2 Figure 2: experiment airplane geometry in AVL

Figure 3 :

 3 Figure 3: Moment of inertia estimation experiment

10th

  International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia.

Figure 5 :Figure 6 :

 56 Figure 5: Torque function of advance ratio J

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: States X

Figure 11 :

 11 Figure 11: Observer output Y o provided by IMU mate states X arbitrarily set to:

Figure 12 :Figure 13 :

 1213 Figure 12: Estimation error Xconverge, a bias remains, in particular for the airplane longitudinal airspeed u. That is why the observer is improved in section 4.2, page 6. Figure13shows the estimation error of the final observer defined by equation (29). All the states

10th

  International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia. Kf ω = 6.49343E -06; (33) Kf ξ = -0.000209663; (34) Kq ω = 1.03864E -07; (35) Kq ξ = 5.02077E -07; (36) Kq u = -9.98976E -05; 093 + 1.1063i -1.093 -1.1063i -3.158 + 4.6121i -3.158 -4.6121i -1 8080 -0.8512 1.7957 -1.8668 -0.9899 0.3489 -4.6903 -3.0280 -2.1717 1.2099 0.6040 0.3691 0.0817 2.9449 -10.2553 1.7359 7.4471 -0.0382 -4.6566 0.0428 -0.2710 0.4464 -2.1005 -14.0328   (46) 10th International Micro-Air Vehicles Conference 22nd-23rd November 2018. Melbourne, Australia.