Zsófia Lendek
email: zsofia.lendek@aut.utcluj.ro

Jimmy Lauber
email: jlauber@univ-valenciennes.fr

Thierry Marie Guerra
email: guerra@univ-valenciennes.fr

Paula Raica
email: paula.raica@aut.utcluj.ro

On stabilization of discrete-time periodic TS systems

Keywords: Takagi-Sugeno fuzzy models, controller design, periodic systems, nonquadratic Lyapunov function

In this paper we consider controller design for periodic Takagi-Sugeno fuzzy models. For this, we use a periodic nonquadratic Lyapunov function defined at the time instants when the subsystems switch. Using the proposed conditions we are able to handle periodic Takagi-Sugeno systems where the local models or even the subsystems are unstable or cannot be stabilized. The application of the conditions is illustrated on numerical examples.

I. INTRODUCTION

This paper deals with a particular class of nonlinear models with periodic parameters. Such models can be found in numerous domains such as computer control of industrial processes, or automotive, aeronautic, and aerospace industries.

For instance examples related to computer control and communication systems are provided in [START_REF] Bittanti | Invariant representations of discrete-time periodic systems[END_REF], to estimate the air/fuel ratio in each cylinder on an internal combustion engine a periodic dynamic model is used in [START_REF] Chauvin | Real-time nonlinear individual cylinder air-fuel ratio observer on a diesel engine test bench[END_REF], [START_REF] Gaiani | Discrete-time analysis of HHC schemes for helicopter vibration attenuation[END_REF] considers a periodic model for the rotor blades of helicopter, etc.

The stability of linear periodic systems is characterized by the monodromy transition matrix and by its eigenvalues [START_REF] Farkas | Periodic motions[END_REF]. For the stabilization of linear periodic models, results have been presented in [START_REF] Farges | Robust H 2 performance analysis and synthesis of linear polytopic discrete-time periodic systems via LMIs[END_REF]. For models including time varying delays, [START_REF] Stepan | Stability of time-periodic and delayed systems -a route to act and-wait control[END_REF] proposed methods based on Floquet's transformation, which led to conditions for exponential stability. Results available for linear periodic systems have been extended to polytopic LPV periodic models, where the stability analysis is based on the use of quadratic [START_REF] Farges | Robust H 2 performance analysis and synthesis of linear polytopic discrete-time periodic systems via LMIs[END_REF], [START_REF] Arzelier | Robust analysis and synthesis of linear polytopic discrete-time periodic systems via LMIs[END_REF] or non-quadratic [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF] Lyapunov functions.

In this paper, to represent nonlinear periodic systems, we use Takagi-Sugeno (TS) fuzzy models [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF]. TS models are nonlinear, convex combinations of local linear models, and are able to exactly represent a large class of nonlinear systems [START_REF] Zs | Stability analysis and nonlinear observer design using Takagi-Sugeno fuzzy models[END_REF].

For the stability analysis and controller design of TS systems the direct Lyapunov approach has been used, involving quadratic Lyapunov functions [START_REF] Tanaka | Fuzzy regulators and fuzzy observers: relaxed stability conditions and LMI-based designs[END_REF]- [START_REF] Sala | Perspectives of fuzzy systems and control[END_REF], piecewise continuous Lyapunov functions [START_REF] Johansson | Piecewise quadratic stability of fuzzy systems[END_REF], [START_REF] Feng | Stability analysis of discrete-time fuzzy dynamic systems based on piecewise Lyapunov functions[END_REF], and nonquadratic Lyapunov functions [START_REF] Guerra | LMI-based relaxed nonquadratic stabilization conditions for nonlinear systems in the Takagi-Sugeno's form[END_REF]- [START_REF] Mozelli | Reducing conservativeness in recent stability conditions of TS fuzzy systems[END_REF]. In particular for discrete-time systems, non-quadratic Lyapunov functions have shown a real improve-ment of the design conditions [START_REF] Guerra | LMI-based relaxed nonquadratic stabilization conditions for nonlinear systems in the Takagi-Sugeno's form[END_REF], [START_REF] Ding | Further studies on LMI-based relaxed stabilization conditions for nonlinear systems in Takagi-Sugeno's form[END_REF]- [START_REF] Lee | Approaches to extended nonquadratic stability and stabilization conditions for discrete-time Takagi-Sugeno fuzzy systems[END_REF]. The stability or design conditions are generally derived in the form of linear matrix inequalities (LMIs).

Periodic TS systems [START_REF] Kerkeni | Estimation of individual incynlinder air mass flow via periodic observer in Takagi-Sugeno form[END_REF], [START_REF] Kerkeni | Individual exaust gas mass flow estimation using a periodic observer design[END_REF], often described by continuous dynamics and discrete dynamics as well as their interactions, have been investigated mainly in the continuous case where the stability is based on the use of a quadratic Lyapunov function [START_REF] Tanaka | Switching control of an R/C hovercraft: stabilization and smooth switching[END_REF]- [START_REF] Ohtake | Switching fuzzy controller design based on switching Lyapunov function for a class of nonlinear systems[END_REF] or a piecewise one [START_REF] Feng | Controller synthesis of fuzzy dynamic systems based on piecewise Lyapunov functions[END_REF], [START_REF]H∞ controller design of fuzzy dynamic systems based on piecewise Lyapunov functions[END_REF]. Although results are available for discrete-time linear switching systems [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF], for discrete-time TS models, few results exist [START_REF] Doo | Output-feedback control of discrete-time switching fuzzy system[END_REF], [START_REF] Dong | H∞ controller synthesis via switched PDC scheme for discrete-time T-S fuzzy systems[END_REF].

In this paper, we consider stability analysis and controller design for discrete-time periodic TS models. To develop the conditions we use a periodic non-quadratic Lyapunov function defined in the moments when the switching takes place. The developed conditions are bilinear matrix inequalities (BMIs), which can be solved using either BMI solvers or an iterative procedure. Using the developed conditions, we are able to prove the stability of a periodic TS system having nonstable local models and stabilize a periodic system having noncontrollable local models.

The structure of the paper is as follows. Section II presents the notations and the general form of the TS models used in this paper. In Section III the proposed conditions for stability analysis are developed, discussed, and illustrated on a numerical example. Section IV extends the conditions for controller design, discusses them, and illustrates their use on a numerical example.

II. PRELIMINARIES

In this paper we consider stability analysis and controller design of discrete-time periodic TS systems. For stability analysis, we consider subsystems of the form

x(k + 1) = rj i=1 h j,i (z j (k))A j,i x(k) = A j,z x(k) (1)
and for controller design

x(k + 1) = rj i=1 h j,i (z j (k))(A j,i x(k) + B j,i u(k)) = A j,z x(k) + B j,z u(k) (2)
where j is the index of the subsystem, j = 1, 2, . . . , n s , n s being the number of the subsystems, x denotes the state vector, r j is the number of rules in the jth subsystem, z j is the scheduling vector, h j,i , i = 1, 2, . . . , r j are normalized membership functions, and A j,i and B j,i , i = 1, 2, . . . , r j , j = 1, 2, . . . , n s , are the local models. We assume that none of the subsystems has a finite escape time.

We consider periodic systems, i.e., the subsystems defined above are activated in a sequence , etc., where p i denotes the number of samples for which the ith subsystem is active. In what follows, we will refer to p i as the period of the ith subsystem.

0 and I denote the zero and identity matrices of appropriate dimensions, and a ( * ) denotes the term induced by symmetry. The subscript z + m (as in A 1,z+m ) stands for the scheduling vector being evaluated at the current sample plus mth instant, i.e., z(k + m). An underlined variable j denotes the modulo of the variable, i.e., j = (j mod n s ) + 1.

In what follows, we will make use of the following results: Lemma 1. [START_REF] Skelton | A Unified Approach to Linear Control Design[END_REF] Consider a vector x ∈ R nx and two matrices

Q = Q T ∈ R nx×nx and R ∈ R m×nx such that rank(R) < n x .
The two following expressions are equivalent:

1) x T Qx < 0, x ∈ {x ∈ R nx , x = 0, Rx = 0} 2) ∃M ∈ R nx×m such that Q + M R + R T M T < 0
Controller design for TS models often lead to double-sum negativity problems of the form

x T r i=1 r j=1 h i (z(k))h j (z(k))Γ ij x < 0 (3)
where Γ ij , i, j = 1, 2, . . . , r are matrices of appropriate dimensions.

Lemma 2. [33]

The double-sum

(3) is negative, if Γ ii < 0 Γ ij + Γ ji < 0, i, j = 1, 2, . . . , r, i < j Lemma 3. [34] The double-sum (3) is negative, if Γ ii < 0 2 r -1 Γ ii + Γ ij + Γ ji < 0, i, j = 1, 2, . . . , r, i = j
Property 1. (Congruence) Given a matrix P = P T and a full column rank matrix Q it holds that

P > 0 ⇒ QP Q T > 0 Property 2.
Let A and B be matrices of appropriate dimensions and ranks, with B = B T > 0. Then

(A -B) T B -1 (A -B) ≥ 0 ⇐⇒ A T B -1 A ≥ A + A T -B

III. STABILITY ANALYSIS OF PERIODIC SYSTEMS

A. Stability conditions

In this section, we consider the stability analysis of periodic TS systems of the form (1). In our previous paper [START_REF] Zs | Switching Lyapunov functions for periodic TS systems[END_REF], we have considered a switching Lyapunov function defined only in the instants when a switching takes place in the system, and the stability conditions have been derived based on the requirement that in each switching part, this Lyapunov function decreases. Generalizing the result in [START_REF] Zs | Switching Lyapunov functions for periodic TS systems[END_REF], this can be formulated as follows:

Theorem 1. The periodic TS system (1) with periods p 1 , p 2 , . . . , p ns is asymptotically stable, if there exist P j,i = P T j,i > 0, M j,i , j = 1, 2, . . . , n s , i = 1, 2, . . . , r j , such that the following conditions are satisfied:

0 B B B @ -Pj,z ( * ) . . . ( * ) ( * ) M A1,0 -M1,0 + ( * ) . . . ( * ) ( * ) . . . . . . . . . . . . . . . 0 0 . . . M A1,p j+1 -1 Ωj+1,j+1 1 C C C A < 0 (4)
where

1 M A 1,b = M j+1,z+b A j+1,z , M 1,b = M j+1,z+b , and Ω j+1,j+1 denotes Ω j+1,j+1 = -M j+1,z+pj+1-1 + ( * ) + P j+1,z+pj+1 .
As already noted, the result above is based on the requirement that the Lyapunov function should decrease in each part. However, this is not necessary, as long as from the beginning and until the end of one complete cycle the Lyapunov function decreases, i.e., the Lyapunov function steadily decreases in each cycle. Consequently, the Lyapunov function may increase at one switch. To quantify this increase, let us impose the condition that V (k + p i ) < δ i V (k), where δ i > 0, but not necessarily subunitary. The decrease of the Lyapunov function during the full cycle can be formulated as

V (k + ns i=1 p i ) < δ π V (k), where δ π = ns i=1 δ i , with δ π < 1.
Then, combining it with Theorem 1, we have the following result.

Theorem 2. The periodic TS system (1) with periods p 1 , p 2 , . . . , p ns is asymptotically stable, if there exist P j,i = P T j,i > 0, M j,i , δ i > 0, j = 1, 2, . . . , n s , i = 1, 2, . . . , r j , such that the following conditions are satisfied:

0 B B B @ -δjPj,z ( * ) . . . ( * ) ( * ) M A1,0 -M1,0 + ( * ) . . . ( * ) ( * ) . . . . . . . . . . . . . . . 0 0 . . . M A1,p j+1 -1 Ωj+1,j+1 1 C C C A < 0 Π ns i=1 δi ≤ 1 (5)
where

M A 1,b = M j+1,z+b A j+1,z+b , M 1,b = M j+1,z+b ,
and Ω j+1,j+1 denotes Ω j+1,j+1 = -M j+1,z+pj+1-1 + ( * ) + P j+1,z+pj+1 .

Proof. Consider the switching Lyapunov function, defined only in the instants when a switching takes place in the system:

V (x(k)) = x(k) T P j,z x(k) if the active subsystem was j
Then, the condition V (x(k+p j+1 )) < δ j V (x(k)), for some δ j > 0, can be written as

V (x(k + p j+1 )) -δ j V (x(k)) = x(k) x(k + p j+1 ) T -δ j P j,z 0 0 P j+1,z+pj+1 x(k) x(k + p j+1 )
The system dynamics during the p j+1 samples are

Υ j+1      x(k) x(k + 1)
. . .

x(k + p j+1 )      = 0 with Υ j+1 =      A j+1,z -I . . . 0 0 0 A j+1,z+1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 . . . A j+1,z+pj+1-1 -I      Using Lemma 1, V (x(k + p j+1 )) < δ j V (x(k)), if there exists M such that      -δ j P j,z 0 . . . 0 0 0 . . . 0 . . . . . . . . . . . . 0 0 . . . P j+1,z+p j+1      + M Υ j+1 + ( * ) < 0 A particular choice of M is M =        0 0 . . . 0 M j+1,z 0 . . . 0 0 M j+1,z+1 . . . 0 . . . . . . . . . . . . 0 0 . . . M j+1,z+pj+1-1       
which leads to the sufficient conditions

0 B B B @ -δjPj,z ( * ) . . . ( * ) ( * ) M A1,0 -M1,0 + ( * ) . . . ( * ) ( * ) . . . . . . . . . . . . . . . 0 0 . . . M A1,p j+1 -1 Ωj+1,j+1 1 C C C A < 0 Π ns i=1 δi ≤ 1 where M A 1,b = M j+1,z+b A j+1,z+b , M 1,b = M j+1,z+b ,
and Ω j+1,j+1 denotes Ω j+1,j+1 = -M j+1,z+pj+1-1 + ( * ) + P j+1,z+pj+1 . During one complete cycle, the above conditions amount to

V (x(k + ns i=1 p i )) < δ ns V (x(k + ns-1 i=1 p i )) < < δ ns δ ns-1 V (x(k + ns-2 i=1 p i )) < . . . < ... < Π ns i=1 δ i V (x(k)) Consequently, if ns i=1 δ i ≤ 1,
V is decreasing during a complete cycle, and the periodic system is asymptotically stable.

B. Discussion and example

In this section we discuss the conditions developed above. First of all, it has to be noted that the conditions of Theorem 2 are not related to the stability of the individual subsystems. Indeed, switching between unstable subsystems can lead to a periodic system that is asymptotically stable. The introduction of the constants δ i actually allows for the increase of the Lyapunov function while one or several subsystems are active. However, due to the condition ns i=1 δ i < 1, during a whole cycle (be that from the first to the first or from the second to second, etc. subsystem), the Lyapunov function decreases. Even if this condition is not satisfied, ns i=1 δ i gives an upper bound on the increase of the Lyapunov function during a cycle and therefore a measure of distance from proving stability. The proposed conditions, as they are stated in Theorem 2, are unfortunately BMIs. However, the search for δ i , i = 1, 2, . . . , n s can be done iteratively. Alternatively, one can use e.g., available BMI solvers to solve directly the BMI problem.

It is important to note that the conditions above can easily be changed to verify the instability of the TS model, as follows.

Theorem 3. The equilibrium point x = 0 of the periodic TS system (1) with periods p 1 , p 2 , . . . , p ns is locally unstable, if there exist P j,i = P T j,i > 0, M j,i , σ i > 0, j = 1, 2, . . . , n s , i = 1, 2, . . . , r j , such that the following conditions are satisfied: The proof follows the same lines as of Theorem 2, by imposing the condition that σ j V (x(k + p j+1 )) > V (x(k)), and is therefore not repeated here.

Let us now illustrate the application of Theorem 2 on an example. Not all the local matrices are stable: A 1,2 and A 3,2 are unstable. However, the periodic system is stable, as indicated by the trajectories in Figure 1. For these particular trajectories, the membership functions used were using the relaxation 2, we obtain positive definite P s and the following δs: δ 1 = 0.79216 δ 2 = 0.79088 δ 3 = 0.79061 Although in the example above we used a BMI solver to solve the conditions of Theorem 2, and therefore prove the stability of the periodic TS system, an iterative procedure can also be used as follows. It has to be noted that for large enough, fixed δs the conditions of Theorem 2 are feasible. Once such δs have been found, the conditions become LMIs and can be solved for P and M . Afterward, one can iterate by solving for δ and M and P and M , respectively.

Example 1. Consider a periodic TS model composed of 3 subsystems, each being active for p

i = 2 samples, i = 1, 2,
h 1,1 (z) = e -x 2 1 , h 1,2 (z) = 1 -h 1,1 (z), h 2,1 (z) = cos(x 1 ) 2 , h 2,2 (z) = 1 -h 2,1 (z), h 3,
Another possibility is to consider the α-sample variation [START_REF] Kruszewski | Stabilization of a class of nonlinear model with periodic parameters in the Takagi-Sugeno form[END_REF] of the Lyapunov function. In this way, the necessity of solving BMIs is circumvented, as all the δs are gathered in a single one, ultimately leading to a decrease of the Lyapunov function through a cycle. However, in this way the number of the LMIs increases very much, eventually leading to unfeasibility due to the limitation of the available solvers.

IV. CONTROLLER DESIGN

A. Design conditions

In this section we extend the conditions developed in Section III-A for stability analysis for controller design. For this, we consider the periodic TS model described by (2) and a controller of the form

u(k) = -F j,z H -1 j,z x(k) (7)
if the jth subsystem is active at sample k. The closed-loop system is then given by

x(k + 1) = (A j,z -B j,z F j,z H -1 j,z )x(k) (8)
2 To solve the BMIs, the solver PenBMI [START_REF] Kočvara | Penbmi, version 2.1[END_REF] has been used.

if the jth subsystem is active, this again being a periodic TS system. For (8) the following result can be formulated. Theorem 4. The controller [START_REF] Arzelier | Robust analysis and synthesis of linear polytopic discrete-time periodic systems via LMIs[END_REF] asymptotically stabilizes the periodic TS model [START_REF] Arzelier | Robust analysis and synthesis of linear polytopic discrete-time periodic systems via LMIs[END_REF] if there exist P j,i = P T j,i > 0, H j,i , δ j > 0, j = 1, 2, . . . , n s , i = 1, 2, . . . , r j , so that

0 B B B @ -δ j Ω 1,1 ( * ) . . . ( * ) ( * ) G 1,0 -H j,z+1 + ( * ) . . . ( * ) ( * ) . . . . . . . . . . . . . . . 0 0 . . . G 1,p j+1 -1 -P j+1,z+p j+1 1 C C C A < 0 Π ns i=1 δ i ≤ 1
where

Ω 1,1 = H j,z + H T j,z -P j,z and G 1,b = A j+1,z+b H j+1,z+b -B j+1,z+b F j+1,z+b , j = 1, 2, . . . , n s .
Proof. Similarly to stability analysis, consider a switching Lyapunov function, defined only in the instants when a switching takes place in the system of the form

V (x(k)) = x(k) T P -1 j,z x(k)
if the active subsystem was j Then, the condition V (x(k+p j+1 )) < δ j V (x(k)), for some δ j > 0, can be written as

V (x(k + p j+1 )) -δ j V (x(k)) = x(k) x(k + p j+1 ) T -δ j P -1 j,z 0 0 P -1 j+1,z+pj+1 x(k) x(k + p j+1 )
The system dynamics during the p j+1 samples are

Υ j+1      x(k) x(k + 1)
. . .

x(k + p j+1 )      = 0 with Υ j+1 =      Υ 1,0 -I . . . 0 0 0 Υ 1,1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 . . . Υ 1,pj+1-1 -I      where Υ 1,b denotes Υ 1,b = A j+1,z+b - B j+1,z+b F j+1,z+b H -1 j+1,z+b . Using Lemma 1, V (x(k + p j+1 )) < δ j V (x(k)), if there exists M such that      -δ j P j,
1 C C C A < 0 Π ns i=1 δ i ≤ 1
where Γ 1,1 = -δ j H T j,z P -1 j,z H j,z and G 1,b = A j+1,z+b H j+1,z+b -B j+1,z+b F j+1,z+b , j = 1, 2, . . . , n s . Applying Property 2 we obtain the conditions of Theorem 4. Moreover, during one complete cycle, the above conditions amount to

V (x(k + ns i=1 p i )) < δ ns V (x(k + ns-1 i=1 p i )) < < δ ns δ ns-1 V (x(k + ns-2 i=1 p i )) < . . . < ... < Π ns i=1 δ i V (x(k))
Consequently, if ns i=1 δ i ≤ 1, V is decreasing during a complete cycle, and the periodic system is asymptotically stabilized.

B. Discussion and example

Let us now discuss the conditions developed above. Similarly to the stability conditions, it has to be noted that the conditions of Theorem 4 are not related to the controllability of the individual subsystems. Switching between subsystems that are not stabilizable on their own can lead to a periodic system that is asymptotically stabilized. Again, if the condition ns i=1 δ i < 1 is not satisfied, ns i=1 δ i gives an upper bound on the increase of the Lyapunov function during a cycle and therefore a measure of distance from stabilization, together with a hint of which subsystem is "problematic". The proposed conditions are BMIs, therefore either a BMI solver is needed to solve them or an iterative procedure has to be employed.

Although similarly to the stability analysis, the conditions can easily be changed to prove that a periodic system cannot be stabilized, this proof would unfortunately be valid only for the control structure and Lyapunov function considered. Consequently, one would only be able to prove that there is no control law of the form (7) that stabilizes the periodic system.

To illustrate the control design using the conditions of Theorem 4 consider the following example.

B 1,1 = B 1,2 = 0 1 B 2,1 = B 2,2 = 0 0
Since B 2,1 = B 2,2 = 0, the local matrices of the second subsystem cannot be stabilized. However, the periodic system can be stabilized, as indicated by the trajectories in Figure 2. For these particular trajectories, the membership functions used were h 1,1 (z) = e -x 2 1 , h 1,2 (z) = 1h 1,1 (z), h 2,1 (z) = cos(x 1 ) 2 , h 2,2 (z) = 1h 2,1 (z), and the initial states were x(0) = (5, -4) T . The control input is presented in Figure 3. Since the input gains of the second subsystem are zero, this subsystem is not controlled. Similarly to stability analysis, although we used a BMI solver to solve the conditions of Theorem 4, an iterative procedure can also be used.

Solving the conditions of Theorem 4 we obtain

While we do not present it in this paper, a similar result can be obtained for observer design, under the assumption that the scheduling variables are known. In such a case, it is possible to design an observer that is able to estimate the states of a periodic TS model of whose local models are unobservable.

  1,b = M j+1,z+b A j+1,z+b , M 1,b = M j+1,z+b + ( * ), and Ω denotes Ω = -P j,z + M A 1,0 + ( * ).

1 Fig. 1 .

 11 Fig. 1. Trajectories of the periodic system in Example 1.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Trajectories of the stabilized periodic system in Example 2.

Example 2 .

 2 Consider a periodic TS model composed of 2 subsystems, each being active for p i = 2 samples, i = 1, 2,

	with the local matrices given by
	A 1,1 =	0.42 0.73 0.89 0.57	A 1,2 =	0.04 0.56 0.67 -0.25
	A 2,1 =	-0.37 -1.47 -0.30 -0.23	A 2,2 =	0.11 1.44 0.31 -0.35

  3 the follow-3 All values are truncated to two decimal places.

	ing controller gains:
	H1,1 = 10 4	" -0.86 3.51 5.68 0.67 «		H1,2 = 10 4	-0.34 3.65 " 5.60 0.02	«
	H2,1 = 10 4	" -0.87 5.61	-0.58 1.67	«	H2,2 = 10 4	0.01 " 5.54 -0.26 2.06	«
	F1,1 = 10 4		" 4.93 3.10	«				F1,2 = 10 4	-0.56 " 3.92 «
	with
	P1,1 = 10 4	" -1.04 4.00	-1.04 0.88	«	P1,2 = 10 4	-0.06 " 3.90	0.93 -0.06	«
	P2,1 = 10 4	" -0.22 6.86	-0.22 2.93	«	P2,2 = 10 4	0.06 3.32 " 6.94 0.06 «
	δ1 = δ2 = 0.5

The modulo is used for the ease of notation, as, after the last subsystem, due to the periodicity, follows the first one.

ACKNOWLEDGMENT

This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS -UEFISCDI, project number PN-II-RU-TE-2011-3-0043, contract number 74/05.10.2011, by International Campus on Safety and Intermodality in Transportation the European Community, the Délegation Régionale à la Recherche et à la Technologie, the Ministére de L'Enseignement supérieur et de la Recherche the Region Nord Pas de Calais and the Centre Nationale de la Recherche Scientifique: the authors gratefully acknowledge the support of these institutions.


Benjamin Pillot

Nadeem Al-Kurdi

Carmen Gervet

Laurent Linguet

Optical network ring analogy for effective solar PV site selection using GIS and optimization modelling

Keywords: Energy planning, Spatial Decision support, Spatio-temporal data MOTS-CLÉS : transition énergétique, aide à la décision spatiale, données spatio-temporelles

Energy transition requires a holistic approach, involving land, resources, environmental and economical data and constraints. The core purpose of energy transition, is to migrate power systems towards renewable energy usage (solar, hydro, biomass, wind), in a technical and economical viable manner. In this article we address this challenge as a spatio-temporal analysis problem combined with decision support targeted for policy makers and investors. We show how the selection of potential PV sites, can efficiently be combined with the optimization of the connection to the grid in terms of operational costs. Our contribution lies in bringing forward a modeling analogy with the SONET problem, addressed in fiber network designs. Improved results compared to an existing GIS-optimization PV site placement approach are illustrated on a real case study and data from the French Guiana.

RÉSUMÉ. La transition énergétique nécessite une approche holistique, impliquant la gestion des sols, ressources potentielles, les données et contraintes environnementales et économiques. L'objectif central est de migrer les systèmes électriques vers l'utilisation des énergies renouvelables (solaire, hydroélectrique, biomasse, éolienne), de manière techniquement et économiquement viable. Dans cet article, nous abordons ce défi comme un problème d'analyse spatio-temporelle intégré à un modèle d'aide à la décision conçu pour les décideurs politiques et les investisseurs. Nous montrons

comment la recherche de sites photovoltaïques peut être combinée efficacement avec l'optimisation du raccordement au réseau en termes de coûts. Notre contribution réside dans la mise en place d'une analogie de modélisation avec le problème SONET, défini dans les conceptions de réseaux de fibre optique. La qualité des résultats comparés à une approche existante de placement de parcs photovoltaïques est présentée sur un cas d'étude en Guyane française.

Introduction

Nowadays, the main goal of most countries' energy planning policy [START_REF] Hache | Renewable energy source integration into power networks, research trends and policy implications: A bibliometric and research actors survey analysis[END_REF] in terms of energy transition, relies on the integration of renewablebased generation in power networks. However, increasing the share of renewable energy (RE) sources is still challenging, as a result of their inherent intermittent and geographical dispersion [START_REF] Hache | Renewable energy source integration into power networks, research trends and policy implications: A bibliometric and research actors survey analysis[END_REF]. Specific planning strategies [START_REF] Gorsevski | A group-based spatial decision support system for wind farm site selection in Northwest Ohio[END_REF][START_REF] Camargo | Spatiotemporal modelling for integrated spatial and energy planning[END_REF]) must be developed accordingly, in order to reach energy transition targets without threatening the existing infrastructures.

Recent works propose an integrated model framework, combining Geographic Information Systems and Robust Optimization, such as the GREECE-OPSPV (Geographical REnewable Energy Candidate Extraction -Optimal Planning and Sizing of PV parks) system [START_REF] Al-Kurdi | Towards robust scenarios of spatio-temporal renewable energy planning: A gis-ro approach[END_REF][START_REF] Pillot | An integrated GIS and robust optimization framework for solar PV plant planning scenarios at utility scale[END_REF]. It identifies the optimal sites for solar photovoltaic (PV) generation at utility scale, according to specific energy planning targets. The approach connects both geographical information system (GIS) and robust linear optimization (RO) through the pipe depicted in Figure 1. The GIS component gathers large heterogeneous sets of spatiotemporal data, and allows for location and size of the best solar PV sites to be retrieved with respect to geographical constraints (restricted areas, land use, distance to grid, etc.), spatial dispersion of the resource, hourly global energy demand and generation, predefined planning scenarios, and the degree of risk adversity authorized by the decision maker. Through a set of spatio-temporal data layers and control parameters, the GIS module (GREECE) converts spatial constraints and parcels into items characterized by de-spatialized attributes and solar resource time series. Based on current electricity generation and demand time series as well as projection scenarios, the RO model (OPSPV) then looks for optimal site candidates (location, size and power), with respect to given temporal constraints (parcel size, hourly electricity demand, maximum penetration of intermittent RE power) and objectives (maximize energy generation and minimize total costs). The model eventually gives an estimate of the risk associated with solar PV investment at utility scale to the decision and policy maker, by means of a Pareto approach (cost vs. energy generation) and according to best and worst case scenarios.

In this integrated GIS and decision support model, when a solar PV site (location and size) is proposed, a new substation is considered for connection to the power network. This approach was justified as the model was designed for private investors in RE plants, thus each solar PV site required its own connection to the grid. However, in island networks (i.e. not interconnected), the network manager might have the authority to decide which site should be turned on and off, depending on the share of intermittent power (Notton et al., 2019). Accordingly, some investors have moved on to diversifying their RE investments (hydro, biomass, etc.) or exploiting power plants with added storage capacity [START_REF] Spaes | Inauguration de la centrale solaire avec stockage de Savane des Pères en Guyane[END_REF][START_REF] Tapachès | Mitigation of the variability of a PV fleet via geographical dispersion and energy storage systems on the Reunion island non-interconnected grid[END_REF]. Another solution though is to rely on co-investment, whereby several sites belonging to multiple investors would be connected to the same substation. To our knowledge this approach hasn't been investigated in terms of simulation models. This would be a cost-effective solution as the maximum power per RE site is constrained by the technical features of the power network. In fact, for multiple small-scale solar PV facilities, the final cost would grow with the number of corresponding substations. In this paper we address this problem by deriving a model that would mitigate the cost by aggregating various power plants around one unique substation. This objective comes with new challenges on the GIS and optimization model.

The problem tackled is to maximize PV energy supply and minimizing costs through small-scale power facilities. Our new framework is based on an analogy with the Synchronous Optical NETwork (SONET) [START_REF] Pelleau | SONET network design problems[END_REF][START_REF] Goldschmidt | SONET/SDH ring assignment with capacity constraints[END_REF], which is a popular network design in the field of fiber-optic technology. Typically, we first tailor the constraints of the SONET optimization problem to our own problem and then integrate them into the RO SAGEO'2021 model of the framework. Since this raises the computational complexity of the combinatorial optimization, the data range of those constraints is first pruned by enhancing the GREECE model with respect to the spatial constraints that apply to the potential sites. It ensures keeping the CPU time reasonably low without loss of solution quality. With the proposed architecture, the best solar PV sites are now gathered into rings, i.e. local aggregates of plants, connected to one sole substation. It results in lower energy unit cost (Me/GWh) as well as means for the decision makers and grid managers to manage their risk adversity, compared to the previous architecture.

To describe the resulting benefits, we compare both approaches, with and without considering constraints from the SONET problem, by applying the GREECE-OPSPV framework to the real case of French Guiana. The article is structured as follows: section 2 presents the revisited version of GREECE-OPSPV based on the SONET analogy; section 3 is the experimental section based on a real-world case study and compare both results from the former and the new version of GREECE-OPSPV; section 5 concludes the paper.

Applying the SONET design to RE sources in power systems

Planning small-scale PV sites for safer grid operation and improved risk management

Safer grid operation. When too much intermittent power is injected into the grid, it might be necessary for the network manager to disconnect some of the solar PV sites [START_REF] Tapachès | Mitigation of the variability of a PV fleet via geographical dispersion and energy storage systems on the Reunion island non-interconnected grid[END_REF]Notton et al., n.d.). When there are only large power plants, the loss of one of them could affect grid stability. The use of smaller facilities therefore gives more flexibility when handling the energy supplied from intermittent sources, and eventually ensures safer grid operation once solar PV sites have been commissioned. This idea of targeting smaller solar PV sites is even more relevant in small not interconnected electricity networks, whose low inertia actually implies higher impact from intermittent RE sources on frequency variability and grid stability (Notton et al., n.d.;[START_REF] Tapachès | Mitigation of the variability of a PV fleet via geographical dispersion and energy storage systems on the Reunion island non-interconnected grid[END_REF]. Typically, the GREECE-OPSPV model allows risk adversity to be assessed upstream using Pareto optimisation through best and worst case scenarios [START_REF] Al-Kurdi | Towards robust scenarios of spatio-temporal renewable energy planning: A gis-ro approach[END_REF][START_REF] Pillot | An integrated GIS and robust optimization framework for solar PV plant planning scenarios at utility scale[END_REF]. In this paper, we consider smaller PV sites, and look for handling the risk downstream: whether the riskier solution is chosen or not, the grid manager must remain capable of eventually switching off a plant without jeopardizing the whole network.

Risk adversity management. When there is one unique investor deploying multiple sites, it is also more interesting to invest in various small-scale facilities rather than in few bigger ones. As a matter of fact, the grid manager will not care about the economical risk but rather about the technical risk: he will eventually shut down any solar PV site in case of too much intermittent power injected into the network (Notton et al., n.d.;[START_REF] Tapachès | Mitigation of the variability of a PV fleet via geographical dispersion and energy storage systems on the Reunion island non-interconnected grid[END_REF]. This technical liability might not be related to technical constraints, as this is the case for EDF in French Guiana which disconnects plants in the chronological order of their connection to the grid [START_REF] Ctg | Systèmes énergétiques insulaires Guyane -bilan prévisionnel de l'équilibre offre / demande d'électricité[END_REF][START_REF] Edf | Open Data EDF Guyane[END_REF]. As a result, while the investor may ensure for instance to still have 6 or 7 out of 8 smallscale PV plants running when some sites are eventually disconnected, he could loose half of the production if he only runs two large-scale facilities. This is an economical risk that some investors might not be willing to take.

Clusters of RE sources: towards the SONET analogy. Previous works on the GREECE-OPSPV model, sought each solar PV power plant to be built along its own substation [START_REF] Al-Kurdi | Towards robust scenarios of spatio-temporal renewable energy planning: A gis-ro approach[END_REF][START_REF] Pillot | An integrated GIS and robust optimization framework for solar PV plant planning scenarios at utility scale[END_REF]. As mentioned, this approximation is no longer valid when the number of sites dramatically increases, which is the case if we only allow small-scale plants to be eventually built. In order to keep the solutions economically relevant, we revisit this model so that several power plants might actually be connected to the same substation. We present the underlying methodology in the following section.

Adapting SONET constraints to the RE planning problem

In the field of fiber-optic technology, one of the most popular network designs is the Synchronous Optical NETwork (SONET) [START_REF] Pelleau | SONET network design problems[END_REF][START_REF] Goldschmidt | SONET/SDH ring assignment with capacity constraints[END_REF]. In one possible topology, each customer is connected to exactly one local ring through add-drop multiplexers (ADM), and those local rings are all connected to one federal ring through a digital cross connector (DXC) [START_REF] Pelleau | SONET network design problems[END_REF]. This topology is depicted in Figure 2a. The cost [START_REF] Pelleau | SONET network design problems[END_REF] and power network analogy of DXCs is much higher than that of ADMs, so the number of DXCs must be minimized, that is the number of rings. This is known as the SONET Ring Assignment Problem (SRAP) with capacity constraints. This can be formally depicted as a node-partitioning problem for a given graph G [START_REF] Goldschmidt | SONET/SDH ring assignment with capacity constraints[END_REF]. Nodes of G stand for the customers to be connected and the edge weights represent the traffic demand between sites. The analogy with power networks is depicted in Figure 2b and can be summarized as follows: we may think of the federal ring as the electricity network and of the DXCs as the different substations connected to it. Various power plants can then be connected to each substation/local ring with respect to the available hosting capacity. Those power plants can be of different type, such as dispatchable (biomass, hydroelectricity, geothermal, etc.) and non-dispatchable (solar, wind, etc.) RE sources SAGEO'2021 for example. Each ring may be seen as an aggregate of power plants that are spatially close enough to each other to be connected to the same substation.

We now show how constraints for this problem can be derived and adapted from the SRAP [START_REF] Goldschmidt | SONET/SDH ring assignment with capacity constraints[END_REF], as well as developed to answer specific needs. In order to avoid combinatorial explosion, the spatio-temporal GIS model is used to determine pairwise distance between all solar PV potential sites belonging to the set P S. For each ring k, we will therefore limit the available sites to the set B k , that is the set of site indices for which sites P S j have distance d kj from site P S k below a threshold D (including P S k ):

B k = {j | j ∈ N = {1, . . . , n}, d kj D} ∀k ∈ N = {1, . . . , n} (1)
We then apply the constraints from the SONET problem to revisit the former optimisation problem (see Table 1 for the nomenclature):

> max < b k i h SA ik × P pv h,i (2a) s.t. i∈B k SA ik × P nom C ∀k ∈ N = {1, . . . , n}, (2b) k x ik 1 ∀i ∈ N = {1, . . . , n}, ( 2c
) i∈B k x ik M × y k ∀k ∈ N = {1, . . . , n}, (2d)
M × i∈B k x ik y k ∀k ∈ N = {1, . . . , n}, ( 2e
)
i∈B k k SA ik × P pv h,i + Eint h 0.35 × Dem h ∀h, ( 2f
)
i∈B k k SA ik × P pv h,i + Eint h + Ep h Dem h ∀h, (2g)
SA ik Smax i × x ik , ∀i ∈ B k , ∀k ∈ N = {1, . . . , n}, (2h) x ik × Smin SA ik , ∀i ∈ B k , ∀k ∈ N = {1, . . . , n}, (2i)
x ik , y k ∈ {0, 1}.

Let x ik = 1 if site P S i is selected and belongs to ring k and x ik = 0 otherwise. Let y k = 1 if power plants are connected to ring k (i.e. ring is active) and y k = 0 otherwise. Aggregate of PV power per ring is kept below a threshold C corresponding to the substation maximum hosting capacity (2b). Constraint (2c) ensures that a selected site only belongs to one ring. Constraints (2d) and (2e) guarantee that a ring with one or more selected PV sites is active. Constraint (2f) prevents the amount of intermittent energy from being over 35 % of the total forecast energy demand [START_REF] Ctg | Systèmes énergétiques insulaires Guyane -bilan prévisionnel de l'équilibre offre / demande d'électricité[END_REF]EDF, 2017).

Satisfaction of the forecast energy demand is defined by constraint (2g), using existing resources augmented with new PV generation. Constraints (2h) and (2i) relate the size of the PV sites, lying between Smin and Smax i , to whether they are selected or not. This relationship is required to link both the energy production and the different costs. If a plant size is not null then the site is selected, and conversely if a site is not selected then its size is forced to be null.

Essentially the strategic energy planning over some given time horizon has two main objective functions: 1) to maximize the total hourly energy production over the year through new PV energy generation (2a), 2) to minimize the total costs related to PV installation, connection to the grid, etc. Since these functions are in different units they are not combined into a single weighted function that would not be meaningful, but instead solved by seeking the Pareto frontier, i.e. optimizing each function while constraining the other one. As a result, the optimization problem (2a)-( 2i) is applied for every given constrained cost value in the range of the Pareto.

Minimize costs: Modeling non-linear functions. The total cost Cost corresponds to the sum of the costs in every ring k. Each cost in ring k is the aggregate of the capital cost Cap k of all the PV sites built in the ring, plus the connection cost to the grid Ccon k , the substation cost Csta k and the operational & maintenance costs Cop k :

Cost = k Cap k + Cop k + Ccon k + Csta k (3)

SAGEO'2021

Capital cost Cap k and operational & maintenance costs Cop k are defined as piecewise linear functions [START_REF] Pillot | An integrated GIS and robust optimization framework for solar PV plant planning scenarios at utility scale[END_REF], and depend on the nominal power range of the PV sites implemented in ring k. Regarding connection costs Ccon k , using the centroid of all sites as the connection bridge to the sites would in fact make the problem nonlinear. To keep linearity, we have thus considered in first approximation the maximum distance to the grid among PV sites implemented in ring k:

Ccon k = Clan × max i∈B k (dg i × x ik ) ∀k (4)
where Clan stands for the unit cost of transmission lines (e/km). Finally, following the same idea as for the plant capital cost, we also define the substation cost Csta k as a piecewise function. The piecewise model allows for the problem to remain linear and for economies of scale in building substations to be included into the analysis. The substation cost depends on the aggregated PV nominal power in ring k, which actually determines the final substation hosting capacity:

Csta k =      a4 × P nom × ( i∈B k SA ik ) + y4 if 0MW ≤ ( i∈B k SA ik ) × P nom ≤ 10MW a5 × P nom × ( i∈B k SA ik ) + y5 if 10MW ≤ ( i∈B k SA ik ) × P nom ≤ 50MW a6 × P nom × ( i∈B k SA ik ) + y6 if 50MW ≤ ( i∈B k SA ik ) × P nom (5)
3 Case study: ring vs. site approach in French Guiana

In this section, we compare the SONET-based ring approach with the site approach [START_REF] Al-Kurdi | Towards robust scenarios of spatio-temporal renewable energy planning: A gis-ro approach[END_REF][START_REF] Pillot | An integrated GIS and robust optimization framework for solar PV plant planning scenarios at utility scale[END_REF]. We have applied both methods in French Guiana, according to characteristics and policy targets for the horizon 2030 in the region. In the context of the French Energy Transition Act, it is expected to triple RE source capacities by 2023 (CTG, 2017), first by increasing solar PV, then biomass.

Data and processing

The spatio-temporal GIS GREECE model first extracts the potential sites that will feed the OPSPV module depending on geographical constraints and land management. The model prunes the unrestricted territory in order to finally get suitable parcels with respect to minimum and maximum surface area. We refer the reader to [START_REF] Pillot | An integrated GIS and robust optimization framework for solar PV plant planning scenarios at utility scale[END_REF] for further details about the methodology. It resulted in a set P S of 133 land parcels between 1.5 and 50 ha, depicted in Figure 3. They correspond to all the potential sites where solar PV plants could be built, such that all the constraints hold. We will rely on those plots to compare both ring and site approaches.

Beyond geographic location, each potential site generated by GREECE also comes along with specific features and resource time series (e.g. solar irradiation). Those features are extracted from heterogeneous data layers and remote sensing data. Essentially, GREECE allows for de-spatialization and discretiza- tion of the potential sites in order to feed the combinatorial optimization model with items defined as tuples of digitalized attributes. In the present study, for simplicity purposes, items have been defined as tuples of 3 attributes: land surface area, distance to the grid, and time series of solar irradiation values. Our optimization model then converts those values into maximum nominal power, connection costs, and solar PV plant output power profiles.

Finally, hourly energy demand for the 2030 horizon has been projected according to 2016 records (EDF, 2019) and EDF estimations for worst case (5% annual growth) and best case (2% annual growth).

Results and analysis

In this paper, we compare two main aspects of relevance to the decision maker and network manager (power plant energy investors in French Guiana and EDF) between the GREECE-OPSPV model (site approach) and our proposed SONET-based model (ring approach). These aspects are related to what we previously stated as safer grid operations and better risk adversity management (see section 2.1). The comparison is based on the worst case scenario implemented in the site approach. In all cases, we only consider rings made of small-scale facilities whose nominal power is not greater than 5 MW. We first compare Pareto solutions for both approaches (Figure 4) in two fair distinct cases: (a) only small-scale facilities up to 5 MW can be built in the site approach, while no limitation is set on the ring hosting capacity (i.e maximum power that can be injected into the corresponding substation, sum of all PV sites' power attached to that ring); (b) solar PV sites up to 20 MW can be built in the site approach, and so the same applies to ring hosting capacity in the ring approach, which cannot exceed 20 MW as well. The 20 MW limit corresponds to the specific technical constraints of the French Guiana power network [START_REF] Ctg | Systèmes énergétiques insulaires Guyane -bilan prévisionnel de l'équilibre offre / demande d'électricité[END_REF]EDF, 2017). This comparison investigates the relative contribution of each approach with respect to economical benefits, and risk management in terms of spatial distribution of the sites. When only small-scale solar PV sites can be integrated into a power network (a), the Pareto solutions depicted in Figure 4a show that the ring approach gives better results, with for instance 10 MW of extra installed power and 16 GWh of extra power generation compared to the site approach for the same overall cost (60 Me). When both approaches are subject to the same constraints (b), the resulting Pareto solutions remain similar in both cases regardless of the cost value (see Figure 4b). However, as stated in section 2.1, small-scale power stations allow for more manageable risk adversity once it is operating. Hence the ring approach gives the best energy planning strategy with respect to both the decision maker and the network manager needs.

To describe the resulting spatial difference between both models, we plot in Figure 5 the selected PV sites against the selected rings regarding the second case (Figure 4b) and with respect to the same Pareto solution (cost value = 80 Me). In both examples, the selected PV sites are scattered throughout the Northern shore of French Guiana, along the power network (i.e. cost is lower for sites close to the grid). The total installed power is similar (45.4 MW against 44.9 MW) but is allocated among 3 facilities in the site approach (Figure 5a), and among 3 rings and 11 small-scale facilities in the ring approach (Figure 5b). The characteristics of the PV sites in both examples are given in Table 2 and Table 3.

Table 2. Characteristics of the solar PV sites depicted in Figure 5a site ID Nominal power (MW)

Conclusion

In this paper, we exploited the strengths of spatial decision support, applied to the field of energy transition. In particular, we studied the integration of spatial analysis on spatio-temporal data (solar radiation) as a pre-processing step to an optimization module. This module introduces a novel approach to planning the placement of solar PV sites, given the GIS computation of potential sites, such that the costs are further optimized compared to a more traditional approach. We showed and modelled an analogy with the SONET model, to aggregate SAGEO'2021 potential sites around one shared substation that feeds the power grid. The results show substantial gain in costs per KWh produced, as well as enhanced energy stability and risk management given a larger number of selected sites. Further and ongoing work include the extension of our GIS and optimization model, to assess biomass potentials and its energy power embedded into the grid.

Figure 1 .

 1 Figure 1. Schematic view of the original GREECE-OPSPV integrated model framework (Al-Kurdi et al., 2019)

  Figure2. Fiber optical rings[START_REF] Pelleau | SONET network design problems[END_REF] and power network analogy

Figure 3 .

 3 Figure 3. Set of potential sites derived by the GREECE model in Fr. Guiana.

  Figure 4. Pareto solutions for ring and site approach.

Integrated model framework GREECE module Optimization module

		Base layer
		(county, etc.)
	Geographical data layers	Topography Land use Restricted area layers Road, grid, etc.	Remove restricted areas Compute de-spatialized features Aggregate PV power time series		Partition larger polygons Remove smaller polygons and polygons beyond distance threshold
			Set	of candidates
		Resource
		dataset
	Control parameters	Distance thresholds Land surface area Buffers	Parcel size PV existing + new maximal size 35 % of total power demand Existing production + new PV hourly demand Maximize energy generation Minimize total costs
			Subset of	power stations
			Size		Location
			Capacity		& Features

Table 1 .

 1 Nomenclature for the OPSPV problem formulation

	h ∈ H = {1, . . . , 8760}	Hour (per year)
	i ∈ N = {1, . . . , n}	Site index
	k ∈ N = {1, . . . , n}	Ring index
	n	Number of candidate sites derived from GREECE
	x ik , y k	Boolean decision variables
	B k	Boundary set of nearest PV sites corresponding to ring k
	C	Substation maximum hosting capacity (kW)
	M	Big number to enforce Boolean inference
	Eint h	Current hourly production from intermittent energy sources (kWh)
	Ep h	Current hourly production from dispatchable energy sources (kWh)
	Dem h	Estimated global (forecasted) hourly power demand (kWh)
	P nom	Nominal power per unit area (kW/m 2 )
	SA ik	Surface area of new selected PV site belonging to ring k (m 2 )
	Smin, Smaxi	Minimum and maximum area for each candidate parcel (m 2 )
	P Si	Potential Site
	P pv h,i	Estimated hourly production per PV unit (kWh/m 2 )
	dgi	Minimal distance from the grid to PV site centroid (m)
	dij	Distance between site P Si and site P Sj
	Clan	Transmission line unit cost (e/m)
	Ccap k	Capital cost of implementation of new PV power plants (e)
	Cop k	Annual operational cost per PV power plant (e)
	Ccon k	Connection costs for each ring, transmission lines (e)
	Csta k	Capital cost for new substation per ring (e)

Table 3 .

 3 Characteristics of the rings depicted in Figure 5bRing ID Number of sites Installed power (MW) site capacities (MW)

		site 1	6.87
		site 2	20
		site 3	18.52
		Total	45.39
	Ring 1	2	5.56	{4.72, 0.84}
	Ring 2	4	19.3	{5, 5, 4.74, 4.56}
	Ring 3	5	19.99	{5, 5, 5, 2.8, 2.19}
	Total	11	44.85	-


Gaël Loubet

Alassane Sidibe

Alexandru Takacs

Daniela Dragomirescu

This paper presents a cyber-physical system based on a wireless sensor network dedicated to structural health monitoring of reinforced concretes throughout their lifetime. This cyber-physical system is intended to implement a communicating reinforced concrete. Two types of nodes compose this WSN. The sensing node is fully wireless, can measure various parameters (such as temperature, relative humidity, mechanical strain, or resistivity), is battery-free, and is wirelessly and remotely powered and controlled via a radiative electromagnetic power transfer system by the second type of nodes, the communicating node. The communicating node connect the WSN to the digital world.

INTRODUCTION

Wireless Sensor Networks (WSN) can be employed to deploy Cyber-Physical Systems (CPS), both by monitoring and/or controlling the physical world, and by connecting the physical and digital worlds. They are suitable for Structural Health Monitoring (SHM) applications in various fields, such as the civil engineering industry. Currently, one of the main limitations to the deployment of WSN is their energy autonomy. To overcome this bottleneck, ambient energy harvesting and Wireless Power Transfer (WPT) solutions are investigated [START_REF] Peruzzi | A review of energy harvesting techniques for Low Power Wide Area Networks (LPWANs)[END_REF]. In this context, the McBIM project (Material communicating with the Building Information Modelling (BIM)) [START_REF] Derigent | Materials communicating with the BIM: aims and first results of the McBIM project[END_REF], funded by the French National Research Agency (ANR), design and implement a CPS dedicated to the SHM for civil engineering industry. This one is composed of a WSN directly embedded in the reinforced concrete.

In this paper, the design and the implementation of the WSN embedded into the reinforced concrete will be presented. Its energy autonomy, its sensing capabilities, as well as its integration into reinforced concrete will be detailed. _____________

DESIGN AND IMPLEMENTATION OF COMMUNICATING REINFORCED CONCRETE

In order to make communicating the reinforced concretes, the choice has been made to embed a WSN with the ability to sense various parameters, to process, store and share wirelessly the collected data both with the physical and digital worlds. The aim of making it useful throughout its life -expressed in decades-and for each step of its lifecycle, the issues of energy autonomy, reconfigurability, measurement ability and accessibility have been considered.

Architecture of the Cyber-Physical System and its Wireless Sensor Network

The proposed WSN is part of a full CPS and has two kinds of nodes: the communicating nodes (CN) and the sensing nodes (SN) as shown in Figure 1.

The CN are organized in a mesh network, accessible and powered in a conventional way. They will collect, process, store and share together, as well as with the digital world and its virtual models (e.g. a digital twin) via the Internet, the data measured and sent wirelessly by the SN. Moreover, the CN will power the SN wirelessly over several meters.

The SN are organized in a star topology around a CN, inaccessible because randomly deployed and buried into the reinforced concrete, and autonomous in terms of energy by being battery-free and powered wirelessly by the CN. These must measure some parameters to quantify the health of the reinforced concrete and/or its environment, then pre-process and transmit wirelessly the collected data to the CN. The WSN is based on the long-range LoRa ultra-low power wireless communication technology, but experiments with the Bluetooth Low Energy (BLE) are in progress. The power supply of the SN by the CN is achieved by radiative electromagnetic WPT system operating in the far-field. Further details on the design and implementation of the nodes are available in [START_REF] Loubet | Implementation of a Wireless Sensor Network Designed to be Embedded in Reinforced Concrete[END_REF][START_REF] Loubet | LoRaWAN Battery-Free Wireless Sensors Network Designed for Structural Health Monitoring in the Construction Domain[END_REF].

Long-Term Energy Autonomy of the Sensing Nodes

With the aim of making the SN energy autonomous in the long term (that says for decades), and because these are inaccessible once deployed (because buried into reinforced concrete), these are designed to be wirelessly powered in the far-field (up to several meters), battery-free, and able to cold-start.

Wired solutions have been ruled out in favor of wireless ones. The former require manual installation, are expensive and heavy, and are difficult to deploy on a large scale.

Due to slow evolution of monitor parameters (e.g. temperature, humidity, mechanical strain, corrosion rate, etc.), the frequency of measurement can be low (e.g. on an hourly, daily, or weekly base). That is why an energy "harvest and store, then use" strategy has been chosen. Due to the use of a Power Management Unit (PMU) commercially available, it is possible to manage only by hardware the power harvesting and the energy storage and use. Using the cold start capability of this PMU allows to start from an empty energy buffer. The SN can be kept functional for a long time, even if not used and supplied during a long period of time.

The choice of a (super)capacitor has been made for the energy buffer. Regardless the chosen wireless communication technology (LoRa, BLE, etc) and the sensing technology, a low to moderate capacity is needed to supply the SN in the targeted operation mode: "harvest and store, then use", with "use" consisting solely of the collection and transmission of data. Low power losses are sought. Long lifetime with a lot of charge and discharge cycles, as well as full discharge ability are needed for a longterm use. So, the batteries do not meet all these requirements: they have too high capacities, limited lifetime, as well as limited number of charge and discharge cycles and cannot be completely empty. The rechargeable ones require a recharge system, but today only wired or short-range wireless ones are available that forces the operators to know their location to recharge these, which is difficult, time-consuming and expensive.

The proposed choice is to use a far-field radiative electromagnetic WPT system to provide the required energy to the SN. This system is independent of the environment of deployment and its unavailable or fluctuating ambient energy sources [START_REF] Matiko | Review of the application of energy harvesting in buildings[END_REF] and of the targeted structure and application. Moreover, it has a sufficient range of use, and it is scalable in all civil engineering use cases. It provides a generic, long-range and longterm solution. The available near-field electromagnetic WPT solutions are limited in terms of range up to tens of centimeters.

Thus, due to its design based on battery-free SN supplied by a far-field WPT system, the proposed system can be used for decades, especially for SHM applications in civil engineering, whatever the environment and the targeted structure.

The CN are still not energy autonomous. But, it could be possible to use batteries to power these due to their accessibility. Few propositions have already been made to power with photovoltaic cells, both LoRaWAN gateways [START_REF] Peruzzi | A review of energy harvesting techniques for Low Power Wide Area Networks (LPWANs)[END_REF] and proof-of-concept radiative electromagnetic power sources [START_REF] Georgiadis | Solar powered class-E active antenna oscillator for wireless power transmission[END_REF].

Reconfigurability of the Sensing Nodes

The SN can be reconfigured on demand in terms of periodicity of measurement and data wireless transmission. The CN can remotely control and reconfigure the periodicity of functioning of the SN by tuning the WPT system, and thus by controlling the among of power available on each of them. This is a purely hardware-based solution that does not require any modification of the SN. This function can be improved by designing and implementing a beam-forming solution for the source of the WPT system. The global energy efficiency could be enhanced.

Measurement Capabilities of the Sensing Nodes

The SN have been designed in order to provide a generic platform on which various kinds of sensors can be connected. In order to limit the needs in terms of energy and computing capacity, only direct and temporally punctual measurements for nondestructive testing (NDT) methods have been considered. Methods requiring signal processing and/or energy-consuming equipment have been discarded. In the case of the monitoring of reinforced concretes during all their lifespan, several sensors have already been implemented and tested. These allow to get estimations of the temperature, relative humidity, mechanical strain and corrosion rate and thus, allow to monitor the curing process via the maturity method or the exploitation phase thanks to various measurements. By adding an identifier and by storing the data in the material, the traceability could be covered.

TEMPERATURE AND RELATIVE HUMIDITY SENSOR

Temperature and relative humidity sensors have first been implemented. The offthe-shelf Texas Instruments HDC20*0 sensors have been chosen for a temperature range from -40 °C to +125 °C and a relative humidity range from 0 % to 100 %, respectively with an accuracy of ±0.2 °C and 2 %.

Thermodiodes designed and manufactured by the University of Cambridge have been used as temperature sensors [START_REF] De Luca | Experimental, analytical and numerical investigation of non-linearity of SOI diode temperature sensors at extreme temperatures[END_REF]. These are highly linear, accurate and very lowpower. The subsidiary circuits (zero-temperature coefficient current source and scaling buffer circuit) have been designed with components that do not achieve the lowest possible power consumption. The proposed architecture provides a temperature range larger than from -30 °C to +125 °C, with an accuracy lower than ±0.16 °C, measurable on a voltage from 0 V to 1 V.

STRAIN GAUGE

As a simple and low-cost proof-of-concept, a naked linear strain gauge has been implemented and tested in a quarter Wheatstone bridge configuration to monitor the mechanical deformation. To increase the linearity and to decrease the power consumption, the bridge has been powered by a current source, whilst an amplifier and a scaling circuit has been used to get a measurement voltage between 0 V and 1 V. This cheap configuration must be improved to move beyond the proof-of-concept stage and to provide a more energy efficient and accurate system. This can be achieved by designing more precise electronic and sensor: multidirectional and/or in half-or fullbridge configuration and/or with a dedicated mechanical specimen to get a load cell.

RESISTIVITY SENSOR

In order to be able to estimate the corrosion rate, a resistivity sensor designed by civil engineering researchers from the LMDC lab has been implemented, as presented in Fig. 2. This is fully buried in the reinforced concrete sample under test and is based on a direct current driven Wenner probe, with four identical and equally spaced probes. A buffer circuit allows the measurement of a voltage. This first implementation must be improved to move beyond the proof-of-concept stage and to provide a more energy efficient and accurate system. This can be achieved by designing a more precise electronic and by considering other kinds of configuration, potentially with a larger number of probes. The periodicity of measurement has been assumed to be greater than two hours, and as the duration of the measurement is short, the polarization of the concrete has not been considered. However, to avoid this effect, an alternative current source can be employed at the cost of a more complex electronic.

Integration of the sensing nodes into the reinforced concrete

Currently, the choice has been made to cast all the SN into the reinforced concrete. This solution provides the highest security level, but reduces the wireless communicating range and available power, and results in longer (re)charging times. There are two alternatives: casting only the sensor in the core of the concrete and the rest on the surface; or casting the SN in the core of the concrete and its antenna on the surface. For the first, the SN are accessible, thus replaceable and updatable. For the second, the SN are inaccessible and the constraints of the electromagnetic propagation are released. In the two cases, measurements are performed at the most relevant locations, but the use of wires between the different parts can create mechanical weaknesses and access to the surface is also an access for pollutants and contaminants.

EXPERIMENTAL RESULTS

Several qualitative and quantitative tests in several configurations have been carried out, as the one presented in Fig. 2 [START_REF] Loubet | Implementation of a Wireless Sensor Network Designed to be Embedded in Reinforced Concrete[END_REF][START_REF] Loubet | LoRaWAN Battery-Free Wireless Sensors Network Designed for Structural Health Monitoring in the Construction Domain[END_REF]. The SN is embedded in three 15 cm x 15 cm x 15 cm cavities of a reinforced concrete beam. The SN is powered and controlled wirelessly up to several meters by the CN.

According to Fig. 3, the SN need about 1.6 seconds to complete a full process: power-up, initialization, measurement, pre-processing, LoRaWAN transmission, regardless of the sensor and transmission power. The energy required is strongly correlated to the transmission power (about 92 mJ at +4 dBm and 172 mJ at +16 dBm), and little to the sensor (up to +7.5 mJ when using two thermodiodes simultaneously).

Regarding the WPT, the duration of the first charge from an empty energy buffer and of the recharges is function of the electromagnetic energy available at the input of the SN, as shown in Fig. 4. The radiated electromagnetic power is essentially function of the distance with the power source, the CN. In this case: 868 MHz and +33 dBm or 2 W EIRP, according to local regulations [START_REF] Etsi | Electromagnetic compatibility and radio spectrum matters (ERM); short range devices (SRD) intended for operation in the bands 865 MHz to 868 MHz and 915 MHz to 921 MHz; guidelines for the installation and commissioning of radio frequency identification (RFID) equipment at UHF[END_REF]. These times have been characterized between -14 dBm and +15 dBm, and are respectively about 34 h and 36 min for first charge and 1 min for recharge, and about 14 h and 38 min for first change and 30 s for the recharge. Fig. 5 presents a medium-term characterization of a concrete sample embedding a resistivity sensor driven by a SN and located in the neighborhood of another SN equipped with a temperature and relative humidity sensor. Even if this test has been carried out in the air and with a conducted electromagnetic power to be independent of the environment fluctuations, a consistent behavior can be observed during the drying process of the concrete, and a correlation between the environment parameters and the resistivity is highlighted (especially when a fan heater is used near the setup, but also with the day and night cycle or after a rainy period).

IMPROVEMENTS AND PERSPECTIVES

A full proof-of-concept of the targeted CPS has been achieved. Some parameters have still to be improved: the range of use that is limited by the WPT system and the regulations to several meters, currently about 11 meters; the global energy efficiency; the volume and the packaging of the SN.

Both the range of use and the energy efficiency can be increase, especially by improving the energy harvesting system by using a more efficient rectenna, a higher gain antenna, a PMU requiring less power, a low loss energy buffer and by using a power source generating optimized electromagnetic waves.

The energy efficiency can be also improved by reducing the energy needed by optimizing the sensor, the pre-processing and the wireless transmission protocol (for LoRaWAN: the data-rate, the payload).

The antenna is the main dimensioning element and it depends of the targeted performances (highest performances, bigger the antenna) and the environment (e.g. into or on reinforced concretes). A packaging to cast the WSN into reinforced concrete must be designed to protect it without altering its functioning.

CONCLUSION

This paper proposes an innovative SHM system, dedicated to the monitoring of reinforced concretes throughout their life, which can be easily adapted to various applications. This CPS is based on a battery-free WSN buried into reinforced concrete composed of low-cost wireless SN that can drive different types of sensors and are completely autonomous for decades. The Sensor Nodes are powered wirelessly by farfield radiative electromagnetic waves, and based on ultra-low power wireless communication technology. They are remotely reconfigurable in terms of periodicity of use (measurement and wireless communication) through the tuning of the power source used by the WPT system. It has to be noted that this is one of the first implementations of a resistimeter fully embedded into concrete [START_REF] Corva | Miniature resistance measurement device for structural health monitoring of reinforced concrete infrastructure[END_REF], and the first fully wireless and battery-free.

Figure 1 .

 1 Figure1. Architecture of the proposed cyber-physical system and its wireless sensor network.

Figure 2 .

 2 Figure 2. Photographs of: a sensing node with a resistivity sensor (left) and the experimental setup with sensing nodes embedded in a reinforced concrete beam (right).

Figure 3 .

 3 Figure 3. Power consumption of the sensing nodes against the sensor and the transmission power.

Figure 4 .

 4 Figure 4. Durations of the first charge and recharges, and rectifier output voltage against the electromagnetic power applied at the input of the sensing node for a frequency of 868 MHz.

Figure 5 .

 5 Figure 5. Evolutions of the computed resistivity and of the measured output voltage of the resistivity sensor, and of the temperature and relative humidity, for several days of concrete drying.

ACKNOWLEDGMENT

This research has been funded by the French National Research Agency (ANR) under the McBIM project, grant number ANR-17-CE10-0014. The authors would like to thank all their French partners in the McBIM project: CRAN (Nancy), LIB (Dijon), and FINAO SAS / 360SmartConnect (Trans-en-Provence).