Francesco Ferrante
email: francesco.ferrante@gipsa-lab.fr

Chiara Ravazzi
email: chiara.ravazzi@ieiit.cnr.it.

Fabrizio Dabbene
email: fabrizio.dabbene@ieiit.cnr.it

Fabrizio

An LMI Approach for Structured H_∞ State Feedback Control

Keywords: Distributed/decentralized control, Linear Matrix inequalities, Linear systems, Structured control design, State feedback optimal control methods

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

INTRODUCTION

Coordination and control of large-scale systems, composed by interconnected independent subsystems, is a relevant problem in many applications, such as transportation networks [START_REF] Ataşlar | A decentralized control approach for transportation networks[END_REF], flow networks [START_REF] Bauso | The linear saturated decentralized strategy for constrained flow control is asymptotically optimal[END_REF][START_REF] Bauso | Optimization of long-run average-flow cost in networks with time-varying unknown demand[END_REF], communication networks [START_REF] Ephremides | Control and optimization methods in communication network problems[END_REF][START_REF] Moreno Banos | A linear programming approach to large-scale linear optimal control problems[END_REF], and smart-grids [START_REF] Ayar | A distributed control approach for enhancing smart grid transient stability and resilience[END_REF], to mention just a few. Dealing with these large, complex, and distributed systems poses new challenges. In fact, the implementation of centralized controllers is in general not a viable approach for this kind of systems, and controllers must be designed in a decentralized/distributed way in order to be robust and scalable [START_REF] Ikeda | An inclusion principle for dynamic system[END_REF]. Using only information about a subset of components to produce control laws imposes several restrictions to the structure of the overall control schemes, and the complexity of the controller depends both on the number of the states and on the particular architecture [START_REF] Ikeda | An inclusion principle for dynamic system[END_REF]. In [START_REF] Delvenne | The price of distributed design in optimal control[END_REF] the performance degradation using only local information is quantified and compared with the global optimum obtained when the global information of the system is available in synthesis. In [START_REF] Tanaka | The bounded real lemma for internally positive systems and H-infinity structured static state feedback[END_REF] distributed and scalable control methods are proposed for positive systems but the optimality is no longer guaranteed.

The quadratic invariance principle [START_REF] Lessard | Convexity of decentralized controller synthesis[END_REF] provides a powerful tool to identify control problem with structural constraints that are convex. However, also ⋆ The authors' research has been partially supported by the PRIN 2017 Project, Prot. 2017S559BB. in this case, the problem of finding the optimal distributed control can often be intractable.

Besides decentralized/distributed control where the controller is required to have a block-diagonal structure, more complicated constraints arise in sparse optimal control [START_REF] Polyak | Sparse feedback in linear control systems[END_REF], where the constraints are in the form of sparsity requirements for the feedback matrix [START_REF] Lin | Augmented Lagrangian approach to design of structured optimal state feedback gains[END_REF] and the optimization problem becomes nonconvex. A large portion of literature proposes to solve directly the nonconvex problem. For instance, in [START_REF] Lin | Augmented Lagrangian approach to design of structured optimal state feedback gains[END_REF] the authors propose to tackle the optimization via augmented Lagrangian and alternating direction method of multipliers (ADMM), or sequential convex-programming [START_REF] Fardad | On the design of optimal structured and sparse feedback gains via sequential convex programming[END_REF]. Although these methods are very flexible and general, the convergence to the optimal solution can not be guaranteed.

In this paper, we study structured H ∞ optimal control design, namely we aim at designing a state feedback obeying to specified structural constraints, so that the system is stabilized while minimizing the gain of the disturbance-tooutput channel. This problem has attracted a large interest in the past years, and several results are available in the literature. However, it should be remarked that most of the approaches to this topic are restricted to very specific classes of systems and problems. These assumptions allow one to derive closed-form solution of the control law, but substantially restrict their practical applicability. For instance, closed-form expressions have been derived be for linear time-invariant discrete time systems with symmetric and Schur state matrix [START_REF] Lidström | Hinfinity optimal distributed control in discrete time[END_REF][START_REF] Rantzer | Structure preserving control[END_REF], or for internally-positive systems, for which the closed-loop state space matrix needs is Met-zler [START_REF] Tanaka | The bounded real lemma for internally positive systems and H-infinity structured static state feedback[END_REF]. Other approaches, as already discussed, directly tackle the ensuing nonconvex problem [START_REF] Lian | Sparsity-constrained mixed H 2 /H ∞ control[END_REF], or resort to dynamic feedback, see e.g. [START_REF] Langbort | Distributed control design for systems interconnected over an arbitrary graph[END_REF] and references therein. Specifically, in this work we assume the feedback matrix is constrained to belong to a given linear space. We remark that the formulation of the structural constraints in terms of linear space is very general, and, as extensively discussed in our previous work [START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF], it captures important problems such as network decentralized/distributed control, overlapping control, and sparse control as particular cases. Our methodology is mainly based on a linear matrix inequalities (LMI) formulation of the problem. In this setting, the design of stabilizing structured static feedback controllers is recast into the solutions to some linear matrix inequalities with structured decision variables. Those matrix inequalities naturally stem from Lyapunov/dissipation inequalities. One of the challenging features of this problem is that stability conditions generally take the form of bilinear matrix inequalities in the controller gain and the matrix defining the underlying Lyapunov/storage functions. Hence, those cannot be directly employed for controller design. An approach to overcome this problem consists of performing linearizing change of variables preserving the structure of the resulting controller gain. This is done, e.g., in [START_REF] Blanchini | Network-decentralized control strategies for stabilization[END_REF], by considering a specific structure of the underlying Lyapunov function. However, imposting a predetermined structure to Lyapunov functions is in general a major source of conservatism. In our recent paper, [START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF], we show how by relying on an approach that is reminiscent of the S-variable approach [START_REF] Ebihara | Svariable approach to LMI-based robust control[END_REF], one can obtain less conservative conditions for structured feedback control design. The key idea consists of resorting to structured "dilated" stability conditions [START_REF] Pipeleers | Extended LMI characterizations for stability and performance of linear systems[END_REF] to remove the need of imposing a specific structure to the underlying Lyapunov function. The idea of using dilated conditions to release contraints on the Lyapunov matrix when the feedback gain has a certain structure was originally illustrated in [START_REF] De Oliveira | Extended H 2 and H ∞ norm characterizations and controller parametrizations for discrete-time systems[END_REF] in some specific examples for the case of discrete-time systems. The main novelty in [START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF] consists of providing a systematic and computationally affordable approach to design stabilizing controllers fulfilling general geometric constraints. This paper extends the results in [START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF] to H ∞ static feedback control design. The methodology we propose relies on the solution to some linear matrix inequalities coupled to a line search on a scalar parameter and turns out to be less conservative than previous methods presented in literature. Through numerical examples we demonstrate the simplicity of the method and performance of the optimal control law.

Outline

We formally present the problem of structured H ∞ state feedback control in Section 2 and review some preliminary results that are useful for subsequent developments in Section 3. The main theoretical contribution of this paper is presented and discussed in Section 4. Then, the effectiveness of the proposed approach is shown in exam-ples retrieved from literature on Decentralized Control and Overlapping Control (see Section 5). Finally, Section 6 collects some concluding remarks and discussions on future directions of research.

Notation

We denote by R and C the sets of real and complex numbers, respectively. The symbols R n×m and S n p represent the set of n × m real matrices and the set of real n × n symmetric positive definite matrices, respectively. The symbol R n stands for the set of n × n nonsingular real matrices. Given A ∈ R m×n , we denote its transpose by A T , and, when n = m, we define He(A) = A + A T . The notation I denotes the identity matrix whose size is determined from context. For a symmetric matrix A, negative and positive definiteness are denoted, respectively, by A ≺ 0 and A ≻ 0. The symbol ⋆ stands for symmetric block in symmetric matrices. The symbol ⊗ is used for Kronecker product between matrices, respectively. Given a linear time-invariant continuous-time system with state x, input w, and output z of the form:

Σ : ẋ = Ax + Ew z = Cx + Dw. we denote by Σ ∞ the H ∞ norm of the transfer function of Σ from w to z.

PROBLEM STATEMENT

We consider a linear continuous-time plant of the form:

ẋ = Ax + Bu + Ew z = Cx + Dw (1)
where x ∈ R nx is the plant state, u ∈ R nu is the control input, w ∈ R nw is an exogenous perturbation, and z ∈ R nz is a regulated output. We assume that matrices

A ∈ R nx×nx , B ∈ R nx×nu , C ∈ R nz×nx , E ∈ R nx×nw ,
and D ∈ R nz×nw are constant and given. The problem we consider in this paper consists of designing a structured state-feedback control law u = Kx for (1) ensuring closedloop stability and an H ∞ performance from the input w to the output z. More precisely, let S ⊂ R nu×nx be a linear space and

Σ : ẋ = (A + BK)x + Ew z = Cx + Dw. (2)
be the closed-loop system. The problem we solve can be formalized as follows: 2) is asymptotically stable and the following bound holds: Σ ∞ ≤ γ. ⋄

Problem 1. (Structured H ∞ State Feedback). Given γ > 0, design K ∈ S such that the closed-loop system (

PRELIMINARY RESULTS

For a given controller gain, consider the closed-loop system (2). With the objective of streamlining the design of the controller, as in [START_REF] Ebihara | Svariable approach to LMI-based robust control[END_REF], we consider the dual of closed-loop system (2), which reads as follows

Σ D : χ = (A + BK) T χ + C T v f = E T χ + D T v. (3)
that has state χ ∈ R nx , input v ∈ R nz , and output f ∈ R nw . In the result given next, we recall that Σ ∞ = Σ D ∞ . Lemma 1. (Dual Representation). Let Σ and Σ D be defined respectively in (2) and (3). Then, the following identity holds:

Σ ∞ = Σ D ∞ .
Proof. Let G Σ and G ΣD be the transfer functions of Σ and Σ D , respectively. Then, the result follows directly from the following identity:

G ΣD (s) = G Σ (s) T ∀s ∈ C
It is well-know from H ∞ theory (see, e.g., (Boyd et al., 1994, Chapter 2, page 26)) that Σ D ∞ < γ if and only if there exists P ∈ S np p such that the following matrix inequality holds:

  He((A + BK)P ) P C T E ⋆ -γI D ⋆ ⋆ -γI   ≺ 0. ( 4
)
Building upon that, we state a "dilated" version of ( 4), which will be used to formulate the main result of this paper. This condition is based on the general results in [START_REF] Pipeleers | Extended LMI characterizations for stability and performance of linear systems[END_REF][START_REF] Ebihara | Svariable approach to LMI-based robust control[END_REF].

Theorem 1. Let P ∈ S np p and K ∈ R nu×np be given. The following items are equivalent.

(i) Inequality (4) holds;

(ii) There exists X ∈ R np×np and α > 0 such that

He           (A + BK)X α(A + BK)X + P 0 E -X -αX 0 0 CX αCX - γ 2 I D 0 0 0 - γ 2 I           ≺ 0
(5)

Proof. Following the results in [START_REF] Pipeleers | Extended LMI characterizations for stability and performance of linear systems[END_REF], Section 5, Extension IV), it follows that (4) holds if and only if there exists X and α > 0 such that:

He (A + BK)X α(A + BK)X 0 -X -αX 0 CX αCX 0 +   γ -1 EE T P γ -1 ED T ⋆ 0 0 ⋆ ⋆ γ -1 DD T -γI   ≺ 0 (6) (7)
Thus, using the fact that γ > 0, via Schur's complement, it follows that (6) is equivalent to (5). This concludes the proof.

Building on the above result, the following "design oriented" equivalent condition can be obtained. The proof of the result follows directly from Theorem 1, hence it is omitted. Corollary 1. Assume that there exist a nonsingular matrix

X ∈ R np×np , P ∈ S np p , R ∈ R nu×np , and α > 0 such that He           AX + BR α(AX + BR) + P 0 E -X -αX 0 0 CX αCX - γ 2 I D 0 0 0 - γ 2 I           ≺ 0 (8) RX -1 ∈ S. (9)
Then, K = RX -1 solves Problem 1.

MAIN RESULT

We are now in the position of stating the main result of this paper. Theorem 2. Let {S 1 , S 2 , . . . , S k } be a basis of S and

L := [ S 1 | S 2 | . . . | S k ] Define the following structure set Υ := {Q ∈ R np×np : ∃Λ ∈ R k s.t. L(I k ⊗Q) = L(Λ⊗I np )}.
Assume that there exist P ∈ S np p , R ∈ span{S 1 , S 2 , . . . , S k }, α > 0, and X ∈ R np×np such that:

X ∈ Υ, (10a) He 
         AX + BR α(AX + BR) + P 0 E -X -αX 0 0 CX αCX - γ 2 I D 0 0 0 - γ 2 I           ≺ 0. (10b)
Then, X is nonsingular and K = RX -1 solves Problem 1.

Sketch of the proof. Nonsingularity of X follows directly from (10b). Using inequality (10b) and applying Corollary 1, we deduce that K = RX -1 is such that Σ D < γ. The proof can be completed by observing that, as shown in [START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF], proof of Theorem 2), R ∈ span{S 1 , S 2 , . . . , S k } and (10a) imply that RX -1 ∈ S.

APPLICATIONS AND EXAMPLES

In this section, we compare the performance show the effectiveness of the proposed approach in two examples. The first example is retrieved from the literature of Network Decentralized Control. The second example pertains to the so-called Overlapping Control problem.

Network Decentralized Control

We analyze the setting proposed in [START_REF] Blanchini | Network-decentralized control strategies for stabilization[END_REF], in which a collection of N decoupled linear systems is interconnected through a set of local controller nodes, each having access only to local state information. We revisit the example in [START_REF] Blanchini | Network-decentralized control strategies for stabilization[END_REF] where the networked system is composed of N = 5 subsystems and the overall dynamics can be written as in (1) with

A = diag(A 1 , A 2 , . . . , A 5 ), B =      B u -B d 0 0 0 0 0 B u -B d 0 0 -B d 0 0 B d -B u 0 0 0 0 0 B d B u 0 0 0 0 0 B u B d      , E = I
where for all i = 1, 2, . . . , 5

A i = -ξ i β i 0 ξ i -β i 0 0 1 0 , B d = 0 1 0 , B u = 1 0 0
In particular, for all i ∈ {1, 2, 4} we select β i = 0, ξ 1 = 15, ξ 2 = 20, ξ 3 = 16, ξ 4 = 16.7, ξ 5 = 14, β 3 = 12, and β 5 = 22. In addition, we assume C = I and D = 0.

As illustrated in [START_REF] Blanchini | Network-decentralized control strategies for stabilization[END_REF][START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF], in this setting, the proposed control strategy can be performed by considering a standard feedback controller with a block-structured feedback gain. In particular, the feedback gain needs to have the same zero-block structure of B T . It can be shown1 that the space of matrices with the same zero-blocks of B T is a linear space of R 6×15 with dimension equal to 33. We use Theorem 2 to design a controller fulfilling the prescribed structural constraints and ensuring a guaranteed level of H ∞ performance. Let Υ be defined as in Theorem 2. Then, simple manipulations in this case show that any matrix X ∈ Υ has the following structure:

X = diag (X 1 , X 2 , X 3 , X 4 , X 5 ) (12)
where X i ∈ R 3×3 , i = 1, . . . , 5, are nonsingular arbitrary matrices.

Therefore, Theorem 2 can be applied with X structured as in ( 12) and R with the same zero-blocks of B T , that is:

R =        r 1,1 r 1,2 r 1,3 0 0 0 0 0 0 0 0 0 0 0 0 r 2,1 r 2,2
r 2,3 r 2,4 r 2,5 r 2,6 0 0 0 0 0 0 0 0 0 0 0 0 r 3,4 r 3,5 r 3,6 r 3,7 r 3,8 r 3,9 0 0 0 0 0 0 0 0 0 0 0 0 r 4,7 r 4,8 r 4,9 r 4,10 r 4,11 r 4,12 0 0 0 0 0 0 0 0 0 0 0 0 r 5,10 r 5,11 r 5,12 r 5,13 r 5,14 r 5,15 0 0 0 r 6,4 r 6,5 r 6,6 0 0 0 0 0 0 r 6,13 r 6,14 r 6,15

      
We solve Problem 1 for a minimal value of2 γ: inf P,X,R,α,γ γ subject to (10a), (10b), P ≻ 0, α > 0.

Once the structure of the set Υ is made explicit, the above optimization problem can be solved via standard SDP solvers, with the caveat of performing a line search on the positive scalar α. Numerical experiments show that the selection of α is not critical to ensure the feasibility of the problem, but it has an impact on the resulting solution.

To show the influence of the variable α on the design of the controller, in Fig. 1 we report the norm of the gain K against the value of α. With the objective of selecting a good tradeoff between control effort and disturbance rejection, in this example we select α = 0.1554, which leads to γ = 1.7887, and K and P as in (11). From (11a) one can understand the main novelty of the proposed approach: we observe that K has indeed the desired zero-block structure, while the matrix P is completely free and does not present any zero-block structure; see (11b). This degree of freedom in the choice of P renders the proposed design much less conservative with respect to other approaches.

Overlapping Control

We consider the control with overlapping information structure constraints problem in [START_REF] Zecevic | Control of complex systems: Structural constraints and uncertainty[END_REF], Chapter 2) and we revisit (Zecevic and Siljak, 2010, Example 2.16). Let

A = 1 4 0 1 2 2 0 -2 3 , B = 1 0 0 0 0 1 , E = 0 0 1 , D = 0, C = 0 1 0 and S = K ∈ R 2×3 : K = a b 0 0 c d | (a, b, c, d) ∈ R 4 . ( 13
)
setup analyzed in this example is schematically represented in Fig. 2. In particular, our example falls in the so-called Type II problem defined in (Zecevic and Siljak, 2010, Fig. 2.6), since the control inputs do not directly affect the "shared" state x 2 , which however directly influenced by the disturbance w. This is recognized to be a critical situation, considerably more challenging than the Type I overlapping control, in which the shared state is directly influenced by the control inputs. Notice that as a performance output we take y = x 2 , i.e., the shared state. Now we show how our approach hinging upon Theorem 2 can be successfully adopted in this case to solve Problem 1. It can be easily shown that the set S in ( 13) is a linear space of R 2×3 . In particular, as illustrated in [START_REF] Ferrante | On the design of structured stabilizers for LTI systems[END_REF], in this case the set Υ in Theorem 2 can be written as:

Υ = Q ∈ R 3 : Q = q 11 q 12 0 0 q 22 0 0 q 32 q 33 , q ij ∈ R .
As such, Problem 1 can be solved by solving (10b) with R = r 11 r 12 0 0 r 22 r 23 , X =

x 11 x 12 0 0 x 22 0 0 x 23 x 33 Also in this case, with the objective to minimize the effect of the disturbance w onto the performance output y, we recast Problem 1 as the solution to (13). In particular, solving in (13) in Matlab R using the YALMIP package [START_REF] Lofberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] combined with the solver MOSEK (Andersen and [START_REF] Andersen | The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm[END_REF], one gets: To further emphasize the benefit of the proposed design strategy in achieving disturbance rejection despite the structural contraints on the control gain K, we compare the response of the closed-loop system for the optimal gain K here above and for the (nonoptimal) gain K s = -3.831 -5.744 0 0 -5.059 -8.922

proposed in (Ferrante et al., 2020, Section IV.b) to achieve closed-loop asymptotic stability for the same system under the same overlapping control constraint. For the sake of comparison, in Fig. 3 we report the response of the plant output to an energy bounded disturbance and of the control energy for the two controllers, i.e., the "optimal controller" u = Kx and the stabilizing control law u = K s x. Simulations clearly show the benefit of the proposed design strategy. In particular, it is interesting to notice that, despite the constrained structure of the control gain, the proposed optimal design leads to improved disturbance rejection and this with reduced control effort with respect to the nonoptimal feedback law u = K s x.

P =                    
                    (

CONCLUSIONS

This paper proposes a very general approach for the design of state feedback H ∞ controllers obeying to structural constraints. The constraints are expressed in terms of LMIs depending on a design matrix R which is constrained to belong to a pre-defined structure, while the Lyapunov matrix P is not forced to obey any structural constraint but it is instead left completely free. This approach is suitable for a large range of possible applications, as also shown by our numerical simulations. We are now working into extending the proposed structured design approach for cases even more realistic, involving for instance nonlinear terms and/or different costs, as H 2 or guaranteed-cost design.

Fig. 1 .

 1 Fig. 1. Variation of γ (left) and K (right) versus α. The red asterisk denotes the solution selected in this example.

Fig. 2 .

 2 Fig. 2. Topology of the overlapping control problem in Example 5.2.

Fig. 3 .

 3 Fig.3. Plant response (left) and associated control energy (right) for optimal and nonoptimal feedback.

Code available at https://github.com/f-ferrante/IFAC2020_ FerranteDabbeneRavazzi

To avoid the occurrence of overly large control gains, a specific lower bound on (10b) is considered in the solution to the optimization problem.


Keywords: Observer, reduced order, parameter estimation, tunneling current, STM application

In the context of experimental modeling for a Scanning Tunneling Microscope (STM) prototype, a reduced order observer is proposed for parameter adaptation in the tunneling current model. The underlying dynamical model is first recalled, and the proposed observer is then presented, and discussed. The approach is finally illustrated with real data taken from the prototype.

Towards observer-based tunneling current calibration in an experimental STM device.

Gildas Besançon * Alina Voda * Andrei Popescu *

INTRODUCTION

Nanosciences are a source of challenging problems in estimation and control (Voda, 2010;Clévy et al., 2011;[START_REF] Eleftheriou | Control Techno. for Emerging Micro & Nanoscale Systems[END_REF][START_REF] Fleming | Design, Modeling and Control of Nanopositioning Systems[END_REF]. In the specific area of Scanning Probe Microscopy, the socalled Scanning Tunneling Microscope (STM) is the one which first gave birth to the whole field, in the early eighties [START_REF] Binnig | Scanning tunneling microscopy[END_REF]. As meant by its name, it is based on the tunneling current phenomenon, and is the first instrument that allowed surface imaging at atomic resolution.

Even though now pretty well known, this approach still presents challenges in control, as addressed in [START_REF] Bonnail | Variable structure control of a piezoelectric actuator for a scanning tunneling microscope[END_REF] or some former studies of ours [START_REF] Ahmad | Controller design for a closed-loop scanning tunneling microscope[END_REF][START_REF] Ahmad | Robust digital control approach for high performance tunneling current measurement system[END_REF][START_REF] Ryba | Modelling and control of 3d stm-like scanning device with application to surface reconstruction[END_REF] for instance, up to more recent works of [START_REF] Tajaddodianfar | Stability analysis of a scanning tunneling microscope control system[END_REF][START_REF] Tajaddodianfar | On the effect of local barrier height in scanning tunneling microscopy: Measurement methods and control implications[END_REF]. In particular, it requires specific calibrations for the device to be properly operated, and in the present paper, a simple observer approach is proposed to identify a parameter that enters into the dynamical description of tunneling effect, thus making it instrumental in the use of an STM for surface imaging (see e.g. [START_REF] Besançon | Experimental nanosurface imaging with a scanning-tunneling-like device under 3D control[END_REF][START_REF] Popescu | Control-observer technique for surface imaging with an experimental platform of scanning-tunnelingmicroscope type[END_REF]). This work is based on an experimental device which was developed in Gipsa-lab for the purpose of such studies [START_REF] Blanvillain | Pullin control during nanometric positioning by near field position sensing[END_REF][START_REF] Blanvillain | Subnanometer positioning and drift compensation with tunneling current[END_REF].

The formal statement of the estimation problem under study is first given in section 2, together with the presentation of the underlying device. The observer solution, finally given under the form of a reduced order one, is then discussed in section 3, and related application results based on real data from our experimental device are presented in section 4. Some conclusions and considerations for future works end the paper in section 5.

1 Institute of Engineering Univ. Grenoble Alpes

PROBLEM STATEMENT

STM device under consideration

Standard STM operation (Julian [START_REF] Chen | Introduction to Scanning Tunneling Microscopy[END_REF]) in short means approaching a tip to the surface under study, close enough (less than 1nm) so that an applied voltage between both (bias voltage) results in a (quantum) current effect, called tunneling current. The control task is then to regulate this distance.

In the considered Gipsa STM device (shown in figure 1), the vertical motion is driven by a piezo actuator, fed by a voltage generated numerically (via a PC) and appropriately amplified (via a Voltage Amplifier), while the obtained current is in its turn amplified and measured (via a Current Sensor), and brought back to the numerical equipment (PC). For a mathematical description of the whole dynamical behaviour, let us recall modeling basics for each element (see again figure 2 for notations).

First, the piezoactuator dynamics can typically be represented by a second order model (eg as in [START_REF] Ryba | Modelling and control of 3d stm-like scanning device with application to surface reconstruction[END_REF]), that is:

Z p (s) = G p w 2 p s 2 + 2ζ p w p s + w 2 p U a (s) (1)
where s stands for Laplace variable, Z p (s) and U a (s) Laplace transforms of piezo displacement and its applied (amplified) voltage, respectively, and G p , w p , ζ p the piezo DC gain, natural pulsation, and damping coefficient respectively.

Similarly, the dynamics of both voltage amplifier and current sensor can be captured by simple 1st order models, that is:

U a (s) = G v w v s + w v U (s) and U c (s) = G i w i s + w i I t (s) (2)
respectively, where U, U c and I t refer to Laplace transforms of the applied voltage, the voltage read by the current sensor and the tunneling current respectively, while G v , G i denote voltage and current gains respectively, and w v , w i the corresponding bandwidths.

In practice, the piezo actuator bandwidth in our device being about ten times larger than those of the other elements, its dynamics will be neglected in the present study. This means that piezo model (1) will be restricted to its DC gain G p . Finally, the tunneling current can be usually expressed according to the tip motion z p with respect to the sample position z s , and initial tip-sample surface d 0 , as depicted by figure 3. This results in the following expression:

i t = gV b exp(-k(d 0 + z p -z s ))
(3) where g is a constant depending on the materials, V b is the applied voltage, and k the typical exponential rate of the current variation with respect to tip-sample distance.

This means that we can finally end up with a state space representation of the following form:

ẋ1 (t) = -w i x 1 (t) + G i w i gV b exp(-k(d 0 + x 2 (t) -z s (t))) ẋ2 (t) = -w v x 2 (t) + G p G v w v u(t) y(t) = x 1 (t) (4
) where x 1 stands for u c , x 2 for z p and u for the applied voltage.

Estimation problem under consideration

Typically in the above model, z s will vary in some unknown way during a scanning operation mode, and from the knowledge of everything else, its variation can be recovered, giving an image of the surface (Julian [START_REF] Chen | Introduction to Scanning Tunneling Microscopy[END_REF]). In the present paper instead, we will focus on the characterization of the tunneling parameters, and in particular coefficient g, which may be inaccurately known in practice. This step thus comes prior to any scanning, and can be also an opportunity to calibrate the operation with respect to values of d 0 and initial sample level z s for instance. In other words, we consider here the problem of estimating parameter gV b exp(-k(d 0 + z s )) when z s is supposed to be constant (no scanning). Of course if d 0 and z s are known, then g can be obtained, and conversely, whenever gV b is known, z s + d 0 can be estimated. Notice that those parameters cannot be separately identified.

Setting θ := gV b exp(-k(d 0 + z s )) then, we are brought to an observer problem for the following system:

ẋ1 (t) = -w i x 1 (t) + θG i w i exp(-kx 2 (t)) ẋ2 (t) = -w v x 2 (t) + G p G v w v u(t) y(t) = x 1 (t) (5)
where y denotes the available measurement, and u the known applied voltage, while θ (and state variable x 2 ) are to be estimated.

OBSERVER SOLUTION

High gain approach

It is clear that model ( 5) can be extended with unknown parameter θ as an additional state variable in the following way:

ẋ1 (t) = -w i x 1 (t) + θG i w i exp(-kx 2 (t)) ẋ2 (t) = -w v x 2 (t) + G p G v w v u(t) θ = 0 y(t) = x 1 (t) (6)
and that an observer can be built for this extended model. So called uniform observability property [START_REF] Gauthier | Observability for any u(t) of a class of nonlinear systems[END_REF] can indeed here be checked, and a diffeomorphism can be found to turn the model into an appropriate form for a well-known high gain observer design [START_REF] Gauthier | A simple observer for nonlinear systems-applications to bioreactors[END_REF]:

Proposition 1. The system below:

ξ = Aξ -Λ(λ)K(Cξ -y) + ϕ(ξ, u) (7)   x1 x2 θ   =      ξ 1 ξ 3 kw v ξ 2 ξ 2 G i w i exp( ξ 3 w v ξ 2 )      (8) with A = 0 1 0 0 0 1 0 0 0 , C = (1 0 0), Λ(λ) =   λ 0 0 0 λ 2 0 0 0 λ 3   and ϕ(ξ, u) =    -w i ξ 1 -kG p G v w v ξ 2 u -w v ξ 3 + kG p G v w 2 v ξ 2 u -kG p G v w v ξ 3 u + ξ 2 3 ξ 2   
provides an observer for system (6), in the sense that xi (t) -x i (t) for i = 1, 2, and θ(t) -θ, can asymptotically decay to zero for any initial condition, whenever K is chosen such that A -KC is Hurwitz, ϕ is saturated according to some bounded domain of x 1 , x 2 , θ, and λ is large enough.

The formal verification follows classical lines of high gain observers, and is omitted here, since it will not be the approach finally chosen. The usual advantage indeed in such an approach is the guarantee of convergence even in the presence of nonlinearities, as well as the single tuning coefficient λ (high gain).

But a classical counterpart is the possible amplification of the measurement noise effect (see [START_REF] Astolfi | Sensitivity to high-frequency measurement noise of nonlinear high-gain observers[END_REF] for instance), which is quite significant in our application, and makes this high gain technique indeed too sensitive to the noise. One could think of some possible improvements (e.g. as in [START_REF] Ahrens | High-gain observers in the presence of measurement noise: A switched-gain approach[END_REF][START_REF] Boizot | An adaptive high-gain observer for nonlinear systems[END_REF], but we instead propose here to consider a more simple approach, presented in next section.

Reduced order approach

As an alternative approach, owing to the property that trajectories in x 2 all converge to each other, we can propose a more simple reduced order observer (e.g. in the spirit of [START_REF] Besançon | Remarks on nonlinear adaptive observer design[END_REF]) as follows: Proposition 2. For system (5) with u such that x 1 and x 2 remain within positive bounded intervals, the system below provides asymptotic estimates x2 and θ for x 2 and θ respectively, for any γ > 0:

ẋ2 (t) = -w v x2 (t) + G p G v w v u(t) ζ(t) = -γG i w i exp(-kx 2 (t))[ ζ(t) + γy(t)] + γw i y(t) θ(t) = ζ(t) + γy(t) (9)
Proof. The verification goes as follows: first set ζ := θ-γx 1 and notice that the parameter estimation error θ -θ then coincides with e ζ := ζ -ζ. This means that we are brought to the study of the latter. Then, one can check that:

ėζ = -γG i w i exp(-kx 2 )e ζ -γG i w i [exp(-kx 2 ) -exp(-kx 2 )] × [ ζ + γy] =: -γG i w i exp(-kx 2 )e ζ + δ
(10) where δ depends on the difference between x2 and x 2 , as well as ζ and y. It is clear, from the system and observer equations, that e 2 := x2 -x 2 satisfies: ė2 = -w v e 2 meaning that x2 exponentially approaches x 2 . In addition, from equations ( 9) and boundness of y, it results that ζ is also bounded. This means that δ is exponentially vanishing. Now since x 2 remains bounded, e ζ is clearly ISS w.r.t. δ, and since δ goes to zero, so does e ζ (in addition, it is also clear that one can tune the rate of decay via the choice of parameter γ).

Beyond its simplicity (second order only, and still with a single tuning parameter, γ), this approach also ensures some properties with respect to the noise effect (similarly to the adaptive observer approach of [START_REF] Zhang | Adaptive observer for mimo linear time varying systems[END_REF] for instance): in short, the error is guaranteed to converge in mean value, and the tuning parameter γ can even be chosen according to a filtering purpose. If indeed the measurement takes the form:

y = x 1 + ν (11)
where ν stands for some measurement noise assumed to be a gaussian white noise, the estimation error e ζ with the above notations, now satisfies:

ėζ = -γG i w i exp(-kx 2 )e ζ -γ 2 G i w i exp(-kx 2 )ν + γw i ν -γG i w i [exp(-kx 2 ) -exp(-kx 2 )] × [ ζ + γx 1 + γν] =: -γG i w i exp(-kx 2 )e ζ -γ 2 G i w i exp(-kx 2 )ν + γw i ν +δ 1 (12)
where δ 1 is still vanishing. It is here clear that x 2 and x2 are independent of ν, so that the mean value of e ζ satisfies:

dē ζ dt = -γG i w i exp(-kx 2 )ē ζ + δ1 (13)
where . refers to the mean. From this, and the fact that δ 1 vanishes, we get that ēζ also asymptotically decays to zero.

In a similar way, it can be checked that the variance of e ζ asymptotically behaves as the solution

Q of: Q = -2γG i w i exp(-kx 2 )Q + γ 2 w 2 i (γG i exp(-kx 2 ) + 1) 2 V (14
) where V stands for the noise intensity. From this, the smaller γ is, the smaller this variance is in turn.

ESTIMATION RESULTS WITH REAL DATA

The above approach has been tested with data recorded from our Gipsa-lab STM prototype [START_REF] Blanvillain | Subnanometer positioning and drift compensation with tunneling current[END_REF] The corresponding numerical values for model ( 5) are summarized in table 1 below.

Gv 15 V /V fv 4 kHz G i 1e9 V /A f i 13 kHz k 16.5nm -1 Gp 1.2 nm/V
Table 1. Numerical data For this system, measurements are obtained when operating the STM in a fixed position, under simple PI regulation of the distance between tip and sample (via output direct regulation).

The corresponding output measurement and available control input are presented in figures 5 and 6 respectively, over 5 seconds (under a setpoint control during 3 seconds first, then followed by a square wave tracking configuration). The estimation for coefficient gV b with the above reduced order observer is then tested with those data, and considering some rough estimate for d 0 + z s here set to 52.13 nm.

Noting yet, from the control profile, that after about 1.5 second, a drift-like disturbance appears (as often with piezoactuation in such devices [START_REF] Rakotondrabe | Complete open loop control of hysteretic, creeped, and oscillating piezoelectric cantilevers[END_REF]), the estimation of gV b is here considered on earlier times only, say up to time 1.4 second. This means that the convergence should be fast enough (γ large enough), but with enough noise attenuation (γ not too large).

A first estimation result is displayed on fig. 7, obtained with γ = 1e10: the convergence can be seen to be here very fast (notice that the parameter estimate is initialized at 0), but the estimate is very noisy (it seems to reach a value of about 1e -3, but with variations up to 7 times this value). A second result with a lower γ (1e8) is then presented in figure 8, illustrating the noise effect reduction (as expected from the discussion after proposition 2), and emphasizing a final value around 0.0011. This is confirmed by the superposed result of an additionally filtered version of this estimate, displayed in dahsed line on the same figure (just using a low-pass filter). Notice then that the estimate of gV b can be frozen after time t = 1.4 second for instance, and going on with the same observer allows now to estimate the drift, in terms of tip-sample distance.

The corresponding estimation result is shown in figure 9. Notice finally that combining this drift estimate with the frozen estimate of gV b and that of x 2 , the actual tunneling current can be estimated in its turn, and the corresponding result is presented in figure 10. An additionally filtered estimate is displayed as well, emphasizing the current regulation around a value of about 0.3nA. From the estimated gV b and the available output measurement, this means a distance of about 0.92 nm.

CONCLUSIONS

In this paper an observer approach has been discussed for a problem of parameter estimation in tunneling current modeling. The approach has been formally presented, and illsutrated on the basis of real data taken from a prototype experiment. The focus has been made on a single parameter (or combination of parameters), which can be either related to tunneling current characteristics, or to tip-sample distance.

The extension to the full parametric characterization of the tunneling current will be part of future developments, as well as its possible use in control and imaging purposes.

Fig. 1 .

 1 Fig. 1. STM prototype.This operation is summarized by the schematic view of figure 2 below. Full technical details on this instrumentation can be found in[START_REF] Popescu | Control-observer technique for surface imaging with an experimental platform of scanning-tunnelingmicroscope type[END_REF] for instance.

  Fig. 3. Vertical Z-motion.

TIPFig. 4 .

 4 Fig. 4. Real tip and its reflection.

Fig. 5 .

 5 Fig. 5. Measured output (as delivered by current sensor, under closed loop operation).

Fig. 6 .

 6 Fig. 6. Control input (as fed to the voltage amplifier, during closed loop operation).

Fig. 7 .

 7 Fig. 7. Estimation result for gV b with γ = 1e10.

Fig. 8 .

 8 Fig. 8. Estimation result for gV b with γ = 1e8 (and filtered version in dashed red line).

Fig. 9 .

 9 Fig. 9. Estimation result for distance drift, started after time 1.4 sec (and based on the estimate for gV b obtained before time 1.4 sec).

Fig. 10 .

 10 Fig. 10. Tunneling current estimate (reconstructed from gV b and drift estimates) -top: direct estimate, bottom: filtered estimate.

  (see tip picture in fig. 4 hereafter).

ACKNOWLEDGEMENTS

This work has been partly supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir.


A Aguilera-Gonzalez
email: a.aguilera-gonzalez@estia.fr

R Lopez-Rodriguez

I Vechiu

S Bacha

Model Predictive Control for Optimal Energy Management of an Island Wind Storage Hybrid Power Plant

Keywords: Modeling and simulation of power systems, Control of renewable energy resources, Model predictive control, Smart grids

This paper presents a Model Predictive Control-based Energy Management System for compliance with the day-ahead power dispatching plan of a hybrid power plant connected to the Guadeloupe Island electrical grid. The hybrid power plant combines a wind farm and a Li-ION battery energy storage system. The proposed EMS handles several operation rules, in order to solve the optimization problem while considering the production forecasting data as well as the battery lifespan. A simulation study is implemented via PowerFactory and Matlab.

INTRODUCTION

Dependence on fossil fuels imports is a major problem in island territories. Renewable energy resources such as sun and wind are playing a key role in terms of reduction of petroleum-based fuels reliance as they are free and abundant in those overseas locations. However, the fluctuations caused by their random behavior render the scheduling more difficult while increasing the system operational costs [START_REF] Jin | Dynamic modeling and design of a hybrid compressed air energy storage and wind turbine system for wind power fluctuation reduction[END_REF].

Due to the variability of wind speed, wind turbines (WT) have intermittent power output and, consequently, the difficulties grow to meet day-ahead generation schedules. These issues are amplified in the context of island grids as the French Archipelago of Guadeloupe, which should be managed to enhance grid efficiency without affecting stability and energy quality [START_REF] Rodrígues | Grid code reinforcements for deeper renewable generation in insular energy systems[END_REF].

Hybrid Power Plants (HPP) provide a solution to these issues by associating renewable energy sources (RES) with other technologies of generation or storage into one smart generation facility. The smartness of the hybrid systems lies in its Energy Management System (EMS), which allows power dispatching according to a generation plan, as well as the supply of ancillary services like frequency and voltage regulation.

Several EMSs based on control methods and data acquisition have been developed to operate those HPP [START_REF] Luna | Mixed-integer-linear-programming-based ems for hybrid PV-wind-battery microgrids: Modeling, design, and experimental verification[END_REF][START_REF] Bukar | A review on stand-alone photovoltaic-wind energy system with fuel cell: System optimization and energy management strategy[END_REF]. More recently, Model Predictive Control (MPC) has received increasing attention because it can incorporate both forecasts and updated information to decide the future behavior of the system while handling constraints efficiently. MPC integrates the RES production forecasts in the optimization problem, either in terms of the expected generation or in terms of the primary source, (wind speed, solar irradiance, etc) [START_REF] Abdeltawab | Market-oriented energy management of a hybrid wind-battery energy storage system via mpc with constraint optimizer[END_REF][START_REF] Zhang | Uncertainty-resistant stochastic mpc approach for optimal operation of chp microgrid[END_REF][START_REF] Chen | Quantifying cumulative effects of stochastic forecast errors of renewable energy generation on energy storage soc and application of hybrid-mpc approach to microgrid[END_REF]. However, there is still a large number of difficulties with regard to effective load management, economic dispatch, and especially its high dependency on the accuracy of expected forecast data.

In this paper, an EMS is proposed for an existing windbattery storage HPP injecting power into the Guadeloupe island grid. Predictive control and optimization are combined to guarantee compliance with the day-ahead power dispatching plan while managing the storage system charge/discharge cycles efficiently, in order to accomplish the day-ahead power dispatch plan while considering the battery lifespan. To test the strategy, realistic simulation scenarios were considered with different forecasting errors to prove the ability of the proposed method to compensate them. The complete system was implemented in PowerFactory while the control algorithm is executed in MATLAB. This paper is organized as follows: Section 2 presents a description of the Grid-connected hybrid power plant. Section 3 introduces the HPP PowerFactory model. In Section 4 the MPC strategy proposed is developed and finally in Section 5 the results are shown and discussed.

GRID-CONNECTED HYBRID POWER PLANT

The Guadeloupe archipelago relies heavily on thermal nonrenewable sources (coal, diesel and combustion turbines) for electricity generation. In 2018, 71, 8% of the electricity consumed was provided by the combustion of fossil fuels [START_REF] Belfort | Les chiffres clés de l'énergie en guadeloupe. bilan 2018. Report, Observatoire régional de l'énergie et du climat de la Guadeloupe. Comité de l'observatoire régional de l'énergie et du climat[END_REF]. The electrical grid of the Guadeloupe Islands has a 63kV high voltage transmission system, composed of two big loops of overhead lines (see Fig. 1). The loop covering Basse-Terre is interconnected at the Jarry Nord/Jarry Sud substation with Grande-Terre's loop, which is composed of 14 nodes that correspond to HTB/HTA substations each comprising two 63/21 kV step down transformers connecting loads and, in most cases, also reactive compensation stations. This paper proposes a study to assess and validate the operation of a new hybrid wind/storage power plant at a particular point of the Guadeloupe electric grid (Sainte-Rose).The PowerFactory model of the Guadeloupe island electrical system shown in Fig. 2 was implemented from the data presented by [START_REF] Marin | Study of wind power integration in the island grids[END_REF]. The shaded area corresponds to the point of common coupling (PCC) located within the node Sainte-Rose, where the HPP is connected.

Sainte-Rose

Bouillante

SYSTEM MODELING

In this paper, a PowerFactory/Matlab co-simulation study was carried out in order to assess the impact of a windstorage hybrid power plant, injecting power into one of the HTA busbars of the Sainte-Rose substation. The grid model was validated through simulation with respect to available data.

Wind generation system modeling

The wind generation system is depicted in Fig. 3. This system is composed of four 2MW induction generators containing the model of a DFIG (Doubly Fed Induction Generators) wind turbines and converters. The transformers, busbars and underground lines used to rely the wind turbines to the PCC are also part of the system. The WT model is based on the steady-state characteristics of the turbine presented by [START_REF] Rueda | Implementation of Simplified Models of DFIG-Based Wind Turbines for RMS-Type Simulation in DIgSILENT PowerFactory[END_REF]. Such model allows calculate the power extracted from the WT and considers electrical, mechanical, aerodynamic and control aspects.

Energy storage system modeling

The battery energy storage system (BESS) is based on the commercial storage solution Intensium Max 20M of Safe. For the project, four 1 MW 580 kWh Li-ion storage were considered. The modified Shepherd's model proposed by [START_REF] Tremblay | Experimental validation of a battery dynamic model for ev applications[END_REF] was employed to describe the output voltage of a cell as a function of its capacity in ampere-hours (Ah). The battery cell terminal charge and discharge voltages are calculated as follows:

V ch = E 0 + R • i * -K Q m i t -0.1Q m i * -K Q m Q m -i t i t + Θ V dch = E 0 + R • i * -K Q m Q m -i t i * -K Q m Q m -i t i t + Θ (1)
where

Θ = Ā • e -B•i t Qmax
is the term representing the exponential voltage. E 0 represents the battery constant voltage (V), R is the battery internal resistant (Ω), and K is the polarization constant (V/Ah). Meanwhile Q m represents the maximum battery cell capacity (Ah), i * is the filtered current, i t is the dynamic current, Ā is the exponential zone amplitude (V), and B is the exponential zone time constant inverse (Ah -1 ). In this work, to determine the parameters E 0 , R, Ā, B and K, an identification method from the manufacturing discharge curve was employed. This method applies optimization techniques to find the values of the X vector, subject to X min and X max as follows:

X = E 0 Ā B R K X min = E min 0 Āmin Bmin R min K min X max = E max 0 Āmax Bmax R max K max (2)
to minimize the following cost function:

f obj (X) = i=n i=1 Q=Qnom Q=0 (V mes (Q, I i ) -V (Q, I i )) 2 (3)
where I i = I 1 , ..., I n represents the discharge currents, which allows that the quadratic error among the manufacturer's voltage V mes and the model estimations V , is minimized. For that, the amperes-hour delivered Q and the discharge current I are varied.

The four storage units comprising the BESS implemented in PowerFactory is shown in Fig. 4. It was also added the respective converters, two 2-winding transformers, one 3winding, three busbars, and three underground lines. 2018). For this, a current input was imposed on both system representations (Matlab and PowerFactory) to observe the dynamic of battery states during the charging and discharging process.

DC BUS 4 LV BUS 4 DC BUS 3 LV BUS 3 DC BUS 2 LV BUS 2 DC BUS 1 LV BUS 1 PCC LOCAL LI-ION MV BUS 3-4 MV BUS 2 MV BUS 1 V LI-ION 4 CV 4 V LI-ION 3 CV 3 2-Winding.. transfo c.. V LI-ION 2 CV 2 2-Winding.. transfo c.. V LI-ION 1 CV 1 LINE 3-4 NA2YSY 1x240rm 12/20kV it LINE 2 NA2YSY 1x240rm 12/20kV it LINE 1 NA2YSY 1x240rm 12/20kV it 3-Winding.. 3-Winding..

MPC CONTROLLER

The EMS proposed warrants the continuous power injection from the HPP into the PCC, according to a commitment production profile determined from wind forecast data. For this, a control-oriented linear model of the BESS is required by the MPC controller to generate optimal control actions. The control scheme is presented in Fig. 5. As can be seen, the power reference profile P BESS ref is obtained from the commitment generation schedule P SCHED and the output power of the wind energy system P W ECS , as:

P BESS ref = P SCHED -P W ECS (4)
where the signal P SCHED represents the reference for the total power injected into the main grid. P IN J is calculated as:

P IN J = P W ECS + P BESS (5)

Storage System modeling for MPC

The proposed EMS uses predictive control to optimize the use of the storage system and thereby, the power injection into the utility grid. The prediction is based on the linear model of the storage system. With this purpose, the follow linear time-invariant (LTI) form is used:

x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (6)
where k represents the sampling time, x(k) ∈ R n represents the system states, u(k) ∈ R nu the decision variable or input, and y(k) ∈ R nr , the output controlled. Based on the discrete-time state-space representation obtained from Eq. ( 6), the BEES model is given by the next equations set:

i t (k + 1) = i t (k) + t i BESS (k) i * (k + 1) = (1 -α)i * (k) + α i BESS (k) (7) SOC(k + 1) = 100 1 - i t (k) Q - t i BESS (k) Q
where i t is the dynamic current, i * is the filtered current, SOC is the state-of-charge and Q is the battery capacity.

i BESS represents the battery current as the input system and α is the mitigating factor defined by:

α = t (τ + t ) (8)
t represents the simulation step time and τ is the filter time constant. The model output vector containing the variables to optimize is given by:

y(k) = [i * (k) SOC(k)] T (9)
From the measured voltage v BESS coming for the BESS (see Eq. ( 1)) and the filtered current i * given by the model, the battery power can be computed as:

P BESS = v BESS • i * ( 10
)
where P BESS is negative in case of charge (if v BESS = V ch ) and positive in case of discharge (if v BESS = V dch ).

MPC Algorithm

The EMS proposed combines an MPC controller with an optimization strategy based on the representation given by Eq (6). A multi-step ahead map is then required to apply the predictive methodology, which uses a state trajectory vector containing estimates for the next N p sampling instants (i.e. the length of the optimization window is of N p samples). This vector is described by:

x(k + 1) =   

x(k + 1) . . .

x(k + N p )    (11)
In the same way, the N c elements conforming the future control sequence, (N c is called the control horizon) are contained in the following vector:

ũ(k + 1) =     u(k) . . . u(k + N c -2) u(k + N c -1)     (12)
where x(k|ũ(k)) denotes the fact that the state trajectory x(k) calculation takes place at the current time step k and results from the control actions sequence ũ(k) [START_REF] Alamir | A Pragmatic Story of Model Predictive Control: Self-Contained Algorithms and Case-Studies[END_REF]. Applying the one-step-ahead prediction map, the instant k + 2 is calculated as:

x(k + 2) = Ax(k + 1) + Bu(k + 1) = A[Ax(k) + Bu(k)] + Bu(k + 1) = A 2 x(k) + ABu(k) + Bu(k + 1) (13)
For simplicity N p = N c it is considered. Then, rewritten Eq. ( 13) in a more generally form, for any i ∈ 1, ..., N p , the state-space representation is given by:

x(k + i) = A i x(k) + [A i-1 B, ..., AB, B]     u(k) . . . u(k + i -2) u(k + i -1)     (14)
Then, it is possible to define a matrix Π (nu,Np) i allowing the selection of the i th vector of dimension n u from the N p elements of vector u (Alamir, 2013):

    u(k) . . . u(k + N p -2) u(k + N p -1)     =       Π (nu,Np) 1 . . . Π (nu,Np) i-1 Π (nu,Np) i       (15)
where each element of the Π (nu,Np) i matrix is given by:

Π (nu,Np) i = I nui O nui×(Np-inu) (16)
so, the Eq. ( 14) can then be rewritten as:

x(k + i) = A i x(k) + [A i-1 B, ..., AB, B] • Π (nu,Np) i ũ(k) (17)
and the prediction of future states can be written in a compact form as:

x(k + i) = Φ i x(k) + Ψ i • Π (nu,Np) i ũ(k) ( 18
)
where

Φ i = A i and Ψ i = [A i-1 B, ..., AB, B] • Π (nu,Np) i
are constant matrices depending on state and input matrices, A and B respectively.

Power dispatch optimization

The aim of the EMS is ensuring that the power supplied to the main grid respect the operating rules while maximizing the plant's income. The disrespect of the injection band leads to penalties reducing the power plant profits. The power injections with excursions of 60 consecutive seconds outside the limits are penalized with non-payment of the power supplied to the grid for the next 10 minutes. The plant revenue (PR) is a key indicator performance (KPI), which was established for monitoring the HPP profits. During the operation time, it can be calculated considering the energy selling price (SP in /c per kW h) as:

P R = sim.time t=0 P IN J (t) • SP (t) • ∂(t) (19)
where ∂(t) = 0 when a penalty condition is active 1 in other cases

The optimization objectives and the constraints to which the cost function is subject are explained in the following paragraphs:

Objectives

• Power injection band respect: A tolerance region (injection band) is established on the basis of the injection schedule. This objective consists in minimizing the difference among P IN J and P SCHED , which is equivalent to attracting P IN J towards the center of the tolerance region (see Fig. 6).

Fig. 6. Power injection band.

The upper and lower bounds are given by:

P uplim = P SCHED + 0.25 • P M AX P lowlim = P SCHED -0.25 • P M AX ( 20
)
were P M AX is the plant generation capacity (8M W ).

• State-of-charge maximization: This objective consists in minimizing the difference among the current SOC and reference established at 90%.

Constraints

• Maximal power injection: When the wind-generated power is close to maximal of the WP capacity and the BESS is being discharged, P max IN J is fixed as follow:

P max IN J = P SCHED + tol • P M AX (21)
conversely, its minimum value is zero when the WT is not generating and the BESS is not being discharged.

In order to avoid penalties, rather than setting constraints on the injection band limits, the constraint P max IN J is placed on the band ceiling as:

P IN J ≤ P max IN J (22)
In this way, instead of allowing injections greater than the upper limit, extra available power can be used to charge the BESS. This constraint was defined through the current i * , which is a controlled output.

• Rate of change of power injected: the rate of change limits of the power injected were defined by the objectives. This can be put in inequality form, giving:

P M AX 300 M W/s ≤ dP IN J dt ≤ P M AX 30 M W/s (23) -P M AX 600 M W/s ≤ dP IN J dt ≤ -P M AX 60 M W/s
which are constraints for upward and downward steps of dP IN J /dt. In other words, these constraints limit the rate of change of P IN J so as to avoid big changes in the power transferred towards the grid.

• State-of-charge: The BESS must be operated respecting the recommendations of the manufacturer in terms of depth-of-discharge and charging rates as to avoid its replacement before 15 years. To achieve this, a SOC operation region is established as follows:

20% ≤ SOC ≤ 80% (24)
• BESS maximum charge current: The maximum continuous charge and discharge currents (respectively 3280A and -6400A) are defined for the BESS. Hence, the battery current is constrained as:

-6400 ≤ i * ≤ 3280 (25)

Prediction and optimization strategy

The aim of the optimization strategy is that of generating a control signal that optimizes a cost function over the prediction horizon. The following quadratic cost function is then defined:

Γ(k) = N p i=1 λ 1 P IN J -P IN Jref 2 + (26) λ 2 SOC -SOC ref 2 + Π (nu,N p) i ũ T Q u Π (nu,N p) i ũ
where the control objectives are expressed in a quadratic form. λ 1 and λ 2 are weights allowing to increase/decrease the importance of an objective with respect to the other. ũ is the vector containing the N p elements conforming the control sequence going from u(k) to u(k An appropriate optimization algorithm is used to reduce future errors based on such a quadratic function and in the presence of linear constraints [START_REF] Alamir | A Pragmatic Story of Model Predictive Control: Self-Contained Algorithms and Case-Studies[END_REF]. The MPC prediction algorithm is defined based on the BESS model to feed the cost function at every time step (see [START_REF] Aguilera-González | Mpc ems for a grid-connected renewable energy/battery hybrid power plant[END_REF]).

+ N p -1). Q u ∈ R nu×nu

RESULTS

To test the robustness of our strategy, a 24h profile of wind power production from Guadeloupe Island was analyzed. Two scenarios representing 3-hours commitment profiles were chosen (as can be seen in Fig. 7). These results deal with the electrical grid of Guadeloupe as well as the HPP modeled under PowerFactory environment, while the MPC and the optimization strategy was executed in Matlab. • Case 1: In Fig. 8 a 3-hours commitment profile with steps of 30 minutes is shown, a prediction horizon (N p ) of 15 seconds was defined. In this first case, greater importance was given to respect of the injection band over the SOC maximization. Fig. 8b shows that the injected power P IN J stays inside the tolerance region of 2 MW above and bellow of the commitment profile, (P SCHED ). Between t = 3.5 hours and t = 3.8 hours the power transferred to the grid is less than the commitment. This is because the WECS output smaller than expected that caused the storage system to discharge from its initial SOC of 30 % until reaching 20 % (the lower SOC level allowed), as depicted in Fig. 9. In this second case, the Fig. 10b shows that the injected power P IN J is greater than commitment (WECS overproduction between 15.5-16 hours), but this injection stays always inside the tolerance region. During this time, the wind being stronger than expected makes to charge the battery until the upper limit, 80% SOC (Fig. 11). As shown, the proposed controller is able to efficiently manage the BESS with respect to the constraints defined.

CONCLUSION

In this work, an energy management strategy for the power dispatch of a wind-BESS HPP injecting power into the electrical grid of Guadeloupe Island is proposed. For that, a model predictive control and optimization strategy were developed to drive the storage system with respect to two optimization objectives: respect of a defined region for the power injection, and maximization of the storage system SOC. The strategy uses the model of the BESS to find optimal solutions that bring the system's predicted output close to a trajectory of future power injections defined.

The controller performance was verified in co-simulation PowerFactory/Matlab using real wind speed data collected in Guadeloupe Island. The methodology proved its ability to efficiently manage the power output of the hybrid power plant and meet generation schedules with respect to the constraints defined. The prediction algorithm developed here is compatible with a quadratic programming solver tool to minimize a cost function subject to a set of linear constraints. Those constraints represent the physical and operational burdens of the hybrid plant studied. In future work, forecast reliability analysis for the implementation of the MPC control strategy should be considered.

ACKNOWLEDGMENT

  Fig. 1. Electrical substations on the Guadeloupe island.

Fig. 2 .

 2 Fig. 2. Guadeloupe island electrical model.

Fig. 3 .

 3 Fig. 3. Wind generation system model.

Fig. 4 .

 4 Fig. 4. Energy storage system model. This model was validated through comparison with the battery model available on Simulink and presented in Aguilera-González et al.(2018). For this, a current input was imposed on both system representations (Matlab and PowerFactory) to observe the dynamic of battery states during the charging and discharging process.

  Fig. 5. HPP control scheme.

  is a symmetric positive definite matrix for weighting and adjusting the control effort of the inputs, and Π (nu,N p) i is an identity matrix concatenated with a zero matrix whose sizes change with i, i.e.: Π nu,N p i = I nui×nui O nui×(N p-1)nu (27) Here, Γ can be interpreted as the cost obtained by the evaluation of the future control actions function ũ using the future references P IN Jref and SOC ref and the present system state x(k) within the time horizon [k, k + N p ].

Fig. 7 .

 7 Fig. 7. Wind and power profiles during 24h.

  Wind production and the day-ahead dispatch power profiles.

Fig. 8 .

 8 Fig. 8. Profile and Injection band at morning (2h-5h).

Fig. 10 .

 10 Fig. 10. Profile and Injection band at afternoon (15h-18h).

  Fig. 11. Battery SOC in case 2.

  In Fig.10a second 3-hours commitment profile at afternoon is shown. In this case, more importance was given to the objective of the SOC maximization. Initial SOC was fixed at 78%.

	11,00						DIgSILENT
	[MW]
	9,00
	7,00
	5,00
	3,00
	15,0 1,00	15,6	16,2	16,8	17,4	[h]	18,0
	Pwecs: Pwecs
	Psched: Psched
	(a) Wind production and the day-ahead dispatch power profiles.
	11,00						DIgSILENT
	[MW]
	9,00
	7,00
	5,00
	3,00
	15,0 1,00	15,6	16,2	16,8	17,4	[h]	18,0
	Upper bound: Upper bound
	Pinj: Pinj
	Psched: Psched
	Lower bound: Lower bound
		(b) Power injection band.
							38,00
							[%]
							34,00
							30,00
							26,00
							22,00
							2,00 18,00	2,60	3,20	3,80	4,40	[h]	5,00

SOC: SOC

DIgSILENT

Fig.

9

. Battery SOC in case 1.

• Case 2:

The authors gratefully acknowledge the financial support of the Aquitaine region and BPI France in the context of InsulGrid project.


Hugues Mounier ˚silviu-Iulian Niculescu

Controllability of linear delay systems and of their sampled versions

Keywords: Delay systems, Controllability, Sampling, Module theoretic setting, π-freeness

This paper focuses on the controllability preservation through sampling of linear time-delay systems. We make use of a module theoretic framework acting as a unifying one for most of the existing delay system controllability notions. The controllability properties are envisioned through ring theoretic properties. Some illustrative examples complete the presentation.

INTRODUCTION

Delay systems, e.g. systems modeled by differentialdifference equations represent the simplest class of infinitedimensional processes. They naturally arise in technological systems, where the delay is quite often situated in the input or the output. This infinite dimensional character is merely lost through discretization.

Contrarily to the case of finite-dimensional systems, where the controllability notion is unique, there has been numerous, and seemingly unrelated controllability generalizations for delay systems. More specifically, notions have been elaborated either from a functional analytic approach (see, e.g., [START_REF] Bartosiewicz | Approximate controllability of neutral systems with delays in control[END_REF]; [START_REF] Manitius | Function space controllability of retarded systems: a derivation from abstract operator conditions[END_REF]; [START_REF] Yamamoto | Reachability of a class of infinitedimensional linear systems: an external approach with applications to general neutral systems[END_REF]) or from a more algebraic viewpoint (see, e.g., with a similar formalism as used here [START_REF] Bourlès | Linear Time-Varying Systems[END_REF]; [START_REF] Fabiańska | Applications of the quillen-suslin theorem to multidimensional systems theory[END_REF]; [START_REF] Mounier | Time delay systems[END_REF]; [START_REF] Rudolph | Flatness based control of distributed parameter systems: Examples and computer exercises from various technological domains[END_REF]), with systems over rings [START_REF] Loiseau | Algebraic tools for the control and stabilization of time-delay systems[END_REF]; [START_REF] Sontag | Linear systems over commutative rings: a survey[END_REF]), with distributional rings [START_REF] Bourlès | Generalized convolution behaviors and topological algebra[END_REF], [START_REF] Vettori | Module theoretic approach to controllability of convolutional systems[END_REF]), with dual trajectory modules over an operator ring (see, for instance, [START_REF] Gluesing-Luerssen | Linear Delay-Differential Systems with Commensurate Delays: An Alebraic Approach[END_REF]; [START_REF] Pillai | A behavioral approach to control of distributed systems[END_REF]). There has thus also been many bridging notions between continuous delay systems and sampled data systems (Fridman (2014a,b); [START_REF] Louisell | Delay differential systems with timevarying delay: new directions for stability theory[END_REF]).

In [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF], an algebraic setting was introduced, based on module theoretic properties, which yielded a unifying framework within which most of theses notions could be compared. Having the previous facts in mind, it is natural to ask which controllability property can be preserved under sampling.

The main features of our algebraic framework can be summarized as follows:

' The intrinsic character of the various notions and definitions.

' Varying the coefficients ring and the module properties yield a great flexibility in defining structural notions. ' A framework enabling the unification of most of the existing structural notions of the delay system's literature (see [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF]). ' The constructive character of the structural notions, as demonstrated by the numerous computer algebra implementations (see, e.g., [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF]; [START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF]). ' The independence of the obtained structural results from often tricky existence and uniqueness prerequisites, as is customary for infinite dimensional systems through functional analytic techniques. ' The same framework can be used both for finite-and infinite-dimensional linear systems.

The sampling we shall use are classical ones, and we will not explicitly study the effect of sampled data control on a continuous system.

The aim of this contribution is to investigate whether controllability notions can be preserved through sampling.

The remaining of the paper is organized as follows: In Section 2, general algebraic definitions are given, and the case of finite dimensional and delay systems are given in Section 3. Discretization and interplay between controllability and sampling are studied in Sections 4 and 5. Finally, some concluding remarks end the paper.

PREREQUISITES AND R-SYSTEMS

Most of the commutative algebra we use may be found in standard textbooks, such as [START_REF] Eisenbud | Commutative Algebra with a view toward Algebraic Geometry[END_REF]; Lang (71); [START_REF] Rotman | An Introduction to Homological Algebra[END_REF].

In particular, the basic definitions of module, ideal, free module, . . . , can be found for example in pages 11 to 17 of [START_REF] Eisenbud | Commutative Algebra with a view toward Algebraic Geometry[END_REF]. The notion of torsion freeness is in (Lang,71,p. 388) or (Rotman, 1979, p. 224). For projectivity, see, e.g., (Eisenbud, 1995, A3.2 p. 615) or (Rotman, 1979, Th. 3.11 and 3.14, p. 62-63).

Preliminary definitions

All the rings and algebras we shall consider are commutative, with unity, and without zero divisors. Definition 2.1. Let R be a ring. An R-linear system, or an R-system, or a (linear) system over R, Σ is a finitely generated R-module. l

In other words, an R-system is an R-module with finite free presentation (see, e.g., [START_REF] Eisenbud | Commutative Algebra with a view toward Algebraic Geometry[END_REF], p.17), (Rotman, 1979, p. 90)):

E φ ÝÑ W ÝÑ Σ ÝÑ 0. A presentation matrix P Σ of Σ is one associated with φ.
Notation. We denote by rξs the R-submodule spanned by a subset ξ of Σ. Definition 2.2. An R-linear dynamics, or an R-dynamics, or a (linear) dynamics over R, is an R-system Σ equipped with an input u " pu 1 , . . . , u m q, i.e., a finite subset of Σ, such that the quotient module Σ{rus is torsion. The input u is said to be independent if rus is a free R-module of rank m. l Definition 2.3. An output y " py 1 , . . . , y p q is a finite subset of Σ. An R-system equipped with an input and an output is called an input-output R-system. l

General controllability notions

Varying the rings under consideration through suitable tensor products yields various controllability notions. Definition 2.4. Let A be an R-algebra. The R-system Σ is said to be A-torsion free controllable (resp. A-reflexive controllable, A-projective controllable, A-free controllable) if the A-module A b R Σ is torsion free (resp. projective, free). l

Elementary homological algebra [START_REF] Rotman | An Introduction to Homological Algebra[END_REF] yields: Proposition 2.1. The A-free controllability implies the Aprojective controllability, which implies the A-reflexive controllability, which implies the A-torsion free controllability. l

Localisation (e.g. inversion of some elements) preserves module properties in the following sense (see, e.g., (Rotman, 1979, 4.81, p. 198

)): Proposition 2.2. If S Ď R is multiplicative, then localisa- tion Σ Þ Ñ S ´1Σ " S ´1R b R Σ defines an exact functor Mod R Ñ Mod S ´1R . Thus if Σ is free (resp. projective), then S ´1Σ is free (resp.

projective). l

The following result is borrowed from [START_REF] Rowen | Ring Theory[END_REF], Proposition 2.12.17, p. 233): Proposition 2.3. Let Σ be an R-system, A an R-algebra, and S a multiplicative part of A such that Σ is S ´1R-free controllable. Then, there exists an element π in S such that Σ is Rrπ ´1s-free controllable. l Definition 2.5. Let Σ be an R-system, A an R-algebra, and S a multiplicative part of A containing an element π such that Σ is Rrπ ´1s-free controllable. Such a system Σ will be called π-free. l

We have the following direct consequence: Corollary 2.1. Let Σ be an R-torsion free controllable Rsystem and S a multiplicative part of R such that S ´1R is a principal ideal ring. Then, there exists π P S such that Σ is Rrπ ´1s-free controllable and Σ is π-free. l Remark 2.1. Note that the previous definitions allow for an extreme flexibility in the choice of the structural notions to be considered. First, three distinct algebraic notions are underlined (torsion freeness, projectivity and freeness) and last, but not least, the base change, through extension of scalars, yields a vast number of notions. l

Let us end this subsection by another notion, the spectral controllability, very useful in the infinite-dimensional case.

The following definition (given in [START_REF] Woittennek | Controllability of networks of spatially one-dimensional second order PDE -an algebraic approach[END_REF]) extends previous ones (see, e.g., [START_REF] Kirillova | On the problem of controllability of linear systems with aftereffect[END_REF], [START_REF] Asmykovich | Spectrum control in systems with delay[END_REF], [START_REF] Bartosiewicz | Approximate controllability of neutral systems with delays in control[END_REF], [START_REF] Rocha | Controllability for delay-differential systems[END_REF]) in our context. Definition 2.6. Let R be any ring that is isomorphic to a subring of the ring O of entire functions with pointwise defined multiplication. Denote the embedding R Ñ O by L . A finitely presented R-system with presentation matrix P is said to be spectrally controllable if the Omatrix P " L pP q satisfies Dk P N : @σ P C : rk C P pσq " k

The following proposition establishes a first algebraic characterisation for spectral controllability Proposition 2.4. (See [START_REF] Woittennek | Controllability of networks of spatially one-dimensional second order PDE -an algebraic approach[END_REF]).

Let R be any ring that is isomorphic to a subring of the ring O of entire functions with pointwise defined multiplication. A finitely presented R-system is spectrally controllable if, and only if, the module Σ O " O b R Σ is torsion free.

The next proposition strongly links spectral controllability and R-freeness, which yields a fairly constructive character.

Proposition 2.5. (See [START_REF] Woittennek | Controllability of networks of spatially one-dimensional second order PDE -an algebraic approach[END_REF]).

Let R be any Bézout domain that is isomorphic to a subring of O with the embedding R Ñ O denoted by L . Then the notions of spectral controllability and Rfree controllability are equivalent if and only if L maps non-units in R to non-units in O.

Decomposition

The classical decomposition of a module over a PID into its torsion and torsion free parts generalizes to Bézout domains (see, e.g. (Hazewinkel, 1995, p. 381)). Proposition 2.6. Let R be a Bézout domain, and Σ an Rmodule. One has the following decomposition:

Σ " tΣ ' Σ{tΣ where tΣ is the torsion submodule of Σ. l

FINITE DIMENSIONAL AND DELAY SYSTEMS

Finite dimensional systems

Two common instances of R-systems are given when R is the ring of differential operators Rr d dt s, which yields continuous time linear finite-dimensional systems [START_REF] Fliess | Some basic structural properties of generalized linear systems[END_REF], and when R is Rrδ h s, with δ h the shift operator of delay h, h being the sampling interval, with δ 0 being the identity operator. This yields linear discrete time finitedimensional systems [START_REF] Fliess | Reversible linear and nonlinear discretetime dynamics[END_REF].

On such rings (being common examples of principal ideal rings), all the notions of torsion freeness, projectivity and freeness coincide, resulting in the well known Kalman controllability notion. See [START_REF] Fliess | Some basic structural properties of generalized linear systems[END_REF]; [START_REF] Fliess | Finite poles and zeros of linear systems: an intrinsic approach[END_REF] for several fruitful consequences of this framework.

Various rings and delay systems

We shall consider the following rings in the sequel:

' R p " Rr d dt , δ L s, the polynomial ring; a ring on which delay systems are defined. In the incommensurate delay case, it becomes R p " Rr d dt , δ L s, where δ L " pδ L1 , δ L2 , . . . , δ Lr q and the dimension of the In the incommensurate delay case, it becomes R s " Rrδ h s. The latter ring will not be used here, since in computer implementations, only commensurate delays are used.

Q vector space rL 1 , . . . , L r s Q is r. ' R e "

Delay systems

The definition we shall adopt is the following. Definition 3.1. A pointwise linear delay system is a finitely generated R p -module. A pointwise linear commensurate delay system is a finitely generated R p -module. l Definition 3.2. Let Λ be a pointwise linear delay system. We shall call the exponential version of Λ, denoted by Λ e the R e -system obtained by the functor isomorphism ψ

DISCRETIZATION

We shall define the discretization in the following way. Definition 4.1. Let Λ be a delay system, let h be a real number dividing L (i.e. such that L " N h with N P N). Consider then the ring morphism φ D p,q,h , where p, q P R s φ D p,q,h : R p Ñ pqq ´1R s d dt Þ Ñ ppδ h q qpδ h q this morphism corresponding to an identity approximation, i.e., p and q being such that lim

hÞ Ñ0 ppe ´h d dt q qpe ´h d dt q " d dt
A discretizer of Λ with sampling period h is a functor D p,q,h : Mod Rp Ñ Mod pqq ´1Rs built in the following way: ' Backward difference discretization. In this case, the morphism φ D p,q,h is φ D p,q,h p d dt q "

r Λ " Rr d dt , δ h , qpδ h q ´1s b Rp Λ pextnq D p,q,h Λ " φ Dp,q,h
1 ´δh h i.e. p " p1 ´δh q{h, and q " 1. We shall denote D B,h this discretization. ' Forward difference discretization. In this case, the morphism φ D p,q,h is

φ D p,q,h p d dt q " δ ´1 h ´1 h "
1 ´δh hδ h i.e. p " p1 ´δh q, and q " hδ h . We shall denote D F,h this discretization. ' Tustin discretization. In this case, the morphism φ D p,q,h is

φ D p,q,h p d dt q "
2p1 ´δh q hp1 `δh q i.e. p " 2p1 ´δh q, and q " hp1 `δh q. We shall denote D T,h this discretization. l

INTERPLAY BETWEEN CONTROLLABILITY AND SAMPLING

Ring and dimension level analyis

On a gross level, one has the following proposition, whose proof follows directly from the results in Section 2. Proposition 5.1. Let Λ be a delay system. We have the following:

(1) For Λ: R p -free ô R p -projective ñ ö R p -torsion free

(2) For Λ e : R e -free ô R e -projective ñ ö R e -torsion free (3) For Λ R d : R d -free ô R d -projective ô R d -torsion free (4) For D ˚Λ: R s -free ô R s -projective ô R s -torsion free, where D ˚" D p,q,h , with p, q P R s . Remark 5.1. Thus, from these general module properties implications and equivalences, one is tempted to conclude that from a controllability viewpoint, to work with Λ R d in a continuous setting or to work on D ˚Λ in a discrete setting is equivalent. We shall see in the sequel that this is not true. l

A Disturbing situation

On the ring and dimension level analysis, things appear to be easy to handle. But these levels are quite coarse, and when one looks deeper into a finer, Fitting ideal analysis, the situation appears to be more tricky, as the following examples highlight.

Example 5.1. Consider the following example 9 xptq " s 0 xptq `p1 ´hs 0 q N uptq ´upt ´Lq whose associated module is Λ. The associated presentation matrix is

´d dt ´s0 p1 ´hs 0 q N ´e´L d dt ¯" ´m1 p d dt , e ´L d dt q m 2 p d dt , e ´L d dt q
Hence, the system is spectrally controllable if s 0 ‰ 0 and not spectrally controllable if s 0 " 0. Through backward difference discretization, we get ˆ1 ´δh h ´s0 p1 ´hs 0 q N ´δN h ˙" and the system is not spectrally controllable. The discretization D p,q,h has the following presentation matrix ˆppδ h q qpδ h q ´s0 ppδ h q N ẇhich has no common factor. Hence, the discretized system is controllable. l

p1 ´hs 0 ´δh q ˆ1 h rpδ h q Ṫhus D B,

CONCLUSION

We have exhibited some discretizer notion of for delay systems envisioned through an algebraic module theoretic setting. A ring theoretic and dimension analysis has been conducted, under which strong similarities occur between properties of the continuous system on a distributed delay ring and of the discrete system on the difference operator ring. Several simple examples showed that this analysis is not sufficient and one needs to investigate further through finer tools, such as the Fitting ideal. This is the topic of further research, in progress.

e

  from Mod Rp to Mod Re (where Mod R is the category of R-modules, for a ring R), with ψ e the ring isomorphism: ψ e : R p ÝÑ R e δ L Þ ÝÑ e ´L d dt l Definition 3.3. Let Λ be a pointwise linear delay system. The extension of scalars R a b Re Λ e is called the advanced extension, or R a -extension of Λ, and is denoted as Λ a . The extension of scalars R d b Re Λ e is called the distributed extension, or R d -extension of Λ, and is denoted as Λ d . l

  r Λ prestq where (extn) is an extension of scalars, (rest) a restriction of scalars, and φ Dp,q,h is the restriction of scalars functor associated to φ D p,q,h . l Remark 4.1. Let us remark that Rrδ L s " Rrδ N h s Ă Rrδ h s " R s Let us give three examples of this general definition. Examples 4.1.

  h Λ e is non controllable, whatever the value of s 0 is. Let p and s 0 such that pps 0 q " 0. Consider the system 9 xptq " s 0 xptq `ppe

						l
	Example 5.2. Another example is
		h 9 xptq " xptq `upt ´Lq
	whose associated presentation matrix is
		´1	´h d dt	e	dt ´L d	ānd
	thus, the system is spectrally controllable. Its back-
	ward difference discretization has the presentation matrix
	ˆ1	´h 1 ´δh h	δ N h	h ˙" `δh δ N	Hence,
	the discretized system is non controllable.	l
	Example 5.3. ´L d dt qpuqptq
	The associated presentation matrix is
		p	d dt	´s0 ppe	´L d dt qq

ACKNOWLEDGMENTS

The authors thank an anonymous reviewer for most valuable and helpful comments that helped the authors to improve the overall quality of the paper.