Raphaël Lallement
email: raphael.lallement@laas.fr

Lavindra De Silva
email: lavindra.desilva@nottingham.ac.uk

Rachid Alami
email: rachid.alami@laas.fr

HATP: Hierarchical Agent-based Task Planner Demonstration

Keywords: Multi-robot planning, HTN planning, Geometric planning

INTRODUCTION

Hierarchical Task Network (HTN) planning is a proven approach to solving complex, real world planning problems more efficiently than planning from first principles when "standard operating procedures" (or "recipes") can be supplied by the user [START_REF] Nau | SHOP2: An HTN Planning System[END_REF][START_REF] Nau | SHOP: Simple Hierarchical Ordered Planner[END_REF]. By planning for tasks in the same order that they are later executed, total-order HTN planners always know the complete state of the world at each planning step. This enables writing more expressive planning domains than what is possible in partial-order HTN planning, such as preconditions with calls to external procedures [START_REF] Nau | SHOP: Simple Hierarchical Ordered Planner[END_REF]. Such features have facilitated the use of total-order HTN planners in agent systems and seen them excel in AI games [START_REF] De | Planning on Demand in BDI Systems[END_REF][START_REF] Dix | IM-PACTing SHOP: Putting an AI Planner Into a Multi-Agent Environment[END_REF][START_REF] Nau | Applications of SHOP and SHOP2[END_REF][START_REF] Stephen | Success in Spades: Using AI Planning Techniques to Win the World Championship of Computer Bridge[END_REF].

This paper describes the Hierarchical Agent-based Task Planner (HATP), a total-order HTN planner. Since its first implementation [START_REF] Montreuil | Planning Human Centered Robot Activities[END_REF], HATP has had various extensions and integrations over the years, such as support for splitting a solution into multiple streams and assigning them to the agents in the domain [START_REF] Alili | A Task Planner for an Autonomous Social Robot[END_REF]; modelling their beliefs as distinct world states [START_REF] Warnier | When the Robot Puts Itself in Your Shoes. Managing and Exploiting Human and Robot Beliefs[END_REF]; allowing "social rules" to be included by the user to define what kind of agent behaviour is appropriate [START_REF] Alili | A Task Planner for an Autonomous Social Robot[END_REF][START_REF] Montreuil | Planning Human Centered Robot Activities[END_REF]; allowing tasks to be planned by taking the human's safety and comfort into account [START_REF] Akin Sisbot | Supervision and Motion Planning for a Mobile Manipulator Interacting with Humans[END_REF]; and to interleave HTN and geometric planning [START_REF] Alili | Interleaving Symbolic and Geometric Reasoning for a Robotic Assistant[END_REF][START_REF] Lavindra De Silva | A New Approach to Combined Symbolic-Geometric Backtracking in the Context of Human-Robot Interaction[END_REF][START_REF] Lavindra De Silva | An Interface for Interleaved Symbolic-Geometric Planning and Backtracking[END_REF].

Since many of these implementations have remained prototypes, we have significantly enhanced them as well as HATP itself, and integrated them into a stand-alone HATP distribution, which is now available as open source software (under a BSD 2-Clause License). This paper presents some of our recent improvements to HATP, and gives an overview of its user-friendly language, which treats agents as distinct entities; its mechanisms for effective control over decomposition; and its integration into our robotics framework.

LANGUAGE FEATURES

The HATP syntax is a product of lessons learnt over years of practical integration with real agent/robotic systems and other applications. Like in standard HTN planning, an HATP domain consists of a set of user supplied HTN methods ("recipes") and operators. For defining these, however, HATP employs some useful structured programming concepts such as structures, conditionals, and typing. Moreover, the agent types and object types in an HATP domain/problem are defined separately as a collection of structures called entities. An entity-structure has a set of attributes, each of which represents either a data value, or a relationship between the entity and other entities. An attribute can be one of the standard data types boolean, integer, or string. An attribute can also correspond to a set of objects, which can be manipulated in preconditions and effects via standard set operations. For example, an entity representing a robot-agent may have the attribute carry of type Container declared as a set, denoting that the robot can carry multiple objects of type Container.

An HATP initial state is an instantiation of the defined entities, with value assignments to their attributes; e.g. the robot-agent above could be instantiated with "robot = new Agent;" and its containers added via the statements "robot.carry add con1; robot.carry add con2;" where con1 and con2 are instantiated container objects.

To show that HATP's syntax is sound with respect to traditional HTN languages, we have developed in [START_REF] Lavindra De Silva | The HATP Hierarchical Planner: Formalisation and an Initial Study of its Usability and Practicality[END_REF] a formal mapping from HATP's syntax into that of SHOP-the dominant HTN planner in the literature. Intuitively, attributes of entities correspond to predicate symbols in SHOP, and the entity names and attribute values to the parameters of predicates. For example, if variables Rob, Loc1, and Loc2 represent robots and locations, then, the assignment "Rob.at = Loc2" will map to the two predicates "at(Rob,Loc2)" and "¬at(Rob,Loc1)" in SHOP (where Loc1 is the robot's previous location), and "Loc2.empty = false" to "¬empty(Loc2)".

HATP supports developing agent-oriented domains by making a distinction between agents and other types of objects. Agents are treated as "first class citizens", and it is possible to define different types of agents. By virtue of making this distinction, HATP is able to split a solution (sequence of action nodes) found into multiple solution "streams", one per agent in the solution, and to add causal links between streams for synchronisation [START_REF] Alili | A Task Planner for an Autonomous Social Robot[END_REF]. The various streams may then be executed concurrently in a multi-agent system by synchronising where necessary.

CONTROL OVER SEARCH

HATP provides ways to control how variables occurring in methods are bound at runtime, via the SELECT (which we abbreviate to SEL), SEL-ORDERED, and SEL-ONCE constructs. The SEL construct binds the given variable in the usual way. In essence, the construct amounts to a "backtrack point" that allows all values of the associated variable-and thus all ground instances of the method-to be considered. SEL-ORDERED binds the variable in some given order, governed by a user-supplied ordering relation; the variable can be bound in ascending or descending order with respect to the relation. Finally, SEL-ONCE binds the variable only once-any remaining bindings are disregarded. This offers a reduction in the branching factor at the expense of completeness, as some of the ignored bindings may also yield HATP solutions. As an example of the value of SEL-ONCE, consider the setting where if one robot cannot navigate to a destination then none of the others will be able to either. Here, there is no need to consider all possible robot-object bindings: trying a single binding will suffice.

Recently, we have enabled the encoding of domain-dependent backtrack preferences by a developer, allowing the more "promising" backtrack points to be explored first, resulting in faster search. These preferences can be conditional, and the next backtrack point to explore can be at a different "level" in the search tree. For example, if it is impossible to plan a particular route to a restaurant (e.g. due to a road closure), it might be faster in some settings, before trying the other (possibly numerous) routes, to "jump" to an earlier backtrack point in the search tree and pick a different restaurant.

SYMBOLIC-GEOMETRIC PLANNING

While an HTN hierarchy allows one to intuitively reason about high-level tasks in terms of more specific tasks, and eventually in terms of basic actions, the latter are still "abstractions" of the lowest possible level of detail, and they make certain assumptions about the world. For example, an HATP operator move(Robot,From,To) might assume that as long as location To is adjacent to location From, that the robot at From will be able to navigate to location To. In reality, however, this will not work when certain geometrical characteristics of the robot and the route make the move physically impossible. Combining HATP with geometric planning algorithms used in robotics is therefore crucial to be able to obtain primitive solutions that are viable in the real world.

To this end, we have developed a tighter integration with a geometric planner (GP) by taking inspiration from the prototype integrations in [START_REF] Lavindra De Silva | A New Approach to Combined Symbolic-Geometric Backtracking in the Context of Human-Robot Interaction[END_REF][START_REF] Lavindra De Silva | An Interface for Interleaved Symbolic-Geometric Planning and Backtracking[END_REF]. Interleaved HTN and geometric planning is now a built-in feature of HATP. The main improvements are the following new syntactic constructs (which concretise the concepts discussed in [START_REF] Gharbi | Combining Symbolic and Geometric Planning to Synthesize Human-Aware Plans: Toward More Efficient Combined Search[END_REF]) for action definitions, which can be used to: (i) consult the GP, during HTN planning, to determine against the current 3D world state whether there exists a viable trajectory for an action; (ii) apply such a trajectory to the 3D world; (iii) compute the symbolic facts resulting from the trajectory, including any relevant side effects such as distanceMoved(robot1,15); (iv) transparently add such facts into HATP's symbolic world state; and (v) constrain the geometric goal states (positions and orientations) considered.

The listing below illustrates our new constructs. The precondition adds constraints to ensure that all robots in the same travel group as Rob are visible to it after the move; the constraints are considered by the GP when building goal states corresponding to To. The precondition then simulates a move in the 3D world to check whether a trajectory exists for one such goal state. If so, the "projects" line applies the trajectory to the 3D world, and the "effects" line computes and adds the resulting symbolic facts into HATP's state. The "cost" line retrieves the actual moved distance, computed by the GP. To isIn From.adjacent; checkMove(Rob, From, To) == true; // simulate to check that a move is possible }; projects { applyMove(Rob, From, To) }; // do the move within the 3D world effects { Rob.at = To; }; // any side effect is added transparently cost { moveCost(From, To) }; // used for finding the optimal plan duration { moveDur(From, To) }; // used to estimate a plan's execution period We evaluated our tighter integration of symbolic and geometric planning with a new case study, which required extending the GP to be able to model the dynamics of quadrotor UAVs, in order to simulate the assembly of complex structures by multiple cooperating UAVs. The screenshots below illustrate the assembly of coloured bars by a single UAV with an attached robotic arm. Each bar is moved from a holder (top left) and inserted into a matching slot, or into two vertical bars. The high-level tasks (e.g. 'pick' and 'insert') are planned by HATP, and their trajectories by the GP. Any correct high-level plan involves assembling the long blue bar last; otherwise, there would not exist a trajectory (due to a lack of space under the blue bar) for the UAV to assemble the long green bar.

Our second scenario is collaborative manipulation involving a robot and human [START_REF] Gharbi | Combining Symbolic and Geometric Planning to Synthesize Human-Aware Plans: Toward More Efficient Combined Search[END_REF]. At one level of abstraction, the robot must find and pick the objects on a table and move them onto another; at the lower level of abstraction, the robot must check that objects are graspable, and visible and reachable to the human after the move.

LAAS ROBOTICS ARCHITECTURE (LRA)

While HATP can be used as a stand-alone planner, it is also part of the LRA [START_REF] Fleury | GenoM: A Tool for the Specification and the Implementation of Operating Modules in a Distributed Robot Architecture[END_REF]: an interconnected set of components responsible for different aspects of the overall system. An important component is the Move3D [START_REF] Thierry Simeon | Move3D: A Generic Platform for Path Planning[END_REF] GP, which represents the real world in 3D and uses this as input for geometric planning. The robot uses various sensors to update its 3D state in real-time, e.g. a tag-based stereovision system is used for object identification and localisation, and a Kinect (Microsoft) sensor is used for tracking humans. The execution controller is the Procedural Reasoning System [START_REF] Félix Ingrand | PRS: A High Level Supervision and Control Language for Autonomous Mobile Robots[END_REF], responsible for invoking HATP when a task needs to be planned, and for executing the solution returned, by invoking the corresponding actuators on the (possibly simulated) robot. More recently, we have made HATP accessible to the robotics community by making each of HATP's constituent modules a ROS node. This enables standard and safe interoperability between HATP modules (e.g. via strongly typed data structures), and between HATP modules and other ROS nodes, e.g. the ones that control robots in the LRA.


Jocelyne Erhel
email: jocelyne.erhel@inria.fr

Tangi Migot
email: tangi.migot@insa-rennes.fr

Characterizations of Solutions in Geochemistry: Existence, Uniqueness and Precipitation Diagram

Keywords: geochemistry, existence, uniqueness, precipitation diagram, modeling, polynomial system, convex optimization, precipitation-dissolution reaction

In this paper, we study the properties of a geochemical model involving aqueous and precipitationdissolution reactions at a local equilibrium in a diluted solution. This model can be derived from the minimization of the free Gibbs energy subject to linear constraints. By using logarithmic variables, we define another minimization problem subject to different linear constraints with reduced size. The new objective function is strictly convex, so that uniqueness is straightforward. Moreover, existence conditions are directly related to the totals, which are the parameters in the mass balance equation. These results allow us to define a partition of the totals into mineral states, where a given subset of minerals are present. This precipitation diagram is inspired from thermodynamic diagrams where a phase depends on physical parameters. Using the polynomial structure of the problem, we provide characterizations and an algorithm to compute the precipitation diagram. Numerical computations on some examples illustrate this approach.

Introduction

Many environmental studies rely on modeling geochemical reactions such as aquifer contamination, underground waste disposal, remediation. In some of these applications the reactions can be assumed to be fast compared with their environment and so the reactions are considered in a local equilibrium. In this latter case, the model involves non-linear algebraic equations [START_REF] Bethke | Geochemical Reaction Modeling : Concepts and Applications[END_REF]. Here, we focus on aqueous and precipitation-dissolution reactions so only liquid and mineral species are involved. Many applications consider geochemical reactions coupled with hydrodynamic processes as in reactive transport model [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF]. It is therefore of great interest to get a deep knowledge of the geochemical properties and to solve efficiently the algebraic equations.

We consider here a diluted solution, as in some reactive transport models such as [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Kräutle | The semismooth Newton method for multicomponent reactive transport with minerals[END_REF] or some exercises such as the MoMaS benchmark [START_REF] Carrayrou | Reactive transport benchmark of MoMaS[END_REF][START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problemsthe MoMaS benchmark case[END_REF]. We use an approach based on the definition of primary and secondary aqueous species, which are related by mass action laws [START_REF] Zeleznik | Calculation of complex chemical equilibria[END_REF]. Regarding the minerals, their appearance or disappearance is governed by saturation thresholds. The mathematical model is completed by mass balance equations, with so-called totals which are given data of the system, as for example in [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF][START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF].

In order to solve the non-linear equations resulting from the geochemical model, the main challenge is to find out which minerals are precipitated and which ones are dissolved, in other words to check whether saturation thresholds are active or not. Some numerical algorithms are based on a combinatorial approach [START_REF] Migot | Gas phase appearance and disappearance as a problem with complementarity con-straints[END_REF][START_REF] Bethke | Geochemical and biogeochemical reaction modeling[END_REF][START_REF] Parkhurst | Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF]. Because the system can be formulated as a non-linear complementarity problem [START_REF] Cottle | The Linear Complementarity Problem[END_REF], numerical algorithms are frequently based on semi-smooth Newton methods [START_REF] Georget | A speciation solver for cement paste modeling and the semismooth Newton method[END_REF][START_REF] Kräutle | The semismooth Newton method for multicomponent reactive transport with minerals[END_REF] or interior-point methods [START_REF] Leal | Computational methods for reactive transport modeling: a Gibbs energy minimization approach for multiphase equilibrium calculations[END_REF][START_REF] Saaf | Computing general chemical equilibria with an interior-point method[END_REF][START_REF] Saaf | A study of reactive transport phenomena in porous media[END_REF].

Complementarity problems have been extremely popular in the field of optimization theory and have been extensively studied, [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF] and references therein. Many first-order optimality conditions of a minimization problem with inequality constraints can be cast as a complementarity problem. In the other way round, the geochemical model can be viewed as first-order optimality conditions.

This approach was previously used connected with the minimization of the Gibbs energy, see for instance in [START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF][START_REF] Smith | Chemical Reaction Equilibrium Analysis: Theory and Algorithms[END_REF]. More recent theoretical studies also took this path in the context of reactive transport in [START_REF] Brunner | Multiphase multicomponent flow in porous media with general reactions: efficient problem formulations, conservative discretizations, and convergence analysis[END_REF][START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF]. The Gibbs energy is a function of all the species (primary and secondary aqueous species, minerals). Linear constraints are the mass balance equations (dimension the number of aqueous components) and positive quantities of minerals (dimension the number of minerals). Minimizing the Gibbs energy with these constraints is a convex problem, ensuring a global minimum. However, it requires additional assumptions or transformations to get existence or uniqueness, see for instance [START_REF] Brunner | Multiphase multicomponent flow in porous media with general reactions: efficient problem formulations, conservative discretizations, and convergence analysis[END_REF][START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF].

In this paper, we propose a different formulation, which is new as far as we know. The objective function is based on the mass balance equation and the mass action laws, whereas linear constraints are conditions on the saturation thresholds. Under mild assumptions on the stoichiometry, we can prove that the first-order optimality conditions are nothing else than the geochemical model considered. This minimization problem has several merits compared to the aforementioned one. The variables of the objective function are only the aqueous components and the number of constraints is only the number of minerals. Moreover, the objective function is strictly convex, so that it is easy to assess the uniqueness of a solution. The analysis of existence is directly related to the totals, which are input data of the model. Finally, the strict convexity as well as the reduced dimensions would allow interesting numerical performances.

In order to do so, we have to use logarithmic variables for the aqueous species, thus to assume that the variables are strictly positive. This usage of logarithms is quite common in the literature, see for example [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of momas with specy[END_REF][START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMaS benchmark[END_REF][START_REF] Hoffmann | A parallel global-implicit 2-d solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF][START_REF] Mayer | Solution of the momas reactive transport benchmark with min3pmodel formulation and simulation results[END_REF][START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. In this paper, we carefully analyze the mathematical model and prove that it is equivalent to a reduced system with only strictly positive aqueous variables, completed by an independent trivial system with non negative aqueous variables. To the best of our knowledge, our theoretical justification of logarithms is different from the literature, as in [START_REF] Brunner | Multiphase multicomponent flow in porous media with general reactions: efficient problem formulations, conservative discretizations, and convergence analysis[END_REF][START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF].

As we already pointed out, one motivation is to characterize the set of precipitated minerals and the complement set of dissolved minerals. First, we charac-terize the set of data T such that there exists a solution. Therefore, we can associate the totals in T to a unique geochemical system. Then, we introduce the concept of precipitation diagram. More precisely, we build a partition of the set T such that each subset is associated to a geochemical system with a given set of precipitated minerals. This partition defines a diagram, in the same spirit as thermodynamic diagrams (or stability diagrams), see for instance [START_REF] Connolly | An automated strategy for calculation of phase diagram sections and retrieval of rock properties as a function of physical conditions[END_REF][START_REF] De Capitani | The computation of equilibrium assemblage diagrams with Theriak/Domino software[END_REF]. However, in our case, the axis are the total concentrations and not some physical quantities like the pressure or the temperature.

We propose to build precipitation diagrams thanks to symbolic computations. This approach provides a graphical view of the conditions under which some minerals can precipitate. With numerical values of given totals, the symbolic approach can compute the solution, although this is not the main objective. Moreover, the diagram can be used in the context of reactive transport models, for example to provide an initial guess of non-linear algorithms or to find out in which subset are the totals. Then the complementary problem could be replaced by an easier system of equations, which can be solved by a classical Newton's method.

We developed a software in Maple and used it for several test cases. These experiments show that precipitation diagrams are quite complex, even for simple examples. Thus, it would not be easy to compute them numerically. We also point out that the CPU and memory requirements are reasonable.

The paper is organized as follows: in Section 2, we describe the geochemical model and show that it can be reduced to equations with strictly positive variables, justifying the use of logarithmic unknowns. We introduce an elementary example of geochemical systems as a common thread throughout the paper. Then, in Section 3, we define the equivalent optimization problem and prove uniqueness of the solution along with some results of existence. We provide an example where existence is conditioned by the totals. Based on these theoretical outputs, we introduce in Section 4 our concept of a precipitation diagram. Our numerical experiments of Section 5 illustrate the complexity of these diagrams and the information they provide. Finally, we discuss some future work in Section 6.

A Model for Geochemistry

In this section, we describe a mathematical model for geochemistry, which has already been used in the context of reactive transport, see [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMaS benchmark[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF][START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. We consider a species/components model frequently used in geochemistry modeling as in [START_REF] Bethke | Geochemical Reaction Modeling : Concepts and Applications[END_REF] or [START_REF] Steefel | Approaches to modeling of reactive transport in porous media[END_REF].

We consider here aqueous and precipitation-dissolution reactions, so that species can be in solution or minerals. For the sake of simplicity, we do not include reactions with fixed species other than minerals.

We first describe the mathematical model, then we show that we can consider aqueous species with positive concentrations. Therefore we can introduce an equivalent model with logarithms.

Mass action laws and mass balance equations

We consider here a diluted aqueous solution, so that the activity of water is equal to 1 and the activity of an aqueous species is equal to its molar concentration. We also assume that a precipitation-dissolution reaction involves only one mineral.

A model of this type is proposed in the MoMaS benchmark [START_REF] Carrayrou | Reactive transport benchmark of MoMaS[END_REF][START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problemsthe MoMaS benchmark case[END_REF]. These assumptions are also done in some numerical models of reactive transport [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Kräutle | The semismooth Newton method for multicomponent reactive transport with minerals[END_REF].

Aqueous species can be split through a Morel table [START_REF] Morel | Principles and application of aquatic chemistry[END_REF]. Variables of the problem are then the concentrations of the primary aqueous components c ∈ R Nc + , those of the secondary aqueous species α ∈ R Nα + , and the quantities of minerals p ∈ R Np + . All of them are non-negative quantities.

Without loss of generality, the stoichiometric coefficients can be chosen as integers. Let S ∈ Z Nα×Nc be the stoichiometric matrix for aqueous reactions associated with the positive constants of reaction κ α ∈ R Nα , κ α > 0.

The law of mass action is directly used to model chemical equilibrium [START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF]. Here, because aqueous species can vanish, we write the law of mass action in a product form [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF][START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF]:

α i (c) = κ αi k∈{1,...,Nc:S ik =0} c S ik k , i = 1, ..., N α .
Remark 1 Since for all k, c k ≥ 0, then for all i, α i ≥ 0. Moreover, if, for some i S ik < 0, then the law of mass action requires c k > 0.

Precipitation-dissolution reactions cannot be treated in the same way, since the activity of a mineral is constant. These are heterogeneous reactions, which are represented with a matrix E ∈ Z Np×Nc of integer stoichiometric coefficients and a vector of solubility products denoted by κ p ∈ R Np with κ p > 0. Appearance and disappearance of the mineral j, j = 1, . . . , N p , are governed by a saturation threshold γ j (c) defined by

γ j (c) = k∈{1,...,Nc:E jk =0} c E jk k , j = 1, ..., N p .
So, there are two possible cases, either the fluid is undersaturated and there is no mineral, p j = 0, κ pj -γ j (c) > 0, j = 1, ..., N p , or the fluid is saturated and the mineral is present,

p j ≥ 0, κ pj -γ j (c) = 0, j = 1, ..., N p .
By convention we use a strict inequality in the undersaturated case to avoid any confusion. This kind of equilibrium problem is well-known in optimization theory under the name of complementarity problem, [START_REF] Cottle | The Linear Complementarity Problem[END_REF]. This problem may be reformulated as

p ≥ 0, κ p -γ(c) ≥ 0, p T (κ p -γ(c)) = 0.
The model is completed by a mass balance equation, written for each primary aqueous component. Let T ∈ R Nc be the vector of total concentrations, which is given data in the context of a closed geochemical system. The mass balance equation is then written

T = c + S T α(c) + E T p.
Finally, combining these three equations gives the geochemical model:

       c + S T α(c) + E T p -T = 0, p T (κ p -γ(c)) = 0, p ≥ 0, κ p ≥ γ(c), c ≥ 0, and c k > 0 ∀k ∈ Z, (1)
with Z = {k : 1 ≤ k ≤ N c , ∃i S ik < 0 or ∃j E jk < 0}.
This system is a closed system of (N c + N p ) equations with inequality constraints. With our assumptions the coefficients are integer, so that all the expressions are polynomials. This observation will be used in Section 4.3.

The Gibbs energy G(c, α, p) is defined by chemical potentials for c ≥ 0, α ≥ 0, p ≥ 0. The usual gradient of G uses logarithms and requires c > 0, α > 0 [START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF]. If G is defined in the domain c > 0, α > 0, then equations [START_REF] Migot | Gas phase appearance and disappearance as a problem with complementarity con-straints[END_REF], where the constraint c ≥ 0 is replaced by c > 0, are nothing else than the KKT conditions applied to the minimization of the Gibbs energy, G(c, α, p), subject to the mass conservation law and the constraint p ≥ 0 [START_REF] Leal | Computational methods for reactive transport modeling: a Gibbs energy minimization approach for multiphase equilibrium calculations[END_REF]. Since G is convex and the constraints are linear, the KKT conditions are necessary and sufficient provided that the constraints are feasible.

A logarithmic model

In order to use log(c), we have to study systems where an aqueous species could vanish. We first get rid of the particular case of a tracer, which appears for example in the MoMaS benchmark.

When a stoichiometric coefficient S ik is zero, the component c k does not participate in the reaction involving the secondary species α i (c). Similarly, it does not participate in the reaction with the mineral p j if E jk = 0. A particular case arises when all the coefficients associated to a component are zero, which does not react with any other species. Such a component is called a tracer. We already noticed in Remark 1 and in ( 1) that if at least one stoichiometric coefficient S ik or E jk is negative, then c k must be positive.

We now consider the case where all coefficients S ik and E ik of a component c k are non-negative. If the total T k of a component c k is zero, then this component is absent, as well as the secondary species and minerals associated. Such cases appear in the MoMaS benchmark.

Proposition 2 We assume that ∀i, S ik ≥ 0, ∀j, E jk ≥ 0 and that the k th mass balance equation is satisfied. Then, T k ≥ 0. Moreover, if T k = 0, then c k = 0, and S ik α i (c) = 0, i = 1, . . . , N α , E jk p j = 0, j = 1, . . . , N p .

Proof By definition, the k th mass balance equation is written

T k = c k + i∈{1,...,Nc:S ik =0} S ik α i (c) + j∈{1,...,Nc:E jk =0} E jk p j .
The conclusions follow as the right-hand side is a sum of non-negative terms.

Remark 2 With the assumptions of Proposition 2, the model has no solution if T k < 0.

Proposition 3 We assume that ∀i, S ik ≥ 0, ∀j, E jk ≥ 0. We also assume that (c, p) is solution of (1). If T k > 0, then c k > 0.
Proof Note that (c, p) solution of (1) implies that the k th mass balance equation is satisfied and so T k ≥ 0 by the previous proposition.

We prove that if c k = 0, then T k = 0. First, we notice that c k = 0 implies ∀i, S ik α i (c) = 0 and ∀j, E jk γ j (c) = 0. Now, if E jk > 0, then γ j (c) = 0 < κ pj thus p j = 0 by the complementarity condition in [START_REF] Migot | Gas phase appearance and disappearance as a problem with complementarity con-straints[END_REF]. Finally, ∀j, E jk p j = 0, thus T k = 0. Therefore, any tracer and any component associated to non-negative coefficients with a null total can be removed from the system. Moreover, the secondary species or minerals associated with this component can also be removed from the model. Indeed, the reactions involving a primary component with a zero concentration can be removed from the system.

All in all, it is possible to define a reduced model where all the concentrations c k of the primary components are positive.

Proposition 4 Model (1) is equivalent to a reduced model and a trivial model, where in the reduced model:

-There is no tracer: ∀k, ∃i, S ik = 0 or E ik = 0.

-All the components are reactive: ∀k, c k > 0.

-All the secondary species are reactive: ∀i, α i (c) > 0.

-All the saturation thresholds are positive: ∀i, γ i (c) > 0.

and, in the trivial model:

-For a tracer k: c k = T k . -For a component k with ∀i, S ik ≥ 0, ∀j, E jk ≥ 0 and T k = 0: c k = 0 and ∀i, S ik α i (c) = 0,∀j, E jk p j = 0.
Proof The result follows by combining Remark 1 with propositions 1, 2 and 3.

From now on, we consider such a reduced model, where the number of primary components is still denoted by N c and the positivity constraint is now c > 0:

   c + S T α(c) + E T p -T = 0, p T (κ p -γ(c)) = 0, c > 0, p ≥ 0, κ p ≥ γ(c).
(2) Finally, we can introduce logarithmic variables and define an equivalent model. In this result and through the rest of the paper the logarithm and exponential of a vector are taken componentwise, that is log(z) = (log(z i )) 1≤i≤n and exp(z) = (exp(z i )) 1≤i≤n for all z ∈ R n .

Proposition 5 Model (2) is equivalent to the comple- mentarity problem    exp(x) + S T exp(log(κ α ) + Sx) + E T p -T = 0, p T (Ex -log(κ p )) = 0, p ≥ 0, Ex ≤ log(κ p ), (3)
where x = log(c).

It can be noted that some conditions on T are necessary to ensure the existence of a solution.

Corollary 1 If a primary component k has only non negative stoichiometric coefficients and if T k ≤ 0, then system (2) has no solution.

Proof If the stoichiometric coefficients associated to the primary component k are non-negative, the corresponding total mass (exp(x) + S T exp(log(κ α ) + Sx) + E T p) k is strictly positive. Thus, if there exists a solution, the total T k satisfies T k > 0.

Example with two salts: sodium chloride and potassium chloride

We conclude this section by giving a first example of a geochemical system in order to illustrate and highlight some characteristics. Chemical reactions are only precipitation-dissolution in order to simplify the presentation. Therefore, there are as many reactions as minerals and the system is defined by the stoichiometric matrix E and the solubility products κ p .

This first example, with two salts, considers precipitation of sodium chloride N aCl and potassium chloride KCl, involving the three aqueous components, sodium ion N a + , potassium ion K + and chloride ion Cl -. The two reactions, along with the solubility products, are summed up in the following table:

Mineral Reaction κ p N aCl N aCl N a + + Cl -37.5837 KCl KCl K + + Cl - 7.6208
The stoichiometric matrix of this system is:

E 1 0 1 0 1 1 , which is of full rank N p .
Since all the coefficients are positive, the data T must satisfy T > 0.

Equivalence with an Optimization Problem

We already noticed that system (1) is the first-order optimality condition for minimizing the Gibbs energy subject to the mass balance law and non-negativity of the mineral quantities.

In this section, we prove that the geochemical system (3) is the first-order optimality condition of another minimization problem, with an objective function and inequality constraints of reduced dimensions. This result is the keystone for studying existence and uniqueness of solutions of the geochemical model. Furthermore, this study is directly related to the totals T , which are the parameters of the mass balance law. Another approach could use variational analysis tools, for instance from [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF].

Objective function and inequality constraints

We first define the objective function of the optimization problem for x ∈ R Nc by f (x) = e T Nc exp(x) + e T Nα exp(log(κ α ) + Sx) -T T x, (4) where e n is the vector of dimension n whose components are all ones.

The function f is C ∞ in R Nc with a gradient given by ∇f (x) = exp(x) + S T exp(log(κ α ) + Sx) -T, and a Hessian matrix given by [START_REF] Brunner | Multiphase multicomponent flow in porous media with general reactions: efficient problem formulations, conservative discretizations, and convergence analysis[END_REF] where D(y) denotes a diagonal matrix with the elements of y on the diagonal.

∇ 2 f (x) = D (exp(x)) + S T D (exp(log(κ α ) + Sx)) S,
Then, we define inequality constraints by denoting the feasible set:

F = {x ∈ R Nc , g(x) ≤ 0},
where g(x) = Ex -log(κ p ). Clearly, the feasible set is convex.

A typical assumption that we will use in the sequel to guarantee non-emptiness of F is to require E to be a full rank matrix.

Lemma 1 Assume that E is of full rank. Then, F is non-empty.

Proof For E a full rank matrix, there exists a solution of the equation Ex = log(κ p ).

Finally, we introduce the convex minimization problem min

x∈F f (x). ( 6
)
Remark 3 The objective function f is based on mass balance and constraints are saturation limits, so that problem [START_REF] Buchberger | Ein algorithmus zum auffinden der basiselemente des restklassenrings nach einem nulldimensionalen polynomideal[END_REF] has N c unknowns and N p linear constraints.

On the other hand, minimizing the Gibbs energy G and satisfying mass balance and positivity of minerals involves generally (N c +N α +N p ) variables and (N c +N p ) linear constraints. Reduced dimensions would be useful in numerical solving algorithms.

Uniqueness of a solution

We are now able to prove the equivalence between the model (3) and the minimization problem [START_REF] Buchberger | Ein algorithmus zum auffinden der basiselemente des restklassenrings nach einem nulldimensionalen polynomideal[END_REF]. Indeed, the complementarity problem is nothing else than the Karush-Kuhn-Tucker (KKT) optimality conditions [START_REF] Bonnans | Numerical optimization: theoretical and practical aspects[END_REF]. Uniqueness of a solution immediately follows.

We first show that the objective function is strictly convex:

Lemma 2 The Hessian matrix of f , denoted ∇ 2 f (x), is symmetric positive definite for any x and the function f is strictly convex in R Nc .

Proof The elements of the diagonal matrices in the Hessian matrix given in (5) are positive. Any matrix of the form S T DS, where D is a diagonal matrix with positive elements, is symmetric positive semi-definite, while the diagonal matrix exp(x) is symmetric positive definite (spd). Thus, for all x, the matrix ∇ 2 f (x) is also spd.

The strict convexity of f follows as a consequence.

Remark 4 It can be shown, with some assumptions or a change of variables, that the Gibbs energy is also strictly convex, see for example [START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF][START_REF] Brunner | Multiphase multicomponent flow in porous media with general reactions: efficient problem formulations, conservative discretizations, and convergence analysis[END_REF][START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF]. Here, the strict convexity of f is straightforward.

From this lemma, and using the linearity of the constraints, it becomes easy to prove that the equations (3) are the KKT conditions of the minimization problem [START_REF] Buchberger | Ein algorithmus zum auffinden der basiselemente des restklassenrings nach einem nulldimensionalen polynomideal[END_REF] and to prove the uniqueness of a solution.

Theorem 1 Assume that the feasible set is non-empty, i.e. F = ∅. Then, equations (3) are equivalent to solving the minimization problem [START_REF] Buchberger | Ein algorithmus zum auffinden der basiselemente des restklassenrings nach einem nulldimensionalen polynomideal[END_REF]. Moreover, if the minimization problem has a solution, it is unique.

Proof We first note that any x ∈ F satisfies a linear constraint qualification. Since the objective function f is convex, the feasible set is non-empty convex and a linear constraint qualification is satisfied at any feasible point, the KKT conditions are necessary and sufficient. Using the gradient of f and the Jacobian E T of g, it appears that these KKT conditions are exactly equations (3). Thus, x is a solution of (6) if and only if there exists a Lagrange multiplier p ∈ R Np such that x and p are solutions of (3).

Uniqueness of x is ensured by the strict convexity of f , [START_REF] Kyparisis | On uniqueness of Kuhn-Tucker multipliers in nonlinear programming[END_REF], proved in Lemma 2.

If we further assume that the matrix E is of rank N p , the Lagrange multiplier is also unique.

Corollary 2 Assume that the matrix E is of rank N p . If the minimization problem has a solution, then it is unique and the Lagrange multiplier p is unique.

Proof Note that the feasible set is non-empty here thanks to Lemma 1.

The gradients of the inequality constraints are the lines of the matrix E and are linearly independent if E is of rank N p . Thus, any x ∈ F satisfies a linear independence constraint qualification (LICQ), so that if there is a solution, then the Lagrange multiplier p is unique. [START_REF] Kyparisis | On uniqueness of Kuhn-Tucker multipliers in nonlinear programming[END_REF].

Consequently, Corollary 2 guarantees uniqueness of the solution of the geochemical system (3), whenever such solution exists.

Existence of a solution

In Corollary 1, we have shown that, with non-negative stoichiometric coefficients, the existence of a solution requires that T > 0. Now, we prove that, for any full rank stoichiometric matrix E, there exists a solution of (6), provided that T > 0.

Theorem 2 Assume that the matrix E is of rank N p , and that T > 0. Then, the geochemical model ( 2) has a unique solution.

Proof Let F (x) = e T Nc exp(x) -T T x, x ∈ R Nc , clearly ∀x ∈ R Nc , f (x) > F (x). Now, ∀k = 1, . . . , N c , lim x k →+∞ (exp(x k ) -T k x k ) = +∞; since, by assumption T k > 0, we also have ∀k = 1, . . . , N c , lim x k →-∞ (exp(x k ) -T k x k ) = +∞.
Thus the function F is coercive, in other words lim x →+∞ F (x) = +∞, and it follows that the function f is also coercive.

This property ensures the existence of a solution of the minimization problem, since ( 6) is a problem of minimizing a continuous coercive function over a closed set. The result follows using Corollary 2, which ensures uniqueness of the solution of (2).

Remark 5 Some assumptions are also required to prove the existence of a solution which minimizes the Gibbs energy. Here, these assumptions are directly related to the totals T .

We denote by T the set of totals T such that the problem has a solution. This set is then easy to define when all the stoichiometric coefficients are positive.

Corollary 3 Assume that S ≥ 0 and the matrix E is of rank N p with E ≥ 0. Then, the geochemical model (2) has a solution if and only if T > 0:

T = {T ∈ R Nc , T > 0}
Proof The result follows from Corollary1 and Theorem 2.

However, in the general case, it seems quite challenging to find the set T , as illustrated by the following example.

Example with two minerals: calcite and gypsum

We consider here an example where one of the stoichiometric coefficient is negative. The two reactions involve two minerals, calcite CaCO 3 and gypsum CaSO 4 , with three aqueous components, hydrogen ion H + , calcium ion Ca 2+ and sulfate ion SO 2- 4 . We simplify the mathematical model with the following assumptions: the activity of carbon dioxide CO 2 is equal to 1, and the solubility product of calcite is equal to 1.

Mineral Reaction

κp CaCO 3 CaCO 3 + 2H + Ca 2+ + CO 2 + H 2 0 1 CaSO 4 CaSO 4 Ca 2+ + SO 2- 4 4 10 -5
The stoichiometric matrix E is given by

E = -2 1 0 0 1 1 , which is of rank N p .
The following result provides necessary and sufficient conditions on T for the existence of a solution.

Proposition 6

The set T is given by

T = {(T 1 , T 2 , T 3 ) ∈ R 3 , T 2 > 0, T 3 > 0, T 1 + 2T 2 > 0}.
Proof Mass balance equations and positivity constraints imply the necessary conditions T 2 > 0, T 3 > 0 and

T 1 + 2T 2 > 0.
Now, we show that they are sufficient. Recall that f (x) = exp(x 1 ) + exp(x 2 ) + exp(x 3 ) -

T 1 x 1 -T 2 x 2 -T 3 x 3 , and that g 1 (x) = -2x 1 + x 2 - log(κ 1 ), g 2 (x) = x 2 + x 3 -log(κ 2 ).
For k = 2 and k = 3, the functions exp

(x k ) -T k x k are coercive since T k > 0. If x 1 is finite or x 1 → +∞, then lim |x2|→+∞ f (x) = +∞. Now, since g 1 (x) ≤ 0 we have -x 2 ≥ -2x 1 -log(κ 1 ) thus, since T 2 > 0, also -T 2 x 2 ≥ -2T 2 x 1 -T 2 log(κ 1 ) and finally f (x) ≥ exp(x 1 )+exp(x 2 )+exp(x 3 )-T 3 x 3 - (T 1 + 2T 2 )x 1 -T 2 log(κ 1 ). Since T 1 + 2T 2 > 0, for any (x 2 , x 3 ), we get lim x1→-∞ f (x) = +∞.
We conclude that the function f is coercive and that the minimization problem has a solution.

Precipitation Diagram

Throughout this section, we assume that the matrix E is of full rank. Then, for any T ∈ T , there is a unique solution of the chemical model, which is the speciation of the geochemical system associated with the data T . It gives the positive concentrations of the primary components (the unknown c or x = log(c)) and the nonnegative quantities of the minerals (the Lagrange multiplier p). Now, we propose to characterize this speciation by the presence or lack of each mineral. Then, the aim is to find the subset of totals T leading to a speciation with a given subset of precipitated minerals. We build a partition of T into so-called mineral states, associated to subsets of minerals. This partition defines a diagram, in the same spirit as thermodynamic diagrams [START_REF] Connolly | An automated strategy for calculation of phase diagram sections and retrieval of rock properties as a function of physical conditions[END_REF][START_REF] De Capitani | The computation of equilibrium assemblage diagrams with Theriak/Domino software[END_REF].

This diagram can be formally defined by using the notions of strongly active or weakly active constraints. We first give some definitions based on the optimization problem, then propose a procedure to compute the diagram using symbolic calculus.

Definitions

For T ∈ T , we denote by (x(T ), p(T )) the unique solution of the minimization problem [START_REF] Buchberger | Ein algorithmus zum auffinden der basiselemente des restklassenrings nach einem nulldimensionalen polynomideal[END_REF] with the data T . We characterize the geochemical system by the minerals, introducing the subset M (T ) ⊂ {1, 2, . . . , N p } of precipitated minerals defined by

M (T ) = {i, p i (T ) > 0}. ( 7
)
Conversely, let I be a subset of {1, 2, . . . , N p } and Ī its complement set. We call a mineral state associated to I the subset of T such that, for any T in this subset, M (T ) = I.

Definition 1 A mineral state associated to a subset I of {1, 2, . . . , N p } is defined by

M I = {T ∈ T : M (T ) = I},
where M (T ) is defined by [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of momas with specy[END_REF]. In other words, the constraints with indices in the subset I are strongly active:

This definition builds a partition of the space T into at most 2 Np mineral states, which we call a precipitation diagram.

Proposition 7

The non-empty mineral states build a partition of T . We call this partition a precipitation diagram.

Proof By definition of M (T ), it holds ∪ I M I = T . Moreover, for all I, J such that M I and M J are nonempty, then

I = J =⇒ M I ∩ M J = ∅.
The mineral states are interfaced by what we call state boundaries, which are obtained by removing a mineral which is dissolved and still at saturation. Therefore, each mineral state is associated to at most N p boundaries and there are at most N p 2 Np-1 non-empty boundaries (since each of them is defined by one given mineral and a subset of the other N p -1 minerals).

Definition 2 A state boundary of the mineral state associated to I and j ∈ I is defined by

M I\{j} ∩ {T ∈ T : p j (T ) = g j (x(T )) = 0}.
In other words, the constraint of index j is weakly active and those of other indexes in I are strongly active.

The intersection of boundaries is a special state where all minerals are dissolved and still at saturation. We call it the critical state.

Definition 3

The critical state is the subset C defined by

C = {T ∈ T : ∀i ∈ {1, 2, . . . , N p }, p i (T ) = g i (x(T )) = 0}.
In other words, all constraints of the critical state are weakly active.

Example with two salts: sodium chloride and potassium chloride

Let us illustrate these definitions by the example with two salts of Section 2.3. Here, all the totals must be strictly positive to ensure a unique solution.

We further assume that the species satisfy electrical neutrality. Using the charge of each ion, the balance equation can be written c 1 +c 2 -c 3 = 0, thus c 3 = c 1 +c 2 can be eliminated from the equations.

By combining the mass balance equations, we also get the property T 1 + T 2 -T 3 = 0, thus we can restrict the study to only two parameters T 1 , T 2 .

The 2D precipitation diagram, with T 1 > 0 and T 2 > 0, is shown in Figure 1.

The space is divided into four states, with four boundaries. The critical set is the point (T 1,c T 2,c ):

T 1,c = κ 1 / √ κ 1 + κ 2 , T 2,c = κ 2 / √ κ 1 + κ 2 .
The four state boundaries are:

T 1 = T 1,c , T 2 = T 2,c , T 1 (T 1 + T 2 ) = κ 1 , T 2 (T 1 + T 2 ) = κ 2 .
Here, it is easy to compute by hand the solutions and to characterize the diagram. However, for more complex geochemical systems, the task becomes tricky.

Computing a precipitation diagram

We propose to follow a combinatorial approach exploring the 2 Np possible mineral states and using symbolic computations. This framework is tractable for a moderate number of minerals, say a dozen. A by-product of this procedure will be to get a symbolic expression of the solution (c, p), but the main objective is to build the precipitation diagram.

The combinatorial approach should be more robust than fixing a mesh-grid on the domain of T and repeatedly computing numerically the mineral state. Indeed, as can be seen in Figure 1, it might be easy to miss a phase if the mesh size is not large enough. In Figure 3, it might also be difficult make the difference between a zone without solution and a numerical failure. On the other hand, it might become very costly to refine the mesh near the boundaries.

We come back to the non-logarithmic system (2). Because we consider integer stoichiometric coefficients, equations are polynomial equalities and inequalities, which define semi-algebraic sets. Therefore, solving means finding solutions depending on algebraic numbers, which are roots of univariate polynomials.

In order to simplify the symbolic computations, we can take advantage of some invariants, such as the charge balance described for the example with two salts.

Invariants with conservative variables

Since we assume that the matrix E is of full rank N p , the QR factorization of the matrix E T yields an orthogonal matrix (Q 1 , Q 2 ) ∈ R Nc×Nc and a non-singular triangular matrix R ∈ R Np×Np such that

E T = (Q 1 , Q 2 ) R 0 = Q 1 R.
It can be noted that EQ 1 = R T and EQ 2 = 0. It follows that p can be written as

p(c) = R -1 Q T 1 (T -c -S T α(c)). (8)
Also, the variables Q T 2 c are decoupled from minerals thanks to the relation

Q T 2 T = Q T 2 (c + S T α(c)). (9)
Finally, these simplifications lead to a model involving only the N c variables c:

   Q T 2 (c + S T α(c)) = Q T 2 T, p(c) T (κ p -γ(c)) = 0, c > 0, p(c) ≥ 0, κ p ≥ γ(c). (10)
The quantities Q T 2 c are called conservative variables, because they are often related to some balance equation. In many cases, these equations can be used to eliminate some unknowns, in order to simplify the precipitation diagram. A similar transformation was used in [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF] to simplify the system and prove strict convexity of the Gibbs energy.

An example is the charge balance, when the primary components are ions. The corresponding vector q in Q 2 represents the charges of these ions. In order to ensure electrical neutrality, it is common to require that q T T = 0. This invariant can be used to remove a parameter T k and to eliminate the unknown c k .

Symbolic computations

Let M I with I ⊂ {1, 2, . . . , N p } be a given mineral state and let T ∈ T be a parameter of the symbolic procedure. Then T ∈ M I if and only if (c(T ), p(T )) satisfies: [START_REF] Connolly | An automated strategy for calculation of phase diagram sections and retrieval of rock properties as a function of physical conditions[END_REF] where p(c) is given by [START_REF] Carrayrou | Reactive transport benchmark of MoMaS[END_REF].

               Q T 2 (c + S T α(c)) -Q T 2 T = 0, γ i (c) = κ pi , ∀i ∈ I, p i (c) = 0, ∀i ∈ Ī, c > 0, p i (c) > 0, ∀i ∈ I, γ i (c) ≤ κ pi , ∀i ∈ Ī,
We recall that this system can usually be simplified by using conservative variables and balance equations.

In this form, it is a set of N c polynomial equations with N c positivity constraints, and N p polynomial inequalities. We first solve the N c polynomial equations and get solutions in terms of algebraic numbers, which are function of T . We keep only those such that c > 0.

Then, these computations allow to simplify the polynomial inequalities that now depend on the algebraic numbers computed in the previous step. Thus, we get an implicit description of the state boundaries in function of T .

Different approaches exist to obtain a description of the solutions of a set of polynomial equations. The usual way of representing this set of solutions is through a triangular system of equations, sometimes called regular chains. A different approach consider the rational univariate representation of F. Rouillier [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF]. Several methods exist to compute a triangulation of a polynomial system, probably the most used one is to compute Gröbner bases [START_REF] Sterk | Some tapas of computer algebra, volume 4 of Algorithms and Computation in Mathematics[END_REF][START_REF] Faugere | Efficient computation of zero-dimensional Gröbner bases by change of ordering[END_REF]. This tool is in general useful as it also determines if the system is inconsistent and if it has a finite or an infinite number of solutions. The oldest and most well-known method to compute such a basis is the Buchberger's algorithm [START_REF] Buchberger | Ein algorithmus zum auffinden der basiselemente des restklassenrings nach einem nulldimensionalen polynomideal[END_REF], but more recent and efficient methods were developed, for instance Faugère's F4 and F5 algorithms [START_REF] Faugere | A new efficient algorithm for computing Gröbner bases (F 4)[END_REF]. Computation of a Gröbner basis requires in the worst case d 2 O(n) operations where d is an upper bound on the degrees of the polynomials and n the number of variables. Nowadays, most of the computer algebra systems contain routines to do so. So we do not get into the details of these computations.

Numerical Experiments

In this section, we show numerical results with three examples. Precipitation diagrams were computed with a code developed in Maple TM , implementing the symbolic procedure described above.

These diagrams are computed only once for a given geochemical system. Nevertheless, time and memory requirements could be high, regarding the number of mineral states and boundaries. We provide CPU and memory measurements for the last example, with 64 states.

Three salts: sodium chloride, potassium chloride, and carnallite

We consider here an example with three minerals, sodium chloride N aCl, potassium chloride KCl and carnallite KM gCl 3 , reacting with four aqueous components, sodium ion N a + , potassium ion K + , magnesium ion M g 2+ and chloride ion Cl -. The three reactions and the solubility products are given in the following table: The stoichiometric matrix E of this system is given by

E =   1 0 0 1 0 1 0 1 0 1 1 3   .
Since the matrix E is of full rank N p and all the coefficients are non-negative, the set of parameters T with a unique solution is given by

T = {T ∈ R Nc , T > 0}.
The kernel of E has one vector q defined by the charges of ions: q = (1, 1, 2, -1) T . Thus, the charge balance is a conservative variable and electrical neutrality requires that T 1 + T 2 + 2T 3 = T 4 . As described above, we use electrical neutrality to remove the total T 4 and the concentration c 4 of chloride ion. The 3D precipitation diagram, with T i > 0, i = 1, 2, 3, is shown in Figure 2. There is only one critical point (T 1,c , T 2,c , T 3,c ). The state with three precipitated minerals is bounded by the three hyperplanes:

T 1 = T 1,c , T 2 = T 2,c + T 3 -T 3,c , T 3 = T 3,c .
More details can be found in [START_REF] Migot | Analyse mathématique de modèles géochimiques[END_REF].

Two minerals: calcite and gypsum

This example was already described in Section 3.4.

Here too, the matrix Q 2 has one vector which is defined by the charge of ions. Electrical neutrality allows to remove the total T 3 and the concentration c 3 of sulfate ion.

The 2D precipitation diagram is shown in Figure 3. The critical set is the point (T 1,c T 2,c ) where T 1,c is the unique positive solution of κ 1 y 3 (1/2+κ 1 y)-κ 2 = 0 and T 2,c = κ 1 T 2 1,c . The four state boundaries are given by T

1 = T 1,c , T 2 = T 2,c + 1/2(T 1,c -T 1 ), T 2 = κ 1 T 2 1 , T 2 (1/2T 1 + T 2 ) = κ 2 .

Acid mine drain

We now give an example coming from the literature of acid mine drain, [START_REF] Molson | Reactive transport modelling of acid mine drainage within discretely fractured porous media: Plume evolution from a surface source zone[END_REF]. An acid mine drain (AMD) is an acid mineral solution which appears in some kind of mines or in the context of waste storage. The presence of AMD leads to acidification of the water and appearance of sulphuric acid (H 2 SO 4 ). Figure 4 gives the various species in the initial porous media. We are interested in the precipitation-dissolution reactions of Figure 5. The solubility products come from the database MINTEQ/A2.

In order to simplify the chemistry, we consider only the minerals and the primary species. So, we get a stoichiometric matrix E with three independent systems, .

Even if the polynomial equations possess a favorable structure, there are still a high number of states (2 6 = 64) to explore. In the We observe that our symbolic software was able to compute the precipitation diagram quite quickly, with low memory requirements. The time and memory required to load the library of functions used in the computation are here to make a comparison point. Since the generation of the reduced chemical system involves only very simple operations, it requires almost no time and memory. Computing each mineral state with its boundaries out of 64 subsets is the most expensive step, as expected.

Conclusion and Perspectives

The study presented in this article shows some of the main characteristics of a geochemical system, considering a diluted solution with aqueous and precipitationdissolution reactions. We introduced a new minimization problem, with a strictly convex objective function. It was then straightforward to show that under mild classical assumptions the geochemical system has at most one solution. We also discussed existence of this solution with conditions on the totals in mass balance equations, which are here parameters of the model. We defined a set of totals such that the geochemical system has a unique solution.

Then, we introduced a precipitation diagram which provides a partition of this set into the various possible states of mineral reactions. We gave some characterization of this diagram and described a symbolic procedure to compute it, based on the polynomial structure of the problem. Numerical results using a code developed in Maple illustrated some examples of precipitation diagram and showed that the whole computation is actually reasonable in time and memory.

Reactive transport modeling is one of the main applications using the geochemical model. It is clear that the computation of a precipitation diagram cannot be included in a full code for this application since in this case the non-linear equation has to be solved millions of time. However, it can be very useful in describing the chemical reactions and so it can be used to find a good initial point for a numerical code or to compute the exact solution in a verification process. These extensions to reactive transport and applications of the precipitation diagram are left for further research.

Our results currently apply to diluted solutions, with a simplified model of activities. An important perspective is to extend them in a more general context, where activities are functions of the quantities of species.

Proposition 1

 1 Let c k ≥ 0 be an aqueous tracer: ∀i, S ik = 0 and ∀j, E jk = 0. If the k th mass balance equation is satisfied, then T k ≥ 0 and c k = T k . Proof The k th mass balance equation is reduced to T k = c k and c k ≥ 0.

Fig. 1 :

 1 Fig. 1: Precipitation diagram with the two salts N aCl, KCl.

Fig. 2 :

 2 Fig. 2: Precipitation diagram with the three salts N aCl, KCl and KM gCl 3 .

Fig. 3 :

 3 Fig. 3: Precipitation diagram with calcite and gypsum (simplified system).

Fig. 4 :

 4 Fig. 4: Chemical species in the Acid Mine Drain test case from [29].

Fig. 5 : 4 c

 54 Fig. 5: Mineral reactions and constants in the Acid Mine Drain test case from [29].

  Table below, we provide the time and memory required for the different steps of computations.

	Steps	Time (in sec.)	Memory (in MB)
	Loading function library	0.96	33.93
	Generation of the reduced	< 1	< 1
	chemical system
	Computation of every phase	2.14	16.00
	Total	3.18	49.93


Claire Thérial

Francoise S Lucas

Ludivine Lesage

Vincent Rocher

Alexandre Goncalves

Pierre Servais

Jean-Marie Mouchel

Faecal indicator bacteria in combined sewer overflows and in the Seine River

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

Introduction Material and methods

Results

Conclusions

Faecal indicator bacteria in combined sewer overflows and in the Seine River.

Claire Thérial 1 , Françoise Lucas 1 , Ludivine Lesage 2 , Vincent Rocher 3 , Alexandre Goncalves 3 , Pierre Servais 4 and Jean-Marie Mouchel 2 .

Claire Therial: claire.therial@u-pec.fr During dry periods there was a strong daily cycle for turbidity, TSS, NH 4 + (data not shown) and IE concentrations in raw water, a weaker one for EC concentrations, and no visible daily cycle for conductivity (data not shown). There was also seasonal variations in the raw water between summer and autumn for NH 4 + and FIB (data not shown).

For wet periods, dilution phenomenon with runoff water explained the variation of FIB densities, NH 4 + concentration and conductivity in the CSOs, but not the variation of TSS and turbidity (data not shown). Additional sources of FIB seemed to occur during storm events since there was no significant correlation between FIB and NH 4 + concentrations (data not shown). This suggests that FIB could be brought by the rain runoff and/or eroded from sewer network. Furthermore, the CSO had a significant impact on the river water quality concerning NH 4 + and FIB concentrations. However elimination processes of FIB, such as dilution and sedimentation, contributed to lower this impact.

This study should be used to develop a model to predict FIB densities during storm events in surface waters. This model should take into account all the CSOs in the Paris area in order to evaluate at the annual impact of CSOs on the Seine River microbial quality [START_REF] Poulin | Modélisation de la contamination fécale en Seine : impact des rejets de temps de pluie[END_REF]. Impact of CSOs on surface waters:

• CSO impact was significant at Bougival for NH 4 + concentration, conductivity and FIB densities (fig. 56).

• FIB concentrations in surface water increased during rain events, despite the dilution of CSOs in the river water and the sedimentation of TSS and FIB.

• There was no correlation between FIB concentrations and turbidity or TSS concentration (data not shown).

Figure 2 :

 2 Figure 2: Map of Ile-de-France region.Seine centre waste water treatment plant raw water, . sewer outlets Clichy and La Briche CSO,,,,,,, Bougival, surface water.

Figure 1 :Figure 5 :

 15 Figure 1: The River Seine at Bougival The annual volume of sewage water passing through Paris drainage network is 900 millions m 3 of which 25 millions flows without treatment in the Seine River during wet weather. A quarter of these combined sewer overflows (CSO) occurs at Clichy and La Briche outfalls. The importance of these CSOs can degrade the microbiological quality of the Seine River and create sanitary problems.The goal of this study was to measure the concentrations of faecal indicator bacteria (FIB) in the Seine water after CSOs and to estimate their impact on surface water quality. Escherichia coli (EC) and intestinal enterococci (IE), were measured during dry and wet conditions inside the sewer network and in the Seine River at Bougival (fig.1-2).

Figure 6 :

 6 Figure 6: FIB concentration in surface water during dry and wet periods.

1 .

 1 Laboratoire Eau Environnement Systèmes Urbains (LEESU), Université Paris-est Créteil, 2. UMR Sisyphe, Université Pierre et Marie Curie, 3. Département du développement et de la prospective, SIAAP, 4. Ecologie des Systèmes Aquatiques (ESA), Université Libre de Bruxelles. Bougival Seine Centre MUG/EC for E. coli (NF-EN ISO 9308-3) and MUD/SF for enterococci (NF-EN ISO 7899-1)

Figure 3 :

 3 Figure 3: FIB concentration in raw water/CSOs. A: E. coli, B: Intestinal enterococci. Black circles for wet weather, colour circles for dry weather. Black line represents sinusoidal adjustment with the first harmonic. Replicated data until 48h for wet weather.

Dry weather

Wet weather (n= 6) Surface water (n=11) CSO(n=11)

Sampling strategy

FIB concentration using Most Probable Number estimation [START_REF] Jarvis | Reconsideration of the derivation of Most Probable Numbers, their standard deviations, confidence bounds and rarity values[END_REF] CSOs:

• Daily cycles (fig. 3) were observed in raw waters during dry weather for both FIB concentrations depending on human activity.

• During dry weather, the average conductivity of raw water was 107±4 mS/m and lowered during wet weather (fig. 4). Depending on the rain event, FIB concentrations in CSOs also decreased from 2 to 16 fold during wet periods, suggesting a dilution of raw water by runoff water containing low FIB concentrations (fig. 3 and4).

• FIB concentrations were highly variable in both raw water and CSOs (fig. 4).

Analysis Raw water (n=8)

Hour


Aurélie Clodic
email: aclodic@laas.fr

Rachid Alami
email: rachid@laas.fr

Bilge Mutlu

Cordula Vesper
email: vesperc@ceu.hu

Elisabeth Pacherie
email: elisabeth.pacherie@ens.fr

Julie A Shah
email: julie_a_shah@csail.mit.edu

FJA@HRI15 -Towards a Framework for Joint Action

Keywords: I.2.9 [Robotics]: Commercial robots and applications, Manipulators, H.1.2 [User/Machine systems]: Human factors, J.4 [Social and behavioral sciences]: Psychology Experimentation, Design, Human Factors Human-Robot Interaction, Action Theory, Joint Action Theory, Teamwork

The HRI 2015 Workshop "Towards a Framework for Joint Action" is a full-day workshop held in conjunction with the 10th ACM/IEEE International Conference on Human-Robot Interaction, in Portland (USA) on March 2nd 2015. The first edition of the workshop took place at RO-MAN 2014. This workshop aims at bringing together researchers from several disciplines to discuss the development of frameworks for analyzing and designing human-robot joint action. It is meant to provide the opportunity to researchers interested in joint action, roboticists but also philosophers and psychologists, to discuss in depth the topic and to contribute to the elaboration of a framework for human-robot joint action. To achieve this goal, we propose to the community to tackle a COMMON EXAMPLE (as it is sometimes done in robotics planning competition) with the goal to identify the capacities and skills needed for the successful performance of the joint action. This should enable us to build upon each other's experience to further develop ongoing work. The proposed example is described on the workshop website: fja.sciencesconf.org

WEBSITE

fja.sciencesconf.org

OVERVIEW

For more than a decade, the field of human-robot interaction has generated many valuable contributions of interest to the robotics community at large. The field is vast, going all the way from perception (e.g., tactile or visual) to action (e.g., manipulation, navigation) and decision (e.g., interaction, human-aware planning). However, when it comes to the development of future robot assistants or robotic team-mates in mixed human-robot teams, there is a need for a deeper understanding of human-robot joint action that could provide a framework for the different contributions and studies.

It is interesting to observe, from a roboticist point of view, that research on human joint action is a topic of intense research in cognitive psychology and philosophy. This observation lead us to start a multi-disciplinary initiative, successfully launched during a first workshop at the RO-MAN 2014 conference (fja2014.sciencesconf.org) with more than 30 participants. With this second edition, we would like to continue to analyze the fundamental assumptions as well as detailed empirical findings from cognitive psychology and philosophy and connect them with various ongoing research activities in robotics, from the design of Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact Owner/Author. control architectures to human-robot interaction. More specifically, our goal is to bring various aspects of existing work together and examine how they can help us define the kind of integrative framework needed for the design of an autonomous robot that can engage in long-term interaction with a human partner. This framework should be able to serve two complementary purposes. On the one hand, it should help us define with precision what a robot needs to understand about the human it interacts with for the interaction to be successful and thus what capacities the robot should be equipped with to ensure it can build this understanding. On the other hand, the robot also needs to be understood by its human partner and this framework should help us clarify how this understanding operates and what is needed to enable the robot to behave appropriately and in a way that manifests what it is doing to the human partner.

The aim of this workshop series is to bring together researchers in robotics, psychologists and philosophers. This will create a unique opportunity for scientific exchange between these disciplines. In particular, psychologists and philosophers will be given the opportunity to present recent developments in joint action research, while roboticists will be able to discuss the challenges they face with regard to human-robot interaction and more precisely human-robot joint activity. We imagine and hope that this workshop series will create a context in which interested researchers can come together and discuss the development of an integrative framework for research on human-robot joint action like it has been the case for the first edition.

TOPICS

• joint goal / joint action description and modeling

• joint goal establishment and negotiation

• planning when the goal is a joint goal

ORGANIZERS --

Aurélie Clodic (aurelie.clodic@laas.fr)

--Cordula Vesper (vesperc@ceu.hu)

--Rachid Alami (rachid.alami@laas.fr)

--Elisabeth Pacherie (elisabeth.pacherie@ens.fr)

--Bilge Mutlu (bilge@cs.wisc.edu)

--Julie Shah (julie_a_shah@csail.mit.edu)

ACKNOWLEDGMENTS

This workshop is partly supported by the European Commission under contract FP7-ICT-600877 (SPENCER) and ANR-Project ROBOERGOSUM (DECISION ANR-12-CORD-0030).


E Camby
email: ecamby@ulb.ac.be

G Caporossi
email: gilles.caporossi@hec.ca

M H M Paiva

M E V Segatto
email: segatto@ele.ufes.br

Expected distance based on random walks

Keywords: random walks, distance, topological index, graph Mathematics Subject Classification (2000) MSC 05C81, MSC 05C12, MSC 05C21, MSC 05C38

By considering a graph as a network of resistances, Klein and Randić [14] proposed the definition of a distance measure. Indeed, if each edge of the graph represents a resistance of 1Ω, the equivalent resistance of the graph between each pair of vertices may be used as a distance. Based upon random walks in graphs, Stephenson and Zelen [17] built a computational model to find the probability that each edge is used. From a mathematical point of view, both articles are based upon exactly the same model and the link between random walks and the electrical representation was established by Newman [16] when defining an alternative to Freeman's betweenness centrality [9, 10] based upon random walks.

In the present paper, the similitude between these two processes is exploited to propose a new random walks based distance measure that may be defined as the expected length of a walk between any pair of vertices. We call it the expected distance and we prove that it is actually a distance. From this new definition, the RW Index is proposed that sums the expected walks lengths between pairs of vertices exactly in the same way as the Wiener index sums the shortest paths distances or the Kirchhoff index sums the equivalent

Introduction

Let G = (V, E) be a graph without loops or multiple edges on n = |V | vertices and m = |E| edges. All basic definitions can be found in [START_REF] Diestel | Graph theory[END_REF]. A graph G is said connected if there exists a sequence of edges joining any pair of vertices, i.e. there exists a geodesic path between any pair of vertices. Connected graphs with minimum number of edges are called trees. Then, a tree on n vertices has n -1 edges. We consider in this paper only connected graphs. The Wiener index [START_REF] Wiener | Structural determination of paraffin boiling points[END_REF], likely the oldest and most popular topological index, was originally defined for trees as:

(s,t)∈E n s × n t , (1)
where n s is the number of vertices of the connected component containing s if the edge (s, t) was removed, and similarly for n t . In trees, the Wiener index corresponds to the sum of distances between pairs of vertices in G:

W (G) = 1 2
s∈V t∈V

d st = s<t∈V d st , (2)
where d st denotes the distance between s and t, i.e. the number of edges in a geodesic path joining s and t. The definition given in (2) easily extends to general graphs and is commonly used as the definition of the Wiener index.

The concept of distance is a key feature in the design of molecular descriptors. However, the distance based upon the geodesic paths has only integer values and inevitably has the same value for a large number of pairs of vertices. To overcome this weakness, Klein and Randić [START_REF] Klein | Resistance distance[END_REF] proposed an alternative way to define the distance between pairs of vertices. Considering the graph as a network of resistors and each edge being associated to a resistor of 1Ω, they proposed to consider the equivalent resistance r st between s and t as the definition of the distance between s and t. Then, analogously to the Wiener index, the sum of the equivalent resistances among all pairs of vertices corresponds to the Kirchhoff index. In practice, these equivalent resistances can be computed from the flow conservation used for electric circuits. The latter could be described as a set of equations characterized by the Laplacian matrix L = D -A where D is a diagonal matrix of degrees whereas A is the adjacency matrix.

Studying random walks in social networks, Stephenson and Zelen [17] also used a matrix that turns out to be related to the Laplacian matrix. The analogy between random walks and electric circuits was later demonstrated formally by Newman [START_REF] Newman | A measure of betweenness centrality based on random walks[END_REF]. Indeed, considering a unit flow of 1 ampere (1A) between two vertices s and t, the intensity on a given edge e of the network corresponds to the probability that e is used by a random walk from s to t. A general reference to random walks on graphs is given by Doyle and Snell [START_REF] Doyle | Random walks and electric networks[END_REF].

In this paper, we introduce the expected distance between any pair of vertices as the expected length of a path between them based upon the so-defined probabilities in random walks. Instead of only considering the shortest paths between pairs of vertices as the Wiener index does, from this distance, we define a new index based on Random Walks, called RW index.

In Section 2, we present required background on random walks. In Section 3, we define the expected distance and the RW index while we compare the latter to the Kirchhoff and the Wiener indices in Section 4. In Section 5, we propound a decomposition of the RW index based upon edge or vertex centrality measures and the relation with a similar decomposition of the Wiener [START_REF] Caporossi | Centrality and betweenness: vertex and edge decomposition of the Wiener index[END_REF] and Kirchhoff indices. The values of the RW index for some families of graphs, conjectures on their bounds and partial results are presented in Section 6 while Section 7 concludes.

Background

As explained in the introduction, the Wiener index of a graph G is usually defined by the sum of distances between pairs of vertices:

W (G) = 1 2
s∈V t∈V

d st = s<t∈V d st (3)
This index gives an indication on how pairs of vertices are close. Another similar index is the Kirchhoff index. Before defining this index, we give some basic tools.

The mathematical aspects involved in the present paper are related to the use of a graph as representation of an electrical circuit. Suppose a graph G that represents a network of resistors, each edge corresponding to a unit resistor. If an electrical flow of 1A is sent through the network from vertex s to vertex t, this flow will be distributed among the circuit in a way that respects both the flow conservation equations and the electrical properties of the resistors. Indeed, according to the flow conservation rule, the sum of flows from s is 1 as well as the sum of flows to t and the net flow, i.e. the difference between inflow and outflow, associated to each other vertex is always 0. These rules are described by the following equations, where i uv st is the intensity of the flow from u to v when s is the source and t is the sink:

k:(s,k)∈E i sk st = 1 (4) k:(t,k)∈E i tk st = -1 (5)
k:(q,k)∈E

i qk st = 0 ∀q ∈ {1 . . . n} \ {s, t} (6)
The physical representation of the resistor is such that if an intensity i uv st flows through the edge (u, v) from u to v, the potentials d

(u) st and d (v)
st in these vertices will respect the following general equation:

Potential Difference = Resistance × Intensity. ( 7
)
In our case, with the appropriated notation, it means that

d (u) st -d (v) st = i uv st , (8)
since the edge (u, v) has a resistance of 1Ω. A consequence is that

-1 d (u) st -d (v) st 1 ∀(u, v) ∈ E. (9)
The potential d

(u)
st associated to each vertex u could be computed by the resolution of the system of equations ( 4)-( 6) and [START_REF] Feng | Sharp upper bounds for the number of spanning trees of a graph[END_REF]. Actually, this system of equations are represented by a matrix equation, where appears the Laplacian matrix L:

LV = S, (10)
where

V = d (v) st
is the vector of potentials and S is the net flow vector. Accordingly, S has 0 entries at each position except s and t where they are 1 and -1 respectively. The problem of computing the flow on each edge consists in finding the potential vector V in [START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF]. From a technical point of view, since only differences between potentials are considered, adding a constant to each potential will not affect the realizability of the system, which has actually an infinity of solutions. Indeed, 0 is always an eigenvalue of L and an associated eigenvector is the vector V 1 whose all of its components are 1. Accordingly, the matrix L is not invertible. A way to handle this problem and to find a solution in [START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF] is to consider a similar system by adding 1 to each entry of L, as used by Klein and Randić [START_REF] Klein | Resistance distance[END_REF]. Other methods can be used [START_REF] Gutman | Generalized inverse of the laplacian matrix and some applications[END_REF][START_REF] Newman | A measure of betweenness centrality based on random walks[END_REF] to solve [START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF]. The obtained system will then have a unique solution that will also be a solution to the original problem, since the sum of the lines of matrices L and S are both equal to 0.

From a physical point of view, the equivalent resistance r st between s and t is obtained by [START_REF] Entringer | Distance in graphs[END_REF]:

d (s) st -d (t) st = r st × i st = r st (11)
since the global intensity i st from s to t is 1A.

Now, since we have all tools, we can define the Kirchhoff index [START_REF] Klein | Resistance distance[END_REF]: it is the sum of equivalent resistances among all pairs of vertices in a graph G, i.e.

Kf (G) = 1 2 s∈V t∈V r st = s<t∈V r st . (12)
Notice that the equivalent resistance may also be computed using the power dissipation law, which is described by the following general equation:

Dissipated Power = Potential Difference × Intensity. ( 13
)
The dissipated power P uv st on an edge (u, v) between s and t is rewritten:

P uv st = (d (u) st -d (v) st ) × i uv st by (8)
= (d

(u) st -d (v) st ) 2 . ( 14
)
Since the global dissipated power P st between s and t is the sum of dissipated powers on each edge, this global dissipated power is then

P st = (u,v)∈E (d (u) st -d (v) st ) 2 . ( 15
)
From a high point of view, the latter is also, by [START_REF] Kirchhoff | Ueber die auflösung der gleichungen, auf welche man bei der untersuchung der linearen vertheilung galvanischer ströme geführt wird[END_REF],

P st = (d (s) st -d (t) st ) × i st ist=1 = d (s) st -d (t) st by (11)
= r st ,

thus the equivalent resistance r st between s and t is defined by d

(s) st -d (t)
st , as well as

(u,v)∈E (d (u) st -d (v) st ) 2 .

The expected distance and the RW index

Based upon the same model as the electric one, when viewing the network from the random walks point of view, we define the expected length dst of the walk from s to t as the sum of the weighted lengths of paths from s to t. Observe that the weight of a path is its probability, i.e. the absolute value of its intensity. Moreover, the expected length dst can be decomposed by the edges, i.e. dst =

(u,v)∈E |d (u) st -d (v) st |. ( 17
)
As expected, the expected distance defined by the expected length between any pair of vertices is actually a distance, as established in the following theorem.

Theorem 1 Let G = (V, E) be a graph. The expected distance is a distance on V .

Proof By definition in [START_REF] Stephenson | Rethinking centrality: Methods and examples[END_REF], the expected distance is well-defined. Reversing the source and the sink in a flow is imposed to take opposite potentials. Since the expected distance is defined by means of the absolute value of difference of potentials, it is symmetric. Separating distance is ensured because the only case where all potentials are equal is when the source are the sink.

The last property to check is the triangle inequality. Let s, r, t be three vertices in G. We need to prove that dst dsr + drt . By the electric property, we know that if f 1 , respectively f 2 , is a flow through the network G from s to r with potentials d (.) sr on V , respectively from r to t with potentials d

(.) rt , then f = f 1 + f 2 is a flow from s to t with potentials d (.) st = d (.) sr + d (.) rt . (18)
Accordingly, we have dst =

(u,v)∈E |d (u) st -d (v) st | (18) = (u,v)∈E |(d (u) sr + d (u) rt ) -(d (v) sr + d (v) rt )| = (u,v)∈E |(d (u) sr -d (v) sr ) + (d (u) rt -d (v) rt )| (u,v)∈E |d (u) sr -d (v) sr | + (u,v)∈E |d (u) rt -d (v) rt | = dsr + drt . 2
Similarly to the Wiener and the Kirchhoff indices, we define the RW index of a graph G as:

RW (G) = 1 2 s∈V t∈V dst = s<t∈V dst . ( 19
)
4 Comparison of the Wiener, Kirchhoff, and RW indices

The following theorem establishes a comparison of the Wiener, Kirchhoff, and RW indices for all graphs. In addition, it yields that these three indices are the same only in the class of trees.

Theorem 2 Let G = (V, E) be a graph, then we have:

Kf (G) W (G) RW (G). ( 20
)
Moreover, the following statements are equivalent:

(i) Kf (G) = W (G), (ii) W (G) = RW (G), (iii) Kf (G) = RW (G), (iv) G is a tree.
Proof To be self-contained, we proof the first inequality. Let s and t be two vertices of G. From an electric point of view, the equivalent resistance r st is always the length d st of a shortest path from s to t when there is only one path between s and t, otherwise the equivalent resistance is strictly smaller than this value d st . Then, r st d st for every pair of vertices s, t ∈ V . Accordingly, Kf (G) W (G). Now, we show the second inequality. Since the expected length dst is by definition a convex combination of lengths of paths between s and t, dst is always greater than the smallest length of these paths, i.e. dst d st . Therefore, W (G) RW (G).

We focus on the equivalences. First, we assume that G is a tree. By definition, there exists only one path between every pair of vertices. Thus, r st = d st = dst for every pair of vertices s and t. Accordingly,

Kf (G) = W (G) = RW (G).
The last case is when G contains a cycle. Let C be an arbitrary cycle in G with the ordered vertices v 1 , . . . , v k . According to the Kirchhoff's circuit laws, assuming that

v k+1 = v 1 , k i=1 i vivi+1 v1v2 = 0, i.e. i v1v2 v1v2 = k i=2 i vivi+1 v1v2 .
Moreover, i v1v2 v1v2 = 0 since the edge (v 1 , v 2 ) is the shortest path between v 1 and v 2 . Accordingly,

k i=2 i vivi+1 v1v2
= 0, i.e. there exists an edge (u, v) such that i uv v1v2 = 0. It is necessary that this edge is from a path from v 1 to v 2 with non-zero probability. Notice that the length of this path must be greater than d v1v2 = 1. Thus, there are 2 paths from v 1 to v 2 with a non-zero probability and with different length, so d v1v2 < dv1v2 , which implies that W (G) < RW (G), since d st dst ∀s, t ∈ V . Moreover, by the same argument, we have that r v1v2 < d v1v2 . Because for every pair of vertices s, t ∈ V , r st d st , we conclude that Kf (G) < W (G) and then Kf (G) < RW (G). 2

5 Vertex and edge decompositions of Wiener, Kirchhoff, and RW indices

Decompositions of the Wiener index

The computation of the Wiener index may be achieved through two distinct, but related, concepts, that are the transmission and the edge betweenness [START_REF] Anthonisse | The rush in a directed graph[END_REF][START_REF] Freeman | A set of measures of centrality based on betweenness[END_REF][START_REF] Girvan | Community structure in social and biological networks[END_REF]. Recall that the edge betweenness of an edge e is defined [START_REF] Girvan | Community structure in social and biological networks[END_REF] as the number of shortest paths between pairs of vertices that use e. On the one hand, the transmission T v of the vertex v is defined as:

T v = u∈V d uv . (21)
The vertex decomposition of the Wiener index is naturally:

W (G) = 1 2 v∈V T v . (22)
On the other hand, if b uv denotes the edge betweenness associated to the edge (u, v), then Caporossi et al. [START_REF] Caporossi | Centrality and betweenness: vertex and edge decomposition of the Wiener index[END_REF] proved the edge decomposition of the Wiener index, i.e.

W (G) = 1 2 u∈V v:(u,v)∈E b uv = (u,v)∈E:u<v b uv . (23)

Decompositions of the Kirchhoff index

Similarly to the decompositions of the Wiener index through transmission and edge betweenness, one can also decompose the Kirchhoff index through the related concepts of equivalent resistance and power dissipation. As stated by Klein and Randić [START_REF] Klein | Resistance distance[END_REF], the equivalent resistance between any pair of vertices may be considered as a distance measure. It is then possible to compute the resistance transmission RT v of a vertex v as follows:

RT v = u∈V r uv . (24)
The Kirchhoff index then writes by the following vertex decomposition:

Kf (G) = 1 2 v∈V RT v (25)
If P uv the power dissipation associated to each edge (u, v), independently of the source and the sink, is defined by:

P uv = s,t∈V d (u) st -d (v) st 2 , (26)
then the edge decomposition of the Kirchhoff index of a graph G is:

Kf (G) = 1 2 u∈V v:(u,v)∈E P uv = (u,v)∈E:u<v P uv . (27)

Decompositions of the RW index

In the same way as the Wiener and the Kirchhoff indices, the RW index can be decomposed according to vertices and edges. As the transmission is defined by the use of shortest paths, the random walks transmission RW T v of a vertex v is defined as:

RW T v = u∈V duv . ( 28
)
The vertex decomposition of the RW index of a graph G is then

RW (G) = v∈V RW T v . (29)
Moreover, the edge decomposition of the Wiener and the Kirchhoff indices can be also extended to the RW index. Indeed, while the edge betweenness represents the number of shortest paths through pairs of vertices that uses a given edge, the random edge betweenness rb uv of the edge (u, v) represents the sum of the probabilities that the edge (u, v) is used over all pairs of vertices as follows:

rb uv = s,t∈V |d (u) st -d (v) st |. (30)
Accordingly, the RW index of a graph G has the following edge decomposition:

RW (G) = 1 2 u∈V v:(u,v)∈E rb uv = (u,v)∈E:u<v rb uv . (31)

Further remarks

As the adjusted 1 vertex betweenness [START_REF] Caporossi | Centrality and betweenness: vertex and edge decomposition of the Wiener index[END_REF], a vertex related value, was defined by

b v = u:(u,v)∈E b uv , (32)
1 We are interested by the adjusted vertex betweenness, a centrality measure computing, for a given vertex u, the number of pairs of vertices such that u belongs to the shortest path between this pair of vertices. This definition allows that u can be an extremity of the shortest path, at the opposite of the vertex betweeness defined by Freeman [START_REF] Freeman | A set of measures of centrality based on betweenness[END_REF][START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF].

we can define, in the same way, the vertex power dissipation p v and the vertex random betweenness rb v of a vertex v by the following equations:

p v = u:(u,v)∈E p uv , (33)
rb v = u:(u,v)∈E rb uv . ( 34
)
6 Bounds and values of the RW index 6.1 Value of RW index for some families of graphs

In this section, we study the RW index for some families of graphs: paths, cycles, complete graphs and stars. Recall that a star is a tree on n vertices with a vertex of degree n -1.

Theorem 3 We denote by

-C n a cycle on n vertices, -P n a path on n vertices, -S n a star on n vertices, -K n a complete graphs on n vertices.

Then, the following results yield:

(i) RW (C n ) = RW (P n ) = n 6 (n 2 -1), (ii) RW (S n ) = RW (K n ) = (n -1) 2 .
Proof As mentionned in Theorem 2, in the class of trees, the three indices match. Moreover, Entringer, Jackson and Snyder [START_REF] Entringer | Distance in graphs[END_REF] proved that W (P n ) = n 6 (n 2 -1) and W (S n ) = (n -1) 2 , proving two of the four equalities.

We consider the cycle C n on n vertices. Since every vertex has the same role in the graph, we can fix one of them, say s, as a source. Let k be the distance between the source and the sink. By the property of the cycle, there are only two paths between s and t: one of length k and another one of length n -k. So, the potential difference all along each path is constant and inversely proportional to its length, i.e. the shortest path has a probability of n-k n and the longest one k n . Accordingly,

RW (C n ) = 1 2 s∈V t∈V dst = n 2 n-1 k=1 k n -k n + (n -k) k n = n 2 n-1 k=1 2 n (k(n -k)) = n-1 k=1 (k(n -k)) = n 6 (n 2 -1).
The last case is the complete graph K n . We fix arbitrarily a vertex s as a source and another vertex t as a sink, since every pair of vertices has the same feature. Paths of length strictly greater than 2 receive 0 as probability because theses paths contain at least two successive vertices, distinct from the source and the sink, and because these two vertices have an identical role, and so the same potential. Thus, paths with non-zero probability have length 1 or 2. By a similar argument of the previous case, there are one path of length 1 with probability 2 n and n -2 paths of length 2 with probability 1 n . Therefore, We observe that in the class of trees, stars and paths are extremal graphs of the RW index, since it is the case for the Wiener index [START_REF] Entringer | Distance in graphs[END_REF]. More generally, we present the two following conjectures for extremal graphs, found by AutoGraphiX [START_REF] Caporossi | Variable neighborhood search for extremal vertices: The system AutoGraphiX-III[END_REF][START_REF] Caporossi | Variable neighborhood search for extremal graphs: 1 The autographix system[END_REF].

RW (K n ) = 1
Conjecture 1 The graph on n vertices that minimizes the RW index is either the star S n or the complete graph K n .

Conjecture 2

The graph on n vertices that maximizes the RW index is either the path P n or the cycle C n .

If previous conjectures are true, we deduce the following conjecture by Theorem 3.

Conjecture 3 Let G be a graph. Then, (n -1) 2 ≤ RW (G) ≤ n 6 (n 2 -1).

Since the intensity of the edges are also related to the number of spanning trees [START_REF] Kirchhoff | Ueber die auflösung der gleichungen, auf welche man bei der untersuchung der linearen vertheilung galvanischer ströme geführt wird[END_REF][START_REF] Lyons | Probability on Trees and Networks[END_REF], it is interesting to observe that the graphs on n vertices minimizing the number of spanning trees [START_REF] Feng | Sharp upper bounds for the number of spanning trees of a graph[END_REF] are exactly the ones minimizing the RW index, according to the Conjecture 1, and that is not the case for the Kirchhoff index.

Concluding remarks and extensions

A possible extension of the RW index could be on graphs with weights on the edges. Indeed, in this paper, we consider the case when each edge is associated to a resistor of 1Ω. In the case of weighted graphs, we could set a resistor of w e Ω on each edge e with the weight w e .

2 6. 2

 22 Bounds on RW index

This work was partially supported by a travel grant from the F.R.S.-FNRS, by a post-doc grant "Bourse d'Excellence WBI.WORLD" from Fédération Wallonie-Bruxelles (Belgium) and also by NSERC, Foundation HEC Montréal (Canada), CAPES, and CNPq (Brazil).