Didier Henrion

Daniel Liberzon

Alexandre Dolgui

Claude Jard

Stephane Lafortune

Jean-Jacques J Loiseau

Jose-Luis Santi

Alice Pascale

Ma Camille

Emilienne Mère

PhD Santi Esteva

J Aguilar-Martin

J L De La Rosa

J Colomer

J C Hennet

E Garcia

J Melendez

V Puig

Jose-Luis L Villa
email: jvilla@emn.fr

M Morari

K E Arzen

M Duque
email: maduque@uniandes.edu.co

A Gauthier
email: agauthie@uniandes.edu.co

N Rakoto
email: naly.rakoto@mines-nantes.fr

Eduardo Mojica
email: ea.mojica70@uniandes.edu.co

P Caines

Nicanor Quijano

P Riedinger

N Rakoto -Supervision

M Quijano

C Ocampo-Martinez

H Gueguen

Jean-Sebastien Besse

Naly Rakoto-Ravalontsalama

Eduardo Mojica-Nava

Germán Obando

Modelling, Control and Supervision for a Class of Hybrid Systems

Keywords: method of moment, optimal control, switched systems Distributed optimisation, resource allocation, multi-agent systems

dynamique.

Introduction

The HDR (Habilitation à Diriger des Recherches) is a French Degree that you get some years after the PhD. It allows the candidate to apply for some University Professor positions and/or to apply for a Research Director position at CNRS. Instead of explaining it in details, the selection phases process after an HDR is summarized with a Petri Net model in Figure 1. After obtaining the HDR Thesis Degree, the candidate is allowed to apply for a Research Director position at CNRS, after a national selection. On the other hand, in order to apply for some University Professor positions, the candidate should first apply for a National Qualifcation (CNU). Once this qualification obtained, the candidate can then apply to some University Professor positions, with a selection specific to each university. This HDR Thesis is an extended abstract of my research work from my PhD Thesis defense in 1993 until now. This report is organized as follows.

• Chapter 1 is a Curriculum Vitae I had and am having the following administrative responsabilities at Mines Nantes:

• 1997-2000: Last Year's Option AII (Automatique et Informatique Industrielle)

•

Research Activities

My main topics of research are the following:

1. Analysis and control of hybrid and switched systems 2. Supervisory control of discrete-event systems

These will be detailed in Chapter 2 and Chapter 3, respectively. The following other topics of research will not be presented. However, the corresponding papers can be found in the Complete List of Publications.

• Resource Allocation

• Holonic Systems

• Inventory Control -Invited Session, Knowledge Based Systems, IEEE ISIC 1999, Cambridge, MA, USA, Sep. 1999 (jointly organized and chaired with Karl-Erik Årzèn).

Funded and Submitted Projects

-Workshop on G2 Expert System, LAAS-CNRS, Toulouse, France, Oct. 1995 (jointly organized and chaired with Joseph Aguilar-Martin).

Complete List of Publications

A summary of the papers, classified per year, from 1994 to 2017, is given in the following

Chapter 2

Analysis and Control of Hybrid and Switched Systems

Modeling and Control of MLD systems

Piecewise affine (PWA) systems have been receiving increasing interest, as a particular class of hybrid system, see e.g. [2], [13], [11], [16], [14], [12] and references therein. PWA systems arise as an approximation of smooth nonlinear systems [15] and they are also equivalent to some classes of hybrid systems e.g. Linear complementarity systems [9]. On the other hand Mixed Logical and Dynamical (MLD) systems have been introduced by Bemporad and Morari as a suitable representation for hybrid dynamical systems [3]. MLD models are obtained originally from PWA system, where propositional logic relations are transformed into mixed-integer inequalities involving integer and continuous variables. Then mixed-integer optimization techniques are applied to the MLD system in order to stabilize MLD system on desired reference trajectories under some constraints. Equivalences between PWA systems and MLD models have been established in [9]. More precisely, every well-posed PWA system can be rewritten as an MLD system assumung that the set of feasible states and inputs is bounded and a completely well-posed MLD system can be rewritten as a PWA system [9]. Conversion methods from MLD systems to equivalent PWA models have been proposed in [4], [5], [6] and [?]. Vice versa, translation methods from PWA to MLD systems have been studied in [3] (the original one), and then in [8], [?]. A tool that deals with both MLD and PWA systems is HYSDEL [17].

The motivations for studying new methods of conversion from PWA systems into their equivalent MLD models are the following. Firstly the original motivation of obtaining MLD models is to rewrite a PWA system into a model that allows the designer to use existing optimization algorithms such as mixed integer quadratic programming (MIQP) or mixed integer linear programmimg (MILP). Secondly there is no unique formulation of PWA systems. We can always address some particular cases that introduce some differences in the conversions. Finally, it has been shown that the stability analysis of PWA systems with two polyhedral regions is in general NP-complete or undecidable [7]. The conversion to MLD systems may be another way to tackle this problem.

Piecewise Affine (PWA) Systems

A particular class of hybrid dynamical systems is the system described as follows.

{ ẋ(t) = A i x(t) + a i + B i u(t) y(t) = C i x(t) + c i + D i u(t) (2.1)
where i ∈ I, the set of indexes, x(t) ∈ X i which is a sub-space of the real space R n , and R + is the set of positive real numbers including the zero element. In addition to this equation it is necessary to define the form as the system switches among its several modes. This equation is affine in the state space x and the systems described in this form are called piecewise affine (PWA) systems NR HDR 22 [15], [9]. The discrete-time version of this equation will be used in this work and can be described as follows.

{ x(k + 1) = A i x(k) + b i + B i u(k) y(k) = C i x(k) + d i + D i u(k) (2.2)
where i ∈ I is a set of indexes, X i is a sub-space of the real space R n , and R + is the set of positive integer numbers including the zero element, or an homeomorphic set to Z + .

Mixed Logical Dynamical (MLD) Systems

The idea in the MLD framework is to represent logical propositions with the equivalent mixed integer expressions. MLD form is obtained in three steps [3], [4]. The first step is to associate a binary variable δ ∈ {0, 1} with a proposition S, that may be true or false. δ is equal to 1 if and only if proposition S is true. A composed proposition of elementary propositions S 1 , . . . , S q combined using the boolean operators like AND, OR, NOT may be expressed with integer inequalities over corresponding binary variables δ i , i = 1, . . . q. The second step is to replace the products of linear functions and logic variables by a new auxiliary variable z = δa T x where a T is a constant vector.

The variable z is obtained by mixed linear inequalities evaluation. The third step is to describe the dynamical system, binary variables and auxiliary variables in a linear time invariant system. An hybrid system described in MLD form is represented by Equations (2.3-2.5).

x(k + 1) = Ax(k)

+ B 1 u(k) + B 2 δ(k) + B 3 z(k) (2.
3)

y(k) = Cx(k) + D 1 u(k) + D 2 δ(k) + D 3 z(k) (2.4) E 2 δ(k) + E 3 z(k) ≤ E 1 u(k) + E 4 x(k) + E 5 (2.5)
where x = [x T c x T l ] ∈ R nc × {0, 1} n l are the continuous and binary states, respectively, u = [u T C u T l ] ∈ R mc × {0, 1} m l are the inputs, y = [y T c y T l ] ∈ R pc × {0, 1} p l the outputs, and δ ∈ {0, 1} r l , z ∈ R rc , represent the binary and continuous auxiliary variables, respectively. The constraints over state, input, output, z and δ variables are included in (2.5).

Converting PWA into MLD Systems

In this subsection two algorithms for converting PWA systems into MLD systems are given. The first case consists of several sub-affine systems with switching regions are explained in detail. The second case deals with several sub-affine systems, each of them belongs to a region which is described by linear inequalities is a variation of the first case. Each case is applied to an example in order to show the validity of the algorithm.

A. Case I

The PWA system is represented by the following equations:

   x(k + 1) = A i x(k) + B i u(k) + f i y(k) = C i x(k) + D i u(k) + g i S ij = {x, u|k T 1ij x + k T 2ij u + k 3ij ≤ 0} (2.6)
where i ∈ I = {1, . . . , n}. The case with jumps can be included in this representation considering each jump as a discrete affine behavior valid during only one sample time. The switching region S ij is a convex polytope which volume, or hypervolume, can be infinite, and the sub-scripts denotes the switching from mode i to mode j. For this purpose we introduce a binary variable δ i for each index of the set I and a binary variable δ i,j for each switching region S ij . In order to gain insight in the following equations, we consider the hybrid the partition and the corresponding automaton is depicted in Figure 2.1. Introductory material on hybrid automata can be found in [1] and [10]. The δ ij variables are not dynamical and, when the elements k in S ij are vectors, the binary variable can be evaluated by the next mixed integer inequality

(δ ij = 1) ⇔ (k T 1ij x + k T 2ij u + k 3ij ≤ 0) (2.7)
which is equivalent to:

{ k 1ij x + k 2ij u + k 3ij -M (1 -δ ij ) ≤ 0 -k 1ij x -k 2ij u -k 3ij + ϵ + (m -ϵ)δ ij ≤ 0 (2.8)
When the elements k in S ij are matrices, it is necessary to introduce some auxiliary binary variables for each row describing a sub-constraint in S ij in the next form:

δ k = 1(⇔ k 1,k x + k 2,k u + k 3,k ≤ 0) δ ij = ∧ k δ k (2.9)
which is equivalent to:

       k 1ij,k x + k 2ij,k u + k 3ij,k -M (1 -δ ij,k ) ≤ 0 -k 1ij,k x -k 2ij,k u -k 3ij,k + ϵ + (m -ϵ)δ ij,k ≤ 0 δ ij -δ ij,k ≤ 0 ∑ k (δ ij,k -1) -δ ij ≤ -1
(2.10)

The binary vector x δ = [δ 1 δ 2 . . . δ n ] T is such that its dynamics is given by:

x δi (k + 1) = (x δi (k) ∧ ∧ j̸ =i ¬δij) ∨ ∨ j̸ =i (x δj (k) ∧ δji) (2.11)
where k is an index of time, and ∧, ∨, and ¬, are standard for the logical operations AND, OR, NOT, respectively. This equation can be explained as follows: The mode of the system in the next time is i if the current mode is mode i and any switching region is enabled in this time, or, the current mode of the system is j different to i and a switching region that enables the system to go into mode i is enabled. Considering that the PWA system is well posed, i.e. for a given initial state [x T i T ] T 0 and a given input u 0,τ there exists only one possible trajectory [x T i T ] T 0,x . That is equivalent to the following conditions:

∑ i∈I x δi = 1, ∏ i∈I
x δi = 0 (2.12)

The dynamical equations for x δ vector are equivalent to the next integer inequalities:

           x δj (k) + δ ji -x δi (k + 1) ≤ 1, ∀i, j ∈ I, i ̸ = j x δi (k) - ∑ j̸ =i δ ij -x δi (k + 1) ≤ 0, ∀i, j ∈ I, i ̸ = j -x δi (k) - ∑ j̸ =i δ ji -x δi (k + 1) ≤ 0, ∀i, j ∈ I, i ̸ = j (2.13)
The first inequality states that the next mode of the system should be mode i if the current mode is j different to i and a switching region for going from mode j to mode i is enabled. The second inequality means that the next mode of the system should be mode i if the current mode is i and any switching region for going from mode i into mode j different to i is enabled. And the third equation states that the system cannot be in mode i in the next time if the current mode of the system is not mode i and any switching region for going from mode i, (j different to i), into mode i is enabled.

NR HDR 24

This form for finding x δ (k + 1) causes a problem in the final model because it cannot be represented by a linear equation in function of x, u, δ and Z. For this reason, x δ (k + 1) is aggregated to the δ general vector of binary variables, and finally assigned directly to x δ (k + 1). The dynamics and outputs of the system can be represented by the next equations:

{ x(k + 1) = Ax(k) + Bu(k) + ∑ i∈I (A i x(k) + B i u(k) + f i ) × x δi (k) y(k) = Cx(k) + Du(k) + ∑ i∈I (C i x(k) + D i u(k) + g i ) × x δi (k) (2.14)
If we introduce some auxiliary variables:

{ Z 1i (k) = (A i x(k) + B i u(k) + f i ) × x δi (k) Z 2i (k) = (C i x(k) + D i u(k) + g i ) × x δi (k) (2.15)
which are equivalent to:

       Z 1i ≤ M x δi (k) -Z 1i ≤ -mx δi (k) Z 1i ≤ A i x(k) + B i u(k) + f i -m(1 -x δi (k)) -Z 1i ≤ -A i x(k) -B i u(k) -f i + M (1 -x δi (k)) (2.16)        Z 2i ≤ M x δi (k) -Z 2i ≤ -mx δi (k) Z 2i ≤ C i x(k) + D i u(k) + g i -m(1 -x δi (k)) -Z 2i ≤ -C i x(k) -D i u(k) -g i + M (1 -x δi (k)) (2.17)
where M and m are vectors representing the maximum and minimum values, respectively, of the variables Z, these values can be arbitrary large. Using the previous equivalences, the PWA system ( 2.2) can be rewritten in an equivalent MLD model as follows:

                       x(k + 1) = A rr x(k) + A br x δ (k) + B 1r u(k) + B 2r δ + B 3r ∑ i∈I Z 1i (k) x δ (k + 1) = A rb x(k) + A bb x δ (k) + B 1b u(k) + B 2b δ + B 3b ∑ i∈I Z 1i (k) y r (k) = C rr x(k) + C br x δ (k) + D 1r u(k) + D 2r δ + D 3r ∑ i∈I Z 2i (k) y δ (k) = C rb x(k) + C bb x δ (k) + D 1b u(k) + D 2b δ + D 3b ∑ i∈I Z 2i (k) (2.18) s.t. E 2   x δ (k + 1) δ ij δ k   + E 3 Z(k) ≤ E 4   x(k) δ ij δ k   + E 1 u(k) + E 5 (2.19)
Using this algorithm, most part of the matrices are zero, because x and y are defined by Z, and x δ is defined by δ. This situation can be avoided by defining the next matrices at the beginning of the procedure:

       A = 1 n (A 1 + . . . + A n ), A i = A i -A, ∀i ∈ I B = 1 n (B 1 + . . . + B n ), B i = B i -B, ∀i ∈ I C = 1 n (C 1 + . . . + C n ), C i = C i -C, ∀i ∈ I D = 1 n (D 1 + . . . + D n ), D i = D i -D, ∀i ∈ I (2.20)
Finally, the equality matrices in (2.18) and (2.19) can be chosen as follows:

       A rr = A, A br = 0 nc×n , B 1r = B, B 2r = 0 nc×(n+m+tk) , B 3r = [I nc×nc 0 nc×pc I nc×nc 0 nc×pc . . . I nc×nc 0 nc×pc ] nc×n×(nc+pc) A rb = 0 n×nc , A bb = 0 x×n , B 1b = 0 n×mc , B 2b = [I n×n 0 n×m 0 n×tk ] n×n×(n+m+tk) , B 3b = 0 nc×n×(nc+pc) (2.21) NR HDR 25        C rr = C, C br = 0 pc×n , D 1r = D, D 2r = 0 pc×(n+m+tk) , D 3r = [0 pc×nc I pc×pc 0 pc×nc I pc×pc . . . 0 pc×nc I pc×pc ] pc×n×(nc+pc) C rb = 0 n×nc , C bb = 0 x×n , D 1b = 0 n×mc , D 2b = [I n×n 0 n×m 0 n×tk ] n×n×(n+m) , D 3b = 0 n×n×(nc+pc) (2.22)
where n C is the number of continuous state variables, m C the number of continuous input variables, p C the number of continuous output variables, n the number of affine sub-systems, m the number of switching regions and tk the number of auxiliary binary variables. The algorithm for converting a PWA system in the form of (2.1) into its equivalent MLD system can be summarized as follows:

B. Algorithm 1

C. Example 1

Consider the system whose behavior is defined by the following PWA model:

   x(k + 1) = A i x(k), i ∈ {1, 2} S 1,2 = {(x 1 , x 2 )|(x 1 ≤ 1.3x 2 ) ∧ (0.7x 2 ≤ x 1 ) ∧ (x 2 > 0) S 2,1 = {(x 1 , x 2 )|(x 1 ≤ 0.7x 2 ) ∧ (1.3x 2 ≤ x 1 ) ∧ (x 2 < 0)
where

A 1 = [ 0.9802 0.0987 -0.1974 0.9802 ] , A 2 = [ 0.9876 -0.0989 0.0495 0.9876 ]
The behavior of the system is presented in Figure 2.2. The initial points are (x 10 , x 20 ) = (1, 0.8). We can see that the system switches between the two behaviors, from A 1 to A 2 in the switching region S 1,2 , and from A 2 to A 1 in the switching region S 2,1 , alternatively. The switched system is stable.

D. Case 2

Consider now the system whose behavior is defined by the following PWA model:

{ x(k + 1) = A i x(k) + b i + B i u(k), i ∈ I, x(k) ∈ X i y(k) = C i x(k) + d i + D i u(k), i ∈ I, x(k) ∈ X i (2.23) with conditions X i ∩ X j̸ =i = ∅, ∀i, j ∈ I, ∪ i∈I X i = X,
where X is the admissible space for the PWA system, and variables and can be represented using the appropriate δ variables instead of x δ (k) variables in the definition of Z in (2.16) and (2.17). However, note that the conditions X i ∩ X j̸ =i = ∅ ∀i, j ∈ I and ∪ i∈I X i = X require a careful definition in the sub-spaces X i in order to avoid a violation to these conditions in its bounds. On the other hand, the MLD representation uses non-strict inequalities in its representation and the ε value in (2.8) and (2.9) should be chosen appropriately. Another way to overcome this situation and to insure an appropriated representation is the use of the following conditions in the bounds of the sub-spaces X i :

X i = {x, u|k 1i x + k 2i u + k 3i ≤ 0}
δ ij = δ i ⊗ δ j
which is equivalent to:

{ δ i + δ j -1 ≤ 0 1 -δ i -δ j ≤ 0 or more generally        ∑ i∈I δ i -1 ≤ 0 1 - ∑ i∈I δ i ≤ 0 (2.24)
We now modify Equations (2.8), (2.10), (2.16), (2.17), (2.21), and (2.22) as follows:

{ k 1i x + k 2i u + k 3i -M (1 -δ i ) ≤ 0 -k 1i x -k 2i u -k 3i + ϵ + (m -ϵ)δ i ) ≤ 0 (2.25) NR HDR 27        k 1i,k x + k 2i,k u + k 3i,k -M (1 -δ i,k ) ≤ 0 -k 1i,k x -k 2i,k u -k 3i,k + ϵ + (m -ϵ)δ i,k ≤ 0 δ i -δ i,k ≤ 0 ∑ k (δ i,k -1) -δ i ≤ -1 (2.26)
The auxiliary variables Z 1i become:

       Z 1i ≤ M δ i (k) -Z 1i ≤ -mδ i (k) Z 1i ≤ A i x(k) + B i u(k) + f i -m(1 -δ i (k)) -Z 1i ≤ -A i x(k) -B i u(k) -f i + M (1 -δ i (k)) (2.27)
where the matrices A i and B i are those previously defined in Equation (2.20).The auxiliary variable Z 2i is now modified according to the following equations:

       Z 2i ≤ M δ i (k) -Z 2i ≤ -mδ i (k) Z 2i ≤ C i x(k) + D i u(k) + g i -m(1 -δ i (k)) -Z 2i ≤ -C i x(k) -D i u(k) -g i + M (1 -δ i (k)) (2.28)
where the matrices C i and D i are those that have been defined in Equation (2.20). Finally the matrices from Equation (2.18) can be chosen as follows:

       A rr = A, A br = 0 nc×n , B 1rr = B, B 2rb = 0 nc×(n+tk) , B 3rr = [I nc×nc 0 nc×pc I nc×nc 0 nc×pc . . . I nc×nc 0 nc×pc ] nc×n×(nc+pc) C rr = C, C br = 0 pc×n , D 1rr = D, D 1rb = [ ], D 2rb = 0 pc×(n+tk) , D 3rr = [0 pc×nc I pc×pc 0 pc×nc I pc×pc . . . 0 pc×nc I pc×pc ] pc×n×(nc+pc) (2.29)
We give now an algorithm that converts a PWA system in the form of (2.23) into its equivalent MLD system.

E. Algorithm 2

1. Compute matrices A, B, C, D and A i , B i , C i and D i using (2.20).

Initialize

E 1 , E 2 , E 3 , E 4 , E 5 matrices.
3. For i = 1 to n include the inequalities using (2.25) or (2.26) that represent the behavior on the n affine regions of the PWA system.

4. For all affine regions include the inequalities in (2.24).

5. For i = 1 to n generate the nc-dimensional Z 1i vector and p c -dimensional Z 2i vector of auxiliary variables Z.

6. For each Z 1i vector introduce the inequalities defined in (2.27). M and m are n c -dimensional vectors of maximum and minimum values of x, respectively.

7. For each Z 2i vector introduce the inequalities defined in (2.28). M and m are p c -dimensional vectors of maximum and minimum values of x, respectively (This completes the inequality matrices).

8. Compute the matrices defined in (2.29) where the binary state variables are removed. 9. End.

F. Example 2

Consider the system whose behavior is defined by the following PWA model:

   x(k + 1) = A i x(k), i ∈ {1, 2} i = 1 if x 1 x 2 ≥ 0 i = 2 if x 1 x 2 < 0 NR HDR 28
where The behavior of the system is presented in Figure 2.4. The PWA system with linear constraints has 4 sub-affine systems. Algorithm 2 produces an MLD system with 12 binary variables (4 variables for the affine sub-system, and 8 auxiliary variables), 16 auxiliary variables Z and 94 constraints. The behavior of the equivalent MLD system is shown in Figure 2.5. We can notice that the behavior of the MLD system is exactly the same as the original PWA model.

A

Stability of Switched Systems

A polynomial approach to deal with the stability analysis of switched non-linear systems under arbitrary switching using dissipation inequalities is presented. It is shown that a representation of the original switched problem into a continuous polynomial system allows us to use dissipation inequalities for the stability analysis of polynomial systems. With this method and from a theoretical point of view, we provide an alternative way to search for a common Lyapunov function for switched non-linear systems. We deal with the stability analysis of switched non-linear systems, i.e., continuous systems with switching signals under arbitrary switching. Most of the efforts in switched systems research have been typically focused on the analysis of dynamical behavior with respect to switching signals. Several methods have been proposed for stability analysis (see [START_REF] Liberzon | Switching in Systems and Control[END_REF], [19], and references therein), but most of them have been focused on switched linear systems. Stability analysis under arbitrary switching is a fundamental problem in the analysis and design of switched systems. For this problem, it is necessary that all the subsystems be asymptotically stable. However, in general, the above stability condition is not sufficient to guarantee stability for the switched system under arbitrary switching. It is well known that if there exists a common Lyapunov function for all the subsystems, then the stability of the switched system is guaranteed under arbitrary switching. Previous attempts for general constructions of a common Lyapunov function for switched non-linear systems have been presented in [20], [21] using converse Lyapunov theorems. Also in [22], a construction of a common Lyapunov function is presented for a particular case when the individual systems are handled sequentially rather than simultaneously for a family of pairwise commuting systems. These methodologies are presented in a very general framework, and even though they are mathematically sound, they are too restrictive from a computational point of view, mainly because it is usually hard to check for the set of necessary conditions for a common function over all the subsystems (it could not exist). Also, these constructions are usually iterative, which involves running backwards in time for all possible switching signals, being prohibitive when the number of modes increases.

The main contribution of this topic of stability of switched systems is twofold. First, we present a reformulation of the switched system as an ordinary differential equation on a constraint manifold. This representation opens several possibilities of analysis and design of switched systems in a consistent way, and also with numerical efficiency [C.39], [C.38], which is possible thanks to some tools developed in the last decade for polynomial differential-algebraic equations analysis [8,10]. The second contribution is an efficient numerical method to search for a common Lyapunov function for switched systems using results of stability analysis of polynomial systems based on dissipativity theory [23], [C.39]. We propose a methodology to construct common Lyapunov functions that provides a less conservative test for proving stability under arbitrary switching. It has been mentioned in [26] that the sum of squares decomposition, presented only for switched polynomial systems, can sometimes be made for a system with a non-polynomial vector fields. However, those cases are restricted to subsystems that preserve the same dimension after a recasting process.

Optimal Control of Switched Systems

Switched Linear Systems

A polynomial approach to solve the optimal control problem of switched systems is presented. It is shown that the representation of the original switched problem into a continuous polynomial systems allow us to use the method of moments. With this method and from a theoretical point of view, we provide necessary and sufficient conditions for the existence of minimizer by using particular features of the minimizer of its relaxed, convex formulation. Even in the absence of classical minimizers of the switched system, the solution of its relaxed formulation provide minimizers.

We consider the optimal control problem of switched systems, i.e., continuous systems with switch-NR HDR 30

ing signals. Recent efforts in switched systems research have been typically focused on the analysis of dynamic behaviors, such as stability, controllability and observability, etc. (e.g., [19], [START_REF] Liberzon | Switching in Systems and Control[END_REF]). Although there are several studies facing the problem of optimal control of switched systems (both from theoretical and from computational point of view [37], [36], [27], [39], there are still some problems not tackled, especially in issues where the switching mechanism is a design variable. There, we see how these difficulties arise, and how tools from non-smooth calculus and optimal control can be combined to solve optimal control problems. Previously, the approach based on convex analysis have been treated in [36], and further developed in [27], considering an optimal control problem for a switched system, these approaches do not take into account assumptions about the number of switches nor about the mode sequence, because they are given by the solution of the problem. The authors use a switched system that is embedded into a larger family of systems and the optimal control problem is formulated for this family. When the necessary conditions indicate a bang-bangtype of solution, they obtain a solution to the original problem. However, in the cases when a bang-bang type solution does not exist, the solution to the embedded optimal control problem can be approximated by the trajectory of the switched system generated by an appropriate switching control. On the other hand, in [36] and [34] the authors determine the appropriated control law by finding the singular trajectory along some time with non null measure.

Switched Nonlinear Systems

The nonlinear, non-convex form of the control variable, prevents us from using the Hamilton equations of the maximum principle and nonlinear mathematical programming techniques on them. Both approaches would entail severe difficulties, either in the integration of the Hamilton equations or in the search method of any numerical optimization algorithm. Consequently, we propose to convexify the control variable by using the method of moments in the polynomial expression in order to deal with this kind of problems. In this paper we present a method for solving optimal control for an autonomous switched systems problem based on the method of moments developed in for optimal control, and in [28], [29], [30] and [32] for global optimization. We propose an alternative approach for computing effectively the solution of nonlinear, optimal control problems. This method works properly when the control variable (i.e., the switching signal) can be expressed as polynomials. The essential of this paper is the transformation of a nonlinear, non-convex optimal control problem (i.e., the switched system) into an equivalent optimal control problem with linear and convex structure, which allows us to obtain an equivalent convex formulation more appropriate to be solved by high performance numerical computing. To this end, first of all, it is necessary to transform the original switched system into a continuous non-switched system for which the theory of moments is able to work. Namely, we relate with a given controllable switched system, a controllable continuous non-switched polynomial system. Optimal control problems for switched nonlinear systems are investigated. We propose an alternative approach for solving the optimal control problem for a nonlinear switched system based on the theory of moments. The essence of this method is the transformation of a nonlinear, non-convex optimal control problem, that is, the switched system, into an equivalent optimal control problem with linear and convex structure, which allows us to obtain an equivalent convex formulation more appropriate to be solved by high-performance numerical computing. Consequently, we propose to convexify the control variables by means of the method of moments obtaining semidefinite programs. The paper dealing with this approach is given in the Appendix 2, paper [J.5].

Chapter 3

Supervisory Control of Discrete-Event Systems

Multi-Agent Based Supervisory Control

Supervisory control initiated by Ramadge and Wonham [START_REF] Ramadge | Supervisory control of a class of discrete-event processes[END_REF] provides a systematic approach for the control of discrete event system (DES) plant. The discrete event system plant be is modeled by a finite state automaton [START_REF] Hopcroft | Introduction to Automata Theory, Languages, and Computation[END_REF], [43]:

Definition 1 (Finite-state automaton). A finite-state automaton is defined as a 5-tuple

G = (Q, Σ, δ, q 0 , Q m , C)
where

• Q is the finite set of states,
• Σ is the finite set of events,

• δ : Q × Σ → Q is the partial transition function, • q 0 ⊆ Q is the initial state, • Q m ⊆ Q is the set of marked states (final states),
Let Σ * be the set of all finite strings of elements in Σ including the empty string ε. The transition function δ can be generalized to δ : Σ * × Q → Q in the following recursive manner:

δ(ε, q) = q δ(ωσ, q) = δ(σ, δ(ω, q)) for ω ∈ Σ *
The notation δ(σ, q)! for any σ ∈ Σ * and q ∈ Q denotes that δ(σ, q) is defined. Let L(G) ⊆ Σ * be the language generated by G, that is,

L(G) = {σ ∈ Σ * |δ(σ, q 0 )!}
Let K ⊆ Σ * be a language. The set of all prefixes of strings in K is denoted by pr(K) with

pr(K) = {σ ∈ Σ * |∃ t ∈ Σ * ; σt ∈ K}. A language K is said to be prefix closed if K = pr(K).
The event set Σ is decomposed into two subsets Σ c and Σ uc of controllable and uncontrollable events, respectively, where Σ c ∩ Σ uc = ∅. A controller, called a supervisor, controls the plant by dynamically disabling some of the controllable events.

A sequence σ 1 σ 2 . . . σ n ∈ Σ * is called a trace or a word in term of language. We call a valid trace a path from the initial state to a marked state (δ(ω, q 0 ) = q m where ω ∈ Σ * and q m ∈ Q m ).

NR HDR 32

In this section we will focus on the Multi-Agent Based Supervisory Control, introduced by Hubbard and Caines [START_REF] Hubbard | Initial investigations on hierachical supervisory control of multi-agent systems[END_REF]; and the modified approach proposed by Takai and Ushio [START_REF] Takai | Supervisory control of a class of concurrent discrete-event systems[END_REF]. The two approaches have been applied to the supervisory control of the EMN Experimental Manufacturing Cell. This cell is composed of two robotized workstations connected to a central conveyor belt. Then, three new semi-automated workstations have been added in order to increase the flexibility aspects of the cell. Indeed, each semi-automated workstation can perform either manual of robotized tasks. These two aspects correspond to the two different approaches of multi-agent product of subsystems, for supervisory control purpose. The results can be found in [C.25].

Switched Discrete-Event Systems

The notion of switched discrete-event systems corresponds to a class of DES where each automaton is the composition of two basic automata, but with different composition operators. A switching occurs when there is a change of the composition operator, but keeping the same two basic automata. A mode behavior, or mode for short, is defined to be by the DES behavior for a given composition operator. Composition operators are supposed to change more than once so that each mode is visited more than once. This new class of DES includes the DES in the context of fault diagnosis where different modes such as e.g., normal, degenerated, emergency modes can be found. The studied situations are the ones where the DES switch between different normal modes, and not necessary the degenerated and the emergency ones.

The most common composition operators used in supervisory control theory are the product and the parallel composition [43], [START_REF] Wonham | Notes on Discrete Event Systems[END_REF] However many different types of composition operators have been defined, e.g., the prioritized synchronous composition [49], the biased synchronous composition [START_REF] Lafortune | The infimal closed controllable superlanguage and its application to supervisory control[END_REF], see [START_REF] Wenck | On composition oriented perspective on controllability of large DES[END_REF] for a review of most of the composition operators. Multi-Agent composition operator [START_REF] Romanovski | On the supervisory control of multi-agent product systems[END_REF], [START_REF] Romanovski | Multi-agent product system: Controllability and non-blocking properties[END_REF] is another kind of operator, which differs from the synchronous product in the aspects of simultaneity and synchronization.

The new class of DES that we define in this paper includes the class of DES in the context of fault diagnosis, with different operating modes. Furthermore this new class addresses especially the DES for which the system can switch from a given normal mode, to another normal mode. More

precisely this new class of DES is an automaton which is the composition of two basic automata, but with different composition operators. A switching corresponds to the change of composition operator, but the two basic automata remains the same. A mode behavior (or mode for short) is defined to be the DES situation for a given composition operator. Composition operators are supposed to change more than once so that each mode is visited more than once.

We give here below some examples of switched DES:

• Manufacturing systems where the operating modes are changing (e.g. from normal mode to degenerated mode)

• Discrete event systems after an emergency signal (from normal to safety mode)

• Complex systems changing from normal mode to recovery mode (or from safety mode to normal mode).

We can distinguish, like for the switched continuous-time systems, the notion of autonomous switching where no external action is performed and the notion of controlled switching, where the switching is forced. The results for this section can be found in [START_REF] Rakoto-Ravalontsalama | Supervisory control of switched discrete-event systems[END_REF].

Switchable Languages of DES

The notion of switchable languages has been defined by Kumar, Takai, Fabian and Ushio in [Kumaret-al. 2005]. It deals with switching supervisory control, where switching means switching between two specifi-cations. In this paper, we first extend the notion of switchable languages to n languages, (n ≥ 3). Then we consider a discrete-event system modeled with weighted automata. The switching supervisory control strategy is based on the cost associated to each event, and it allows us to synthesize an optimal supervisory controller. Finally the proposed methodology is applied to a simple example. We now give the main results of this paper. First, we define a triplet of switchable languages. Second we derive a necessary and sufficient condition for the transitivity of switchable languages (n = 3). Third we generalize this definition to a n-uplet of switchable languages, with n > 3. And fourth we derive a necessary and sufficient condition for the transitivity of switchable languages for n > 3.

Triplet of Switchable Languages

We extend the notion of pair of switchable languages, defined in [START_REF] Kumar | Maximally Permissive Mutually and Globally Nonblocking Supervision with Application to Switching Control[END_REF], to a triplet of switchable languages.

Definition 2 (Triplet of switchable languages). A triplet of languages

(K 1 , K 2 , K 3 ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, 2,
3} are said to be a triplet of switchable languages if they are pairwise switchable languages, that is,

SW (K 1 , K 2 , K 3 ) := SW (K i , K j ), i ̸ = j, i, j = {1, 2, 3}.
Another expression of the triplet of switchable languages is given by the following lemma.

Lemma 1 (Triplet of switchable languages). A triplet of languages

(K 1 , K 2 , K 3 ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, 2,
3} are said to be a triplet of switchable languages if the following holds: The following theorem gives a necessary and sufficient condition for the transitivity of switchable languages.

SW (K 1 , K 2 , K 3 ) = {(H 1 , H 2 , H 3 ) | H i ⊆ K i ∩ pr(H j ), i ̸ = j,
Theorem 1 (Transitivity of switchable languages, n = 3) . Given 3 specifications The proof can be found in [42].

(K 1 , K 2 , K 3 ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, 2, 3} such that SW (K 1 , K 2 ) and SW (K 2 , K 3 ). (K 1 , K 3 )

N-uplet of Switchable Languages

We now extend the notion of switchable languages, to a n-uplet of switchable languages, with (n > 3).

Definition 3 (N-uplet of switchable languages, n > 3). A n-uplet of languages

(K 1 , ..., K n ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, ..., n}, n > 2
, is said to be a n-uplet of switchable languages if the languages are pairwise switchable that is,

SW (K 1 , ..., K n ) := SW (K i , K j ), i ̸ = j, i, j = {1, ..., n}, n > 2.
As for the triplet of switchable languages, an alternative expression of the n-uplet of switchable languages is given by the following lemma.

Lemma 2 (N-uplet of switchable languages, n > 3). A n-uplet of languages

(K 1 , . . . , K n ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, .
.., n}, n > 3 are said to be a n-uplet of switchable languages if the following holds:

SW (K 1 , ..., K n ) = {(H 1 , ..., H n ) | H i ⊆ K i ∩ pr(H j ), i ̸ = j
, and H i controllable}.

Transitivity of Switchable Languages (n > 3)

We are now able to derive the following theorem that gives a necessary and sufficient condition for the transitivity of n switchable languages.

Theorem 2 (Transitivity of n switchable languages, n > 3) . Given n specifications (K 1 , ..., K n ),

K i ⊆ L m (G) with H i ⊆ K i , i = {1, ..., n}.
Moreover, assume that each language K i is at least switchable with another language K j , i

̸ = j. A pair of languages (K k , K l ) is switchable i.e. SW (K k , K l ), if and only if 1. H k ∩ pr(H l ) = H k , and 2. H l ∩ pr(H k ) = H l .
The proof is similar to the proof of Theorem 6 and can be found in [42]. It is to be noted that the assumption that each of the n languages be at least switchable with another language is important, in order to derive the above result.

Perspective 1: Control of Smart Grids

According to the US Department of Energy's Electricity Advisory Committee, "A Smart Grid brings the power of networked, interactive technologies into an electricity system, giving utilities and consumers unprecedented control over energy use, improving power grid operations and ultimately reducing costs to consumers."

The transformation from traditional electric network, with centralized energy production to complex and interconnected network will lead to a smart grid. The five main triggers of Smart grid, according to a major industrial point of view, are 1) Smart energy generation, 2) Flexible distribution, 3) Active energy efficiency, 4) Electric vehicles, and 5) Demand response. From a control point of view, a smart grid is a system of interconnected micro-grids. A microgrid is a power distribution network where generators and users interact. Generators technologies include renewable energy such as wind turbines or photovoltaic cells. The objective of this project is to simulate and control a simplified model of a micro-grid that is a part of a Smart Grid. After a literature review, a simplified model for control will be chosen. Different realistic scenarios will be tested in simulation with MATLAB. Finally different NR

Perspective 2: Simulation with Stochastic Petri Nets

The Air France CDG Airport Hub in Paris-Roissy is dealing daily with 40,000 transfer luggages and 30,000 local luggages (leaving from or arriving at CDG Airport). For this purpose Air France is exploiting the Sorting Infrastructure of Paris Aeroport, and has to propose a Logistical Scheme Allocation for each luggage in order to optimize the sorting and to minimize the number of failed luggages. By failed luggages, we mean a luggage that does not arrive in time for the assigned flight.

The KPI Objective for 2017 is to have less than 20 failed luggages out of 1000 passengers. [2] arise as a suitable representation for Hybrid Dynamical Systems (HDS), in particular for solving control-oriented problems. MLD models can be used for solving a model predictive control (MPC) problem of a particular class of HDS and it is proved that MLD models are equivalent to PieceWise Affine Models in [6]. In the paper by Heemels and coworkers, the equivalencies among PieceWise Affine (PWA) Systems, Mixed Logical and Dynamical (MLD) systems, Linear Complementarity (LC) systems, Extended Linear Complementarity (ELC) systems and Max-Min-Plus-Scaling (MMPS) systems are proved, these relations are transcribed here in Fig. 1.

I. INTRODUCTION Mixed and Logical Dynamical (MLD) models introduced by Bemporad and Morari in

This equivalences are based on some propositions (see [6] for details) A more formal proof can be found in [3], where an efficient technique for obtaining a PWA representation of a MLD model is proposed.

The technique in [3] describes a methodology for obtaining, in an efficient form, a partition of the state-input space. The algorithm in [3] uses some tools from polytopes theory in order to avoid the enumeration of the all possible combinations of the integer variables contained in the MLD model. However, the technique does not describe the form to obtain a suitable choice of the PWA model, even though this part is introduced in the implementation provided by the author in [4]. The objective of this paper is to propose an algorithm of the suitable choice of the PWA description and use the PWA description for obtaining some analysis and control of Hybrid Dynamical Systems.

II. MLD SYSTEMS AND PWA SYSTEMS

A. Mixed and Logical Dynamical (MLD) Systems

The idea in the MLD framework is to represent logical propositions as equivalent integer expressions. MLD form is obtained by three basic steps [5]. The first step is to associate a binary variable δ ∈{0,1} with a proposition S, that may be true or false. δ is 1 if and only if proposition S is true. A composed proposition of elementary propositions S 1 ,…,S q combined using the boolean operators like AND(^), OR (∨), NOT(~) may be expressed like integer inequalities over corresponding binary variables δ i , i=1,…,q.

The second step is to replace the products of linear functions and logic variables by a new auxiliary variable z = δa T x where a T is a constant vector. The z value is obtained by mixed linear inequalities evaluation.

The third step is to describe the dynamical system, binary variables and auxiliary variables in a linear time invariant (LTI) system.

A hybrid system MLD described in general form is represented by (1).

1 2 3 1 2 3 2 3 1 4 5 ( 1 ) ( ) ( ) ( ) ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) x k Ax k B u k B k B z k y k Cx k D u k D k D z k E k E z k E u k E x k E δ δ δ   + = + + +     = + + +     + ≤ + +    ) (1)
where l are the continuous and

[ ] { 0 ,1 } c T T n n C l x x x = ∈ × R
binary states, u u are the inputs, the outputs, and , , represent the binary and continuous auxiliary variables, respectively. The constraints over state, input, output, z and δ variables are included in the third term in (1).

[ ] { 0 ,1 } c T T m m C l u = ∈ × R { 0 ,1 } c l p × , , i i i i i i i x a Bu x u X i x c Du = + + × ∈ ∈ = + + ) ( ) , ( ) i i i i i i i b B u k x u d D u k + + × + l )  [ ] T T p C l y y y = ∈ R x A y C        ( 1 ) ( ( ) ( )
x k Ax k y k C x k  + =     = +   {0,1} l r δ ∈ , t + ∈ Ι , X i k ∈ ∈I c r z ∈ R , + ∈Z

B. PieceWise Affine Systems

A particular class of hybrid dynamical systems is the system described as follows,

(

) 2
where I is a set of indexes, X i is a sub-space of the real space R n , and R + is the set of positive real numbers including the zero element.

In addition to this equation it is necessary to define the form as the system switches among its several modes. This equation is affine in the state space x and the systems described in this form are called PieceWise Affine Systems (PWA). In the literature of hybrid dynamical systems the systems described by the autonomous version of this representation are called Switched Systems.

If the system vanishes when x brings near to zero, i.e. a i and b i are zero, then the representation is called PieceWise Linear (PWL) system.

The discrete-time version of this equation will be used in this work and can be described as follows,

where I is a set of indexes, X i is a sub-space of the real space R n .

III. MLD SYSTEMS INTO PWA SYSTEMS

The MLD framework is a powerful structure for representing hybrid systems in an integrated form. Although E 1 , E 2 , E 3 , E 4 and E 5 matrices are, in general, large matrices, they can be obtained automatically. An example is the HYSDEL compiler [10].

However, some analysis of the system with the MLD representation are computationally more expensive with respect to some tools developed for PWA representations. Exploiting the MLD and PWA equivalencies, it is possible to obtain analysis and control of a system using this equivalent representations. Nevertheless, as it is underlined in [3], this procedure is more complex with respect to the PWA into MLD conversion, and there exist more assumptions. To our knowledge, the only previous approach has been proposed by Bemporad [3]. We propose then a new approach of translating MLD into PWA systems.

The MLD structure can be rewritten as follows,

1 1 2 3 1 1 2 3 2 3 2 1 1 5 ( ) ( 1 ) ( ) ( ) ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) c c l l c c l l c c l l u k x k Ax k B B B k B z k u k u k y k Cx k D D D k D z k u k u k E k E z k E x k E E E u k δ δ δ          + = + + +                         = + + +                        + ≤ + +               (4)
Here, the binary inputs are distinguished from the continuous inputs, because they induce switching modes in the system, in general. Supposing that the system is well posed, z(k) has only one possible value for a given x(k) and u(k), and can be rewritten as:

1 2 3 ( ) ( ) ( ) | [ , ] T T T c z k k x k k u k k m x u b = + + ≤ (5)
Replacing this value in the original equations the system can be represented as,

3 1 1 3 2 3 3 3 1 1 3 2 3 3 4 31 1 32 5 2 33 ( 1 ) ( ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( ) ( ) c c x k A B k x k B B k u k B k y k C D k x k D D k u k D k E E k x E E k u E E E k δ   + = + + + +     + = + + + +     -+ +-+ ≤ - -    (6)
If an enumeration technique is used for generating all the feasible binary states of the [u l T δ T ] T vector, the first problem is to find a value of [x T u T ] T feasible for the problem, that can be obtained solving the linear programming problem,

1 3 4 5 2 1 min [ ] .
.

T T T T c c l l X u z x s t E u E z E x E E E u δ  =      - + - ≤ - +    (7)
The solution is a feasible value [x *T u *T ] T . The next problem is to find k 1 , k 2 and k 3 .

The inequalities can be rewritten as,

3 4 1 1 2 5 41 1 2 5 c c l l c c E z E x E u E u E E E k x E k u E k δ ≤ + + - + = + + 3 (8)
where 5 E includes every constant in the problem, i.e. u l and δ. On the other hand, the E 3 matrix reflects the interaction among the z variables, and we can write:

1 2 F z k x k u k 3 × ≤ + + (9)
The matrix F represents the interaction among the z variables, if the system is well posed F -1 should exist.

With this last equation, for finding 3 k the next linear programming problem is solved,

3 * * 3 3 5 2 1 max . . l l k s t E k E E E u δ       ≤ - +    (10)
The solution to this problem is 3 k , in this case we assume that all components in 5 E are the maximum and minimum values of z and the only solution for the problem is 3

k . With 3 k we can obtain the other matrices.

For obtaining 1 k it is necessary to solve n x , i.e. the length of the state vector, linear programming problems,

* * 3 5 2 1 max . . i i l k s t E k E E E u E δ       ≤ - + +    4 l i (11)
where E 4i represents the column i of the E 4 matrix and

1 i i k k = -3
k is the column i of the matrix 1 k .

For obtaining 2 k it is necessary to solve n u , i.e. the length of the continuous input vector, linear programming problems,

* * 3 5 2 1 max . . i i l k s t E k E E E u E δ       ≤ - + +    1 l c i ( 12
)
where E 1ci represents the column i of the E 1c matrix and

2 i i k k k = -3 is the column i of the matrix 2
k . The matrix F should be found solving n z , i.e. the length of the z vector, linear programming problems,

* * 3 5 2 1 max . . i i l k s t E k E E E u E δ       ≤ - + +    3 l i ( 13
)
where E 3i represents the column i of the E 3 matrix and

3 i i F k k = -is the column i of the matrix .
F Finally, k 1 , k 2 , and k 3 , can be computed as,

1 1 1 2 2 1 3 3 k F k k F k k F k - - -   =     =     =    1 (14)
With these equations, the algorithm for translating the MLD model into PWA model is given as follows, Algorithm 1 1. Find a feasible point for the binary vector, composed by the binary inputs and binary auxiliary variables. 2. Compute 3 k using Eq. ( 10). 3. Compute 1 k , 2 k and F using Eq. ( 11), ( 12) and (13).

4. Compute k 1 , k 2 , and k 3 using Eq. ( 14). 5. Using Eq. ( 6), compute A i , B i , f i , C i , D i and g i and the valid region for this representation. 6. If there exists another feasible point go to step 1. 7. End. Some gains in the algorithm performance can be obtained if the vector z is evaluated after step one, using a linear program for finding the maximum and the minimum in z, if the z min and z max solutions are the same, it is not necessary to calculate steps 3, and 4, and z = z min = z max can be assigned directly.

IV. EXAMPLES

A. The Three-Tank Benchmark Problem

The three-tank benchmark problem has been proposed as an interesting hybrid dynamical system. This Benchmark was proposed in [7] and [8]. See [13] and references there in for some control results using MLD framework in this system. The algorithm described in the last section is used for obtaining a PWA representation of this system.

This system has three tanks each of them interconnected with another as depicted in Fig. 2.

Fig. 2. Three Tank System

The model is written using binary variables (δ i ) and relational expressions,

01 1 1 1 1 2 2 0 2 2 3 3 03 3 3 ( ) 1 1 ( 1 ( ) v v v v v v Z h h h h h h Z h h h h Z h h δ δ δ δ δ δ     = - = ↔ >         = ↔ > = -         = ↔ > = -       2 ) 1 2 1 0 1 0 3 2 0 2 0 3 13 1 3 13 23 2 3 23 ( ) ( ) ( ) ( )
Z Z Z V Z Z Z V Z h h V Z h h V  = -     = -      = -     = -    1 1 1 3 1 1 1 1 1 1 31 1 ( ) ( ) ( ) ( ) ( 1 ) ( ) *( s L q k h k Z k Z k h k h k T C RC R C RC + = + - - - 1 ) 2 2 3 2 2 2 2 2 3 2 2 ( ) ( ) ( ) ( 1 ) ( ) *( q k Z k Z k h k h k Ts C R C RC + = - - - 2 3 1 3 2 3 1 2 3 3 3 1 3 3 2 3 2 1 3 2 ( ) ( ) ( ) ( ) ( ) ( 1 ) ( ) *( N h k Z k Z k Z k Z k h k h k Ts R C R C R C RC R C + = + - + + + + 3
)

)
The simulation of the system using the MLD framework and a Mixed Integer Quadratic Programming MIQP algorithm running in an Intel Celeron 2GHz processor and 256MB of RAM was 592.2s, using the PWA representation the same simulation was 1.33s. The time for obtaining the PWA model using the technique described in this work is 72.90s and the algorithm found 128 regions. Using the algorithm in [4] the computation time of the PWA form was 93.88s and the total regions found was 100 and the simulation took 5.89s. These results are summarized in Table I.

Where Computation Time is the time taken by the computer for computing the PWA model based in the MLD model, and Simulation Time is the time taken by the computer for computing a trajectory given a model, an initial state and an input sequence. The simulation results with MLD model and the error between PWA simulation results and MLD simulation results, for the same input are shown in Fig. 3 In this case, at t=30s, the simulation with the PWA system in the Figure 3.b produces a switching to an invalid operation mode.

B. Car with Robotized Manual Gear Shift

The example of a Hybrid Model of a Car with Robotized Manual Gear Shift was reported in [9] and is used in [3] as example. The car dynamics is driven by the following equation,

e b mx F F x β = -- ( 15
)
where m is the car mass, x and x is the car speed and acceleration, respectively, F e is the traction force, F b is the brake force and β is the friction coefficient. The Transmission Kinematics are given by,

( ) ( ) g s g e s R i x k R i F M k ω = =
where ω is the engine speed, M the engine torque and i is the gear position.

The engine torque M is restricted to belongs between the minimum engine torque C and the maximum engine torque C .

( ) e ω - ( ) e ω +
The model has two continuous states, position and velocity of car, two continuous inputs, engine torque and breaking force, and six binary inputs, the gear positions. The MLD model was obtained using the HYSDEL tool.

The translation of the MLD model took 155.73 s and the PWA model found 30 sub-models, using the algorithm proposed in this work, and the PWA model using the algorithm proposed in [3] took 115.52 s and contains 18 sub-models. The simulation time with MLD model and a MIQP algorithm for 250 iterations took 296.25s, using the PWA model obtained with the algorithm proposed here took 0.17s, and using the PWA model obtained using the algorithm in [4] the simulation took 0.35s. These results are summarized in Table II

C. The Drinking Water Treatment Plant

The example of a Drinking Water Treatment Plant has been reported in [11] and [12]. This plant was modeled using identification techniques for hybrid dynamical systems, and its behavior includes autonomous jumps.

The plant modeled is based in the current operation of drinking water plant Francisco Wiesner situated at the periphery of Bogotá D.C. city (Colombia), which treats on average 12m 3 /s. The volume of water produced by this plant is near to 60% of consumption by the Colombian capital. In this plant, there exist two water sources: Chingaza and San Rafael reservoirs which can provide till 22m 3 /s of water.

The process mixes inlet water with a chemical solution in order to generate aggregated particles that can be caught in a filter. The dynamic of the filter is governed by the differential pressure across the filter and the outlet water turbidity. An automaton associated to the filter executes a back-washing operation when the filter performance is degraded. Because of process non-linearity, the behavior of the system is different with two water sources, that is the case for the particular plant modeled.

The model for each water source includes a dynamic for the aggregation particle process which dynamical variable is called Streaming Current (SC) and is modeled using two state variables, a dynamic for the differential pressure called Head Loss (HL) with only one state variable, a dynamic for the outlet turbidity (T o ) with two state variables.

The identified model consists of four affine models, two for each water source in normal operation, one model in maintenance operation, one model representing the jump produced at the end of the maintenance operation.

i i i i i i x k Ax k B u k f y k C x k D u k g i  + = + +     = + +   ∈ to normal operation}
where water source is an input variables, maintenance operation is executed if outlet turbidity (T o ) is greater than a predefined threshold, or, Head Loss (HL) is greater than a predefined threshold, or, operation time is greater than a predefined threshold.

The MLD model has 7 continuous states (including two variables for two timers in the automaton), 4 continuous inputs (dosage, water flow, inlet turbidity and pH), 3 binary inputs (water source, back-washing operation and normal operation), 8 auxiliary binary variables, and 51 auxiliary variables. The complete model can be obtained by mail from the corresponding authors.

The The simulation results for the same input are shown in Fig. 5, (a) MLD Model (b) Error between MLD and PWA [4] (c) Error between MLD and PWA-This Work Fig. 5. Simulation results for a water plant model.

In this case, at t=168min, the simulation with the PWA system in the Figure 5.b is not valid because there exist no mode in the PWA representation that belongs to the stateinput vector reached in this point. Some other results can be found in [14].

V. CONCLUSIONS This work presents new algorithm for obtaining a suitable choice of the PWA description from a MLD representation. The results are applied to the three-tank benchmark problem, to a car with robotized gear shift and to a drinking water plant, the three examples have been reported in the literature as examples of hybrid dynamical systems modeled with MLD formalism. The simulation results show that the PWA models obtained have the same behavior with respect to the MLD models. However in some cases the obtained PWA model does not have a valid solution for some state-input sub-spaces. As a consequence of the enumeration procedure, our PWA models have more submodels/regions than the algorithm in [3], however we show that the procedure does not spent much more computation time because of the simplicity in its formulation, and it ensures the covering of all regions included in the original MLD model.

Ongoing work concerns the analysis of MLD Systems with some results from PWA systems.

INTRODUCTION

Switched nonlinear control systems are characterized by a set of several continuous nonlinear state dynamics with a logic-based controller, which determines simultaneously a sequence of switching times and a sequence of modes. As performance and efficiency are key issues in modern technological system such as automobiles, robots, chemical processes, power systems among others, the design of optimal logic-based controllers, covering all those functionalities while satisfying physical and operational constraints, plays a fundamental role. In the last years, several researchers have considered the optimal control of switched systems. An early work on the problem is presented in [1], where a class of hybrid-state continuous-time dynamic system is investigated. Later, a generalization of the optimal control problem and algorithms of hybrid systems is presented [2]. The particular case of the optimal control problem of switched systems is presented in [3] and [4]. However, most of the efforts have been typically focused on linear subsystems [5]. In general, the optimal control problem of switched system is often computationally hard as it encloses both elements of optimal control as well as combinatorial optimization [6]. In particular, necessary optimality conditions for hybrid systems have been derived using general versions of the Maximum Principle [7,8] and more recently in [9]. In the case of switching systems [4] and [6], the switched system has been embedded into a larger family of systems, and the optimization problem is formulated. For general hybrid systems, with nonlinear dynamics in each location and with autonomous and controlled switching, necessary optimality conditions have recently been presented in [10]; and using these conditions, algorithms based on the hybrid Maximum Principle have been derived. Focusing on real-time applications, an optimal control problem for switched dynamical systems is considered, where the objective is to minimize a cost functional defined on the state, and where the control variable consists of the switching times [11]. It is widely perceived that the best numerical methods available for hybrid optimal control problems involve mixed integer programming (MIP) [12,13]. Even though great progress has been made in recent years in improving these methods, the MIP is an NP-hard problem, so scalability is problematic. One solution for this problem is to use the traditional nonlinear programming techniques such as sequential quadratic programming, which reduces dramatically the computational complexity over existing approaches [6].

The main contribution of this paper is an alternative approach to solve effectively the optimal control problem for an autonomous nonlinear switched system based on the probability measures introduced in [14], and later used in [15] and [16] to establish existence conditions for an infinitedimensional linear program over a space of measure. Then, we apply the theory of moments, a method previously introduced for global optimization with polynomials in [17,18], and later extended to nonlinear 0 1 programs using an explicit equivalent positive semidefinite program in [19]. We also use some results recently introduced for optimal control problems with the control variable expressed as polynomials [20][21][22]. The first approach relating switched systems and polynomial representations can be found in [23]. The moment approach for global polynomial optimization based on semidefinite programming (SDP) is consistent, as it simplifies and/or has better convergence properties when solving convex problems. This approach works properly when the control variable (i.e., the switching signal) can be expressed as a polynomial. Essentially, this method transforms a nonlinear, nonconvex optimal control problem (i.e., the switched system) into an equivalent optimal control problem with linear and convex structure, which allows us to obtain an equivalent convex formulation more appropriate to be solved by high-performance numerical computing. In other words, we transform a given controllable switched nonlinear system into a controllable continuous system with a linear and convex structure in the control variable.

This paper is organized as follows. In Section 2, we present some definitions and preliminaries. A semidefinite relaxation using the moment approach is developed in Section 3. An algorithm is developed on the basis of the semidefinite approach in Section 4 with a numerical example to illustrate our approach, and finally in Section 5, some conclusions are drawn.

THE SWITCHED OPTIMAL CONTROL PROBLEM

Switched systems

The switched system adopted in this work has a general mathematical model described by

P x.t/ D f .t/ .x.t //, ( 1 )
where x.t/ is the state, f i W R n 7 ! R n is the i th vector field, x.t 0 / D x 0 are fixed initial values, and W OEt 0 , t f 7 ! Q 2 ¹0, 1, 2, ..., qº is a piecewise constant function of time, with t 0 and t f as the initial and final times, respectively. Every mode of operation corresponds to a specific subsystem P x.t/ D f i .x.t //, for some i 2 Q, and the switching signal determines which subsystem is followed at each point of time into the interval OEt 0 , t f . The control input is a measurable function.

In addition, we consider a non-Zeno behavior, that is, we exclude an infinite switching accumulation points in time. Finally, we assume that the state does not have jump discontinuities. Moreover, for the interval OEt 0 , t f , the control functions must be chosen so that the initial and final conditions are satisfied.

Definition 1

A control for the switched system in (1) is a duplet consisting of (i) a finite sequence of modes, and (ii) a finite sequence of switching times such that, t 0 < t 1 < < t q D t f .

Switched optimal control problem

Let us define the optimization functional in Bolza form to be minimized as

J D '.x.t f // C Z t f t 0 L .t/ .t , x.t//dt , ( 2 )
where '.x.t f // is a real-valued function, and the running switched costs

L .t/ W R C R n 7 ! R are continuously differentiable for each 2 Q.
A switched optimal control problem (SOCP) can be stated in a general form as follows.

Definition 2 Given the switched system in (1) and a Bolza cost functional J as in (2), the SOCP is given by min

.t/2Q J.t 0 , t f , x.t 0 /, x.t f /, x.t/, .t// (3)
subject to the state x. / satisfying Equation (1).

The SOCP can have the usual variations of fixed or free initial or terminal state, free terminal time, and so forth.

A Polynomial representation

The starting point is to rewrite (1) as a continuous non-switched control system as it has been shown in [24]. The polynomial expression in the control variable able to mimic the behavior of the switched system is developed using a variable v, which works as a control variable.

A polynomial expression in the new control variable v.t/ can be obtained through Lagrange polynomial interpolation and a constraint polynomial as follows. First, let the Lagrange polynomial interpolation quotients be defined as [25],

l k .v/ D q Y i D0 i ¤k .v i/ .k i/ . ( 4
)
The control variable is restricted by the set D ¹v 2 R jg.v/ D 0 º, where g.v/ is defined by

g.v/ D q Y kD0 .v k/. ( 5 )
General conditions for the subsystems functions should be satisfied.

Assumption 3

The nonlinear switched system satisfies growth, Lipschitz continuity, and coercivity qualifications concerning the mappings

f i W R n 7 ! R n L i W R n 7 ! R
to ensure existence of solutions of (1).

The solution of this system may be interpreted as an explicit ODE on the manifold . A related continuous polynomial system of the switched system (1) is constructed in the following proposition [24].

Proposition 4 Consider a switched system of the form given in (1). There exists a unique continuous state system with polynomial dependence in the control variable v, F.x, v/ of degree q in v, with v 2 as follows:

P x D F.x, v/ D q X kD0 f k .x/l k .v/. ( 6 )
Then, this polynomial system is an equivalent polynomial representation of the switched system (1).

Similarly, we define a polynomial equivalent representation for the running cost L .t/ by using the Lagrange's quotients as follows.

Proposition 5

Consider a switched running cost of the form given in (2). There exists a unique polynomial running cost equation L.x, v/ of degree q in v, with v 2 as follows:

L.x, v/ D q X kD0 L k .x/l k .v/ (7)
with l k .v/ defined in (4). Then, this polynomial system is an equivalent polynomial representation of the switched running cost in (2).

The equivalent optimal control problem (EOCP), which is based on the equivalent polynomial representation is described next.

The functional using Equation ( 7) is defined by

J D '.x.t f // C Z t f t 0 L.x, v/dt , ( 8 )
subject to the system defined in (6), with x 2 R n , v 2 , and x.t 0 / D x 0 , where l k .v/, , and L are defined earlier. Note that this control problem is a continuous polynomial system with the input constrained by a polynomial g.v/. This polynomial constraint is nonconvex with a disjoint feasible set, and traditional optimization solvers perform poorly on such equations, as the necessary constraint qualification is violated. This makes this problem intractable directly by traditional nonlinear optimization solvers. Next, we propose a convexification of the EOCP using the special structure of the control variable v, which improves the optimization process.

SEMIDEFINITE RELAXATION USING A MOMENTS APPROACH

Relaxation of the optimal control problem

We describe the relaxation of the polynomial optimal control problem, for which, regardless of convexity assumptions, existence of optimal solutions can be achieved. Classical relaxation results establish, under some technical assumptions, that the infimum of any functional does not change when we replace the integrand by its convexification. In the previous section, a continuous representation of the switched system has been presented. This representation has a polynomial form in the control variable, which implies that this system is nonlinear and nonconvex with a disjoint feasible set. Thus, traditional optimization solvers have a disadvantaged performance, either by means of the direct methods (i.e., nonlinear programming) or indirect methods (i.e., Maximum Principle). We propose then, an alternative approach to deal with this problem. The main idea of this approach is to convexify the control variable in polynomial form by means of the method of moments. This method has been recently developed for optimization problems in polynomial form (see [17,18], among others). Therefore, a linear and convex relaxation of the polynomial problem ( 8) is presented next. The relaxed version of the problem is formulated in terms of probability measures associated with sequences of admissible controls [15].

Let be the set of admissible controls v.t/. The set of probability measures associated to the admissible controls in is

ƒ D ® D ¹ t º t 2OEt 0 ,t f W supp. t / , a.e., t 2 OEt 0 , t f ¯,
where is a probability measure supported in . The functional J.x, v/ defined on ƒ is now given by

J.x, v/ D '.x.t f // C Z t f t 0 Z L.x.t /, v/d t .v/dt ,
where x.t/ is the solution of

P x.t/ D Z F.x, v/d t .v/, x.t 0 / D x 0 .
We have obtained a reformulation of the problem that is an infinite dimensional linear program and thus not tractable as it stands. However, the polynomial dependence in the control variable allows us to obtain a semidefinite program or linear matrix inequality relaxation, with finitely many constraints and variables. By means of moments variables, an equivalent convex formulation more appropriate to be solved by numerical computing can be rendered. The method of moments takes a proper formulation in probability measures of a nonconvex optimization problem ( [18,23], and references therein). Thus, when the problem can be stated in terms of polynomial expressions in the control variable, we can transform the measures into algebraic moments to obtain a new convex program defined in a new set of variables that represent the moments of every measure [17,18,22]. We define the space of moments as

D ² m D ¹m k º W m k D Z v k d .v/, 2 P . / ³ ,
where P . / is the convex set of all probability measures supported in . In addition, a sequence m D ¹m k º has a representing measure supported in only if these moments are restricted to be entries on positive semidefinite moments and localizing matrices [17,19]. For this particular case, when the control variable is of dimension one, the moment matrix is a Hankel matrix with m 0 D 1, that is, for a moment matrix of degree d , we have

M d .m/ D 2 6 6 6 4
m 0 m 1 m d m 1 m 2 m d C1 . . . . . . . . . m d m d C1 m 2d 3 7 7 7 5 .
The localizing matrix is defined on the basis of corresponding moment matrix, whose positivity is directly related to the existence of a representing measure with support in as follows. Consider the set defined by the polynomial ˇ.v/ D ˇ0 C ˇ1v C ˇd v Á . It can be represented in moment variables as ˇ.m/ D ˇ0 C ˇ1m 1 C ˇÁm Á , or in compact form as ˇ.m/ D P Á D0 ˇ m . Suppose that the entries of the corresponding moment matrix are m , with 2 OE0, 1, : : : , 2d . Thus, every entry of the localizing matrix is defined as l D P d D0 ˇ m C . Note that the localizing matrix has the same dimension of the moment matrix, that is, if d D 1 and the polynomial ˇD v C 2v 2 , then the moment and localizing matrices are

M 1 .m/ D Ä 1 m 1 m 1 m 2 , M 1 .ˇm/ D Ä m 1 C 2m 2 m 2 C 2m 3 m 2 C 2m 3 m 3 C 2m 4 .
More details on the method of moments can be found in [19,26].

Because J is a polynomial in v of degree q, the criterion R Ld involves only the moments of up to order q and is linear in the moment variables. Hence, we replace with the finite sequence m D ¹m k º of all its moments up to order q. We can then express the linear combination of the functional J and the space of moments as follows

min v2 J.x, v/ ! min 2P . / Z J.x, v/d .v/ D min m k 2 Z t f t 0 X i X k L i .x/˛i k m k , ( 9 )
where ˛ik are the coefficients resulting of the factorization of Equation ( 4). Similarly, we obtain the convexification of the state equation

P x.t/ D Z F.x, v/d .v/ D X i X k f i .x/˛i k m k . ( 10
)
We have now a problem in moment variables, which can be solved by efficient computational tools as it is shown in the next section.

Semidefinite programs for the EOCP

We can use the functional and the state equation with moment structure to rewrite the relaxed formulation as a SDP. First, we need to redefine the control set to be coherent with the definitions of localizing matrix and representation results. We treat the polynomial g.v/ as two opposite inequalities, that is, g 1 .v/ D g.v/ > 0 and g 2 .v/ D g.v/ > 0, and we redefine the compact set to be D ¹g i .v/ > 0, i D 1, 2º. Also, we define also a prefixed order of relaxation, which is directly related to the number of subsystems.

Let w be the degree of the polynomial g.v/, which is equivalent to the degree of the polynomials g 1 and g 2 . Considering its parity, we have that if w is even (odd) then r D w=2 (r D .w C 1/=2). In this case, r corresponds to the prefixed order of relaxation. We use a direct transcription method to obtain an SDP to be solved through a nonlinear programming (NLP) algorithm [27]. Using a discretization method, the first step is to split the time interval OEt 0 , t f into N subintervals as t 0 < t 1 < t 2 < : : : < t N D t f , with a time step h predefined by the user. The integral term in the functional is implicitly represented as an additional state variable, transforming the original problem in Bolza form into a problem in Mayer form, which is a standard transformation [27]. Therefore, we obtain a set of discrete equations in moment variables. In this particular case, we have used a trapezoidal discretization, but we could have used a more elaborated discretization scheme. Thus, the optimal control problem can be formulated as an SDP.

Consider a fixed t in the time interval OEt 0 , t f and let Assumption 3 holds. We can state the following SDP of relaxation order r (SDP r ).

Semidefinite program-SDP r : For every j D ¹1, 2, : : : , N º, a semidefinite program SDP r can be described by

J r D min m.t j / h 2 N 1 X j D0 L.x.t j /, m.t j // s.t. x.t j C1 / D x.t j / C h X i X k f i .x.t j //˛i k m k .t j /, x.t 0 / D x 0 , (11)
M r .m.t j // 0, M 0 .g 1 m.t j // 0, M 0 .g 2 m.t j // 0.

Notice that in this case, the localizing matrices are linear. Let us consider the two subsystems case, that is, we have g D v 2 v that leads to polynomials

g 1 D v 2 v and g 2 D v v 2 , thus w D deg g D 2. The localizing matrices are M 0 .g 1 m/ D m 2 m 1 , so M 0 .g 2 m/ D m 1 m 2 .
This happens because we are using the minimum order of relaxation, r D w=2 or r D .w C 1/=2 depending on its parity. It is also known that the optimum J r is not always an optimal solution. However, in this case, a suboptimal solution is obtained, which corresponds to a lower bound on the global optimum J of the original problem. If we are interested in searching for an optimal solution, we can use a higher order of relaxation, that is, r > w=2, but the number of moment variables will increase, which can make the problem numerically inefficient. However, in many cases, low order relaxations will provide the optimal value J as shown in the next section, where we use a criterion to test whether the SDP r relaxation achieves the optimal value J for a fixed time. Still, suboptimal solutions of the original problem are obtained in the iteration that can be used. In order to solve a traditional NLP, we use the characteristic form of the moment and localizing matrices. We know that the moment matrices, and so the localizing matrices, are symmetric positive definite, which implies that every principal subdeterminant is positive [21]. Then, we use the set of subdeterminants of each matrix as algebraic constraints.

Analysis of solutions

Once a solution has been obtained in a subinterval OEt j 1 , t j , we obtain a vector of moments m .t j / D OEm 1 .t j /, m 2 .t j /, : : : , m r .t j /. Then, we need to verify if we have attained an optimal solution. On the basis of a rank condition of the moment matrix [26], we can test if we have obtained a global optimum at a relaxation order r. Also, on the basis of the same rank condition, we can check whether the optimal solution is unique or if it is a convex combination of several minimizers. The next result is based on an important result presented in [26] and used in [19] for optimization of 0 1 problems.

Proposition 6

For a fixed time t j in the interval OEt 0 , t f , the SDP r (11) is solved with an optimal vector solution m .t j /, if

r D rank M r m .t j / D rank M 0 m .t j / , (12)
then the global optimum has been reached and the problem for the fixed time t j has r optimal solutions.

Note that the rank condition ( 12) is a sufficient condition, which implies that the global optimum could be reached at some relaxation of order r and still the rank M r > rank M 0 . It should be noted that for the particular case of minimum order of relaxation, the rank condition yields r D rank M r .m.t j // D rank M 0 .m.t j // D 1, because M 0 D 1. Then, the rank M 0 D 1, which implies that when r > 1, that is, several solutions arise. In this case, we obtain a suboptimal switching solution.

Using the previous result, we can state some relations between solutions that can be used to obtain the switching signal in every t j . First, we state the following result valid for the unique solution case.

Theorem 7

If Problem (11) is solved for a fixed t j 2 OEt 0 , t f and the rank condition in ( 12) is verified with r D rank M r .m .t j // D 1, then the vector of moments m .t j / has attained a unique optimal global solution; and therefore, the optimal switching signal of the switched problem (3) for the fixed time t j is obtained as

.t j / D m 1 .t j /, (13)
where m 1 .t j / is the first moment of the vector of moments m .t j /.

Proof Suppose the problem (11) has been solved for a fixed t j , and a solution has been obtained. Let m .t j / be the solution obtained and the rank condition (12) has been verified. From a result presented in [19], it follows that min

2P . / Z J.x, v/d .v/ D min m k 2 Z t f t 0 X i X k L i .x/˛i k m k ,
where m .t j / D m 1 , : : : , m r is the vector of moments of some measure m . But then, as m is supported on , it also follows that m .t j / is an optimal solution and because of rank M r m .t j / D 1, this solution is unique and it is the solution of the polynomial problem (8). Then, we know that every optimal solution v corresponds to

m .t j / D v .t j /, v

Remark 8

Switched linear systems case. When we have a switched linear system, that is, when each subsystem is defined by a linear system, results presented in Theorem ( 7) can be directly applied, because Assumption (3) is satisfied for linear systems because the Lipschitz condition is satisfied globally [28]. Also, we can notice that if the switched linear system has one and only one switching solution, it corresponds to the first moment solution of the SDP r program for all t 2 OEt 0 , t f , that is, m 1 .t j / D .t j /, for all t j 2 OEt 0 , t f . This can be verified by means of the rank condition (12), which should be r D 1, for all t 2 OEt 0 , t f . This result states a correspondence between the minimizer of the original switched problem and the minimizer of the SDP r , and it can be used to obtain a switching signal directly from the solution of the SDP r . However, it is not always the case. Sometimes, we obtain a non-optimal solution that arises when the rank condition is not satisfied, that is, r > 1. But, we still can use information from the solution to obtain a switching suboptimal solution. In [29], a sum up rounding strategy is presented to obtain a suboptimal switched solution from a relaxed solution in the case of mixedinteger optimal control. We use a similar idea but extended to the case when the relaxed solution is any integer instead of the binary case.

Consider the first moment m 1 . / W OEt 0 , t f 7 ! OE0, q, which is a relaxed solution of the NLP problem for t j when the rank condition is not satisfied. We can state a correspondence between the relaxed solution and a suboptimal switching solution, which is close to the relaxed solution in average and is given by

.t j / D 8 < : dm 1 .t j /e if Z t j t 0 m 1 . /d ıt j 1 X kD0 .t k / > 0.5ıt bm 1 .t j /c otherwise (14)
where d e and b c are the ceiling and floor functions, respectively.

A SWITCHED OPTIMIZATION ALGORITHM

The ideas presented earlier are summarized in the following algorithm, which is implemented in Section 4.2 on a simple numerical example presented as a benchmark in [30]. The core of the algorithm is the inter-relationship of three main ideas:

(i) The equivalent optimal control problem The EOCP is formulated as in Section 2, where the equivalent representation of the switched system and the running cost are used to obtain a polynomial continuous system.

(ii) The relaxation of the EOCP -the theory of moments

The EOCP is now transformed into an SDP of order of relaxation r, which can be solved numerically efficiently. We obtain an equivalent linear convex formulation in the control variable. (iii)The relationship between the solutions of the original switched problem and the SDP solutions The solutions of the SDP r for each t j 2 OEt 0 , t f are obtained; and through an extracting algorithm, the solutions of the original problem are obtained.

Algorithm SDP r -SOCP

The optimal control pseudo-code algorithm for the switched systems is shown in Algorithm 1.

In the next section, we present a numerical example to illustrate the results presented in this work.

Numerical example: Lotka-Volterra problem

We present an illustrative example of a switched nonlinear optimal control problem reformulated as a polynomial optimal control problem. Then, this reformulation allows us to apply the semidefinite relaxation based on the theory of moments. We illustrate an efficient computational treatment to study the optimal control problem of switched systems reformulated as a polynomial expression.

We deal with the Lotka-Volterra fishing problem. Basically, the idea is to find an optimal strategy on a fixed time horizon to bring the biomass of both predator as prey fish to a prescribed steadystate. The system has two operation modes and a switching signal as a control variable. The optimal integer control shows chattering behavior, which makes this problem a benchmark to test different types of algorithms ‡ .

The Lotka-Volterra model, also known as the predator-prey model, is a coupled nonlinear differential equations where the biomasses of two fish species are the differential states x 1 and x 2 , the binary control is the operation of a fishing fleet, and the objective is to penalize deviation from a steady-state. The optimal control problem is described as follows:

min u Z t f t 0 .x 1 1/ 2 C .x 2 1/ 2 dt s.t. P x 1 D x 1 x 1 x 2 0.4x 1 u P x 2 D x 2 C x 1 x 2 0.2x 2 u
x.0/ D .0.5, 0.7/ > , u.t / 2 ¹0, 1º, t 2 OE0, 12.

(

The problem can be represented by the approach described earlier. Consider a subsystem f 0 when the control variable takes value 0, and a subsystem f 1 when the control variable takes value 1. This leads to a two operation modes and a switching control variable . / W OE0, 12 7 ! ¹0, 1º. Thus, by means of the algorithm SDP r -EOCP, an SDP program can be stated. First, we define the order of relaxation as r D w=2 D 1; the constraint control set as

D ¹g i .v/ > 0, g 1 .v/ D v 2 v, g 2 D v v 2 º
; the moment matrix with r D 1, M 1 .m/; and the localizing matrices, M 0 .g 1 m/ and M 0 .g 2 m/. Using the set and the moment and localizing matrices, we set the problem in moment variables obtaining the positive semidefinite program .SDP r /. Solving the SDP r program for each t 2 OE0, 12, with a step time h, we obtain an optimal trajectory, and the moment sequence allows us to calculate the switching signal.

Figure 1 shows the trajectories, the relaxed moment solution, and the switching signal obtained for an order of relaxation r D 1. It can be appreciated that when the relaxed solution has a unique optimal solution, that is, when the rank condition is satisfied, the relaxed solution has an exact unique solution that is integer and corresponds to the switching signal, which shows the validity of Theorem 7. Also, it is shown that when the rank condition is not satisfied, the algorithm proposed gives a suitable solution, that in average is close to the relaxed solution. The algorithm has shown that even if there is no global optimal solution, a local suboptimal solution is found. Furthermore, for the intervals where there is no optimal solution, a suboptimal solution has been found using the relaxed solution. In comparison with traditional algorithms, where a global suboptimal solution based on a relaxation is found, the proposed algorithm is able to detect whether an optimal solution is found in a time interval, which implies that if the system is composed by convex functions, a global optimal solution is found. The computational efficiency is based on the semidefinite methods of solutions.

CONCLUSIONS AND FUTURE WORK

In this paper, we have developed a new method for solving the optimal control problem of switched nonlinear systems based on a polynomial approach. First, we transform the original problem into a polynomial system, which is able to mimic the switching behavior with a continuous polynomial representation. Next, we transform the polynomial problem into a relaxed convex problem using the method of moments. From a theoretical point of view, we have provided sufficient conditions for the existence of the minimizer by using particular features of the relaxed, convex formulation. Even in the absence of classical minimizers of the switched system, the solution of its relaxed formulation provides minimizers. We have introduced the moment approach as a computational useful tool to solve this problem, which has been illustrated by means of a classical example used in switched systems. As a future work, the algorithm can be extended to the case when an external control input and the switching signal should be obtained.

Introduction

The ). These methods exploit the communication capabilities of the agents to coordinate their decisions based on the information received from their neighbours. Fully decentralised methodologies have important advantages, among which we highlight the increase of the autonomy and resilience of the whole system since the dependence on a central authority is avoided.

In this paper, we propose a distributed resource allocation algorithm that does not require a central coordinator. An important characteristic of our method is the capability of handling lower bounds on the decision variables. This feature is crucial in a large number of practical applications, e.g. in [START_REF] Conrad | Resource economics[END_REF], Pantoja and Quijano (2012), and Lee et al. (2016), where it is required to capture the non-negativity of the resource allocated to each entity. We use a Lyapunov-based analysis in order to prove that the proposed algorithm asymptotically converges to the optimal solution under some mild assumptions related to the convexity of the cost function, and the connectivity of the graph that represents the communication topology. In order to illustrate our theoretical results, we perform some simulations and compare our method with other techniques reported in the literature. Finally, we present two engineering applications of the proposed algorithm. The first one seeks to improve the energy efficiency in large-scale air-conditioning systems. The second one is related to the distributed computation of the Euclidean projection onto a given set.

Our approach is based on a continuous time version of the centre-free algorithm presented in Xiao and Boyd (2006). The key difference is that the method in Xiao and Boyd (2006) does not allow the explicit inclusion of lower bounds on the decision variables, unless they are added by means of barrier functions (either logarithmic or exact; [START_REF] Cherukuri | Distributed generator coordination for initialization and anytime optimization in economic dispatch[END_REF]. The problem of using barrier functions is that they can adversely affect the convergence time (in the case of using exact barrier functions) and the accuracy of the solution (in the case of using classic logarithmic barrier functions), especially for large-scale problems [START_REF] Jensen | Operations research models and methods[END_REF]. There are other methods that consider lower bound constraints in the problem formulation. For instance, Dominguez-Garcia, Cady, and Hadjicostis ( 2012 2009) use consensus steps to refine an estimation of the system state, while in our approach, consensus is used to equalise a quantity that depends on both the marginal cost perceived by each agent in the network and the Karush-Kuhn-Tucker (KKT) multiplier related to the corresponding resource's lower bound. In this regard, it is worth noting that the method studied in this paper requires less computational capability than the methods mentioned above. Finally, there are other techniques based on game theory and mechanism design (Kakhbod & Teneketzis, 2012; Sharma & Teneketzis, 2009) that decompose and solve resource allocation problems. Nonetheless, those techniques need that each agent broadcasts a variable to all the other agents, i.e. a communication topology given by a complete graph is required. In contrast, the method developed in this paper only uses a communication topology given by a connected graph, which generally requires lower infrastructure.

The remainder of this paper is organised as follows. Section 2 shows preliminary concepts related to graph theory. In Section 3, the resource allocation problem is stated. Then, in Section 4, we present our distributed algorithm and the main results on convergence and optimality. A comparison with other techniques reported in the literature is performed in Section 5. In Section 6, we describe two applications of the proposed method: (i) the optimal chiller loading problem in large-scale airconditioning systems, and (ii) the distributed computation of Euclidean projections. Finally, in Sections 7 and 8, arguments and conclusions of the developed work are presented.

Preliminaries

First, we describe the notation used throughout the paper and presents some preliminary results on graph theory that are used in the proofs of our main contributions.

In the multi-agent framework considered in this article, we use a graph to model the communication network that allows the agents to coordinate their decisions. A graph is mathematically represented by the pair G = (V, E ), where V = {1, . . . , n} is the set of nodes, and E ⊆ V × V is the set of edges connecting the nodes. G is also characterised by its adjacency matrix A = [a i j ]. The adjacency matrix A is an n × n non-negative matrix that satisfies: a ij = 1 if and only if (i, j) ∈ E, and a ij = 0 if and only if (i, j) / ∈ E. Each node of the graph corresponds to an agent of the multi-agent system, and the edges represent the available communication channels (i.e. (i, j) ∈ E if and only if agents i and j can share information). We assume that there is no edges connecting a node with itself, i.e. a ii = 0, for all i ∈ V; and that the communication channels are bidirectional, i.e. a ij = a ji . The last assumption implies that G is undirected. Additionally, we denote by N i = { j ∈ V : (i, j) ∈ E}, the set of neighbours of node i, i.e. the set of nodes that are able to receive/send information from/to node i.

Let us define the n × n matrix L(G) = [l i j ], known as the graph Laplacian of G, as follows:

l i j = ⎧ ⎨ ⎩ j∈V a i j if i = j -a i j if i ̸ = j. (1)
Properties of L(G) are related to connectivity characteristics of G as shown in the following theorem. We remark that a graph G is said to be connected if there exists a path connecting any pair of nodes. From Equation ( 1), it can be verified that L(G)1 = 0, where 1 = [1, . . . , 1] ⊤ , 0 = [0, . . . , 0] ⊤ . A consequence of this fact is that L(G) is a singular matrix. However, we can modify L(G) to obtain a nonsingular matrix as shown in the following lemma.

Lemma 2.1:

Let L k r (G) ∈ R (n-1
)×n be the submatrix obtained by removing the kth row of the graph Laplacian L(G), and let L k (G) ∈ R (n-1)×(n-1) be the submatrix obtained by removing the kth column of

L k r (G). If G is connected, then L k (G) is positive definite. Furthermore, the inverse matrix of L k (G) satisfies (L k (G)) -1 l k r k = -1, where l k r
k is the kth column of the matrix L k r (G). Proof: First, notice that L(G) is a symmetric matrix because G is an undirected graph. Moreover, notice that according to Equation (1), L(G) is diagonally dominant with non-negative diagonal entries. The same holds for L k (G) since this is a sub-matrix obtained by removing the kth row and column of L(G). Thus, to show that L k (G) is positive definite, it is sufficient to prove that L k (G) is nonsingular.

According to Theorem 2.1, since G is connected, L(G) has exactly n -1 linearly independent columns (resp. rows). Let us show that the kth column (resp. row) of L(G) can be obtained by a linear combination of the other columns (resp. rows), i.e. the kth column (resp. row) is not linearly independent of the rest of the columns (resp. rows).

Since L(G)1 = 0, notice that l ik = -j∈V, j̸ =k l i j , for all i ∈ V, i.e. the kth column can be obtained by a linear combination of the rest of the columns. Furthermore, since L(G) is a symmetric matrix, the same occurs with the kth row. Therefore, the submatrix L k (G) is nonsingular since its n -1 columns (resp. rows) are linearly independent. Now, let us prove that (L k (G)) -1 l k r k = -1. To do so, we use the fact that (L k (G)) -1 L k (G) = I, where I is the identity matrix. Hence, by the definition of matrix multiplication, we have that

n-1 m=1 lk im l k m j = 1 if i = j 0 if i ̸ = j , (2)
where l k i j and lk i j are the elements located in the ith row and jth column of the matrices L(G) and L k (G)

-1 , respectively. Thus,

n-1 m=1 lk im l k mi = 1, for all i = 1, . . . n -1. ( 3
)
Let l k r k m be the mth entry of the vector l k r k . Notice that, according to the definition of L k (G) and since L(G)1 = 0, l k mi = -n-1 j=1, j̸ =i l k m j -l k r k m . Replacing this value in Equation (3), we obtain

- n-1 j=1, j̸ =i n-1 m=1 lk im l k m j - n-1 m=1 lk im l kr km = 1, for all i = 1, . . . n -1.
According to Equation ( 2),

n-1 j=1, j̸ =i n-1 m=1 lk im l k m j = 0. This implies that n-1 m=1 lk im l k r k m = -1, for all i = 1, … , n -1. Therefore, L k (G) -1 l k r k = -1. Theorem 2.
1 and Lemma 2.1 will be used in the analysis of the method proposed in this paper.

Problem statement

In general terms, a resource allocation problem can be formulated as follows (Patriksson, 2008;Patriksson & Strömberg, 2015):

min x φ(x) := n i=1 φ i (x i ) (4a) subject to n i=1 x i = X (4b) x i ≥ x i , for all i = 1, . . . , n, (4c)
where x i ∈ R is the resource allocated to the ith zone;

x = [x 1 , … , x n ] ; φ i : R → R is a strictly convex and differentiable cost function; X is the available resource; and x i , is the lower bound of x i , i.e. the minimum amount of resource that has to be allocated in the ith zone.

Given the fact that we are interested in distributed algorithms to solve the problem stated in Equation ( 4), we consider a multi-agent network, where the ith agent is responsible for managing the resource allocated to the ith zone. Moreover, we assume that the agents have limited communication capabilities, so they can only share information with their neighbours. This constraint can be represented by a graph G = {V, E} as it was explained in Section 2.

Avoiding the individual inequality constraints (4c), KKT conditions establish that at the optimal solution x * = [x * 1 , . . . , x * n ] ⊤ of the problem given in Equation (4a-4b), the marginal costs φ ′ i (x i ) = dφ i dx i must be equal, i.e. φ ′ i (x * i ) = λ, for all i = 1, … , n, where λ ∈ R. Hence, a valid alternative to solve (4a-4b) is the use of consensus methods. For instance, we can adapt the algorithm presented in Xiao and Boyd (2006), which is described as follows:

ẋi = j∈N i φ ′ j (x j ) -φ ′ i (x i ) , for all i ∈ V. ( 5
)
This algorithm has two main properties: (i) at equilibrium,

φ ′ i (x * i ) = φ ′ j (x * j )
if the nodes i and j are connected by a path; (ii) n i=1 x * i = n i=1 x i (0), where x i (0) is the initial condition of x i . Therefore, if the graph G is connected and the initial condition is feasible (i.e. n i=1 x i (0) = X), x asymptotically reaches the optimal solution of (4a-4b) under (5). However, the same method cannot be applied to solve (4) (the problem that considers lower bounds in the resource allocated to each zone) since some feasibility issues related with the constraints (4c) arise.

In the following section, we propose a novel method that extends the algorithm in Equation ( 5) to deal with the individual inequality constraints given in Equation (4c).

Centre-free resource allocation algorithm

Resource allocation among a subset of nodes in a graph

First, we consider the following subproblem: let G = {V, E} be a graph comprised by a subset of active nodes V a and a subset of passive nodes V p , such that V a V p = V. A certain amount of resource X has to be split among those nodes to minimise the cost function φ(x) subject to each passive node is allocated with its corresponding lower bound x i . Mathematically, we formulate this subproblem as:

min x φ(x) (6a) subject to n i=1
x i = X (6b)

x i = x i , for all i ∈ V p . (6c)
Feasibility of ( 6) is guaranteed by making the following assumption.

Assumption 4.1: At least one node is active, i.e. V a ̸ = ∅.

According to KKT conditions, the active nodes have to equalise their marginal costs at the optimal solution. Therefore, a consensus among the active nodes is required to solve (6). Nonetheless, classic consensus algorithms, as the one given in Equation ( 5), cannot be used directly. For instance, if all the nodes of G apply ( 5) and G is connected, the marginal costs of both passive and active nodes are driven to be equal in steady state. This implies that the resource allocated to passive nodes can violate the constraint (6c). Besides, if the resource allocated to passive nodes is forced to satisfy (6c) by setting x * i = x i , for all i ∈ V p , there is no guarantee that the new solution satisfies (6b). Another alternative, is to apply (5) to only active nodes (in this case, the neighbourhood of node i ∈ V a in Equation ( 5) has to be taken as { j ∈ V a : (i, j) ∈ E}, and the initial condition must satisfy i∈V a x i (0) = X -i∈V p x i ). However, the sub-graph formed by the active nodes is not necessarily connected although G is connected. Hence, marginal cost of active nodes are not necessarily equalised at equilibrium, which implies that the obtained solution is sub-optimal. In conclusion, modification of (5) to address ( 6) is not trivial. In order to deal with this problem, we propose the following algorithm:

ẋi = j∈N i (y j -y i ), for all i ∈ V (7a) ẋi = (x i -x i ) + j∈N i (y j -y i ), for all i ∈ V p (7b) y i = φ ′ i (x i ) if i ∈ V a φ ′ i (x i ) + xi if i ∈ V p . (7c)
In the same way as in (5), the variables {x i , i ∈ V} in Equation ( 7) correspond to the resource allocated to both active and passive nodes. Notice that we have added auxiliary variables { xi , i ∈ V p } that allow the passive nodes to interact with their neighbours taking into account the constraint (6c). On the other hand, the term j∈N i (y jy i ), in Equations (7a)-(7b), leads to a consensus among the elements of the vector y = [y 1 , … , y n ] , which are given in Equation (7c). For active nodes, y i only depends on the marginal cost φ ′ i (x i ), while for passive nodes, y i depends on both the marginal cost and the state of the auxiliary variable xi . Therefore, if the ith node is passive, it has to compute both variables x i and xi . Furthermore, it can be seen that, if all the nodes are active, i.e. (V a = V), then the proposed algorithm becomes the one stated in Equation (5).

Notice that the ith node only needs to know y i and the values {y j : j ∈ N i } to compute j∈N i (y j -y i ) in (7a)-(7b). In other words, L(G)y = [ j∈N 1 (y j -y 1 ), . . . , j∈N n (y j -y n )] ⊤ is a distributed map over the graph G [START_REF] Cortés | Distributed algorithms for reaching consensus on general functions[END_REF]. This implies that the dynamics given in Equation ( 7) can be computed by each node using only local information. In fact, the message that the ith node must send to its neighbours is solely composed by the variable y i .

.. Feasibility

Let us prove that, under the multi-agent system proposed in Equation ( 7), x(t) satisfies the first constraint of the problem given by Equation ( 6), for all t ࣙ 0, provided that n i=1 x i (0) = X. Lemma 4.1: The quantity n i=1 x i (t ) is invariant under Equation (7), i.e. if n i=1 x i (0) = X, then n i=1 x i (t ) = X, for all t ࣙ 0.

Proof: It is sufficient to prove that ˙ = 0, where = n i=1 x i . Notice that ˙ = n i=1 ẋi = 1 ⊤ ẋ, where ẋ = [ ẋ1 , . . . , ẋn ] ⊤ . Moreover, according to Equation ( 7),

1 ⊤ ẋ = -1 ⊤ L(G)y. Since G is undirected, 1 ⊤ L(G) = L(G)1 = 0. Therefore, ˙ = 0.
The above lemma does not guarantee that x(t) is always feasible because of the second constraint in Equation ( 6), i.e. x i = x i , for all i ∈ V p . However, it is possible to prove that, at equilibrium, this constraint is properly satisfied.

.. Equilibrium point

The next proposition characterises the equilibrium point of the multi-agent system given in Equation (7).

Proposition 4.1:

If G is connected, the system in Equation (7) has an equilibrium point x * , { x * i , i ∈ V p }, such that:

φ ′ i (x * i ) = λ, for all i ∈ V a , where λ ∈ R is a constant; and x * i = x i , for all i ∈ V p . Moreover, x * i = λ -φ ′ i (x * i ), for all i ∈ V p . Proof: Let x * , { x * i , i ∈ V p }
be the equilibrium point of Equation (7). Since G is connected by assumption, it follows from Equation (7a) that y * i = λ, for all i ∈ V, where λ is a constant. Thus,

y * i = φ ′ i (x * i ) if i ∈ V a , and
y * i = φ ′ i (x * i ) + x * i , if i ∈ V p . Hence, φ ′ i (x * i ) = λ, for all i ∈ V a , and x * i = λ -φ ′ i (x * i )
, for all i ∈ V p . Moreover, given the fact that j∈N i (y * j -y * i ) = 0, it follows from Equation (7b) that x * i = x i , for all i ∈ V p Remark 4.1: Proposition 4.1 states that, at the equilibrium point of ( 7), the active nodes equalise their marginal costs, while each passive node is allocated with an amount of resource equal to its corresponding lower bound. In conclusion, if n i=1 x * i = X, then it follows from Proposition 4.1, that x * minimises the optimisation problem given in Equation (6). Additionally, notice that the values { x * i , i ∈ V p } are equal to the KKT multipliers associated with the constraint (6c).

.. Convergence

Let us prove that the dynamics in Equation ( 7) converge to x * , { x * i , i ∈ V p }, provided that each φ i (x i ) is strictly convex. Proposition 4.2: Assume that φ i (x i ) is a strictly convex cost function, for all i ∈ V. If G is connected, n i=1 x i (0) = X, and Assumption 4.1 holds, then x(t) converges to x * under Equation (7), where x * is the solution of the optimisation problem stated in Equation (6), i.e. x * is the same given in Proposition 4.1. Furthermore, xi converges to x * i , for all i ∈ V p .

Proof: According to Lemma 4.1, since n i=1 x i (0) = X, then x(t) satisfies the first constraint of the problem stated in Equation ( 6), for all t ࣙ 0. Therefore, it is sufficient to prove that the equilibrium point x * , { x * i , i ∈ V p } (which is given in Proposition 4.1) of the system proposed in Equation ( 7) is asymptotically stable (AS). In order to do that, let us express our multi-agent system in error coordinates, as follows:

ė = -L(G)e y ėi = e i -L(G)e y i , for all i ∈ V p e y i = φ ′ i (x i ) -φ ′ i (x * i ) if i ∈ V a φ ′ i (x i ) -φ ′ i (x * i ) + êi if i ∈ V p , (8)
where L(G) is the graph Laplacian of G; e i = x i -x * i , and e y i = y i -y * i , for all i ∈ V; êi = xix * i , for all i ∈ V p ; e = [e 1 , … , e n ] ; e y = [e y 1 , . . . , e y n ] ⊤ ; and L(G)e y i represents the ith element of the vector L(G)e y .

Since Assumption 4.1 holds, V a ̸ = ∅. Let k be an active node, i.e. k ∈ V a , and let e k , e k y be the vectors obtained by removing the kth element from vectors e and e y , respectively. We notice that, according to Lemma 4.1, e k (t) = i ࢠ ν, i ࣔ k e i (t), for all t ࣙ 0. Therefore, Equation ( 8) can be expressed as

ėk = -L k (G)e k y -l k r k e y k e k = -i∈ν,i̸ =k e i ėi = e i -L k (G)e k y + l k r k e y k i , for all i ∈ V p e y i = φ ′ i (x i ) -φ ′ i (x * i ) if i ∈ V a φ ′ i (x i ) -φ ′ i (x * i ) + êi if i ∈ V p , (9)
where L k (G) and l k r k are defined in Lemma 2.1. In order to prove that the origin of the above system is AS, let us define the following Lyapunov function (adapted from Obando, Quijano, & Rakoto-Ravalontsalama, 2014):

V = 1 2 e k ⊤ L k (G) -1 e k + 1 2 i∈V p (e i -êi ) 2 . ( 10
)
The function V is positive definite since G is connected (the reason of this fact is that, according to Lemma 2.1, L k (G) and its inverse are positive definite matrices if G is connected). The derivative of V along the trajectories of the system stated in Equation ( 9) is given by,

V = -e k ⊤ e k y -e k ⊤ (L(G)) -1 l k r k e y k - i∈V p e i (e i -êi ) Taking into account that L k (G) -1 l k r k = -1 (cf. Lemma 2.1), we obtain V = -e k ⊤ e k y + e y k i∈V,i̸ =k e i - i∈V p e i (e i -êi ) = - n i=1 e i (φ ′ i (x i ) -φ ′ i (x * i )) - i∈V p e i êi + i∈V p e i ( êi -e i ) = - n i=1 (x i -x * i )(φ ′ i (x i ) -φ ′ i (x * i )) - i∈V p e 2 i ,
where φ ′ i is strictly increasing given the fact that φ i is strictly convex, for all i ∈ V. Therefore, (x i -

x * i )(φ ′ i (x i ) -φ ′ i (x * i )) ≥ 0, for all i ∈ V, and thus V ≤ 0. Since V does not depend on { êi , i ∈ V p }, it is negative semidefinite. Let S = {{e i , i ∈ V}, { êi , i ∈ V p } : V = 0}, i.e. S = {e i , i ∈ V}, { êi , i ∈ V p } : e i = 0, for all i ∈ V .
Given the fact that G is connected and V ̸ = V p (by Assumption 4.1), then ė = 0 iff e y = 0 (see Equation ( 8)). Therefore, the only solution that stays identically in S is the trivial solution, i.e. e i (t) = 0, for all i ∈ V, êi (t ) = 0, for all i ∈ V p . Hence, we can conclude that the origin is AS by applying the Lasalle's invariance principle.

In summary, we have shown that the algorithm described in Equation ( 7) asymptotically solves the subproblem in Equation ( 6), i.e. (7) guarantees that the resource allocated to each passive node is equal to its corresponding lower bound, while the remaining resource X -i∈V p x i is optimally allocated to active nodes.

Optimal resource allocation with lower bounds

Now, let us consider our original problem stated in Equation (4), i.e. the resource allocation problem that includes lower bound constraints. Let x * = [x * 1 , . . . , x * n ] ⊤ be the optimal solution of this problem. Notice that, if we know in advance which nodes will satisfy the constraint (4c) with strict equality after making the optimal resource allocation process, i.e. I := {i ∈ V : x * i = x i }, we can mark these nodes as passive and reformulate (4) as a subproblem of the form (6). Based on this idea, we propose a solution method for (4), which is divided in two stages: in the first one, the nodes that belong to I are identified and marked as passive; in the second one, the resulting subproblem of the form ( 6) is solved by using (7).

Protocol (7) can be also used in the first stage of the method as follows: in order to identify the nodes that will satisfy (4c) with strict equality at the optimal allocation, we start marking all nodes as active and apply the resource allocation process given by (7). The nodes that are allocated with an amount of resource below their lower bounds at equilibrium are marked as passive, and then ( 7) is newly applied (in this way, passive nodes are forced to meet (4c)). This iterative process is performed until all nodes satisfy their lower bound constraints. Notice that the last iteration of this procedure corresponds to solve a subproblem of the form (6) where the set of passive nodes is equal to the set I. Therefore, this last iteration is equivalent to the second stage of the proposed method.

Summarising, our method relies on an iterative process that uses the continuous-time protocol (7) as a subroutine. The main idea of this methodology is to identify in each step the nodes that have an allocated resource out of their lower bounds. These nodes are marked as passive, so they are forced to satisfy their constraints in subsequent iterations, while active nodes seek to equalise their marginal costs using the remaining resource. In the worst case scenario, the classification between active and passive nodes requires |V| iterations, where |V| is the number of nodes in the network. This fact arises when only one active node becomes passive at each iteration.

The proposed method is formally described in Algorithm 1. Notice that this algorithm is fully decentralised since Steps 4-6 can be computed by each agent using only local information. Step 4 corresponds to solve Equation (7), while Steps 5 and 6 describe the conditions for converting an active node into passive. Let us note that Steps 4-6 have to be performed |V| times since we are considering the worst case scenario. Therefore, each agent needs to know the total number of nodes in the network. This requirement can be computed in a distributed way by using the method proposed in Garin and Schenato (2010, p. 90). We also notice the fact that the agents have to be synchronised (as usual in several distributed algorithms; [START_REF] Cortés | Distributed algorithms for reaching consensus on general functions[END_REF][START_REF] Garin | A survey on distributed estimation and control applications using linear consensus algorithms[END_REF]Xiao & Boyd, 2006) in order to apply the Step 4 of Algorithm 1, i.e. all agents must start solving Equation ( 7) at the same time.

Algorithm 1: Resource allocation with lower bounds

Input:

-Parameters of the problem in Equation ( 4).

-An initial value x (0) , such that n i=1 x (0) i =X. Output: Optimal allocation x * 1 Mark all nodes as active, i.e. Ṽa,0 ← V, Ṽp,0 ← ∅.;

2 xi,0 ← x (0)
i , for all i ∈ V.; 3 for l ← 1 to |V| do 4 xi,l ← x i (t l ), for all i ∈ V, where x i (t l ) is the solution of Equation (7a) at time t l , with initial conditions

x(0) = [ x1,l-1 , . . . , xn,l-1 ] ⊤ , V a = Ṽa,l-1 , V p = Ṽp,l-1 , and { xi (0) = 0, ∀i ∈ V p }.; 5 Ṽp,l ← Ṽp,l-1 {i ∈ Ṽa,l-1 : xi < x i },
and Ṽa,l ← Ṽa,l-1 \{i ∈ Ṽa,l-1 : xi < x i }.;

6 x * ← [ x1,l , . . . , xn,l ] ⊤ .; 7 return x * ;
According to the reasoning described at the beginning of this subsection, we ideally require to know the steady-state solution of Equation ( 7) at each iteration of Algorithm 1 (since we need to identify which nodes are allocated with an amount of resource below their lower bounds in steady state). This implies that the time t l in Step 4 of Algorithm 1 goes to infinity. Under this requirement, each iteration would demand infinite time and the algorithm would not be implementable. Hence, to relax the infinite time condition, we state the following assumption on the time t l . Assumption 4.2: Let x * i,l be the steady state of x i (t) under Equation (7), with initial conditions x(0) = xi,l-1 , V a = Ṽa,l-1 , V p = Ṽp,l-1 , and { xi (0) = 0, ∀i ∈ V p } 1 . For each l = 1, . . . , |V| -1, the time t l satisfies the following condition: x i (t l ) < x i if and only if x * i,l < x i , for all i ∈ V. According to assumption 4.2, for the first |V| -1 iterations, we only need a solution of (7) that is close enough to the steady-state solution. We point out the fact that, if the conditions of Proposition 4.2 are met in the lth iteration of Algorithm 1, then x i (t) asymptotically converges to x * i,l , for all i ∈ V, under Equation (7). Therefore, Assumption 4.2 is satisfied for large values of t 1 , . . . , t |V|-1 .

Taking into account all the previous considerations, the next theorem states our main result regarding the optimality of the output of Algorithm 1.

Theorem 4.1:

Assume that G is a connected graph. Moreover, assume that φ i is a strictly convex function for all i = 1, … , n. If t 1 , . . . , t |V|-1 satisfy Assumption 4.2, and the problem stated in Equation (4) is feasible, then the output of Algorithm 4 tends to the optimal solution of the problem given in Equation (4) as t |V| → ∞.

Proof:

The ith component of the output of Algorithm 1 is equal to xi,|V| = x i (t |V| ), where x i (t |V| ) is the solution of Equation (7a) at time t |V| , with initial conditions [ x1,|V|-1 , . . . , xn,|V|-1 ] ⊤ , V a = Ṽa,|V| , and V p = Ṽp,|V| . Hence, it is sufficient to prove that {x * 1,|V| , . . . , x * n,|V| } solves the problem in Equation ( 4). In order to do that, let us consider the following premises (proof of each premise is written in brackets).

P1: { x1,l , . . . , xn,l } satisfies (4b), for all l = 1, . . . , |V| (this follows from Lemma 4.1, and form the fact that n i=1 xi,0 = X). P2: x * i,l = x i , for all i ∈ Ṽp,l-1 , and for all l = 1, . . . , |V| (this follows directly from Proposition 4.2).

P3: Ṽp,l = Ṽp,l-1 {i ∈ Ṽa,l-1 : x * i,l < x i }, and Ṽa,l = Ṽa,l-1 \{i ∈ Ṽa,l-1 : x * i,l < x i }, for all l = 1, . . . , |V| (this follows from Step 5 of Algorithm 1, and from Assumption 4.2).

P4: If for some l, Ṽp,l = Ṽp,l-1 , then Ṽp,l+ j = Ṽp,l-1 , for all j = 0, . . . , |V| -l (this can be seen from the fact that if the set of passive nodes does not change from one iteration to the next, the steady state of Equation (7a) is the same for both iterations).

P5: Ṽa,l Ṽp,l = V, for all l = 1, . . . , |V| (from P3, we know that Ṽa,l Ṽp,l = Ṽa,l-1 Ṽp,l-1 , for all l = 1, . . . , |V|. Moreover, given the fact that Ṽp,0 = ∅, and Ṽa,0 = V, (see step 1 of Algorithm 1) we can conclude P5). P6: Since the problem in Equation ( 4) is feasible by assumption, then | Ṽp,l | < |V|, for all l = 1, . . . , |V| (the fact that | Ṽp,l | ≤ |V|, for all l = 1, . . . , V, follows directly from P5. Let us prove that | Ṽp,l | ̸ = |V|, for all l = 1, . . . , V. We proceed by contradiction: Assume that there exists some l, such that | Ṽp,l-1 | < |V| and | Ṽp,l | = |V|. Hence, from P2 and P3, we know that x * i,l ≤ x i , for all i ∈ V; moreover, {i ∈ Ṽa,l-1 : x * i,l < x i } ̸ = ∅. Therefore, n i=1 x * i,l < n i=1 x i . According to P1, we know that n i=1 x * i,l = X; thus, X < n i=1 x i , which contradicts the feasibility assumption).

P7: {x * 1,|V| , . . . , x * n,|V| } satisfies the constraints (4c) (in order to prove P7, we proceed by contradiction: assume that {x * 1,|V| , . . . , x * n,|V| } does not satisfy the constraints (4c). Since P2 holds, this assumption implies that {i ∈ Ṽa,|V-1| :

x * i,|V| < x i } ̸ = ∅. Therefore, Ṽp,|V| ̸ = Ṽp,|V|-1 (see P3). Using P4, we can conclude that Ṽp,|V| ̸ = Ṽp,|V|-1 ̸ = • • • ̸ = Ṽp,0 = ∅, i.e. {i ∈ Ṽa,|V|-j : x * i,|V|-j+1 < x i } ̸ = ∅, for all j = 1, . . . , |V|. Thus, according to P3, | Ṽp,|V| | > | Ṽp,|V|-1 | > • • • > | Ṽp,1 | > 0. Hence, | Ṽp,|V| | ≥ |V|, which contradicts P6). P8: i∈ Ṽa,l x * i,l ≥ i∈ Ṽa,l x * i,l+1
(we prove P8 as follows: using P1 and the result in Lemma 4.1, we know that i∈V x * i,l = i∈V x * i,l+1 = X. Moreover, according to P5, V can be expressed as V = Ṽa,l Ṽp,l , where Ṽp,l-1 ⊂ Ṽp,l (see P3). Thus, we have that i∈ Ṽa,l

x * i,l + i∈ Ṽp,l ,i / ∈ Ṽp,l-1 x * i,l + i∈ Ṽp,l-1 x * i,l = i∈ Ṽa,l x * i,l+1 + i∈ Ṽp,l ,i / ∈ Ṽp,l-1 x * i,l+1 + i∈ Ṽp,l-1 x * i,l+1 . Furthermore, since P2 holds, we have that i∈ Ṽa,l x * i,l + i∈ Ṽp,l ,i / ∈ Ṽp,l-1 x * i,l + i∈ Ṽp,l-1 x i = i∈ Ṽa,l x * i,l+1 + i∈ Ṽp,l ,i / ∈ Ṽp,l-1 x i + i∈ Ṽp,l-1 x i . Therefore, i∈ Ṽa,l x * i,l = i∈ Ṽa,l x * i,l+1 + i∈ Ṽp,l ,i / ∈ Ṽp,l-1 (x i -x * i,l
), where x i -x * i,l > 0, for all i ∈ Ṽp,l , i / ∈ Ṽp,l-1 (according to P3). Hence, we can conclude P8).

P9: There exists k, such that k ∈ Ṽa,l , for all l = 1, . . . , |V| (in order to prove P9, we use the fact that, if k ∈ Ṽa,l , then k ∈ Ṽa,lj , for all j = 1, … , l (this follows from P3). Moreover, according to P5 and P6, | Ṽa,|V| | ̸ = 0; hence, there exists k, such that k ∈ Ṽa,|V| . Therefore, P9 holds). P9 guarantees that Assumption 4.1 is satisfied at each iteration.

P10:

φ ′ i (x * i,l ) ≥ φ ′ i (x * i,l+1
), for all i ∈ Ṽa,l (we prove P10 by contradiction: assume that φ

′ i (x * i,l ) < φ ′ i (x * i,l+1
), for some i ∈ Ṽa,l . According to Proposition 4.2, and since P1 and P9 hold, x * i,l has the characteristics given in Proposition 4.1, for all i ∈ V, and for all l = 1, . . . , |V|. Hence, φ ′ i (x * i,l ) has the same value for all i ∈ Ṽa,l-1 , and φ ′ i (x * i,l+1 ) has the same value for all i ∈ Ṽa,l . Moreover, since Ṽa,l ⊂ Ṽa,l-1 (according to P3), we have that

φ ′ i (x * i,l ) < φ ′ i (x * i,l+1
), for all i ∈ Ṽa,l . Thus, x * i,l < x * i,l+1 , for all i ∈ V a,l , because φ ′ i is strictly increasing (this follows from the fact that φ i is strictly convex by assumption). Therefore, i∈ Ṽa,l x * i,l < i∈ Ṽa,l x * i,l+1 , which contradicts P8). Now, let us prove that {x * 1,|V| , . . . , x * n,|V| } solves the Problem in Equation (4). First, the solution {x * 1,|V| , . . . , x * n,|V| } is feasible according to P1 and P7. On the other hand, from P9, it is known that ∃k : k ∈ Ṽa,l , for all l = 1, . . . ,

|V|. Let φ ′ k (x * k,|V| ) = λ, where λ ∈ R. Moreover, let us define V 0 = { j ∈ V : x * i,|V| > x i }, and V 1 = { j ∈ V : x * i,|V| = x i }. If i ∈ V 0 , then i ∈ Ṽa,|V|-1 (given the fact that, if i / ∈ Ṽa,|V-1| ⇒ i ∈ Ṽp,|V-1| ⇒ x * i,|V| = x i ⇒ i / ∈ V 0 ). Hence, φ ′ i (x * i,|V| ) = φ ′ k (x * k,|V| ) = λ (this follows from the fact that φ ′ j (x * j,l
) has the same value for all j ∈ Ṽa,l-1 , which in turn follows directly from step 4 of Algorithm 1, and Proposition 4.2).

If

i ∈ V 1 , then either i ∈ Ṽa,|V|-1 or i ∈ Ṽp,|V|-1 . In the first case, φ ′ i (x * i,|V| ) = φ ′ k (x * k,|V| ) = λ (following the rea- soning used when i ∈ V 0 ). In the second case, ∃l : i ∈ ( Ṽp,l \ Ṽp,l-1 ); hence, φ ′ i (x * i,l ) = φ ′ k (x * k,l
) (this follows from the fact that, if i ∈ ( Ṽp,l \ Ṽp,l-1 ), then i ∈ Ṽa,l-1 ). Furthermore, since i ∈ ( Ṽp,l \ Ṽp,l-1 ), x * i,l < x i (see P3), and given the fact that φ i is strictly increasing, we have that

φ ′ i (x * i,l ) < φ ′ i (x i ). Moreover, according to P10, φ ′ k (x * k,l ) ≥ φ ′ k (x * k,|V| ). Hence, φ ′ i (x i ) > φ ′ k (x * k,|V| ) = λ. In conclusion, if i ∈ V 1 , then φ ′ i (x * i,|V| ) ≥ λ. Thus, we can choose µ i ࣙ 0, for all i ∈ V, such that φ ′ i (x * i,|V| ) -µ i = λ, where µ i = 0 if i ∈ V 0 . Hence, let us note that ∂φ ∂x i | x i =x * i,|V| -µ i -λ = 0, for all i ∈ V, where ∂φ ∂x i | x i =x * i,|V| = φ ′ i (x * i,|V| ). Therefore, {x *
1,|V| , . . . , x * n,|V| , µ 1 , . . . , µ n , -λ} satisfies the KKT conditions for the problem given in Equation ( 4). Furthermore, since φ(x) is a strictly convex function by assumption, then {x * 1,|V| , . . . , x * n,|V| } is the optimal solution to that problem.

Early stopping criterion

Notice that, if the set of passive nodes does not change in the kth iteration of Algorithm 1 because all active nodes satisfy the lower bound constraints (see step 5), then the steady state solutions x * i,k and x * i,k+1 are the same, for all i ∈ V, which implies that the set of passive nodes also does not change in the (k + 1)th iteration. Following the same reasoning, we can conclude that

x * i,k = x * i,k+1 = • • • = x * i,|V|
, for all i ∈ V. Therefore, in this case, {x * 1,k , . . . , x * n,k } is the solution of our resource allocation problem. Practically speaking, this implies that Algorithm 1 does not need to perform more iterations after the kth one. Thus, it is possible to implement a flag z * i (in a distributed way) that alerts the agents if all active nodes satisfy the lower bound constraints after step 4 of Algorithm 1. A way to do that is by applying a min-consensus protocol [START_REF] Cortés | Distributed algorithms for reaching consensus on general functions[END_REF] with initial conditions z i (0) = 0 if the node i is active and does not satisfy its lower bound constraint, and z i (0) = 1 otherwise. Hence, notice that our flag z * i (i.e. the result of the min-consensus protocol) is equal to one, for all i ∈ V, only if all the active nodes satisfy the lower bound constraints, which corresponds to the early stopping criterion described above.

Simulation results and comparison

In this section, we compare the performance of our algorithm with other continuous-time distributed techniques found in the literature. We have selected three techniques that are capable to address nonlinear problems and can handle lower bound constraints: (i) a distributed interior point method (Xiao & Boyd, 2006), (ii) the local replicator equation (Pantoja & Quijano, 2012), and (iii) a distributed interior point method with exact barrier functions [START_REF] Cherukuri | Distributed generator coordination for initialization and anytime optimization in economic dispatch[END_REF]. The first one is a traditional methodology that uses barrier functions; the second one is a novel technique based on population dynamics; and the third one is a recently proposed method that follows the same ideas as the first one, but replaces classic logarithmic barrier functions by exact penalty functions. Below, we briefly describe the aforementioned algorithms.

Distributed interior point (DIP) method

This algorithm is a variation of the one presented in Equation (5) that includes strictly convex barrier functions to prevent the solution to flow outside the feasible region. The barrier functions b i (x i ) are added to the original cost function as follows:

φ b (x) = φ(x) + ϵ n i=1 b i (x i ) b i (x i ) = -ln x i -x i , for all i ∈ V,
where φ b (x) is the new cost function, and ϵ > 0 is a constant that minimises the effect of the barrier function when the solution is far from the boundary of the feasible set. With this modification, the distributed algorithm is described by the following equation:

ẋi = j∈N i φ ′ b j (x j ) -φ ′ b i (x i ) , for all i ∈ V, (11)
where

φ ′ b i (x i ) = dφ i dx i -ϵ db i dx i , i.e. φ ′ b i (x i
) is equal to the marginal cost plus a penalty term induced by the derivative of the corresponding barrier function.

Local replicator equation (LRE)

This methodology is based on the classical replicator dynamics from evolutionary game theory. In the LRE, the growth rate of a population that plays a certain strategy only depends on its own fitness function and on the fitness of its neighbours. Mathematically, the LRE is given by ẋi =

j∈N i (x i -x i )(x j -x j )(v i (x i ) -v j (x j )), v i = -φ ′ i (x i ), for all i ∈ V, (12)
where v i is the fitness perceived by the individuals that play the ith strategy. In this case, the strategies correspond to the nodes of the network, and the fitness functions to the negative marginal costs (the minus appears because replicator dynamics are used to maximise utilities instead of minimise costs). On the other hand, it can be shown that, if the initial condition x(0) is feasible for the problem given in Equation (4), then x(t) remains feasible for all t ࣙ 0, under the LRE.

Distributed interior point method with exact barrier functions (DIPe)

This technique follows the same reasoning of the DIP algorithm. The difference is that DIPe uses exact barrier functions [START_REF] Bertsekas | Necessary and sufficient conditions for a penalty method to be exact[END_REF] to guarantee satisfaction of the lower bound constraints. The exact barrier function for the ith node is given by:

b e i (x i ) = 1 ε [x i -x i ] + , where [•] + = max ( •, 0), 0 < ε < 1 2 max x∈F ∥∇φ(x)∥ ∞ , and F = {x ∈ R n : n i=1 x i = 1, x i ≥ x i }
is the feasible region of x for the problem (4). Using these exact barrier functions, the augmented cost function can be expressed as:

φ e b (x) = φ(x) + n i=1 b e i (x i ).
The DIPe algorithm is given in terms of the augmented cost function and its generalised gradient ∂φ e b (x) = [∂ 1 φ e b (x), . . . , ∂ n φ e b (x)] ⊤ as follows:

ẋi

∈ j∈N i ∂ j φ e b (x) -∂ i φ e b (x) , for all i ∈ V, (13)
where

∂ i φ e b (x) = ⎧ ⎨ ⎩ {φ ′ i (x i ) -1 ε } if x i < x i φ ′ i (x i ) -1 ε , φ ′ i (x i ) if x i = x i {φ ′ i (x i )} if x i > x i
In [START_REF] Cherukuri | Distributed generator coordination for initialization and anytime optimization in economic dispatch[END_REF], the authors show that the differential inclusion (13) converges to the optimal solution of the problem (4), provided that x(0) is feasible.

Comparison

In order to compare the performance of our algorithm with the three methods described above, we use the following simulation scenario: a set of n nodes connected as in Figure 1 (we use this topology to verify the behaviour of the different algorithms in the face of few communication channels since previous studies have shown that algorithms' performance decreases with the number of

(x) = n i=1 e a i (x i -b i ) + e -a i (x i -b i )
, where a i and b i are random numbers that belong to the intervals [1,2] and [- 1 2 , 1 2 ], respectively; a resource constraint X = 1; and a set of lower bounds {x i = 0 : i ∈ V}.

For each n, we generate 50 problems with the characteristics described above. The four distributed methods are implemented in Matlab employing the solver function ode23s. Moreover, we use the solution provided by a centralised technique as reference. The results on the average percentage decrease in the cost function reached with each algorithm and the average computation time (time taken by each algorithm for solving a problem 2 ) are summarised in Table 1. Results of DIPe for 100 and 200 nodes were not computed for practicality since the time required by this algorithm to solve a 100/200-nodes problem is very high.

We notice that the algorithm proposed in this paper always reaches the maximum reduction, regardless of the number of nodes that comprise the network. The same happens with the DIPe algorithm. This is an important advantage of our method compared to other techniques. In contrast, the algorithm based on the LRE performs far from the optimal solution. This unsatisfactory behaviour is due to the small number of links of the considered communication network. In Pantoja and Quijano (2012), the authors prove the optimality of the LRE in problems involving well connected networks; however, they also argue that this technique can converge to suboptimal solutions in other cases. On the other hand, the DIP method provides solutions close to the optimum. Nonetheless, its performance decreases when the number of nodes increases. This tendency is due to the influence of barrier functions on the original problem. Notice that, the larger the number of nodes, the bigger the effect of the barrier functions in Equation (11).

Regarding the computation time, although convergence of the proposed method is slower than the one shown by LRE and DIP, it is faster than the convergence of the method based on exact barrier functions, i.e. DIPe. Therefore, among the methods that guarantee optimality of the solution, our technique shows the best convergence speed. Computation time taken by DIPe is affected by the use of penalty terms that generate strong changes in the value of the cost function near to the boundaries of the feasible set. The drastic variations of the generalised gradient of exact barrier functions produces oscillations of numerical solvers around the lower bounds (a visual inspection of the results given in Figure 3 of [START_REF] Cherukuri | Distributed generator coordination for initialization and anytime optimization in economic dispatch[END_REF] confirms this claim). These oscillations are the main responsible for the low convergence speed shown by DIPe. On the other hand, LRE and DIP exhibit the fastest convergence. Hence, LRE and DIP are appealing to be implemented in applications that require fast computation and tolerate suboptimal solutions.

Applications

This section describes the use of the approach developed in this paper to solve two engineering problems. First, we present an application for sharing load in multiple chillers plants. Although this is not a large-scale application (multi-chiller plants are typically comprised of less than ten chillers; Yu & Chan, 2007), it aims to illustrate the essence of the proposed method and shows algorithm's performance in small-size problems. One of the reasons to use a distributed approach in small-/medium-size systems is due to the need of enhancing systems resilience in the face of central failures (e.g. in multiple chiller plants, central failures can occur due to cyber-attacks (Manic, Wijayasekara, Amarasinghe, & Rodriguez-Andina, 2016) against building management systems (Yu & Chan, 2007)). The second application deals with the distributed computation of the Euclidean projection of a vector onto a given set. Particularly, we use the proposed algorithm as part of a distributed technique that computes optimal control inputs for plants composed of a large number of sub-systems. This application aims to illustrate the performance of the method proposed in this paper when coping with large-scale problems.

Optimal chiller loading

The optimal chiller loading problem in multiple chiller systems arises in decoupled chilled-water plants, which are widely used in large air-conditioning systems [START_REF] Chang | Optimal chilled water temperature calculation of multiple chiller systems using Hopfield neural network for saving energy[END_REF]). The goal is to distribute the cooling load among the chillers that comprise the plant for minimising the total amount of power used by them. For a better understanding of the problem, below we present a brief description of the system. A decoupled chilled-water plant comprised by n chillers is depicted in Figure 2. The purpose of this plant is to provide a water flow f T at a certain temperature T s to the rest of the air-conditioning system. In order to do this task the plant needs to meet a cooling load C L that is given by the following expression:

C L = m f T (T r -T s ), (14)
where m > 0 is the specific heat of the water, and T r is the temperature of the water returning to the chillers. Since there are multiple chillers, the total cooling load C L is split among them, i.e. C L = n i=1 Q i , where Q i is the cooling power provided by the ith chiller, which, in turn, is given by

Q i = m f i (T r -T i ), (15)
where f i > 0 and T i are, respectively, the flow rate of chilled water and the water supply temperature of the ith chiller. As it is shown in Figure 2, we have that f T = n i=1 f i . In order to meet the corresponding cooling load, the ith chiller consumes a power P i that can be calculated using the following expression [START_REF] Chang | Optimal chilled water temperature calculation of multiple chiller systems using Hopfield neural network for saving energy[END_REF]:

P i = k 0,i + k 1,i m f i T r + k 2,i (m f i T r ) 2 + + k 3,i -k 1,i m f i -k 4,i m f i T r -2k 2,i (m f i ) 2 T r T i + k 5,i + k 6,i m f i + k 2,i (m f i ) 2 T 2 i , (16)
where k j, i , for j = 0, … , 6, are constants related to the ith chiller. If we assume that the flow rate f i of each chiller is constant, then P i is a quadratic function of the temperature T i . The optimal chiller loading problem involves the calculation of the chillers' water supply temperatures that meet the total cooling load given in Equation ( 14), and minimise the total amount of power consumed by the chillers, i.e. n i=1 P i . Moreover, given the fact that each chiller has a maximum cooling capacity, we have to consider the following additional constraints:

m f i (T r -T i ) ≤ Q i for all i = 1, . . . , n, (17)
where Q i is the maximum capacity (rated value) of the ith chiller.

Summarising, the optimal chiller loading problem can be expressed as follows: min

T 1 ,...,T n n i=1 P i (T i ) s.t. n i=1 m f i (T r -T i ) = C L T i ≥ T r -Q i m f i , for all i = 1, . . . , n. (18)
Now, let us consider that we want to solve the aforementioned problem in a distributed way by using a multiagent system, in which each chiller is managed by an agent that decides the value of the water supply temperature. We assume that the ith agent knows (e.g. by measurements) the temperature of the water returning to the chillers, i.e. T r , and the flow rate of chilled water, i.e. f i . Moreover, agents can share their own information with their neighbours through a communication network with a topology given by the graph G. If each P i (T i ) is a convex function, then the problem can be solved by using the method proposed in Algorithm 1 (we take, in this case, x i = f i T i ). The main advantage of this approach is to increase the resilience of the whole system in the face of possible failures, due to the fact that the plant operation does not rely on a single control centre but on multiple individual controllers without the need for a centralised coordinator.

.. Illustrative example

We simulate a chilled-water plant comprised by 7 chillers. 3 The cooling capacity and the water flow rate of each chiller are, respectively, Q i = 1406.8 kW, and f i = 65 kg.s -1 , for i = 1, … , 7; the specific heat of the water is m = 4.19 kW.s.kg -1 . degC -1 ; the supply temperature of the system is T s = 11 degC; and the coefficients k j, i of Equation ( 16) are given in Table 2. We operate the system at two different cooling loads, the first one is 90% of the total capacity, i.e. C L = 0.9 n i=1 Q i , and the second one is 60% of the total capacity, i.e. C L = 0.6 n i=1 Q i . The P i -T i curves are shown in Figure 3(a) for both cases, it can be noticed that all functions are convex. In order to apply 3(a)) are more loaded than the less efficient ones (i.e. chiller 2 and chiller 5). This can be noticed from the fact that their supply temperatures, in steady state, reach the minimum value. Furthermore, the energy consumption is minimised and power saving reaches to 2.6%. The results for the second cooling load, i.e. C L = 5908.6 kW, are shown in Figure 3(c), where it can be noticed a similar performance to that obtained with the first cooling load. However, in this case, it is not necessary that the supply temperatures reach the minimum value to meet the required load. Newly, energy consumption is minimised and power saving reaches to 2.8%. As it is stated in Section 4, convergence and optimality of the method is guaranteed under the conditions given in Theorem 4.1.

i =  i =  i =  i =  i =  i =  i =  k , i 
In both cases we use the early stopping criterion given in Section 4.

Although other techniques have been applied to solve the optimal chiller loading problem, e.g. the ones in Chang and Chen (2009), they require centralised information. In this regard, it is worth noting that the same objective is properly accomplished by using our approach, which is fully distributed.

Distributed computation of the Euclidean projection

Several applications require computing the Euclidean projection of a vector in a distributed way. These applications include matrix updates in quasi-Newton methods, balancing of traffic flows, and decomposition techniques for stochastic optimisation (Patriksson, 2008). The problem of finding the Euclidean projection of the vector ξ onto a given set X is formulated as follows:

min ξ ∥ ξ -ξ ∥ 2 2 s.t. ξ ∈ X , (19)
where ∥ • ∥ 2 is the Euclidean norm. The vector that minimises the above problem, which is denoted by ξ * , is the Euclidean projection. Roughly speaking, ξ * can be seen as the closest vector to ξ that belongs to the set X . In Barreiro-Gomez, Obando, Ocampo-Martinez, and Quijano (2015), the authors use a distributed computation of the Euclidean projection to decouple large-scale control problems. Specifically, they propose a discrete time method to address problems involving plants comprised of a large number of decoupled sub-systems whose control inputs are coupled by a constraint. The control inputs are associated with the power applied to the subsystems, and the constraint limits the total power used to control the whole plant. At each time iteration, local controllers that manage the sub-systems compute optimal control inputs ignoring the coupled constraint (each local controller uses a model predictive control scheme that does not use global information since the sub-systems' dynamics are decoupled). Once this is done, the coupled constraint is addressed by finding the Euclidean projection of the vector of local control inputs (i.e. the vector formed by all the control inputs computed by the local controllers) onto a domain that satisfies the constraint associated with the total power applied to the plant.

For a better explanation of the method, consider a plant comprised of n sub-systems. Let ûi (k) ≥ 0 be the control input computed by the ith local controller at the kth iteration ignoring the coupled constraint (non-negativity of ûi (k) is required since the control signals correspond to an applied power). Let û(k) = [ û1 (k), . . . , ûn (k)] ⊤ be the vector of local control inputs, and let u * (k) be the vector of control signals that are finally applied to the sub-systems. If the maximum allowed power to control the plant is U > 0, the power constraint that couples the control signals is given by n i=1 u * i (k) ≤ U . The vector u * (k) is calculated by using the Euclidean projection of û(k) onto a domain that satisfies the power constraint, i.e. u * (k) is the solution of the following optimisation problem (cf. Equation ( 19)):

min u(k) ∥ û(k) -u(k)∥ 2 2 (20a) s. t. n i=1 u i (k) ≤ U (20b) u i (k) ≥ 0, for all i = 1, . . . , n, (20c)
where u i (k) denotes the ith entry of the vector u(k).

Notice that u * (k) satisfies the power constraint and minimises the Euclidean distance with respect to the control vector ûk that is initially calculated by the local controllers. Computation of u * (k) can be performed by using the approach proposed in this paper because the problem stated in Equation ( 20) is in the standard form given in Equation ( 4) except for the inequality constraint (20b). However, this constraint can be addressed by adding a slack variable.

.. Illustrative example

Consider a plant composed of 100 sub-systems. Assume that, at the kth iteration of the discrete time method presented in Barreiro-Gomez, Obando, Ocampo-Martinez, and Quijano (2015), the control inputs that are initially computed by the local controllers are given by the entries of the vector û(k) = [ û1 (k), . . . , û100 (k)] ⊤ , where ûi (k) is a random number chosen from the interval [0, 1] kW. Furthermore, assume that the maximum allowed power to control the plant is U = 40 kW. To satisfy this constraint, the Euclidean projection described in Equation ( 20) is computed in a distributed way using Algorithm 1 with the early stopping criterion described in Section 4.

The results under a communication network with path topology (see Figure 1) are depicted in Figure 4. The curve at the top of Figure 4 describes the evolution of the Euclidean distance. Notice that the proposed algorithm minimises this distance and reaches the optimum value (dashed line), which has been calculated employing a centralised method. On the other hand, the curves at the bottom of Figure 4 illustrate the evolution of the values 100 i=1 u i (k) (solid line) and min {u i (k)} (dash-dotted line). These curves show that the constraints of the problem stated in Equation (20) performance even considering that the communication graph is sparse and the optimal solution is not in the interior of the feasible domain. As shown in Section 5, this characteristic is an advantage of Algorithm 1 over population dynamics techniques as the one proposed in Barreiro-Gomez, Obando, Ocampo-Martinez, and Quijano (2015) to compute the Euclidean projection in a distributed way.

Discussion

The method developed in this paper solves the problem of resource allocation with lower bounds given in Equation (4). The main advantage of the proposed technique is its distributed nature; indeed, our approach does not need the implementation of a centralised coordinator. This characteristic is appealing, especially in applications where communications are strongly limited. Moreover, fully distributed methodologies increase the autonomy and resilience of the system in the face of possible failures. In Section 5, we show by means of simulations that the performance of the method presented in this paper does not decrease when the number of nodes (which are related to the decision variables of the optimisation problem) is large, or the communication network that allows the nodes to share information has few channels.

In these cases, the behaviour of our approach is better than the behaviour of other techniques found in the literature, such as the DIP method, or the LRE. Moreover, it is worth noting that our technique addresses the constraints as hard. This fact has two important consequences: (i) in all cases, the solution satisfies the imposed constraints, and (ii) the objective function (and therefore the optimum) is not modified (contrary to the DIP method that includes the constraints in the objective function decreasing the quality of the solution as shown in Section 5.4). Other advantage of the method proposed in this paper is that it does not require an initial feasible solution of the resource allocation problem (4). Similarly to the DIPe technique, our method only requires that the starting point satisfies the resource constraint (4b), i.e. we need that n i=1 x i (0) = X. Notice that an initial solution x(0) that satisfies (4b) is not hard to obtain in a distributed manner. For instance, if we assume that only the kth node has the information of the available resource X, we can use (x k (0) = X, {x i (0) = 0 : i ∈ V, i ̸ = k}) as our starting point. Thus, an initialisation phase is not required. In contrast, other distributed methods, such as DIP and LRE needs an initial feasible solution of the problem (4), i.e. a solution that satisfies (4b) and (4c). Finding this starting point is not a trivial problem for systems involving a large number of variables. Therefore, for these methods, it is necessary to employ distributed constraint satisfaction algorithm (as the one described in Domınguez-Garcıa & Hadjicostis, 2011) as a first step.

On the other hand, we notice that to implement the early stopping criterion presented at the end of Section 4, it is required to perform an additional min-consensus step in each iteration. Despite this fact, if the number of nodes is large, this criterion saves computational time, because in most of the cases, all passive nodes are identified during the first iterations of Algorithm 1.

Conclusions

In this paper, we have developed a distributed method that solves a class of resource allocation problems with lower bound constraints. The proposed approach is based on a multi-agent system, where coordination among agents is done by using a consensus protocol. We have proved that convergence and optimality of the method is guaranteed under some mild assumptions, specifically, we require that the cost function is strictly convex and the graph related to the communication network that enables the agents to share information is connected. The main advantage of our technique is that it does not need a centralised coordinator, which makes the method appropriate to be applied in large-scale distributed systems, where the inclusion of centralised agents is undesirable or infeasible. As future work, we propose to use a switched approach in order to eliminate the iterations in Algorithm 1. Moreover, we plan to include upper bound constraints in our original formulation.

Notes

Disclosure statement

No potential conflict of interest was reported by the authors.

On Switchable Languages of Discrete-Event Systems with Weighted Automata

Michael Canu and Naly Rakoto-Ravalontsalama

Abstract-The notion of switchable languages has been defined by Kumar, Takai, Fabian and Ushio in [11]. It deals with switching supervisory control, where switching means switching between two specifications. In this paper, we first extend the notion of switchable languages to n languages, (n ≥ 3). Then we consider a discrete-event system modeled with weighted automata. The use of weighted automata is justified by the fact that it allows us to synthesize a switching supervisory controller based on the cost associated to each event, like the energy for example. Finally the proposed methodology is applied to a simple example.

Keywords: Supervisory control; switching control; weighted automata.

I. INTRODUCTION

Supervisory control initiated by Ramadge and Wonham [15] provides a systematic approach for the control of discrete event system (DES) plant. There has been a considerable work in the DES community since this seminal paper. On the other hand, from the domain of continuous-time system, hybrid and switched systems have received a growing interests [12]. The notion of switching is an important feature that has to be taken into account, not only in the continuous-time domain but for the DES area too.

As for non-blocking property, there exist different approaches. The first one is the non-blocking property defined in [15]. Since then other types of nonblocking properties have been defined. The mutually non-blocking property has been proposed in [5]. Other approaches of mutually and globally nonblocking supervision with application to switching control is proposed in [11]. Robust non-blocking supervisory control has been proposed in [1]. Other types of non-blocking include the generalised non-blocking property studied in [13]. Discrete-event modeling with switching maxplus systems is proposed in [17], an example of mode M. Canu is with Univ. los Andes, Bogota, Colombia, email: m.canu134@uniandes.edu.co N. Rakoto-Ravalontsalama is with IMT Atlantique and LS2N, France, e-mail: naly.rakoto@mines-nantes.fr switching DES is described in [6] and finally a modal supervisory control is considered in [7].

In this paper we will consider the notion of switching supervisory control defined by Kumar and Colleagues in [11] where switching means switching between a pair of specifications. Switching (supervisory) control is in fact an application of some results obtained in the same paper [11] about mutually non blocking properties of languages, mutually nonblocking supervisor existence, supremal controllable, relative-closed and mutually nonblocking languages. All these results led to the definition of a pair of switchable languages [11].

In this paper, we first extend the notion of switchable languages to n languages, (n ≥ 3). Then we consider a discrete-event system modeled with weighted automata. The switching supervisory control strategy is based on the cost associated to each event, and it allows us to synthesize an optimal supervisory controller. Finally the proposed methodology is applied to a simple example.

This paper is organized as follows. In Section II, we recall the notation and some preliminaries. Then in Section III the main results on the extension of n switchable languages (n ≥ 3) are given. An illustrative example of supervisory control of AGVs is proposed in Section IV, and finally a conclusion is given in Section V.

II. NOTATION AND PRELIMINARIES

Let the discrete event system plant be modeled by a finite state automaton [10], [4] to which a cost function is added. Definition 1: (Weighted automaton). A weighted automaton is defined as a sixtuple G = (Q, Σ, δ, q 0 , Q m , C) where • Q is the finite set of states, • Σ is the finite set of events, • δ : Q × Σ → Q is the partial transition function,

• q 0 ⊆ Q is the initial state, • Q m ⊆ Q is the set of marked states (final states),

• C : Σ → N is the cost function.

Let Σ * be the set of all finite strings of elements in Σ including the empty string ε. The transition function δ can be generalized to δ : Σ * × Q → Q in the following recursive manner: δ(ε, q) = q δ(ωσ, q) = δ(σ, δ(ω, q)) for ω ∈ Σ * The notation δ(σ, q)! for any σ ∈ Σ * and q ∈ Q denotes that δ(σ, q) is defined. Let L(G) ⊆ Σ * be the language generated by G, that is, L(G) = {σ ∈ Σ * |δ(σ, q 0 )!} Let K ⊆ Σ * be a language. The set of all prefixes of strings in K is denoted by pr(K) with pr(K) = {σ ∈ Σ * |∃ t ∈ Σ * ; σt ∈ K}. A language K is said to be prefix closed if K = pr(K). The event set Σ is decomposed into two subsets Σ c and Σ uc of controllable and uncontrollable events, respectively, where Σ c ∩ Σ uc = ∅. A controller, called a supervisor, controls the plant by dynamically disabling some of the controllable events.

A sequence σ 1 σ 2 . . . σ n ∈ Σ * is called a trace or a word in term of language. We call a valid trace a path from the initial state to a marked state (δ(ω, q 0 ) = q m where ω ∈ Σ * and q m ∈ Q m ). The cost is by definition non negative. In the same way, the cost function C is generalized to the domain Σ * as follows:

C(ε) = 0 C(ωσ) = C(ω) + C(σ) for ω ∈ Σ * In other words, the cost of a trace is the sum of the costs of each event that composes the trace. Definition 2: (Controllability) [15]. A language K ⊆ L(G) is said to be controllable with respect to (w.r.t.) L(G) and Σ uc if pr(K)Σ uc ∩ L(G) ⊆ pr(K). Definition 3: (Mutually non-blocking supervisor) [5]. a supervisor f : L(G) → 2 Σ-Σu is said to be (K 1 , K 2 )mutually non-blocking if

K i ∩ L m (G f ) ⊆ pr(K j ∩ L m (G f )), for i, j ∈ {1, 2}. (1)
In other words, a supervisor S is said to be mutually non-blocking w.r.t. two specifications K 1 and K 2 if whenever the closed-loop system has completed a task of one language (by completing a marked trace of that language), then it is always able to continue to complete a task of the other language [5].

Definition 4: (Mutually non-blocking language) [5]. A language H ⊆ K 1 ∪K 2 is said to be (K 1 , K 2 )-mutually non-blocking if H ∩K i ⊆ pr(H ∩K j ) for i, j ∈ {1, 2}.

The following theorem gives a necessary and sufficient condition for the existence of a supervisor.

Theorem 1: (Mutually nonblocking supervisor existence) [5]. Given a pair of specifications K 1 , K 2 ⊆ L m (G), there exists a globally and mutually nonblocking supervisor f such that L m (G f ) ⊆ K 1 ∪ K 2 if and only if there exists a nonempty, controllable, relative-closed, and (K 1 , K 2 )-mutually non-blocking sublanguage of

K 1 ∪ K 2 .
The largest possible language (the supremal element) that is controllable and mutually non-blocking exists, as stated by the following theorem.

Theorem 2: (SupMRC(K 1 ∪ K 2 ) existence) [5]. The set of controllable, relative-closed, and mutually nonblocking languages is closed under union, so that the supremal such sublanguage of K 1 ∪ K 2 , denoted supM RC(K 1 ∪ K 2 ) exists.

Recall that a pair of languages K 1 , K 2 are mutually nonconflicting if pr(K 1 ∩ K 2 ) = pr(K 1 ) ∩ pr(K 2 ) [18]. K 1 , K 2 are called mutually weakly nonconflicting if K i , pr(K j ) (i ̸ = j) are mutually nonconflicting [5].

Another useful result from [5] is the following. Given a pair of mutually weakly nonconflicting languages K 1 , K 2 ⊆ L m (G), the following holds ( [5], Lemma 3). If K 1 , K 2 are controllable then K 1 ∩ pr(K 2 ), K 2 ∩ pr(K 1 ) are also controllable.

The following theorem is proposed in [11] and it gives the formula for the supremal controllable, relativeclosed, and mutually nonblocking languages.

Theorem 3: (SupMRC(K 1 ∪ K 2 )) [11]. For relative-closed specifications

K 1 , K 2 ⊆ L m (G), supM RC(K 1 ∪ K 2 ) = supRC(K 1 ∩ K 2 ).
The following theorem, also from [11] gives another expression of the supremal controllable, relative-closed, and mutually nonblocking languages. Theorem 4: [11] Given a pair of controllable, relativeclosed, and mutually weakly nonconflicting languages K 1 , K 2 ⊆ L m (G), it holds that supM RC(K 1 ∪ K 2 ) = (K 1 ∩ K 2 ).

And finally the following theorem gives a third formula of the supremal controllable, relative-closed, and mutually nonblocking languages.

Theorem 5: [11] For specifications

K 1 , K 2 ⊆ L m (G), supM RC(K 1 ∪ K 2 ) = supM C(supRC(K 1 ∩ K 2 )).
In order to allow switching between specifications, a pair of supervisors is considered, such that the supervisor is switched when the specification is switched. The supervisor f i for the specification K i is designed to enforce a certain sublanguage H i ⊆ K i . Suppose a switching in specification from K i to K j is induced at a point when a trace s ∈ H i has been executed in the f i -controlled plant. Then in order to be able to continue with the new specification K j without reconfiguring the plant, the trace s must be a prefix of H j ⊆ K j . In other words, the two supervisors should enforce the languages H i and H j respectively such that H i ⊆ pr(H j ). Hence the set of pairs of such languages are defined to be switchable languages as follows.

Definition 5: (Pair of switchable languages) [11]. A pair of specifications K 1 , K 2 ⊆ L m (G) are said to be switchable languages if SW (K 1 , K 2 ) := {(H 1 , H 2 )|H i ⊆ K i ∩ pr(H j ), i ̸ = j, and H i controllable}.

The supremal pair of switchable languages exists and is given by the following theorem.

Theorem 6: (Supremal pair of switchable languages) [11]. For specifications

K 1 , K 2 ⊆ L m (G), supSW (K 1 , K 2 ) = (supM C(K 1 ∪ K 2 ) ∩ K 1 , supM C(K 1 ∪ K 2 ) ∩ K 2 ).

III. MAIN RESULTS

We now give the main results of this paper. First, we define a triplet of switchable languages. Second we derive a necessary and sufficient condition for the transitivity of switchable languages (n = 3). Third we generalize this definition to a n-uplet of switchable languages, with n > 3. And fourth we derive a necessary and sufficient condition for the transitivity of switchable languages for n > 3.

A. Triplet of Switchable Languages

We extend the notion of pair of switchable languages, defined in [11], to a triplet of switchable languages. Definition 6: (Triplet of switchable languages). A triplet of languages (K 1 , K 2 , K 3 ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, 2, 3} are said to be a triplet of switchable languages if they are pairwise switchable languages, that is, SW (K 1 , K 2 , K 3 ) := SW (K i , K j ), i ̸ = j, i, j = {1, 2, 3}.

Another expression of the triplet of switchable languages is given by the following lemma.

Lemma 1: (Triplet of switchable languages). A triplet of languages (K 1 , K 2 , K 3 ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, 2, 3} are said to be a triplet of switchable languages if the following holds:

SW (K 1 , K 2 , K 3 ) = {(H 1 , H 2 , H 3 ) | H i ⊆ K i ∩
pr(H j ), i ̸ = j, and H i controllable}.

B. Transitivity of Switchable Languages (n = 3)

The following theorem gives a necessary and sufficient condition for the transitivity of switchable languages.

Theorem 7: (Transitivity of switchable languages, n = 3) . Given 3 specifications (K 1 , K 2 , K 3 ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, 2, 3} such that SW (K 1 , K 2 ) and SW (K 2 , K 3 ). (K 1 , K 3 ) is a pair of switchable languages, i.e. SW (K 1 , K 3 ), if and only if 1) H 1 ∩ pr(H 3 ) = H 1 , and 2) H 3 ∩ pr(H 1 ) = H 3 .

Proof: The proof can be found in [3].

C. N-uplet of Switchable Languages

We now extend the notion of switchable languages, to a n-uplet of switchable languages, with (n > 3).

Definition 7: (N-uplet of switchable languages, n > 3). A n-uplet of languages (K 1 , ..., K n ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, ..., n}, n > 2, is said to be a n-uplet of switchable languages if the languages are pairwise switchable that is, SW (K 1 , ..., K n ) := SW (K i , K j ), i ̸ = j, i, j = {1, ..., n}, n > 2.

As for the triplet of switchable languages, an alternative expression of the n-uplet of switchable languages is given by the following lemma.

Lemma 2: (N-uplet of switchable languages, n > 3).

A n-uplet of languages (K 1 , . . . , K n ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, ..., n}, n > 3 are said to be a n-uplet of switchable languages if the following holds:

SW (K 1 , ..., K n ) = {(H 1 , ..., H n ) | H i ⊆ K i ∩
pr(H j ), i ̸ = j, and H i controllable}.

D. Transitivity of Switchable Languages (n > 3)

We are now able to derive the following theorem that gives a necessary and sufficient condition for the transitivity of n switchable languages.

Theorem 8: (Transitivity of n switchable languages, n > 3) . Given n specifications (K 1 , ..., K n ), K i ⊆ L m (G) with H i ⊆ K i , i = {1, ..., n}. Moreover, assume that each language K i is at least switchable with another language K j , i ̸ = j.

A pair of languages (K k , K l ) is switchable i.e. SW (K k , K l ), if and only if 1) H k ∩ pr(H l ) = H k , and 2) H l ∩ pr(H k ) = H l .

Proof: The proof is similar to the proof of Theorem 6 and can be found in [3]. It is to be noted that the assumption that each of the n languages be at least switchable with another language is important, in order to derive the above result.

IV. EXAMPLE: SWITCHING SUPERVISORY CONTROL OF AGVS

The idea of switching supervisory control is now applied to a discrete-event system, modeled with weighted automata. We take as an illustrating example the supervisory control of a fleet of fleet automated guided vehicles (AGVs) that move in a given circuit area.

The example is taken from [9]. A circuit is partitioned into sections and intersections. Each time an AGV moves in a new intersection or a new section, then the automaton will move to a new state in the associated automaton. An example of an area with its associated basic automaton is depicted in Figure 1.

The area to be supervised is the square depicted in Figure 1 (left). The flow direction with the arrows are specified the four intersections {A, B, C, D} and the associated basic automaton are given in Figure 1 (right). The basic automaton is denoted G basic = (Q b , Σ b , δ b , ∅, ∅) where the initial state and the final state are not defined. The initial state is defined according to the physical position of the AGV and the final state is defined according to its mission, that is his position target. A state represents and intersection or a section. Each state corresponding to a section is named XY i where X is the beginning of the section, Y its end and i the number of the AGV. For each section, there are two transitions, the first transition C XY is an input which is controllable and represents the AGV moving on the section from X to Y . The second transition is an output transition U Y which is uncontrollable and represents the AGV arriving to the intersection Y . For example the basic automaton depicted in Figure 1 (right) can be interpreted as follows. If AGV i arrives at section A, then it has two possibilities, either to go to section B with the event C ABi , or the go section D with the event C ADi . If we choose to go to section B, then the next state is AB i . From this state, the uncontrollable event U AB is true so that the following state is B i . And from B i , the only possibility is to exit to Point F with the uncontrollable event exit i . Now consider for example that 2 AGVs are moving in the circuit of Figure 1 (left). Assume AGV 1 is in D and AGV 2 is in AB so that the state is in (D 1 , AB 2 ). AGV 1 is leaving the area when the event exit 1 is true so that the system will be in state (E 1 , AB 2 ). And since AGV 1 is out of the considered area, then the new state will be (E 1 , AB 2 ) = (∅ 1 , AB 2 ) = (AB 2 ) since AGV 1 is out of the area.

We give here below the synthesis algorithm for calculating the supervisor S c as it aws proposed by Girault et Colleagues in [9]. For more details on the synthesis algorithm, the reader is referred to the above paper.

Algorithm 1 -Synthesis algorithm of S C [9] Data: G w,1 , . . . G w,n Result:

Supervisor S C G w ← {G w,1 , . . . G w,n } G u ← {∅} forall G w,i ∈ G w do G u ← G u ∪ U γi (G w,i ) end S C ← S(G u,i ) G u ← G u \{G u,1 } while G u ̸ = ∅ do x ← get(G u ) S C ← S(S C ||x) G u ← G u \{x} end V. CONCLUSIONS
The notion of switchable languages has been defined by Kumar and Colleagues in [11]. It deals with switching supervisory control, where switching means switching between two specifications. In this paper, we have extended the notion of switchable languages to a triplet of languages (n = 3) and we gave a necessary and sufficient condition for the transitivity of two switchable languages. Then we generalized the notion of switchable languages of a n-uplet of languages, n > 3 and we gave also necessary and sufficient condition for the transitivity of two (out of n) switchable languages. Finally the proposed methodology is applied to a simple example for the supervisory control of a fleet of AGVs. Ongoing work deals with a) the calculation of the supremal of n-uplet of switchable languages, and b) the optimal switching supervisory control of DES exploiting the cost of the weighted automata for the synthesis strategy.

2 Figure 1 :

 21 Figure 1: Selection Phases after the HDR

  2001-2004: First and Second Year: Control and Industrial Eng. courses at DAP • 2006-2012: MSc. MLPS (Management of Logistic and Production Systems) • 2012-present: MSc. MOST (Management and Optimization of Supply Chains and Transport) C. Courses given abroad • May 2008: Univ. of Cagliari (Italy): Control of Hybrid Systems (10h) Erasmus • Apr. 2009: Univ. Tec. Bolivar UTB, Cartagena (Colombia): Tutorial on DES (15h) • May 2014: Univ. Tec. Bolivar UTB, Cartagena (Colombia): Intro. to DES (15h) • Dec 2015: ITB Bandung (Indonesia): Simulation with Petri Nets (10h) Erasmus • Apr. 2017: Univ. of Liverpool (UK): Course 1 (10h) Erasmus • May 2017: ITB Bandung (Indonesia): Simulation with Petri Nets (10h) Erasmus

Figure 2 . 1 :

 21 Figure 2.1: Partition and Automaton

26 Figure 2 . 2 :Figure 2 . 3 :

 262223 Figure 2.2: Phase portrait of Example 1 in PWA

Figure 2 . 4 :

 24 Figure 2.4: Phase portrait of Example 2 in PWA

Figure 2 . 5 :

 25 Figure 2.5: Phase portrait of Example 2 in MLD

Figure 3 . 1 :

 31 Figure 3.1: EMN Cell

Figure 4 . 1 :

 41 Figure 4.1: Smart Grid

Figure 4 . 2 : 40 [ 15 ]

 424015 Figure 4.2: CDG Airport Paris-Roissy

Fig. 1 .

 1 Fig. 1. Equivalence relation between hybrid systems Every well-posed PWA system can be re-written as an MLD system assuming that the feasible states and inputs are bounded [6, proposition 4*]. A completely well-posed MLD system can be rewritten as a PWA system [6, proposition 5*].

  , (a) MLD Model (b) Error between MLD and PWA [4] (C) Error between MLD and PWA-This Work Fig. 3. Simulation Results for the Three-Tank System

  Fig. 4. Simulation results for robotized gear shift

  water source1and normal operation, 2, if water source 2 and normal operation, 3, if maintenance operation, 4, change from maintenance operation

Figure 1 .

 1 Figure 1. States and switching signal for the Lotka-Volterra example.

  ) and Tan et al. (2013) have developed a decentralised technique based on broadcasting and consensus to optimally distribute a resource considering capacity constraints on each entity in the network. Nonetheless, compared to our algorithm, the approach in Dominguez-Garcia et al. (2012) and Tan et al. (2013) is only applicable to quadratic cost functions. On the other hand, Pantoja and Quijano (2012) propose a novel methodology based on population dynamics. The main drawback of this technique is that its performance is seriously degraded when the number of communication links decreases. We point out the fact that other distributed optimisation algorithms can be applied to solve resource allocation problems, as those presented in Nedic, Ozdaglar, and Parrilo (2010), Yi, Hong, and Liu (2015), and Johansson and Johansson (2009). Nevertheless, the underlying idea in these methods is different from the one used in our work, i.e. Nedic et al. (2010), Yi et al. (2015), and Johansson and Johansson (

Theorem 2 . 1 (

 21 adapted from Godsil & Royle, 2001): An undirected graph G of order n is connected if and only if rank(L(G)) = n -1.

Figure  .

  Figure . Single path topology for n nodes.

T 1 , f 1 T 2 , f 2 TFigure  .

 1122 Figure . Decoupled chilled-water plant with n chillers.

PFigure

  Figure . (a) P i -T i curves for each chiller, (b) evolution of supply temperatures and total power consumed by the chillers, C L = . kW, (c) evolution of supply temperatures and total power consumed by the chillers, C L = . kW.

Fig. 1 .

 1 Fig. 1. An AGV circuit (left) and its basic automaton (right)

•

  Chapter 2 presents the Analysis and Control of Hybrid and Switched System • Chapter 3 is devoted to Supervisory Control of Discrete-Event Systems • Chapter 4 gives the Conclusion and Future Work

	NR NR				HDR HDR				8 10
	Promo	Course	CM	PC	TD	TP-MP PFE	Resp	Lang.
	A1	Automatique	10h		10h	10h			Fr.
	A2	Optim.			10h	5h			Fr.
	A2 AII	SED		10h	5h	5h			Fr.
	A3 AII	SysHybrides		7.5h		7.5h			Fr.
	MSc. PM3E Control		7.5h		7.5h			Eng.
	MSc. MOST Simulation		5h	5h	5h			Eng.
	MSc. MOST Resp. MSc+UV						90h	Eng.
	A3	PFE superv.					36h		Fr.
	Masters	PFE superv.					36h		Eng.
	Total 1	272h	10h	30h	30h	40h	72h	90h
	Total 2	283up	15up	36up	30up	40up	72up	90up
			Table 1.1: Teaching in 2015-2016
	B. Responsabilities (Option AII, Auto-Prod, MSc MLPS, MSc MOST)

  Diagnosis and Prognosis of Discrete-Event Systems, 48th IEEE CDC Shanghai, China, Dec 2009 (jointly organized and chaired with Shigemasa Takai). -Invited Session, Diagnosis of DES Systems, 1st IFAC DCDS 2007, Paris, France, June 2007 (jointly organized and chaired with Shigemasa Takai). Hybrid Systems, IEEE ISIC 2001, Mexico City, Mexico, Sep. 2001 (jointly organized and chaired with Michael Lemmon).

	NR	HDR	13
	1.6 Organization of Invited Sessions
	-Invited Session, DES and Hybrid Systems, IEICE NOLTA 2006, Bologna, Italy, Sep. 2006
	(jointly organized and chaired with Shigemasa Takai).
	-Invited Session, Supervisory Control, IFAC WODES, Reims, France, Sep. 2004
	(jointly organized and chaired with Toshimitsu Ushio).
	-Invited Session,
	• Participant, French "Contrat Etat-Région" 2000-2006, CER STIC 9 / N.18036,
	J.J. Loiseau PI, Euro 182,940 (US$ 182,940).
	• Co-Principal Investigator (with Ph. Chevrel), Modeling and Simulation of ESP Program,
	Peugeot-Citroen PSA France, Sep. 2000 -Jan 2001, FF 20,000 (US$ 3,000).

• Co-Principal Investigator (with Andi Cakravastia, ITB), LOG-FLOW, PHC NUSANTARA France Indonesia, Project N. 39069ZJ, 2017, Accepted on 31 May 2017. • Participant, "Industrial Validation of Hybrid Systems", France and Colombia ECOS Nord Project N.C07M03, A. Gauthier and J.J. Loiseau PIs, Jan. 2007 to Dec. 2009 (3 years) Euro 12,000. • Co-Principal Investigator (with J. Aguilar-Martin), Control and Supervision of a Distillation Process, Conseil Régional Midi-Pyrénées, France, 1994-1995, FF 200,000 (US$ 30,000). • Participant, European Esprit Project IPCES (Intelligent Process Control by means of Expert Systems), J. Aguilar-Martin PI, 1989-1992, Euro 500,000 (US$ 500,000).

-Invited Session,

Table 1 .

 1 table. Automatic clustering for symbolic evaluation for dynamical system supervision. In Proc. of IEEE American Control Conference (ACC 1992), Chicago, USA, June 1992, vol. 3, pp. 1895-1897.

		Conf. Book Chap. Book Ed. Journal Total
	1994	2
	1995	2	1	1
	1996	1		1
	1997		1
	1998	2
	1999	1
	2000
	2001	4	1
	2002	1
	2003	5	1
	2004	6
	2005	1
	2006	3
	2007	4
	2008	2
	2009	1
	2010			1
	2011
	2012	1
	2013	2
	2014	4		1
	2015	2	1
	2016			4
	2017 1+3*		1*	1+4*

2: Number of published papers per year (as of 30 June 2017) -where (*) means submitted NR HDR 20 [C.1] N. Rakoto-Ravalontsalama and J. Aguilar-Martin,

  1. Compute matrices A, B, C, D and A i , B i , C i and D i using (2.20).2. InitializeE 1 , E 2 , E 3 , E 4 , E 5 matrices.3. For the m switching regions S j,i , include the inequalities defined in (2.8) or (2.10) which define the values of the m auxiliary binary variables δ j,i .4. Generate 2 * nx δi auxiliary binary dynamical variables associated with the n affine models and m auxiliary binary variables δ j,i associated with the m S ij switching regions.

5.

For i = 1 to n include the inequalities using (2.13) representing the behavior on the x δ vector. 6. For i = 1 to n generate the n c -dimensional Z 1i vector and p c -dimensional Z 2i vector of auxiliary variables Z. 7. For each Z 1i vector introduce the inequalities defined in (2.16), by replacing A i , and B i by A i , and B i , computed in Step 1. M and m are n c -dimensional vectors of maximum and minimum values of x, respectively. 8. For each Z 2i vector introduce the inequalities defined in (2.17), by replacing C i , and D i by C i , and D i , computed in Step 1. M and m are p c -dimensional vectors of maximum and minimum values of x, respectively (This completes the inequality matrices). 9. Compute the matrices defined in (2.21) and (2.22) 10. End.

.3.2 Transitivity of Switchable Languages

  and H i controllable}.

	NR	HDR	34
	3

(n = 3)

Chapter 4 Conclusion and Future Work 4.1 Summary of Contributions

  The results can be found in [C.sub1].

	In this HDR Thesis, I have presented a summary of contribution, in Analysis and Control of Hybrid
	Systems, as well as in Supervisory Control of Discrete-event Systems.
	• Analysis and Control of Hybrid and Switched Systems
	-Modeling and Control of MLD Systems
	-Stability of Switched Systems
	-Optimal Control of Switched Systems
	• Supervisory Control of Discrete-Event Systems
	-Multi-Agent Based Supervisory Control
	-Switched Discrete-Event Systems
	-Switchable Languages of DES
	I have chosen to not present some work like

the Distributed Resource Allocation Problem, the Holonic Systems, and the VMI-Inventory Control work. However the references of the corresponding papers are given in the complete list of publications. My perspectives of research in the coming years are threefold: 1) Control of Smart Grids, 2) Simulation with Stochastic Petri Nets and 3) Planning and Inventory Control.

Table I .

 I Computation and Simulation Times

	Representation	Computation	Simulation
		Time (s.)	Time (s.)
	MLD	-	592.20
	PWA-[4]	93.88	5.89
	PWA-This work	72.90	1.33

Table II .

 II , Computation and Simulation Times

	Representation	Computation	Simulation
		Time (s.)	Time (s.)
	MLD	-	296.25
	PWA-[4]	115.52	0.35
	PWA-This work	155.73	0.17

  translation from the MLD model into PWA model took 572.19 s, with the algorithm proposed here, generating 127 sub-models. The translation into PWA model took 137.37s, with the algorithm in[3], generating 14 submodels. The simulation time for 300 iterations with the MLD model and a MIQP algorithm took 4249.301s, the same simulation with the PWA model obtained with the algorithm proposed here took 0.14s, and the same simulation with the PWA model obtained using the algorithm in[4] took 0.31s. These results are summarized in TableIII, Table III. Computation and Simulation Times.

	Representation	Computation	Simulation
		Time (s.)	Time (s.)
	MLD	-	4249.30
	PWA-[4]	137.37	0.31
	PWA-This work	572.20	0.14

  .t j / 2 , : : : , v .t j / m 1 .t j / D v .t j /. Now, using the equivalence stated in Proposition 5, we know that the solutions of the polynomial Problem (8) are solutions of the switching system; and in this case, it is only one. Hence, we obtain .t j / D v .t j /, which implies that .t j / D v .t j / D m 1 .t j /, where m 1 is the first moment of the vector of moments.

	2d	Á	,
	which implies that

Table  .

  Distributed algorithms' performance.

		Percentage decrease, computation time
	Number of nodes	Proposed approach	DIP	LRE	DIPe
	n =  n =  n =  n =  n = 	   % ,  .   s %, . s %, . s %, . s %, . s	   % ,  .   s %, . s %, . s %, . s %, . s	  % ,  .   s    % ,  s %, . s %,  s %, . s %,  s %, . s -%, . s -
	available communication links); a nonlinear cost func-
	tion φ

Table  .

  Chillers' parameters.

  Evolution of the Euclidean distance and constraint satisfaction using the proposed algorithm. Right y-axis corresponds to the dash-dotted line.

	Euclidean distance (kW)	0 1 2 3 4			Euclidean distance evolution Proposed algorithm	Optimal value
		60				Constraints			1
	Power (kW)	40 50					100 i=1 u i (k)	min{u i (k)}	0 0.5
		30	0 0	0.2 0.2	0.4 0.4	0.6 time 0.6	0.8 0.8	1 1	1.2 1.2 -0.5 x 10 x 10	-4 -4
					iteration 1			iteration 2	iteration 3
							are properly satisfied in steady i (k) = 40 kW and min{u * i=1 u * state, i.e. 100 i (k)} = 0 kW. As a final observation, our algorithm exhibits a suitable

Figure .

D Appendix 4 -Paper [C.sub1]:

©  Informa UK Limited, trading as Taylor & Francis Group

Remerciements

VI. ACKNOWLEDMENT This work has been supported in part by « Contrat Etat -Région No STIC 9-18036, 2000-2006 », Nantes, France. The authors are thankful to the reviewers for their valuable comments and suggestions. ACKNOWLEDGEMENTS This study was supported by Proyecto CIFI 2011, Facultad de Ingeniería, Universidad de Los Andes. ACKNOWLEDGMENT Part of this work was carried out when the second author (N.R.) was visiting Prof. Stephane Lafortune at University of Michigan, Ann Arbor, MI, USA, in Sep. 2013. Grant #EMN-DAP-2013-09 is gratefully acknowledged.

Funding G. Obando is supported in part by Convocatoria 528 Colciencias-Colfuturo and in part by OCAD-Fondo de CTel SGR, Colombia (ALTERNAR project, BPIN 20130001000089).