Gabriel Freitas Oliveira

Renato Markele

Ferreira Candido

Vinicius Mariano Gonçalves

Carlos Andrey Maia

Bertrand Cottenceau

Laurent Hardouin
email: laurent.hardouin@univ-angers.fr

Discrete Event System Control in Max-Plus Algebra : Application to Manufacturing Systems*

Keywords: Petri nets, timed event graphs, max-plus algebra

This paper deals with control of industrial systems which can be depicted by timed event graphs. A methodology and software tools are presented to allow engineers to implement automatically controllers in PLC. Three steps are needed and recalled: the modelling of the system in max plus algebra, the design of the controller and the implementation in a Supervisory Control and Data Acquisition (SCADA) system. An example on a real system is given to illustrate these steps.

INTRODUCTION

The industrial manufacturing systems are frequently event-driven systems and their closed-loop control is worth of interest in order to react to disturbances modifying their behaviours. Among these discrete event systems there exists a subclass of system involving synchronization and delay phenomena for which an efficient control theory has been developed during the last decade [START_REF] Schutter | Analysis and control of max-plus linear discreteevent systems: An introduction[END_REF][START_REF] Maia | On the Model Reference Control for Max-Plus Linear Systems[END_REF][START_REF] Lhommeau | Interval analysis in dioid : Application to robust open loop control for timed event graphs[END_REF][START_REF] Shang | An integrated control strategy to solve the disturbance decoupling problem for max-plus linear systems with applications to a high throughput screening system[END_REF][START_REF] Heidergott | Max Plus at Work: Modeling and Analysis of Synchronized Systems : a Course on Max-Plus Algebra and Its Applications[END_REF][START_REF] Cohen | Max-plus algebra and system theory : Where we are and where to go now[END_REF]. Synchronization phenomena appear when meeting between events is needed, e.g. an event starts when all the preceding events are finished. Delay phenomena appear when transportation operation or manufacturing activities need a duration to be achieved, e.g. the event representing a finishing time of a task is equal to its starting time plus the duration of the operation. Although these phenomena lead to a non-linear model in conventional algebra they admit a linear model formulated in some idempotent semirings, the most popular being the max-plus algebra, where the max operation is the sum of the algebraic structure and the classical addition has to be considered as the product. It must be noticed that these systems and their linear models correspond exactly to the Timed Event Graphs (TEGs) which are a subclass of Petri nets where each place has exactly one upstream and one downstream transition and all arcs have weight equal to 1. The time is associated to places and represents the duration a token has to stay in a place before to contribute to the firing of a downstream transition. This graphical model is popular since suitable for engineer depicting the behaviour of a manufacturing system. This paper proposes a methodology and a specific CAPES-COFECUB software to help the engineers to synthesize and implement a closed-loop control of these systems by tacking advantage of the efficient existing control strategies. By starting from a TEG description of the system and a desired behaviour also given as a TEG named reference model, a method is proposed to obtain the algebraic model of the system and of the reference model, then a software is proposed to give automatically the code to implement in a Supervisory Control and Data Acquisition (SCADA) system. The engineer will have to choose the control strategy which will be implemented among different strategies, namely the disturbance decoupling problem, the model matching problem, the observer-based controller. They have for common point to yield an optimal control law according to the just-in-time criterion which aims to delay as much as possible the occurrence of event input while achieving the specific goal of the control strategy. This criterion is quite usual in industry since it leads to produce only the necessary quantity of raw materials and then to reduce as much as possible the useless stock.

In this paper we will recall the observer-based controller strategy in order to illustrate the methodology used in the software. It is organized as follows: First, the algebraic tools necessary to synthesize the control law are recalled, then in section 3 the method to obtain the algebraic model of the system is presented, this section presents an original description which allows the software to obtain automatically an explicit model. The observer-based control strategy which aims to match a reference model is presented in section 4. The software architecture and the methodology adopted are illustrated in section 5, which presents an implementation on a flexible automated system available at University of Angers.

MATHEMATICAL BACKGROUND

An idempotent semiring S is an algebraic structure with two internal operations denoted by ⊕ and ⊗. The operation ⊕ is associative, commutative and idempotent, that is, a ⊕ a = a. The operation ⊗ is associative (but not necessarily commutative) and distributive on the left and on the right with respect to ⊕. The neutral elements of ⊕ and ⊗ are represented by ε and e respectively, and ε is an absorbing element for the law ⊗ (∀a ∈ S, ε ⊗ a = a ⊗ ε = ε). As in classical algebra, the operator ⊗ will be often omitted in the equations, moreover, a i = a ⊗ a i-1 and a 0 = e. In this algebraic structure, a partial order relation is defined by a b ⇔ a, therefore an idempotent semiring S is a partially ordered set (see [START_REF] Baccelli | Synchronization and Linearity : An Algebra for Discrete Event Systems[END_REF][START_REF] Heidergott | Max Plus at Work: Modeling and Analysis of Synchronized Systems : a Course on Max-Plus Algebra and Its Applications[END_REF] for an exhaustive introduction). An idempotent semiring S is said to be complete if it is closed for infinite ⊕-sums and if ⊗ distributes over infinite ⊕-sums. In particular = x∈S x is the greatest element of S ( is called the top element of S).

Theorem 1. [see Baccelli et al. (1992), th. 4.75] The implicit inequality x ax⊕b as well as the equation x = ax⊕ b defined over S, admit x = a * b as the least solution, where a * = i∈N a i (Kleene star operator).

Definition 1. (Residual and residuated mapping). An order preserving mapping f : D → E, where D and E are partially ordered sets, is a residuated mapping if for all y ∈ E there exists a greatest solution for the inequality f (x) y (hereafter denoted f (y)). Obviously, if equality f (x) = y is solvable, f (y) yields the greatest solution. The mapping f is called the residual of f and f (y) is the optimal solution of the inequality.

Example 1. Mappings Λ a : x → a ⊗ x and Ψ a : x → x ⊗ a defined over an idempotent semiring S are both residuated (Baccelli et al. (1992), p. 181). Their residuals are order preserving mappings denoted respectively by Λ a (x) = a• \x and Ψ a (x) = x• /a. This means that a• \b (resp. b• /a) is the greatest solution of the inequality a ⊗ x b (resp. x ⊗ a b).

The set of n×n matrices with entries in S is an idempotent semiring. The sum, the product and the residuation of matrices are defined after the sum, the product and the residuation of scalars in S, i.e.,

(A ⊗ B) ik = j=1...n (a ij ⊗ b jk ) (1) (A ⊕ B) ij = a ij ⊕ b ij , (2)
(A • \B) ij = k=1..n (a ki • \b kj ) , (B• /A) ij = k=1..n (b ik • /a jk ). ( 3
)
The identity matrix of S n×n is the matrix with entries equal to e on the diagonal and to ε elsewhere. This identity matrix will also be denoted e, and the matrix with all its entries equal to ε will also be denoted ε.

Example 2. (max,plus) algebra:

Z max = (Z ∪ {-∞, +∞}, max, +) is a complete idempotent semiring such that a ⊕ b = max(a, b), a ⊗ b = a + b, a ∧ b = min(a, b) with ε = -∞, e = 0
, and = +∞. The order is total and corresponds to the natural order . By extension Z n×n max is a semiring of matrices with entries in Z max . Matrix ε ∈ Z n×m max will be such that all its entries are equal to ε ∈ Z max , matrix e ∈ Z n×n max will be such that all the entries are equal to ε ∈ Z max except the diagonal entries which are equal to e ∈ Z max . being matrices with entries in Z max . As the max-plus multiplication corresponds to the classical addition, its residual corresponds to conventional subtraction, i.e. , 1 ⊗ x 4 admits the solution set X = {x|x 1• \4} with 1• \4 = 4 -1 = 3 being the greatest solution of this set.

Applying the rules of residuation in max-plus algebra to the relation A ⊗ X B results in:

A• \B = 1• \6 ∧ 3• \7 ∧ 5• \8 2• \6 ∧ 4• \7 ∧ ε• \8 = 3 3 Matrix A• \B = [3 3] T is the greatest solution for X which ensures A ⊗ X B. Indeed, A ⊗ (A• \B) = 1 2 3 4 5 ε ⊗ 3 3 = 5 7 8 6 7 8 = B.
Remark 1. Note that residuation achieves equality if a solution exists. [START_REF] Baccelli | Synchronization and Linearity : An Algebra for Discrete Event Systems[END_REF][START_REF] Hardouin | Control and state estimation for max-plus linear systems[END_REF] is formally the quotient dioid of B [[γ, δ]] (the set of formal power series in two commutative variables γ and δ, with Boolean coefficients and with exponents in Z), by the equivalence relation [γ, δ]]. The representative which is minimal with respect to the number of terms is called the minimum representative.

SYSTEM MODELING

3.1 Dioid M ax in [[γ, δ]] Dioid M ax in [[γ, δ]] (see
xRy ↔ γ * (δ -1 ) * x = γ * (δ -1 ) * y. Dioid M ax in [[γ, δ]] is complete. As M ax in [[γ, δ]] is a quotient dioid, an element of M ax in [[γ, δ]] may admit several representatives in B[
A simple geometrical interpretation of the previous equivalence relation is available in the (γ, δ)-plane. Consider a monomial

γ k δ t ∈ B[[γ, δ]], its south-east cone is defined as {(k , t )|k
k and t t}. The south-east cone of a series in B [[γ, δ]] is defined as the union of the south-east cones associated with the monomials composing the considered series. For two elements s 1 and s

2 in B[[γ, δ]], s 1 Rs 2 (i.e., s 1 and s 2 are equal in M ax in [[γ, δ]]
) is equivalent to the equality of their south-east cones. Direct consequences of previous geometrical interpretation are:

• simplification rules in M ax in [[γ, δ]] γ k ⊕ γ t = γ min(k,t) and δ k ⊕ δ t = δ max(k,t)
(4) • a simple formulation of the order relation for monomials γ n δ t γ n δ t ⇔ n n and t t

A simple interpretation of the variable γ and δ for daters is available:

• multiplying a series s by γ is equivalent to shifting the argument of the associated dater function by -1. • multiplying a series s by δ is equivalent to shifting the values of the associated dater function by 1 Example 5. Consider the series s = γδ 2 ⊕ γ 3 δ 3 ⊕ γ 4 δ 1 represented by dots in Fig. 1. The minimum representative of s in M ax in [[γ, δ]] is γδ 2 ⊕ γ 3 δ 3 . This result could be obtained using the simplification rules (4). Besides,

s = k 0 γ k δ -∞ ⊕ k=1,2 γ k δ 2 ⊕ k 3 γ k δ 3
Therefore, the dater d s associated with s is given by

d s (k) =    -∞ if k 0 2 if k = 1, 2 3 if k 3
Fig. 1. s and its south-east cone (hatched)

3.2 Linear state-space representation of TEG in M ax in [[γ, δ]]

From now on, we only consider TEG with at most one place from a transition to another transition. This assumption is not restrictive, as it is always possible to transform any TEG in an equivalent TEG with at most one place from a transition to another transition. The dynamics of a TEG may be captured by associating each transition with a series s ∈ M ax in [[γ, δ]], where d s (k) is defined as the time of firing k of the transition. Therefore, for TEG, γ is a shift operator in the event domain, where an event is interpreted as the firing of the transition, and δ is a shift operator in the time domain.

The transitions of a TEG are divided into three categories:

• state transitions (x 1 , ..., x n ): transitions with at least one input place and one output place. • input transitions (u 1 , ...u p ): transitions with at least one output place, but no input places. • output transitions (y 1 , ...y m ): transitions with at least one input place, but no output places.

Under the earliest functioning rule (i.e., state and output transitions fire as soon as they are enabled), with respect to a place with initially m tokens and holding time t, the influence of its upstream transition on its downstream transition is a positive shift in the time domain of t time units and a negative shift in the event domain of m events. The complete shift operator is coded by the monomial γ m δ t in M ax in [[γ, δ]]. Therefore, consider the place upstream from transition t i and downstream from transition t j , the influence of transition t j on transition t i is coded by the monomial f ij in M ax in [[γ, δ]] defined by f ij = γ mij δ τij where m ij is the initial number of tokens in the place and τ ij is the holding time of the place.

Consequently, a TEG admits a linear state-space representation in M ax in [[γ, δ]].

x = Ax ⊕ Bu ⊕ Rw y = Cx where x ∈ M ax in [[γ, δ]] n is the state, u ∈ M ax in [[γ, δ]] p the input, y ∈ M ax in [[γ, δ]] m the output and w ∈ M ax in [[γ, δ]
] n the additive perturbation of the state. The perturbation w models, for example, unexpected failure, delays or uncertain parameters such as task duration and matrix R how these perturbations affect the inner states, in the sequel R is assumed to be the Identity matrix.

A ∈ M ax in [[γ, δ]] n×n , B ∈ M ax in [[γ, δ]] n×p , C ∈ M ax in [[γ, δ]] m×n and R ∈ M ax in [[γ, δ]
] n×n are matrices with monomial entries describing the influence of transitions on each other.

According to Theorem 1, under the earliest functioning rule, the input-output (resp. perturbation-output) transfer function matrix H (resp. G) of the system is equal to

CA * B (resp. CA * ). y = CA * Bu ⊕ CA * q = Hu ⊕ Gq (5)
Therefore, the condition for holding the maximization (preserving the input-output and perturbation-output behaviors) is rephrased in terms of transfer function matrices. The conditions is now to preserve the input-output and perturbation-output transfer function matrices.

When an element s of M ax in [[γ, δ]] is used to code information concerning a transition of a TEG, then a monomial γ k δ t with k, t 0 may be interpreted as "at most k events occur strictly before time t" (i.e., d s (K) t). An element of s of M ax in [[γ, δ]], used to code a transfer relation between two transitions of a TEG (e.g., an entry of H), is causal (i.e., no anticipation in the time/event domain: all exponents are non-negative) and periodic (i.e., s = p⊕qr * with polynomials p, q and a monomial r = e). For a periodic series s with r = γ ν δ τ , its asymptotic slope σ(s) is defined as ν/τ . Fig. 2. TEG example Example 6. Let us consider a manufacturing system depicted by the TEG given in Fig. 2. The transition labeled u represents the inputs of raw material, it is transported during 2 time units to a machine with 2 treatment spots. Its input is labelled x 1 , the processing time is equal to 5, and the machine's output is labeled x 2 . The processed part is then transported out of the system during 3 time units, the transition y represents when the part is out of the production line. Before accepting a new raw part, the machine must be cleaned, this operation spends 7 time units. The model is then given by

A = ε γ 2 δ 7 δ 5 ε B = δ 2 ε C = ε δ 3
It must be noticed that this system can be realized in a straightforward way in (max, +) or (min, +) form:

x 1 (k) = max(2 + u(k), 7 + x 2 (k -2)) x 2 (k) =5 + x 1 (k) y(k) =3 + x 2 (k), x 1 (t) = min(u(t -2), 2 + x 2 (t -7)) x 2 (t) =x 1 (t -5) y(t) =x 2 (t -3),
where x i (k) represents the firing date of part k and x i (t) represents the number of firing occurred till the time t. These both systems are implicit equations. In order to obtain an explicit model we introduce an original procedure. We propose to split the system in the following way A = A r ⊕ A d ⊕ A g where

• if n ij = 0 and t ij = 0, (A r ) ij = A ij = γ nij δ tij else (A r ) ij = ε • if t ij = 0, (A g ) ij = A ij = γ nij else (A g ) ij = ε • if n ij = 0, (A d ) ij = A ij = δ tij else (A d ) ij = ε
In the present example:

A g = ε ε ε ε A d = ε ε δ 5 ε A r = ε γ 2 δ 7 ε ε
Knowing that x = Ax ⊕ Bu separating the matrix A we have:

x = (A d ⊕ A g ⊕ A r )x ⊕ Bu x = (A d ⊕ A g )x ⊕ A r x ⊕ Bu Using Theorem 1 x = (A d ⊕ A g ) * A r x ⊕ (A d ⊕ A g ) * Bu
Knowing that, we have A = (A d ⊕ A g ) * A r and B = (A d ⊕ A g ) * B and these matrices generate a model in the form

x = Ax ⊕ Bu y = Cx
Which can be realized in an explicit form whether in (max,+) or (min, +)

(A d ⊕ A g ) = ε ε δ 5 ε ⇒ (A d ⊕ A g ) * = e ε δ 5 e A = (A d ⊕ A g ) * A r = ε γ 2 δ 7 ε γ 2 δ 12
and for the input

B = (A d ⊕ A g ) * B = δ 2 δ 7 ,
this leads to the following explicit model

x 1 (k) = max(2 + u(k), 7 + x 2 (k -2)) x 2 (k) = max(7 + u(k), 12 + x 2 (k -2)) y(k) =3 + x 2 (k) x 1 (k) = min(u(t -2), 2 + x 2 (t -7)) x 2 (k) = min(u(t -7), 2 + x 2 (t -12)) y(k) =x 2 (t -3)
These matrices A and B can be programmed into a software without realization problems. The system can be solved by considering Theorem 1, (see Eq. ( 5)) and the transfer matrix y = Hu is given by H = C A * B = CA * B = δ 10 (γ 2 δ 12 ) * . This computation can be easily be performed by using the library MinMaxgd available as a C++ library or alternatively thanks to a web interface [START_REF] Cottenceau | Minmaxgd a library for computation in semiring of formal series[END_REF]; [START_REF] Cândido | Minmaxgdjs : A web toolbox to handle periodic series in minmax[END_REF]).

OBSERVER BASED CONTROLLER

This section presents how to implement an efficient control strategy for dynamical systems considered in the previous section. The control strategy proposed is depicted in Fig. 3. It is inspired from the observer based control for classical linear systems [START_REF] Hardouin | Control and state estimation for max-plus linear systems[END_REF][START_REF] Hardouin | Max-plus Linear Observer: Application to manufacturing Systems[END_REF].

The motivation to control the input of these systems is to decide when the operator should start to achieve an objective, e.g. when do you start the departure of a processing operator in order to achieve the customer demand.

Hence, the aim is to design a controller able to decide when the system should start to work in order to achieve a desired behavior. Classically, a popular production policy is to design a just-in-time policy, that is, to start as late as possible while ensuring the customer demand. It minimizes the internal stock while keeping the performance. The design goal is to get controllers M and P (see Fig. 3) which are matrices ensuring that the control input u = P (v ⊕ M x) be the greatest ( i.e. the one which delays as much as possible the input) in order to achieve a given objective, the reference input v. Signal x is either the real state of the system, (x = x if the state is measurable), or an estimation x observed thanks to an observer, inspired from the Luenberger observer [START_REF] Luenberger | An Introduction to observers[END_REF]). This estimator is given by

x = Ax ⊕ Bu ⊕ L(ŷ ⊕ y). ( 6
)
Where L is an observer matrix to be designed. It is fed by the measured output y of the system and ensures that the real system output be possible to compute the estimator x, especially that disturbance w feeding the system through matrix R. This observer based controller is then a feedback control strategy. The goal is to design P, M, L in order to achieve a desired behavior denoted G ref .

By solving Eq. 6 x is given by

x = Ax ⊕ Bu ⊕ L(Cx ⊕ C x) = (A ⊕ LC) * Bu ⊕ (A ⊕ LC) * LCA * Rw
by repeating in u, the control is:

u = P (v ⊕ M x) = P (M (A ⊕ LC) * BP ) * v ⊕ P M ((A ⊕ LC) * BP M ) * (A ⊕ LC) * LCA * Rw
The development are given in [START_REF] Hardouin | Max-plus Linear Observer: Application to manufacturing Systems[END_REF]) and lead to the optimal design

P opt = (CA * B)• \G ref (7) L opt = ((A * B)• /(CA * B)) ∧ ((A * R)• /(CA * R)) (8) M opt = P opt • \P opt • /(A * BP opt ) (9)

Example

In our sample system, showed in Fig. 2, since it does not have many inner states or inputs and outputs, these calculations can be done by hand and drawn. Using G ref = CA * B, we are able to calculate P opt = (CA * B)• \(CA * B) knowing that (CA * B) = δ 10 (γ 2 δ 12 ) * we obtain

P opt = δ 10 (γ 2 δ 12 ) * • \δ 10 (γ 2 δ 12 ) * = (γ 2 δ 12 ) * the equality holds since a * • \a * = a * (see Baccelli et al. ( 1992
)) which yields that (γ 2 δ 12 ) * • \(γ 2 δ 12 ) * = (γ 2 δ 12 ) * .
With the result of P opt we are able to calculate M opt

M opt = (γ 2 δ 12 ) * • \(γ 2 δ 12 ) * • / δ 2 δ 7 ⊗ (γ 2 δ 12 ) * M opt = δ -2 δ -7 ⊗ (γ 2 δ 12 ) *
Obviously it is not possible to implement a controller that has negative exponents (this controller would be noncausal), hence the solution is to pick only the causal projection. To do this, imagine a Cartesian plane where gamma is the x axis and delta the y axis. Now we put the points according to the desired series, for example

γ -4 δ -1 ⊕ γ -2 δ 2 ⊕ γ 2 δ 3 ⊕ γ 4 δ 4 .
Fig. 4. Gamma-delta plane representation of series γ -4 δ -1 ⊕γ -2 δ 2 ⊕γ 2 δ 3 ⊕γ 4 δ 4 and its causal projection

δ 2 ⊕ γ 2 δ 3 ⊕ γ 4 δ 4 .
It is clear that, using the south-east cone presented in section 2 all the points in the drawing must be inside the north-east quadrant for the series to be realizable. This way, the causal projection is the biggest area possible, in a way that all its corners are inside this quadrant, and is contained in the original area. For our example, it would be δ 2 ⊕ γ 2 δ 3 ⊕ γ 4 δ 4 , as showed in Fig. 4 Using the same reasoning in M opt we get the following causal series:

P r + (M opt ) = γ 2 δ 10 γ 2 δ 5 ⊗ (γ 2 δ 12 ) *
Finally the observer needs to be calculated

L 1 = δ 2 δ 7 ⊗ (γ 2 δ 12 ) * • /δ 10 ⊗ (γ 2 δ 12 ) * = δ -8 δ -3 ⊗ (γ 2 δ 12 ) * L 2 = e γ 2 δ 7 δ 5 e ⊗(γ 2 δ 12 ) * • / δ 8 δ 3 ⊗(γ 2 δ 12 ) * = γ 2 δ 4 δ -3 ⊗(γ 2 δ 12 ) * Then L = L 1 ∧ L 2 : L = γ 2 δ 4 δ -3 ⊗ (γ 2 δ 12 ) * ⇒ P r + (L) = γ 2 δ 4 γ 2 δ 9 ⊗ (γ 2 δ 12 ) * .
Since all the matrices entries are causal, all these controllers are realizable, and thus they can be implemented.

To convert the gamma-delta elements into min-plus or max-plus equations, we propose to split the series s = p ⊕ qr * in the following way : first we define ζ k = qr * thus s = p ⊕ ζ k . By using Theorem 1, we get the relation

ζ k = rζ k ⊕ q.
In our example x = Ax ⊕ Bu ⊕ Ly, for example sake, let us focusing on the first state we obtain, x1 = γ 2 δ 7 x 2 ⊕ δ 2 u ⊕ γ 2 δ 4 ⊗ (γ 2 δ 12 ) * y. Hence in this case p = γ 2 δ 7 x 2 ⊕ δ 2 u, q = γ 2 δ 4 y and r = γ 2 δ 12 . By using the formula we just found, x1 = γ 2 δ 7 x 2 ⊕ δ 2 u ⊕ ζ 1 and

ζ 1 = γ 2 δ 12 ζ 1 ⊕ γ 2 δ 4 y.
Thanks to this expression it is straightforward to put these equations in min-plus or maxplus implicit equations. For simplicity, min-plus will be used because updating the system periodically, in our case each 500ms, is easier than when events occur. The min-plus

equation for x1 is x1 = min(2 + x2 (t -7), u(t -2), ζ 1 (t)) knowing that ζ 1 (t) = min(2 + ζ 1 (t -12), 2 + y(t -4)).
These equations can be easily programmed into a software, since they only require data storage to use the time delays properly.

REAL SYSTEM APPLICATION

In this section, we are interested in applying this method in the real system, depicted in Fig. 5, which is located in Polytech Angers, France.

The system has 2 separated sections, a faster loop and a slower loop as shown in Fig. 5. The slower loop has 6 buttons that do not let the pallets pass while the faster loop has only 4. All the buttons have sensors just before them, as it can be seen in Fig. 6. The pallet's size is such that if they are waiting the button, the sensor will stay active. Each section (between two buttons) has a defined maximum number of pallets. The travel times were measured 10 times, and the time used is the average of them. Each button (labelled B1 to B10), the travel time and the maximum number of pallets are represented in Fig. 7. It is important to notice that there are 3 pallets waiting for B1, 2 waiting for B5 and 1 waiting for B6, as initial conditions of the system. Fig. 7. Buttons, travel times (in sec) and pallet limit

The system is programmed to activate the buttons when there is one pallet waiting for it (the sensor is active), there is at least one space left for the path and at least one control token available. Especially for B3 and B10, there is a forced synchronization, meaning that B3 and B10 will always activate at the same time, requiring 2 control tokens, one free space between B10 and B5 and Now the first step is to put this system in a Petri net model. We will analyze the section between B4 and B1 and then apply the same reasoning to the other sections.

Since each button changes the state of the system, we have the first two inner states, each one with one associated input because of the control tokens as shown in Fig. 8. Then we will need one place for each empty slot in the section, and since 2 pallets cannot be in the same place, these places must have only one token at a time. Besides the tokens, the timings for each place has to be 2 seconds, since the button will only activate 2 seconds after its previous activation, but the sum of all the places between the buttons has to be the total travelling time, using these conditions we have the second step of the model in Fig. 9. It is worth mentioning that the initial pallets are the tokens into (P1), ( P2) and (P3) so the token in (P1) can activate B1 instantaneously.

To summarize the model so far, the token in P1 will activate the B1 transition, meaning that the button will go down and the pallet will enter the path between B1 and B2. After that, the two other pallets will begin moving and wait for the 2 seconds to pass, what explains the 2 seconds travelling time in the model. The last condition that is missing is the total number of pallets for each section.

Since for the section B4 to B1 only 3 pallets are allowed, and all the three are already in it, so we have the final model represented in Fig. 10.

If we replicate this reasoning to all sections combined and knowing that before every button we have a sensor indicating the pallet is there, this means the sensors are the outputs of the system, we get the final model represented in Fig. 11.

We can achieve the system's matrices as proposed in section 3.1, where A 0 will be split into A g , A d and A r to make the model implementable using the same procedure discussed in section 3. The next step is to use these 3 matrices to calculate P opt , L opt and M opt thanks to Eqs. (7-9). All the calculus needed to achieve the final controller were obtained by using the software MinmaxGD [START_REF] Cottenceau | Minmaxgd a library for computation in semiring of formal series[END_REF]. After executing, this software gave us the matrices just like the ones presented in section 4.1.

Transforming those matrices in code is extremely laborious and can lead to unexpected errors. Knowing that, a software has been developed to make the calculus and write a C code directly (see [START_REF] Freitas | State estimation of an automated system in max-plus algebra by using a linear observer[END_REF]), using the method shown at the end of section 4.1. The generated code was then executed into a PC that supervises the real system. The PLCs ensure that the system will comport as specified and provide data for the supervisory system. The code works in a way that it reads the outputs provided by the PLCs, updates the input control tokens consequently updating the input and thus updating the amount of control tokens available for the buttons, then updates the observer for the user to know if the system is comporting adequately every 500 ms. 6. RESULTS

In this section we are going to discuss the effects of having implemented the control into the real system. Since the only place that can occur pallet accumulation is in B4 and B10 and the loop containing B4 is faster, it is the only place that can have accumulation. Our chosen G ref was the system's transfer function, that means, we chose to maintain the original output, but to delay as much as possible the buttons activations in order to reduce stocking. Knowing that, the objective of the controller is to reduce the stocking in B4 without changing the original output. To illustrate the performance of the controller we focus on measurements of x 7 state.

It must be mentioned that the travelling times are not deterministic, hence we focus on the difference between x 7 state with and without control, sometimes this difference is positive and sometimes negative. In Fig. 12 it is possible to see the time evolution of this difference. As expected, the difference is insignificant principally because x 7 is the number of times the state has been activated, which is a number that increases over time. The mean value of the difference is -0.1342, very close to 0, showing that all outputs would have been maintained if the system was perfectly deterministic.

The next topic will be the stocking difference. For this matter the PLCs is programmed to measure x 4 's information, this way the stock is calculated using stock(t) = x 4 (t)x 7 (t). The result can be seen in Fig. 13.

This pattern keeps repeating for the whole simulation. As it is possible to see, the controlled system has at maximum 2 pallets in stock, while the not controlled (natural) system has 3. This result makes sense, since our aim is to maintain the outputs, in this case, for output y 2 to be maintained, the system has to activate B2 at least twice, resulting in a minimum stock of 2. Calculating the mean stock, we get 1.6550 pallets for the natural system, and 1.3175 for the controlled one, meaning a 20.4% stock reduction.

CONCLUSION

In this paper we propose and present a way to implement in an efficient way controllers for max-plus linear systems.

Freely available software tools yield automatically the code for implementation in a SCADA system. In the real example presented, the just-in-time control implemented has reduced the useless internal stocks by 20%, which illustrates the interest of this kind of production policy. The engineers can use this methodology to adapt the control to their specific systems in an easy way. A user-friendly interface will be develop soon in order to help them to select the control strategy they want to implement among the different available possibilities (Observer-based control [START_REF] Hardouin | Control and state estimation for max-plus linear systems[END_REF], Disturbance Decoupling Problem [START_REF] Shang | An integrated control strategy to solve the disturbance decoupling problem for max-plus linear systems with applications to a high throughput screening system[END_REF], Robust control [START_REF] Lhommeau | Interval analysis in dioid : Application to robust open loop control for timed event graphs[END_REF], Stabilization thanks to Feedback Control [START_REF] Cottenceau | On timed event graphs stabilization by output feedback in dioid[END_REF], Model Predictive Control [START_REF] Schutter | Model predictive control for max-plus-linear discrete event systems[END_REF]).

Example 3 .

 3 (min,plus) algebra: Z min = (Z ∪ {-∞, +∞}, min, +) is a complete idempotent semiring such that a ⊕ b = min(a, b), a ⊗ b = a + b, a ∧ b = max(a, b) with ε = +∞, e = 0, and = -∞. The order is total and corresponds to the inverse of the natural order ( i.e. , 2 1. Semiring of matrices Z n×n min is a semiring of matrices with entries in Z min . Example 4. (Matrix operations in Z max ). Given three matrices with entries in Z max ,

  Fig. 3. Observer-based Controller architecture

Fig

  Fig. 5. Slower Loop

FigFig. 10 .

 10 Fig. 8. Modeling 1 Fig. 9. Modeling 2

Fig. 11 .

 11 Fig. 11. Complete Model[START_REF] Gonçalves | Tropical algorithms for linear algebra and linear event-invariant dynamical systems[END_REF]

Fig. 12 .

 12 Fig. 12. Difference between x 7 with (Natural) and without controller (Controlled)


EDF), H. Joumni (CSTB S Nösperger

H Joumni

Mohamed Elmtiri

Olivier Bonin

VEOLIA Oliver Keserue

SETEC Eric Blanc

Armelle Langlois

Designing

VEOLIA), O.Bonin (IFSTTAR -LVMT M El-Mtiri

Designing economic models for the energy transition at the territorial scale

Keywords: economic model, business model, energy transition, cities, externalities

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

BACKGROUND

Energy efficiency investments have by and large been limited to energyintensive industries 'contribution to climate change mitigation [START_REF] Grubb | Planetary Economics, Energy, climate change and the three domains of the sustainable development[END_REF]. In urban projects energy efficiency is les concerted and explicit.

With new large urban projects encouraged to integrate climate change mitigation through high energy performance standards, the impact on project costs is considerable. From the additional conception costs to all the other additions in terms of equipment, facilities and infrastructure, energy efficiency implies a considerable increase in costs when constructing a new housing, office buildings or transport infrastructures. This is even more significant in the case of urban renovation (Guennec et al., 2009, Nösperger et al., 2015). However the poor economic outlook is making it hard to justify this additional investment in the absence of a market for these ambitious projects. Yet, social and environmental co-benefits related to energy efficiency are now apparent (Tirado Herrero et al., 2011;Ürge Vorsatz IEA 2014). How can these contradictions between benefits and costs be solved in order to make projects viable? We shall attempt to solve this paradox by addressing a number of questions:

 When assessing the economic opportunity of an urban project, how can it be ensured that the socio-environmental benefits related to improved energy efficiency are taken into account? How are these benefits qualified from the perspective of different stakeholders and monetised, when relevant?  What types of partnerships, agreements or contracts exist which ensure an equitable distribution or sharing of the economic value stemming from the identified co-benefits? Is this value sufficient to make the energy efficiency investment viable?

Today's energy production facilities are either highly centralised (power plants) or highly decentralised (e.g. photovoltaic panels installed on roofs). What can learn from business models developed for transport (Uber) and internet (Google) that can integrate the multi-scale territorial perspective that characterises the diversity of actors and object in the area of energy efficiency, for example through smart grids. EFFICACITY, the institute for energy transition in the urban environment, has been addressing these questions in a thematic program with the aim of delivering operational tools for use by stakeholders involved in urban projects (local authorities, energy and water suppliers, developers, etc.).

Following a literature review (part 2 of this article), a methodology for a broadened economic assessment was developed by EFFICACITY (presented in part 3) .This methodology was applied in two urban projects to test its validity, which is the subject of sections 4 and 5 of this article. We then end with a concluding section.

LITERATURE REVIEW

Environmental externalities have been studied in numerous projects and research articles and provide a solid basis to help identify broad families of externalities within the framework of an energy efficiency project (Amann, 2006), such as health impacts [START_REF] Loftness | Linking energy to health and productivity in the built environment. Evaluating the Cost-Benefits of High Performance Building and Community Design for Sustainability, Health and Productivity[END_REF], biodiversity protection [START_REF] Herrero | Co-benefits quantified: employment, energy, security and fuel poverty implications of the large-scale, deep retrofitting of the Hungarian building stock[END_REF], and greenhouse gas reduction (IEA, 2014). In the case of environmental externalities, it should be noted that some States set "shadow values" that is to say, the value of the external costs when evaluating public infrastructure projects (Quinet 2008[START_REF] Quinet | Evaluation socio-économique des investissements publics[END_REF], CESD "socio-economic assessment" ...). The positive impact on local economic development (2014Guennec et al IEA. 2009) or productivity [START_REF] Loftness | Linking energy to health and productivity in the built environment. Evaluating the Cost-Benefits of High Performance Building and Community Design for Sustainability, Health and Productivity[END_REF] are also revealed. This literature review has been synthesized in a reference book (Effiacity 2014) and an accompanying database (under development).

The following table provides some examples of the co-benefits of energy efficiency in the case of buildings. In fact, the main difficulty, when evaluating a well-defined energy efficiency project, lies in adapting locally these general externalities. While the enlarged economic evaluation of such urban energy efficiency projects can rely on proven methods (ISO 15686-5, Joumni 2009), or contrary to conventional evaluation methods such as Cost Benefit Analysis, where the State arbitrates between winners and losers, local energy efficiency projects require defining externality values which find endorsement by stakeholders for them to be then be inclined to contribute to a project on a voluntary basis.

The valuation of externalities is thus linked to the creation of economic models able to integrate those actors, traditionally absent or underrepresented, in current economic transactions around building projects.

In fact, beyond the business model, which is a monetary translation, the business model implies (du [START_REF] Tertre | Modèles économiques d'entreprise, dynamique macroéconomique et développement durable[END_REF]Nösperger et al, 2015.) gathering the following prerequisites:



A basis for value creation and distribution (capturing) between the actors; 

The basis and formal structure of relationships between actors (contracts, formal or informal partnerships, "cooperative" society ...);



Sources of mutual investment or financing operations;



The nature of work performed by these players: skills, qualifications, activities produced independently or co-constructed (co-production).

The methodology presented in the next section proposes to address simultaneously the issue of externalities / co-benefits to consider and business models suited to their actual transformation into a monetary or non-monetary value.

METHODOLOGY

EFFICACITY developed an 8-stage methodology intended to assist in identifying relevant externalities of an Energy Efficiency project, determine relevant monetary values, and design suitable partnerships likely to convert them into economic flows (financial or non-financial): Phase 1: Local context identification and definition of alternative solutions (stages 1-3):

-Step 1: Background and overview of the initial situation (nature of the project, scope) and identification of the sets of actors involved in the project.

-

Step 1: context

The planned urban block being subject to evaluation is located in a neighbourhood in the Paris basin shared with new buildings of a university campus, an architectural school, a Regional Express Railway station (RER), college halls of residence and a university canteen. This mixed urban block consists of 4000m² of offices and 13,500 m² student residences. The local authority is also considering building a large nautical centre (2000 m² Basin). The project is in an area with special status, a "new town", managed by a public development agency named EPA (Etablissement Public d'Amenagement).

EPA owns the grounds of the planned urban block and intends to sell it to the property developer who will commission and manage the project and supervise other contracted project managers such as for the halls of residence and the office buildings. This project is the subject of a single contract of Design-Build-Operate & Maintenance of which the incumbent manages the various components having already identified a number of subcontracting partners. EPA intends the cluster located in its district to be competitive internationally against other campuses.

Local energy resources include: an aquifer (120 m deep, 14° C), solar radiation, biomass (subject to demand by the local authority), and fatal heat from the nautical centre. A system composed of 5 wells was installed for the university building, with a discharge well back to the aquifer.

The energy balance of the territory points to:  Net demand for heat from the student canteen in the campus, the building block and nautical centre. The net heat requirement of the building block is due to the predominance of residential housing (and high domestic hot water needs)  Net cold requirements from the university building. This building will therefore significantly use a cooling tower before releasing water above or below the temperature of aquifer. However, the low flow rate of the aquifer will heat up the aquifer in the absence of additional uses for the waste heat.

Steps 2 & 3: Alternative energy solutions

The reference project, for calculating the additional cost of the energy efficiency solution, is a building block powered by a gas boiler and airconditioning from adiabatic air coolers. The reference block is 30% more efficient than new period buildings.

The comparative analysis that follows, with respect to this reference solution, takes an energy efficient starting point (baseline).

In the context of this article, we will develop only the most promising alternatives, especially from the perspective of their intensity of use in a context of basically very efficient buildings. It is noteworthy that all of these alternatives include installing photovoltaic sensors whose production will be sold at a fixed feed-in tariff for 20 years.

The energy concept 1 of the building block is based on combining multiple sources of energy: gas, electricity, solar photovoltaic, and wastewater. Energy recovery from the grey waters is more robust (reliable and cost-effective) than solar thermal sources because it is simpler and requires no maintenance. Energy systems are centred on a thermorefrigeratingpump (TFP) that is not able to supply in a balanced way the hot and cold needs of the block. As a result there is a systematic loss of refrigo-calories (cold) having to be released back to the aquifer. This solution could be quite effective from a functional and energy standpoint in theory; unfortunately the temperature rise of the aquifer requires a search for alternatives. Indeed, while the energy balance of Bienvenüe university building wells, on the one hand, and the building block V1 on the other (which are located at a distance) could compensate "partially", the very low flow rate of the aquifer does not allow this "compensation" to be effective.

The Energy Concept 2 relies on the Mutualisation of the abstraction from the aquifer. The thermorefrigeratingpump on this aquifer uses the same wells used for the university building. This solution has the advantage of partially addressing the problem of excess heat emitted annually into the aquifer by the Bienvenue Building (due to an overall cooling need for the year) that leads to its gradual warming (because of the low rate of flow of the aquifer preventing heat to disperse). The decision to adopt this solution is based on the assumption that coupling the university building and the building block is sufficient to balance the needs and reduce very significantly the heat discharge to the aquifer. The innovative character of these buildings with high performance outer casings accentuates the risk of a mismatch between the simulated and observed behaviour. The coupled system "Building block -University building" may still be out of balance with an annual heat surplus, which would still contribute to the gradual warming of the aquifer.

The Energy Concept 3 relies on the Mutualisation of the abstraction from the aquifer and on the addition of a nautical leisure park. The policy targets set by the local authority requires it to meet 50% of the heat demand of the nautical centre with renewable energy, in this case with solar thermal panels. The surface of these panels is downsized to avoid overheating. Assuming that the "super-block" made of the projected building block and the university building to be in net heat excess in the summer (releasing water at 35° C), this balance could be used by the nautical centre to further reduce its summer minimal heat demand and consequently the surface needed for the solar panels. Before such considerations the surface area of the solar panel stood at 2900 m². Investment savings are therefore foreseeable.

Identifying externalities and co-benefits

We have to reason differentially (instead of reasoning on the comparison of absolute life cycle costs of each alternative) with respect to a high performance reference solution. The alternatives will not be compared with the baseline reference with respect to comfort, health, and mobility, which are independent of benefit from the mutualisation of the aquifer energy source. The real-estate value should also not be impacted. Most of the cobenefits discussed here therefore are focused on the impacts of the mutualised aquifer energy/cooling source. These are treated in Section 4.4.

Estimating the market and non-market values

Cost elements (tradeable goods and services)

The following table shows the cost elements related to the tradeable goods and services as compared with the reference solution.

Monetising non-market benefits

The previous financial balance must be completed with the assessment of monetized non-market benefits.

Fighting climate change

The following table shows the cost elements related to market goods and services, always by difference with the reference solution. The three concepts outperform the reference in terms of GHG balance. This difference is greater for the concepts 2 and 3, with the reliance on the thermorefrigerationpump with a coefficient of performance (COP) greater than 3. We used the shadow price of carbon used in public investment evaluations [START_REF] Quinet | Evaluation socio-économique des investissements publics[END_REF]. This value is 52 € / t in 2015 with an escalation rate of 4.5%, equivalent to the discount rate used in public investments under the law of Hotelling (Quinet, op.cit.).

Environment, biodiversity and semiotics

While the project global warming attenuation capacity is indisputable, there may be an impact on biodiversity and the natural environment from the warming of the aquifer. Without corrective action, because of Darcy's velocity (flow rate) of the aquifer and hardness of the water, we could expect this warming over 15 years to result in two consequences:

 The development of algae and amoebae,  Blockage and instability of the release well from the thermorefrigeratingpump at the university building (already installed). The use of cooling towers in the university building slows the phenomenon but does not suppress it.

Preventing the warming up of the aquifer improves the overall robustness of the system (at a 15 year horizon). The alternatives (concepts 2 and 3) for mutualising and compensating the heating and cooling needs reduce the warming of the aquifer. Furthermore, each alternative would give the territory a "showcase" of an innovative university campus, and fulfil the ambition of international renown mentioned above given the university's focus on energy transition and sustainable development. However, this effect is much more significant in the case which includes the water leisure park (concept 3).

The monetization of impacts related to the aquifer was carried out by the method of costs of damage [START_REF] Pearce | Cost Benefit Analysis and the Environment: Recent Developments[END_REF]. The main physical consequence is the impossibility of using a geothermal heat pump for the university building at a 15-year horizon. The assessment of that damage may be based on the discounted residual value of the equipment, valued at € 24 / m² (relative to the surface area of the planned building block) prior to discounting. The discounted value will be chosen to monetize the efforts for preventing the warming of the aquifer in the case of concepts 2 and 3.

In addition, EPA announced it would cover 50% of the investment cost associated with these three concepts in recognition of its innovative nature and international renown.

Operating savings of the university building

The expected synergies between the building block and the university building should reduce the release of cold effluent into the aquifer, enabling both improved energy efficiency of the cooling towers of the Bienvenüe building and decrease the use intensity. Thus, it is predicted that the overall operation of these towers is expected to fall between 20-60%! This is accompanied by reduced maintenance and operating costs and longer service life of the equipment, a notable reduction of the operating costs. This represents an externality to the extent that it constitutes a cost or a benefit to either party, without it being valued economically. This externality is positive as the action of agent A (the use of the aquifer by the building block) results in improved well-being of another, agent B (a decrease in operating costs of the university building) without the latter having had to pay the price (Efficacity 2015).

In the absence of a partnership linking the two parties, the university building is outside the economic reasoning of the actors of the building block.

Based on the information provided, the operational savings from the cooling towers permitted by mutualisation between the university building and the building block (concepts 2 and 3) were estimated at € 0.1 / m² (based on the surface area of the planned building block).

Solar panel investment savings for the nautical centre

By broadening the mutualisation to include the nautical centre the synergies allow a reduction in the surface area of the solar thermal panels (see above) of approximately 10%.

As above, this benefit may be associated with an externality. Based on the information provided, 10% investment savings in solar panels represents € 180,000, or 10.3 € / m² (based on the surface of the building block).

Step 6: Arbitration of the strategies for optimising the allocation of resources

The table below shows the results of the evaluation in overall cost and extensive global cost (ISO 15686-5) of the different variants relative to the reference solution, performed over a period of 20 years with a discount rate of 3% (as per the EN 16627 standard). A negative value indicates a gain relative to the reference solution. Extended evaluation 8.9 2.9 -4.5

It appears that the energy savings will not even cover the additional maintenance costs and the competitiveness of all three concepts depends on the value of the external effects. However, at a 20-year horizon, only e Concept 3 is economically competitive with respect to the reference solution. Choosing this path requires to engage in negotiations involving numerous partnerships.

Step 7 and 8: Designing of an appropriate economic model

and Evaluation of the deployment of the selected strategies

Up to

Step 6, we were in a cost/benefit logic proceeding by an itemised accounting (costs / market and non-market benefits). It must be formalised by the development of relevant partnerships following a logic of contributing actors and economic flows from one actor to another.

Partnerships can be bi-lateral or multilateral (for instance in the case of a "mutual fund" drawing contributions for different benefits or to supervise the bilateral exchanges). There's an apparent need for "systems of rules that bind the actors together, and are developed and adopted during the project" (Pasquelin, 2015). Beyond this governance, the success of such an economic model hinges on genuine cooperation between actors, who share an interest in each other's activities and how they interaction towards achieving a common goal (du [START_REF] Tertre | Modèles économiques d'entreprise, dynamique macroéconomique et développement durable[END_REF]. This instance illustrates, within acceptable technical and economic conditions, the possibility for an assembly comprising the three parties (the university building, the building block and the nautical centre) in the context of a shared willingness to showcase the exemplarity of the project.

THE METHOD IN THE CASE OF A DISTRICT HEATING NETWORK

The second application of this methodology was carried out on a district heating project (steam provided by an incineration plant). The sensitivity analysis of the sector-specific economic models (taking into account identified externalities) has been carried out for each actor's decisions. The two scenarios under study are the retrofitting of building to make them energy efficient, and the application of household waste minimisation policies. The externalities associated with these policies leads us to imagine new configurations of player games that redistribute earnings (profits) and losses for the most beneficial outcome. This constitutes a new way to build an all-encompassing cost approach by establishing new contract arrangement.

Scenarios and simulations

We are interested in all the actors' up-and down-stream of the two scenarios:

 Users of the heating network (households) recipients of the final energy bill and its share of domestic income;  The incineration plant as heat supplier, and the contribution of heat sales to its operating margin (net of investments);  The operator of the district heating system, who adapts the billing of its domestic customers according to the energy source and thus maintains its operating margin. Three parameters are modified independently and simultaneously to develop various scenarios and assess the economic effects associated with the choice of its actors. These three parameters are:

 The percentage of waste incineration: today several local policies target waste minimisation preferring to maximize recycling to the detriment of incineration.

 The percentage of the renovated buildings on the district heating network: building standards improve efficiency and reduce energy consumption. The renovation policy targets a blanket application of new building standards.

 The percentage of residential housing connected to the district heating network. Recent policy orientations, national and local, have made heating network a key policy tool for attaining energy transition objectives.

Results

The fall in volumes of waste incinerated expected from waste minimisation policies may impact significantly the earning of partner stakeholders in a district heating network. The first concerned is the operator of the incineration plant, who sells less heat and therefore reduces its earnings. The plant operator must compensate by increasing prices to be able to purchase an alternative source of energy (gas) and fulfil its contractual obligations with downstream users of the district heating network. Beyond a certain threshold of waste reduction the operator of the heating network must install a booster station, increasing its operating costs. Its margin is decreased but is not erased, since it passes most of the additional cost to its users.

Arrangement of players games

The resilience of the heating network needs rethinking, without compromising either the environmental or economic objectives pursued in the two policies.

In the case of a fall in the volume of waste incinerated, it seems more difficult to remove the losses for all players, but still we can try to maintain the profitability of the system. In such a case, and to maintain a satisfactory level of renewable energy in the energy mix of the network, the thermal renovation of buildings should be encouraged (assuming it will reduce use and cost of fossil fuels to users). However, such action does not compensate the losses for the heating network and incineration plants. We can, in this case imagine a mutually beneficial compromise to encourage thermal renovation for users and to compensate the district heating operator and incinerator for the loss of business (including downsizing the furnaces to reduce operating costs).

The reduction in waste is not necessarily negative for the system; it is possible to "downsize", through public financing, all elements of the chain. The decrease in the production of household waste results in lower operating costs for the community. These savings could be redirected to finance this upgrade, and then proceed with a gradual decrease in user-fees for managing household waste. The scenarios with variations of multiple parameters showed that:  It is possible to compensate for the effects of one of the parameters by those of the other.

 We need to rethink about how to define the energy strategy of the urban area in an integrated way: you have to consider actions on consumption, production and distribution simultaneously. A player should have the supervisory role in planning and upholding the general interest.

CONCLUSION

The application of this methodology turns out to be effective for identifying and valuing externalities. A key point is the building of new economic models able to convert the identified economic value into financial flows (or equivalent) between stakeholders. These models will rely on a formal basis (innovative partnership, contracting) but especially on informal ties (trust between actors, organizational relevance…). The development of relevant formal and informal resources can be delegated to a trusted third-party. To help the latter achieve this complex work, operational tools still need to be developed. These may include:

• a comprehensive reference document covering externalities (literaturebased review) to assist in the complex dialogue between stakeholders on the value of externalities;

• a global life-cycle costing calculation tool including externalities valuation;

• interviews, support guides for the negotiation between stakeholders (benefiting from positive externalities) and who have the willingness to contribute to the project in order to ensure that it become viable;

• specific contracting and partnership supporting tools in the areas of building renovation and district heating.

The development of such tools is the main objective being pursued by Efficacity in the framework of this project.

  Step 2: Identification of technological and organizational solutions in energy efficiency adapted to the situation (and envisaged in the EFFICACITY relevant programs). Who do they concern and to what extent? -Step 3: Classification and selection of a range of solutions Phase 2: Externalities and benefits identification, selection and monetization (stages 4-6) -Step 4: Identification and selection of externalities related to selected solutions. -Step 5: Estimated market and non-market economic values -Step 6: Arbitration of possible strategies to optimize resource allocation. Phase 3: Identification of relevant partnerships/ contractual relationships and business model design (stages 7-8). -Step 7: Design the business model adapted and sensitivity tests -Step 8: Evaluation of the implementation conditions of the selected strategies: economic and contractual conclusions

Table 1 .

 1 Examples of an integrated analysis of building functions and impacts (Source: Efficacity Externality reference tool andNösperger & al., 2015) DDSS 2016, 13th International Conference on Design and decision support systems in architecture and urban planning, June 27-28, 2016, Eindhoven

	Impact or	Impact Description		Stakeholder
	externality
	category
	Health	Occupant thanks to comfort improvement health improvement	Occupier Occupying company
		Improved health thanks to facility	Occupant
		(cycle boxes, shower facility, pedestrian
		paths…) allowing environment-friendly
		mobility (bicycling, walking)
	Productivity	Occupant	wellness	related	Occupying company
		productivity improvement thanks to
		comfort improvement (temperature,
		natural lighting)
		productivity improvement thanks to	Occupying company
		building-related	positive	corporate
		communication

Table 3 :

 3 Financial results of the three concepts in relation to the reference (costs and market

	benefits
	Extra	Impact	on	Impact	on	Electricity
	investment	the yearly energy	annual		sales from solar

Table 3 :

 3 Financial results of three concepts in relation to the reference (costs and market

	benefits)
		Gas	Electricity	Primary	CO2
		Consumption	consumption.	energy	Emissions
		(final energy)	(final energy)	consumption
	Concept	-18.8	+4.6 kWh/m²	-6.9	-3.7 kg/m²
	1	kWh/m²		kWh/m²
	Concept	-22.7	+6.7 kWh/m²	-5.4	-4.3 kg/m²
	2	kWh/m²		kWh/m²
	Concept	-22.7	+6.7 kWh/m²	-5.4	-4.3 kg/m²
	3	kWh/m²		kWh/m²

Table 4 :

 4 Calculation of total cost (simple and extended) of the different concepts (€ / m²)

		Additional	Impact	on	Impact	on
		investment costs	annual	energy	annual maintenance
			costs		costs
	Initial investment	58.0	69.0		72.0
	Operating costs for	-10.9	-12.4		-12.4
	energy and water
	Maintenance costs	11.5	13.1		13.5
	Residual value for	-1.7	-2.0		-2.1
	subtraction
	Social costs CO2/	-4.2	-4.9		-4.9
	period

DDSS 2016, 13th International Conference on Design and decision support systems in architecture and urban planning, June 27-28, 2016, Eindhoven


Laurence Pougnet

Claire V Hoffmann

Richard Pougnet

Fabien Le Ny

Léa Gaitan

Pierrick Cros

Mathilde Artus

Solène Le Gal

Gilles Nevez
email: gilles.nevez@univ-brest.fr

Pneumocystis exhalation by infants developing Pneumocystis primary infection: putative infectious sources in hospitals and the community

Keywords: Pneumocystis jirovecii, primo-infection, airborne transmission, exhalation, Pediatrics

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

Introduction

Pneumocystis jirovecii is a transmissible and uncultivable fungus, specific for humans, that causes severe pneumonia [i.e. Pneumocystis pneumonia (PCP)], in immunosuppressed adults and children [START_REF] Hughes | Pneumocystis Pneumonia[END_REF]. Pneumocystis primary infection linked to first contact with the fungus occurs early in life in both immunosuppressed or immunocompetent infants, mostly those younger than 12 months [START_REF] Vargas | Search for Primary Infection by Pneumocystis carinii in a Cohort of Normal, Healthy Infants[END_REF][START_REF] Larsen | Primary Pneumocystis infection in infants hospitalized with acute respiratory tract infection[END_REF][START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF]. Pneumocystis primary infection in immunocompetent infants appears to be a worldwide phenomenon [START_REF] Vargas | Search for Primary Infection by Pneumocystis carinii in a Cohort of Normal, Healthy Infants[END_REF][START_REF] Larsen | Primary Pneumocystis infection in infants hospitalized with acute respiratory tract infection[END_REF][START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF], the infection being asymptomatic or symptomatic. In this latter case, it occurs contemporary to a self-limiting respiratory tract infection [START_REF] Larsen | Primary Pneumocystis infection in infants hospitalized with acute respiratory tract infection[END_REF][START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF] in the course of which the fungus can be detected using PCR assays in nasopharyngeal aspirates (NPAs) [START_REF] Vargas | Search for Primary Infection by Pneumocystis carinii in a Cohort of Normal, Healthy Infants[END_REF][START_REF] Larsen | Primary Pneumocystis infection in infants hospitalized with acute respiratory tract infection[END_REF][START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF].

Pneumocystis species are host specific and no environmental ecological niche for any of them has been identified so far. For these reasons, it is assumed that patients infected by P. jirovecii are the main part of the fungus reservoir. Indeed, the aerial interindividual transmission of Pneumocystis murina has been demonstrated in mouse models [START_REF] Dumoulin | Transmission of Pneumocystis carinii disease from immunocompetent contacts of infected hosts to susceptible hosts[END_REF] and was strongly suggested for P. jirovecii in humans considering the occurrence of PCP case clusters in hospitals (see review in [START_REF] Pougnet | Airborne Interindividual Transmission of Pneumocystis jirovecii[END_REF]). Moreover, P. jirovecii DNA has been detected in the air surrounding PCP patients or patients colonized by the fungus, consistent with P. jirovecii exhalation and spread from patients in their environment [START_REF] Choukri | Quantification and Spread of Pneumocystis jirovecii in the Surrounding Air of Patients with Pneumocystis Pneumonia[END_REF][START_REF] Gal | Pneumocystis jirovecii in the air surrounding patients with Pneumocystis pulmonary colonization[END_REF][START_REF] Pougnet | Pneumocystis jirovecii Exhalation in the Course of Pneumocystis Pneumonia Treatment[END_REF]. These patients represent potential infectious sources from which, susceptible other individuals may acquire P. jirovecii.

Considering the very high prevalence of Pneumocystis primary infection in infancy [START_REF] Vargas | Search for Primary Infection by Pneumocystis carinii in a Cohort of Normal, Healthy Infants[END_REF][START_REF] Larsen | Primary Pneumocystis infection in infants hospitalized with acute respiratory tract infection[END_REF][START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF], infants in these circumstances may represent the main part of the human reservoir of the fungus. P. jirovecii detection in the surrounding environment of infants would strengthen their potential role as infectious sources. In this context, we performed the detection and characterization of P. jirovecii DNA in air samples from the surrounding environment of

Patients and Methods

Pneumocystis primary infection was defined as detection of P. jirovecii in NPA samples from infant younger than 12 months, an infant who was a priori immune-naïve to P. jirovecii [START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF].

The detection of P. jirovecii was routinely performed in our diagnostic laboratory in NPAs contemporary to that of viruses within the framework of investigation of pulmonary symptoms or fever. This detection was done using a real-time polymerase chain reaction (qPCR) assay targeting the gene of mitochondrial large subunit ribosomal RNA (mtLSUrRNA) on 7500 Real Time PCR System (Applied Biosystems, Foster City, CA, USA) as described elsewhere [START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF].

Fourteen out of 87 infants younger than 12 months (sex ratio male/female 9/5, median age 2 months 24 days [limits, 1 month 0 day -9 months 19 days]), admitted at the Brest University Hospital, Brest, France between January 23th, 2018 and January 8 th , 2020, and diagnosed with Pneumocystis primary infection, were enrolled in the study. Characteristics of the 14 infants are described in table I. All infants were hospitalized in the same pediatric ward, in a single or a two-bed room.

Fourteen air samples were collected in hospital rooms at one-meter distance from infants' head. One cubic meter of air (airflow, 100 L per min. during 10 min.) was collected using the Coriolis® μ air sampler (Bertin Technologies, Montigny-le-Bretonneux, France) which collect and concentrates aerial biological particles into Phosphate-buffered Saline added with 0.002% Tween80. DNA extraction of air samples was performed using the QIAamp DNA Mini Kit (Qiagen, Courtaboeuf, France) according to the recommendations by Choukri et al. [START_REF] Choukri | Quantification and Spread of Pneumocystis jirovecii in the Surrounding Air of Patients with Pneumocystis Pneumonia[END_REF]. DNA was eluted in 100 μL of elution buffer and stored at -80 °C until amplification using the same aforementioned PCR assay [START_REF] Nevez | Pneumocystis primary infection in infancy: Additional French data and review of the literature[END_REF] but with 50 amplification cycles. All air samples and negative controls were run in triplicate whereas positive controls were run in duplicate. An internal positive control (TaqMan® Exogenous Internal Positive Control controls were used to managed putative cross contamination. Moreover, UDG was used to prevent putative cross contamination. The genotyping of P. jirovecii isolates from the nasopharyngal aspirates and air samples was performed by examining the sequences of mtLSUrRNA gene. Genotypes are determined by single nucleotide polymorphisms at positions 85 and 248 [START_REF] Beard | Genetic variation in Pneumocystis carinii isolates from different geographic regions: implications for transmission[END_REF].

The study was noninterventional, which did not require inform consents and ethical approval according to French laws and regulations (CSP Art L1121e1.1).

Statistical analyses were performed using Prism software (version 7.0, GraphPad Software, San Diego, CA, USA). Quantitative variables were compared using Mann-Whitney test. A two-sided P-value < 0.05 was considered significant.

Results and discussion

The results of P. jirovecii DNA detection air samples collected in the surrounding aerial environment of infants' room are presented in table I. In the absence of inhibitors or contamination, P. jirovecii DNA was detected in three out of 14 air samples (21.4%). These three air samples were collected within one day of NPAs collection from infants I4, I9, and I10. P. jirovecii genotyping was successful in only one positive air sample that was collected near infant I9. The allele 2 (A85C248) was identified in both NPA and air samples, objectivizing a genotype match in this pair of samples. This match of genotype was consistent with the fact that P. jirovecii detected in the air sample was effectively exhaled by the infant.

Infant I9 was hospitalized in a single bedroom.

Only 14 infants were enrolled the study, while 87 infants were diagnosed with Pneumocystis primary infection during the study period. The availability of the air sampling device, the disposal of the operator as well as that of the infants may have also hampered in air sampling.

Although the detection of P. jirovecii in NPAs using PCR is a routine technique in our laboratory, it is performed only twice a week. Pneumocystis primary infection is a benign disease frequently contemporaneously to viral infection and with prompt improvement. Most infants infected with P. jirovecii in this context were already discharged from the hospital when the pediatricians were informed of the positive results of P. jirovecii detection in NPAs.

Consequently, in this context the air sampling in infants' room was not be done. The time interval between the date of NPA sampling and that of P. jirovecii detection in NPAs also explains the time interval between the date of NPA sampling and that of air sampling in infants' room. The median of this latter was 3,5 days (range, 1 -10 days). P. jirovecii DNA was detected using the real-time PCR assay in air samples in the rooms of infants I4, I9, and I10 with Ct values of 40.9, 40.6, and 42 respectively. Considering these values, the fungal burdens were probably low and below the quantification limit we previously described elsewhere [START_REF] Pougnet | Pneumocystis jirovecii Exhalation in the Course of Pneumocystis Pneumonia Treatment[END_REF]. We did not perform the quantification in positive air samples for this reason. P. jirovecii DNA was not detected in air samples in the rooms of the remaining 11 infants (I1, I2, I3, I5, I6, I7, I8, I11, I12, I13, I14). Several uncontradictory hypotheses may explain these results: i) P. jirovecii exhalation by these 11 infants was faint or insignificant due to their low P. jirovecii pulmonary burdens. In favor of this hypothesis, it is noteworthy that the median of the Ct values of P. jirovecii DNA amplification in NPAs from these 11 infants was 39.6 whereas it was 36.2 in the three infants from whom P. jirovecii DNA was detected in air samples (P, 0.07, Mann-Whitney test); ii) the time interval between the date of NPA sampling and that of air sampling in infants' room may also be involved. It was one day for the three infants from whom P. jirovecii DNA was detected in air samples whereas it was 4 days (median) (range, 1-10 days) (P, 0.016, Mann-Whitney test) for the 11 remaining infants from whom P. jirovecii DNA was not detected in air samples. This delay in air sampling was concomitant to infant improvement, P. jirovecii clearance from the lungs, and consequently, the decrease or disruption of P. jirovecii exhalation; iii) although, we used the same real time PCR assay that we previously described elsewhere [START_REF] Gal | Pneumocystis jirovecii in the air surrounding patients with Pneumocystis pulmonary colonization[END_REF], a 1-cubic-meter was collected per air sample (100 L per min. during 10 min.) instead of 1.5-cubic meter (300 L per min. during 5 min.) in order to reduce the sound level of the device from 75 dB to 52 dB close to the infants, and to respect the criteria of a non-interventional study. Thus, the air collection of P. jirovecii organisms may have been less efficient.

A plus/minus real-time PCR assay amplifying the mtLSUrRNA gene has been routinely used to perform diagnoses of P. jirovecii infection in our laboratory for many years. For this reason, the same gene was also chosen as the target of the qPCR assay for detection of P. jirovecii in air samples. This multicopy gene provides a highly sensitive amplification, which is required for detection of low fungal burden in NPAs from infants developing primary infection, and in the related air samples. For the same reasons, choosing this gene for P. jirovecii genotyping based on DNA sequence analysis rendered it easy the comparison of putative genotypes in NPAs and air samples. The genotyping was successful in the air sample from infant I9 and unsuccessful in the two other air samples from infants I4 and I10. The vagaries between the real-time PCR-assay that was used for P. jirovecii detection in air samples and the conventional PCR assay that was used for the first step of the method of mtLSUrRNA gene sequencing in this context of low fungal burdens may explain these unsuccessful results. Additionally, some infants were more or less spreaders because they coughed or not, and they were awake or not. This is a hypothesis to explain why fungal loads in air samples were not strictly correlated with fungal loads in NPAs and why detection and genotyping of P. jirovecii failed.

In the course of the present study, none unoccupied rooms in the pediatrics department were available during the winter bronchiolitis epidemic to get negative controls of air samples.

Nonetheless, it was noteworthy, that 11 out of 14 air samples were negative for P. jirovecii detection. Moreover, in our previous studies taken together [START_REF] Choukri | Quantification and Spread of Pneumocystis jirovecii in the Surrounding Air of Patients with Pneumocystis Pneumonia[END_REF][START_REF] Gal | Pneumocystis jirovecii in the air surrounding patients with Pneumocystis pulmonary colonization[END_REF], 43 control air samples were negative providing data excluding putative environmental contaminations.

To sum up, despite the small number of infants enrolled in the study, we obtained the first data on P. jirovecii exhalation by infants with primary infection. This study provides additional data on P. jirovecii human reservoir. Infants may represent potential active infectious sources of the fungus in the hospital environment as previously shown for PCP patients. This observation may have repercussions in terms of infection control in pediatric units monitoring infants with respiratory infections. Infants with primary Pneumocystis infection may also represent potential active sources of infection in the community. Following these original but limited results, a new study considering the rapid turnover of infants hospitalized in pediatric units, and including the daily use of PCR for the detection of P. jirovecii in air samples is warranted.

Table I :

 I Pneumocystis jirovecii detection in the air surrounding infants developing Pneumocystis primary infection

	Infant code (sex)	Age	Clinical presentation	Number of days between NPA sampling and air sampling	Results of P. jirovecii detection in NPA expressed in Ct	Results of P. jirovecii detection in air (Ct)	mtLSUrRNA genotype identification in air samples
	I1 (M)	2 m, 12 days	Bronchiolitis* and pyelonephritis (E. coli)	5	37.8	Negative	-
	I2 (F)	8 m, 24 days	Rhinitis*	3	40.8	Negative	-
			Bronchiolitis
	I3 (M)	3 m	(rhinovirus/enterovirus, Haemophilus influenzae	3	34.9	Negative	-
			concurrent infection)
	I4 (F)	1 m, 23 days	Cow's milk protein allergy**	1	33.4	Positive (40.9)	Unsuccessful amplification
			Acute obstructive
			respiratory deficiency
	I5 (M)	7 m, 12 days	due to rhinopharyngitis,	10	40	Negative	-
			Pierre Robin
			syndrome**
			Bronchiolitis, asthma
	I6 (M)	8 m, 20 days	West Syndrome, brain	9	34.4	Negative	-
			malformation**
	I7 (F)	1 m	Bronchiolitis (virus parainfluenzae type III)	4	39.2	Negative	-

Funding source

This study was supported in part by the European Union (grant number, ERANet-LAC HID-0254).