Lucas Brivadis
email: lucas.brivadis@centralesupelec.fr

Antoine Chaillet
email: antoine.chaillet@centralesupelec.fr

Jean Auriol
email: jean.auriol@centralesupelec.fr

Adaptive observer and control of spatiotemporal delayed neural fields

Keywords: observers, adaptive control, persistence of excitation, neural fields, delayed systems

An adaptive observer is proposed to estimate the synaptic distribution between neurons asymptotically from the measurement of a part of the neuronal activity and a delayed neural field evolution model. The convergence of the observer is proved under a persistency of excitation condition. Then, the observer is used to derive a feedback law ensuring asymptotic stabilization of the neural fields. Finally, the feedback law is modified to ensure simultaneously practical stabilization of the neural fields and asymptotic convergence of the observer under additional restrictions on the system. Numerical simulations confirm the relevance of the approach.

Introduction

Neural fields are nonlinear integro-differential equations used to model the activity of neuronal populations [START_REF] Bressloff | Spatiotemporal dynamics of continuum neural fields[END_REF][START_REF] Coombes | Neural Fields: Theory and Applications[END_REF]. They constitute a continuum approximation of brain structures motivated by the high density of neurons and synapses. Their infinite-dimensional nature allows for accounting for the spatial heterogeneity of the neurons' activity and the complex synaptic interconnection between them. Their delayed version also allows to take into account the non-instantaneous communication between neurons. Yet, unlike numerical models of interconnected neurons, in which every single neuron is represented by a set of differential equations, neural fields remain amenable to mathematical analysis. A vast range of mathematical tools are now available to predict and influence their behavior, including existence of stationary patterns [START_REF] Brivadis | Existence of an equilibrium for delayed neural fields under output proportional feedback[END_REF][START_REF] Faugeras | Persistent neural states: stationary localized activity patterns in nonlinear continuous n-population, q-dimensional neural networks[END_REF] stability analysis [START_REF] Faye | Some theoretical and numerical results for delayed neural field equations[END_REF], bifurcation analysis [START_REF] Atay | Stability and Bifurcations in Neural Fields with Finite Propagation Speed and General Connectivity[END_REF][START_REF] Veltz | Interplay between synaptic delays and propagation delays in neural field equations[END_REF], and feedback stabilization [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF].

This interesting compromise between biological significance and abstraction explains the wide range of neural fields applications, which cover primary visual cortex [START_REF] Bertalmío | Cortical-Inspired Wilson-Cowan-Type Equations for Orientation-Dependent Contrast Perception Modelling[END_REF][START_REF] Pinotsis | Contrast gain control and horizontal interactions in v1: A dcm study[END_REF], auditory system [START_REF] Boscain | A bio-inspired geometric model for sound reconstruction[END_REF], working memory [START_REF] Laing | Multiple bumps in a neuronal model of working memory[END_REF], sensory cortex [START_REF] Detorakis | Structure of receptive fields in a computational model of area 3b of primary sensory cortex[END_REF], and deep brain structures involved in Parkinson's disease [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF].

The refinement of modern technologies (such as multi-electrode arrays or calcium imagining) allows to measure neuronal activity with higher and higher spatial resolution. Using these measurements to estimate the synaptic distribution between neurons would greatly help decipher the internal organization of particular brain structures. Currently, this is mostly addressed by offline algorithms based on kernel reconstruction techniques [START_REF] Alswaihli | Kernel reconstruction for delayed neural field equations[END_REF], although some recent works propose online observers (see [START_REF] Burghi | Online estimation of biophysical neural networks[END_REF] for conductance-based models or [START_REF] Brivadis | Online estimation of hilbert-schmidt operators and application to kernel reconstruction of neural fields[END_REF] for delayed neural fields).

In turn, estimating this synaptic distribution could be of interest to improving feedback control of neuronal populations. A particularly relevant example is that of deep brain stimulation (DBS), which consists in electrically stimulating deep brain structures of the brain involved in neurological disorders such as Parkinson's disease [START_REF] Limousin | Electrical stimulation of the subthalamic nucleus in advanced Parkinson's disease[END_REF]. Several attempts have been made to adapt the delivered stimulation based on real-time recordings of the brain activity [START_REF] Carron | Closing the loop of deep brain stimulation[END_REF]. Among them, it has been shown that a stimulation proportional to the activity of a brain structure called the subthalamic nucleus is enough to disrupt Parkinsonian brain oscillations [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]. Yet, the value of the proportional gain depends crucially on the synaptic strength between the neurons involved: estimating it would thus allow for more respectful stimulation strategies.

In this paper, we thus develop an online strategy to estimate the synaptic distribution of delayed neural fields. This estimation relies on the assumed knowledge of the activation function of the population, the time constants involved, and the propagation delays between neurons, as well as online measurement of a part of the neuronal activity. It exploits the theory of adaptive observers for nonlinear systems developed in [START_REF] Besançon | Remarks on nonlinear adaptive observer design[END_REF][START_REF] Besançon | On adaptive observers for systems with state and parameter nonlinearities[END_REF][START_REF] Pyrkin | An adaptive observer for uncertain linear time-varying systems with unknown additive perturbations[END_REF] and allows to reconstruct the unmeasured quantities based on real-time measurements. We then exploit this feature to propose a stabilizing feedback strategy that may be of particular interest to disrupt pathological brain oscillations. This control law estimates the synaptic kernel in real-time and adapts the stimulation accordingly, thus resulting in a dynamic output feedback controller.

The delayed neural fields model is presented in Section 2 together with an introduction to the necessary mathematical formalism. The synaptic kernel estimation is presented in Section 3, whereas its use for feedback stabilization is presented in Section 4. Numerical simulations to assess the performance of the proposed estimation and stabilization techniques are presented in Section 5.

Problem statement and mathematical preliminaries 2.1 Delayed neural fields

Given a compact set Ω ⊂ R q (where, typically, q ∈ {1, 2, 3}) representing the physical support of a neuronal population, the evolution of the neuronal activity z(t, r) ∈ R n at time t ∈ R + and position r ∈ Ω is modeled as the following delayed neural fields [START_REF] Bressloff | Spatiotemporal dynamics of continuum neural fields[END_REF][START_REF] Coombes | Neural Fields: Theory and Applications[END_REF]:

τ (r)
∂z ∂t (t, r) = -z(t, r) + u(t, r) + Ω w(r, r )S(z(t -d(r, r ), r ))dr .

n ∈ N represents the number of considered neuronal population types; for instance, imagery techniques often allow for discrimination between an excitatory and an inhibitory population, in which case n = 2. τ (r) is a positive definite diagonal matrix of size n × n, continuous in r, representing the time decay constant of neuronal activity at position r. S : R n → R n is a nonlinear activation function; it is often taken as a monotone function, possibly bounded (for instance, a sigmoid). w(r, r ) ∈ R n×n defines a kernel describing the synaptic strength between locations r and r ; its sign indicates whether the considered presynaptic neurons are excitatory or inhibitory, whereas its absolute value represents the strength of the synaptic coupling between them. d(r, r ) ∈ [0, d], for some d > 0, represents the synaptic delay between the neurons at positions r and r that typically mainly results from the finite propagation speed along the axons. Finally, u(t, r) ∈ R n is an input that may represent either the influence of non-modeled brain structures or an artificial stimulation signal. We assume that the neuronal population can be decomposed into z(t, r) = (z 1 (t, r), z 2 (t, r)) ∈ R n 1 × R n 2 where z 1 corresponds to the measured part of the state and z 2 to the unmeasured part. In the case where all the state is measured, we simply write z = z 1 and n 2 = 0. Such a decomposition is natural when the two considered populations are physically separated, as it happens in the brain structures involved in Parkinson's disease [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]. It can also be relevant for imagery techniques that discriminate among neuron types within a given population. Accordingly, we define τ i , S ij w ij and u i of suitable dimensions for each population i, j ∈ {1, 2} so that

τ i (r) ∂z i ∂t (t, r) = -z i (t, r) + u i (t, r) + 2 j=1 Ω w ij (r, r )S ij (z j (t -d ij (r, r ), r ))dr . (2)

Problem statement

In the present paper, we are interested in the following control and observation problems:

Problem 2.1 (Estimation).
From the knowledge of S ij , w 2j , τ i and d ij and the online measurement of u i (t) and z 1 (t) for all i, j ∈ {1, 2}, estimate online z 2 (t), w 11 and w 12 .

Problem 2.2 (Stabilization). From the knowledge of S ij , w 2j , τ i and d ij for all i, j ∈ {1, 2} and the online measurement of z 1 (t), find u 1 in the form of a dynamic output feedback law that stabilizes z 1 and z 2 at some reference when u 2 = 0.

As already said, Problem 2.1 is motivated by the advances in imagery and recording technologies and the importance of determining synaptic distribution in the understanding of brain functioning. The assumption that the transmission delays are known is practically meaningful, as these delays are typically proportional to the distance |r -r | between the considered neurons via the axonal transmission speed, which is typically known a priori. Similarly, the time constants τ i (r) are usually directly dependent on the conductance properties of the neurons. The precise knowledge of the activation function S ij is probably more debatable, although recent techniques allow to estimate them based on the underlying neuron type [START_REF] Carlu | A mean-field approach to the dynamics of networks of complex neurons, from nonlinear integrate-and-fire to hodgkin-huxley models[END_REF].

Problem 2.2 is motivated by the development of deep brain stimulation (DBS) technologies that allow electrically stimulating some areas of the brain whose pathological oscillations are correlated to Parkinson's disease symptoms. In our context, the neuronal activity measured and actuated by DBS through u 1 is denoted by z 1 , which corresponds to a deep brain region known as the subthalamic nucleus (STN). We refer to [START_REF] Detorakis | Closed-loop stimulation of a delayed neural fields model of parkinsonian STN-GPe network: a theoretical and computational study[END_REF] for more details on feedback techniques for DBS. One hypothesis, defended by [START_REF] Holgado | Conditions for the generation of beta oscillations in the subthalamic nucleus-globus pallidus network[END_REF], is that these pathological oscillations may result from the interaction between STN and a narrow part of the brain, the external globus pallidus (GPe). The neuronal activity in this area is inaccessible to measurements or stimulation in clinical practice, but it is internally stable and corresponds to z 2 in our model.

A strategy relying on a high-gain approach answered Problem 2.2 in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]. The system under consideration was similar, except that the nonlinear activation function was not applied to the delayed neuronal activity but to the resulting synaptic coupling. In [START_REF] Faugeras | Absolute stability and complete synchronization in a class of neural fields models[END_REF], system (1) is referred to as a voltage-based model, while [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF] focused on activity-based models. It is proven in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]Proposition 3] that under a strong dissipativity assumption (equivalent1 to our Assumption 3.1 below), for any positive continuous map γ : Ω → R and any square-integral reference signal z ref : Ω → R n 1 there exists a positive constant α * depending on parameters of the system such that, for all α > α * , system (2) coupled with the output feedback law u 1 (t, r) = -αγ(r)(z 1 (t, r) -z ref (r)) and u 2 (t, r) = 0 is globally asymptotically stable (and even input-to-state stable) at some equilibrium whose existence is proved in [START_REF] Brivadis | Existence of an equilibrium for delayed neural fields under output proportional feedback[END_REF]. However, one of the drawbacks of this result is that α * is proportional to the L 2 -norm of w 11 , which is usually unknown or uncertain. In practice, this implies a high-gain choice in the controller, which may lead to large values of u 1 (t, r) that are incompatible with the safety constraints imposed by DBS techniques. On the contrary, our goal in this paper is to propose an adaptive strategy that does not rely on any prior knowledge of the synaptic distributions w 11 and w 12 , but at the price of more knowledge on other parameters of the system.

In a preliminary work [START_REF] Brivadis | Online estimation of hilbert-schmidt operators and application to kernel reconstruction of neural fields[END_REF], we have shown that an observer may be designed in the delayfree case to estimate z 2 (t), w 11 and w 12 , hence to answer Problem 2.1. However, this work was done in a framework that does not encompass time-delay systems, and Problem 2.2 was not addressed at all. With such an observer, a natural dynamic output feedback stabilization strategy would be to choose

u 1 (t, r) = -α(z 1 (t, r) -z 1,ref (r)) + z 1 (t, r) -r ∈Ω ŵ11 (t, r, r )S 11 (z 1 (t - d 11 (r, r ), r ))dr -r ∈Ω ŵ12 (t, r, r )S 12 (ẑ 2 (t -d 12 (r, r ), r
))dr where α > 0 is a tunable controller gain, z 1,ref is a reference signal, ŵ1j (t) denotes the estimation of w 1j made by the observer at time t and ẑ2 (t) is the estimation of z 2 (t). Doing so, if the observer has converged to the state, i.e., ŵ1j = w 1j and ẑ2 = z 2 , then the remaining dynamics of z 1 would be τ 1 (r) ∂z 1 ∂t (t, r) = -α(z 1 (t, r) -z 1,ref (r)) so that z 1 would tend towards z

Definitions and notations

Let q be a positive integer, Ω be an open subset of R q and X be a Hilbert space endowed with the norm • X and scalar product •, • X . Denote by

L 2 (Ω, (X , • X )) := {f : Ω → X Lebesgue-measurable | Ω f 2 F < +∞} the Hilbert space of X -valued square integrable functions. Denote by W 1,2 (Ω, (X , • X )) := {f ∈ L 2 (Ω, (X , • X )) | f ∈ L 2 (Ω, (X , • X )} and by W m,2 (Ω, (X , • X )) := {f ∈ W m-1,2 (Ω, (X , • X )) | f ∈ W m-1,2 (Ω, (X , • X )} for m > 1 the usual Sobolev spaces.
If Ω is a compact set, the above definitions hold by replacing Ω by its interior, and we denote by µ(Ω) := Ω dr the Lebesgue measure of Ω. If I is an interval of R, the space of k-times continuously differentiable functions from I to X is denoted by C k (I, X ). We endow C 0 (I, X ) with the norm defined by x C 0 (I,X ) := sup t∈I x(t) X for all x ∈ C 0 (I, X ). If x ∈ X , denote by x * ∈ X its adjoint. If X is a Hilbert space and Y is a Banach (resp. Hilbert) space, then the Banach (resp. Hilbert) space of linear bounded operators from X to Y is denoted by L(X , Y). For any W ∈ L(X , Y), denote by Ker W its kernel and Ran W its range. The map X x → W x Y defines a semi-norm on W , that is said to be induced by W . It is a norm if and only if W is injective. Set L(X ) := L(X , X ). Denote by Id X the identity operator over X .

For any positive integers n and m and any matrix w ∈ R n×m , denote by w its transpose, Tr(w) its trace, w its norm induced by the Euclidean norm, and w F = Tr(w w) its Frobenius norm. Recall that these norms are equivalent and w w F . Hence, for any positive integers n and m, L 2 (Ω, (R n×m , • )) and L 2 (Ω, (R n×m , • F )) are equivalent Hilbert spaces and

• L 2 (Ω,(R n×m , • )) • L 2 (Ω,(R n×m , • F )) . For all i, j ∈ {1, 2}, set X z i = L 2 (Ω, R n i ) and X w ij = L 2 (Ω 2 , R n i ×n j ), so that X z i (resp.
X w ij ) will be used as the state space of z i (resp. w ij ). Set also X z = X z 1 × X z 2 and X w = X w 11 × X w 12 . By abuse of notations, we write

• (Xw ij , • ) := • L 2 (Ω 2 ,(R n i ×n j , • )) and • (Xw ij , • F ) := • L 2 (Ω 2 ,(R n i ×n j , • F ))
. For any positive constant d and any Hilbert space

X , if x ∈ C 0 ([-d, +∞), X ),
:= sup R n i |S ij | .
Remark 2.3. To ease the reading, we have chosen to consider that the integral over Ω is a Lebesgue integral, i.e., that Ω is endowed with the Lebesgue measure. However, note that our work remains identical when considering any other measure for which Ω is measurable. In particular, an interesting case is when Ω = ∪ N k=1 {r k } for some finite family (r k ) 1 k N in R q and the measure is the counting measure. In that case, (1) can be rewritten as the usual finite-dimensional Wilson-Cowan equation [START_REF] Wilson | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF]: for all k ∈ {1, . . . , N },

τ (r k ) ∂z ∂t (t, r k ) = -z(t, r k ) + u(t, r k ) + N =1 w(r k , r )S(z(t -d(r k , r ), r ). ( 3
)
This case will be further investigated in Section 4.2.

Preliminaries on Hilbert-Schmidt operators

Let q, n and m be positive integers and Ω be an open subset of R q . To any map w ∈ L 2 (Ω 2 , R n×m ), one can associates a Hilbert-Schmidt (HS) integral operator W :

L 2 (Ω, R m ) → L 2 (Ω, R n ) defined by (W z)(r)
= Ω w(r, r )z(r )dr for all r ∈ Ω. The map w is said to be the kernel of W . Let us recall some basic notions on such operators (see e.g. [START_REF] Gohberg | Hilbert-Schmidt Operators[END_REF] for more details). The space

L 2 (L 2 (Ω, R m ), L 2 (Ω, R n )) of HS integral operators is a subspace of L(L 2 (Ω, R m ), L 2 (Ω, R n ))
, and is a Hilbert space when endowed with the scalar product defined by

W a , W b L 2 (L 2 (Ω,R m ),L 2 (Ω,R n )) := w a , w b (L 2 (Ω 2 ,R n×m , • F ))
for all W a and W b in L 2 (L 2 (Ω, R n×m )) with kernels w a and w b , respectively. For any Hilbert basis (e k ) k∈N of L 2 (Ω, R m ), we have that

W 2 L 2 (L 2 (Ω,R m ),L 2 (Ω,R n )) = k∈N W e k for all w ∈ L 2 (Ω 2 , R n×m ).
Let p be another positive integer. If W and P are two HS integral operators with kernels w and ρ, in

L 2 (L 2 (Ω, R m ), L 2 (Ω, R n )) and L 2 (L 2 (Ω, R n ), L 2 (Ω, R p )
) respectively, the composition W P is also a HS integral operators, in L 2 (L 2 (Ω, R m ), L 2 (Ω, R p )). Moreover, its kernel is denoted by w • ρ and satisfies

(w • ρ)(r, r ) = Ω w(r, r )ρ(r , r )dr for all r, r ∈ Ω. If W ∈ L 2 (L 2 (Ω, R m ), L 2 (Ω, R n ))
has kernel w, then its adjoint W * is also a HS integral operator and its kernel w * satisfies w * (r, r ) = w(r , r) for all r, r ∈ Ω. In particular, W is self-adjoint if and only if n = m and w(r, r ) = w(r , r) for all r, r ∈ Ω, and w is a positivedefinite kernel if and only if so is W . In that case, w induces a norm on L 2 (Ω, R n ), defined by z → W z L 2 (Ω,R n ) , that is weaker than or equivalent to the usual norm

• L 2 (Ω,R n ) .
To answer Problem 2.1, we propose to estimate w 1j 's in the norm L 2 (Ω, (R n×m , • F )), which, by the previous remarks, is equivalent to estimate their associated HS operators. This operator-based approach has been followed in [START_REF] Brivadis | Online estimation of hilbert-schmidt operators and application to kernel reconstruction of neural fields[END_REF] to answer Problem 2.1 in the delay-free case. In the present paper, we focus on estimating the kernels rather than their associated operators. From a practical viewpoint, the use of the Frobenius norm corresponds to a coefficientwise estimation of the matrices w 1j (r, r ).

Properties of the system

Let us recall the well-posedness of the system (2) under consideration (that was proved in [24, Theorem 3.2.1]), as well as a bounded-input bounded-state (BIBS) property (that we prove below).

Assumption 2.4. The set Ω ⊂ R q is compact and, for all i, j ∈ {1, 2},

τ i ∈ C 0 (Ω, D n i ++ ), u i ∈ C 0 (R + , X z i ), d ij ∈ C 0 (Ω 2 , [0, d]) for some d > 0, w ij ∈ X w ij , and S ij ∈ C 0 (R n j , R n i ) is bounded and globally Lipschitz.
These assumptions are standard in neural field analysis (see, e.g., [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]). In particular, the boundedness of S reflects the biological limitations of the maximal activity that the population can reach. Proposition 2.5 (Open-loop well-posedness and BIBS). Suppose that Assumption 2.4 is satisfied. Then, for any initial condition

(z 1,0 , z 2,0 ) ∈ C 0 ([-d, 0], X z 1 ) × C 0 ([-d, 0], X z 2 ), the open-loop system (2) admits a unique corresponding solution (z 1 , z 2 ) ∈ C 1 ([0, +∞), X z 1 × X z 2 )∩C 0 ([-d, +∞), X z 1 ×X z 2 ). Moreover, if u i is bounded for all i ∈ {1, 2}, then all solutions (z 1 , z 2 ) of (2)
are such that z i and dz i dt are also bounded.

Proof. Well-posedness. The only difference with [24, Theorem 3.2.1] is that τ depends on r. However, since τ is assumed to be continuous and positive, the proof remains identical to the one given in [24, Theorem 3.2.1].

BIBS. For all i ∈ {1, 2}, let τ i be the smallest diagonal entry of τ i (r) when r spans Ω (which exists since τ i is continuous and Ω is compact). For all t 0, we have by Young's and Cauchy-Schwartz inequalities that

τ i 2 d dt z i (t) 2 Xz i - Ω |z i (t, r)| 2 dr + Ω z i (t, r) u i (t, r)dr + Ω z i (t, r) 2 j=1 Ω w ij (r, r )S ij (z j (t -d ij (r, r ), r ))dr dr - Ω |z i (t, r)| 2 dr + 1 4 Ω |z i (t, r)| 2 dr + Ω |u i (t, r)| 2 dr + 1 4 Ω |z i (t, r)| 2 dr + 2 j=1 Ω 2 w ij (r, r ) 2 |S ij (z j (t -d ij (r, r ), r ))| 2 dr dr - 1 2 z i (t) 2 Xz i + u i (t) 2 Xu i + 2 j=1 S2 ij w ij 2 (Xw ij , • ) .
Hence, if u i remains bounded, then z i also remains bounded by Grönwall's inequality. Moreover,

τ i ∂z i ∂t (t) Xz i z i (t) Xz i + u i (t) Xz i + Ω 2 j=1 Ω w ij (r, r )S ij (z j (t -d ij (r, r ), r ))dr 2 dr z i (t) Xz i + u i (t) Xz i + 2 j=1 Sij w ij (Xw ij , • ) .
Hence dz i dt is also bounded if u i is bounded.

In the rest of the paper, we always make the Assumption 2.4, so that the well-posedness of the system is always guaranteed.

3 Adaptive observer

Observer design

In order to design an observer, we first make a dissipativity assumption on the unmeasured part z 2 of the state.

(t) -z b 2 (t)
Xz 2 is converging towards 0 as t goes to +∞. (This fact will be proved and explained in Remark 3.14.) Therefore, Assumption 3.1 can be interpreted as a detectability hypothesis: the unknown part of the state has contracting dynamics with respect to some norm.

We also stress that Assumption 3.1 is commonly used in the stability analysis of neural fields [START_REF] Faugeras | Absolute stability and complete synchronization in a class of neural fields models[END_REF] and ensures dissipativity even in the presence of axonal propagation delays [START_REF] Detorakis | Incremental stability of spatiotemporal delayed dynamics and application to neural fields[END_REF].

Inspired by the delay-free case investigated in [START_REF] Brivadis | Online estimation of hilbert-schmidt operators and application to kernel reconstruction of neural fields[END_REF], let us consider the following observer:

τ 1 (r) ∂ ẑ1 ∂t (t, r) = -α(ẑ 1 (t, r) -z 1 (t, r)) -z 1 (t, r) + u 1 (t, r) + Ω ŵ11 (t, r, r )S 11 (z 1 (t -d 11 (r, r ), r ))dr + Ω ŵ12 (t, r, r )S 12 (ẑ 2 (t -d 12 (r, r ), r ))dr τ 2 (r) ∂ ẑ2 ∂t (t, r) = -ẑ2 (t, r) + u 2 (t, r) + Ω w 21 (r, r )S 21 (z 1 (t -d 21 (r, r ), r ))dr + Ω w 22 (r, r )S 22 (ẑ 2 (t -d 22 (r, r ), r ))dr τ 1 (r) ∂ ŵ11 ∂t (t, r, r ) = -(ẑ 1 (t, r) -z 1 (t, r))S 11 (z 1 (t -d 11 (r, r ), r )) τ 1 (r) ∂ ŵ12 ∂t (t, r, r ) = -(ẑ 1 (t, r) -z 1 (t, r))S 12 (ẑ 2 (t -d 12 (r, r ), r )) (4)
where α > 0 is a tunable observer gain, to be selected later.

Note that ẑ2 has the same dynamics as z 2 . Hence the dissipativity Assumption 3.1 shall be employed to prove observer convergence. The correction terms are inspired by [START_REF] Besançon | Remarks on nonlinear adaptive observer design[END_REF] that dealt with the finite-dimensional delay-free context.

The well-posedness of the observer system is a direct adaptation of [24, Theorem 3.2.1]. The main differences are that τ i 's are space-dependent and ŵ1i are solutions of a dynamical system. Proposition 3.2 (Observer well-posedness). Suppose that Assumption 2.4 is satisfied. Then, for any initial condition (z 1,0 , ẑ1,0 , z 2,0 , ẑ2,0 , ŵ11,0 , ŵ12,0 )

∈ C 0 ( [-d, 0], X z 1 ) 2 ×C 0 ([-d, 0], X z 2 ) 2 × X w 11 ×X w 12 , the open-loop system (2)-(4) admits a unique corresponding solution (z 1 , ẑ1 , z 2 , ẑ2 , ŵ11 , ŵ12 ) ∈ C 1 ([0, +∞), X 2 z 1 × X 2 z 2 × X w 11 × X w 12 ) ∩ C 0 ([-d, +∞), X 2 z 1 × X 2 z 2 × X w 11 × X w 12 ). Proof.
The proof is based on [26, Lemma 2.1 and Theorem 2.3], and follows the lines of [START_REF] Faye | Some theoretical and numerical results for delayed neural field equations[END_REF]Lemma 3.1.1]. First, note that (2)-( 4) is a cascade system where the observer (4) is driven by the system's dynamics [START_REF] Atay | Stability and Bifurcations in Neural Fields with Finite Propagation Speed and General Connectivity[END_REF]. The well-posedness of ( 2) is guaranteed by Proposition 2.5. Now, let (z 1 , z 2 ) be a solution of (2) and let us prove the existence and uniqueness of (ẑ 1 , ẑ2 , ŵ11 , ŵ12 ) solution of (4) starting from the given initial condition. Let us consider the map F :

R + × X z 1 × C 0 ([-d, 0], X z 2 ) × X w 11 × X w 12 → X z 1 × X z 2 × X w 11 × X w 12 such that (4) can be rewritten as d dt (ẑ 1 , ẑ2 , ŵ11 , ŵ12 )(t) = F (t, ẑ1 , ẑ2t , ŵ11 , ŵ12
). Since τ i are continuous and positive, S ij are bounded and w ij are square-integrable over Ω 2 , d i are continuous and u i ∈ C 0 (R + , X z i ), the map F is well-defined by the same arguments than [24, Lemma 3.1.1]. Let us show that F is continuous, and globally Lipschitz with respect to (ẑ 1 , ẑ1t , ŵ11 , ŵ12 ), so that we can conclude with [26, Lemma 2.1 and Theorem 2.3]. Define F 1 taking values in X z 1 , F 2 taking values in X z 2 , F 3 taking values in X w 11 and F 4 taking values in X w 12 so that F = (F i ) i∈{1,2,3,4} From the proof of [24, Lemma 3.1.1], F 1 and F 2 are continuous and globally Lipschitz with respect to the last variables. From the boundedness of S, F 3 and F 4 are also continuous and globally Lipschitz with respect to the last variables. This concludes the proof of Proposition 3.2.

Let us define the estimation error (z 1 , z2 , w11 , w12 ) = (ẑ 1 -z 1 , ẑ2 -z 2 , ŵ11 -w 11 , ŵ12 -w 12 ). It is ruled by the following dynamical system:

τ 1 (r) ∂ z1 ∂t (t, r) = -αz 1 (t, r) + Ω w11 (t, r, r )S 11 (z 1 (t -d 11 (r, r ), r ))dr + Ω ŵ12 (t, r, r )S 12 (ẑ 2 (t -d 12 (r, r ), r ))dr - Ω w 12 (r, r )S 12 (z 2 (t -d 12 (r, r ), r ))dr = -αz 1 (t, r) + Ω w11 (t, r, r )S 11 (z 1 (t -d 11 (r, r ), r ))dr + Ω w12 (t, r, r )S 12 (ẑ 2 (t -d 12 (r, r ), r ))dr + Ω w 12 (r, r )(S 12 (ẑ 2 (t -d 12 (r, r ), r )) -S 12 (z 2 (t -d 12 (r, r ), r )))dr τ 2 (r) ∂ z2 ∂t (t, r) = -z 2 (t, r) + Ω w 22 (t, r, r )(S 22 (ẑ 2 (t -d 22 (r, r ), r )) -S 22 (z 2 (t -d 22 (r, r ), r )))dr τ 1 (r) ∂ w11 ∂t (t, r, r ) = -z 1 (t, r)S 11 (z 1 (t -d 11 (r, r ), r )) τ 2 (r) ∂ w12 ∂t (t, r, r ) = -z 1 (t, r)S 12 (ẑ 2 (t -d 12 (r, r ), r )) (5)

Observer convergence

In what follows, we wish to exhibit sufficient conditions for the convergence of the observer towards the state, meaning the convergence of the estimation error (z 1 , z2 , w11 , w12 ) towards 0. To do so, we introduce a notion of persistence of excitation over infinite-dimensional spaces.

Definition 3.3 (Persistence of excitation). Let X be a Hilbert space and Y be a Banach space. A continuous signal g : R + → X is persistently exciting (PE) with respect to a bounded linear operator P ∈ L(X , Y) if there exist positive constants T and κ such that

t+T t | g(τ ), x X | 2 dτ κ P x 2 Y , ∀x ∈ X , ∀t 0. ( 6
)
Remark 3.4. If X = Y is finite-dimensional and P is a self-adjoint positive-definite operator, then Definition (3.3) coincides with the usual notion of persistence of excitation since all norms on X are equivalent. However, if X = Y is infinite-dimensional, then there does not exist any PE signal with respect to the identity operator on X . (Actually, it is a characterization of the infinite dimensionality of X ). Indeed, if P = Id X , then (6) at t = 0 together with the spectral theorem for compact operators implies that T 0 g(τ )g(τ ) * dτ is not a compact operator, which is in contradiction with the fact that the sequence of finite range operators N j=0 g( jT N )g( jT N ) * converges to it as N goes to infinity. This is the reason for which we introduce this new PE condition which is feasible even if X is infinite-dimensional. Indeed, P induces a semi-norm on X that is weaker than or equivalent to • X . Remark 3.5. When X is infinite-dimensional, note that there exist signals that are PE with respect to an operator P inducing a norm on X (weaker than • X ), and not only a semi-norm. For example, consider X = Y = l2 (N, R) the Hilbert space of square summable real sequences. The signal g : R + → X defined by g(τ ) = ( sin(kτ ) k 2 ) k∈N is PE with respect to P : X → X defined by

P (x k ) k∈N = x k k 2 k∈N
with constants T = 2π and κ = π since 2π 0 sin 2 (kτ )dτ = π for all k ∈ N. Remark 3.6. If X = L 2 (Ω, R n ) for some positive integer n and if P is a HS integral operator with kernel ρ, then (6) is equivalent to

t+T t | g(τ ), x X | 2 dτ κ Ω Ω ρ(r, r )v(r )dr 2 dr, ∀x ∈ X , ∀t 0. ( 7
)
As explained in Remark 3.4, the role of ρ is to weaken the norm with respect to which g has to be PE. If one changes the Lebesgue measure for the counting measure over a finite set Ω as suggested in Remark 2.3, a possible choice of ρ is the Dirac mass: ρ(r, r ) = 1 if r = r , 0 otherwise. In that case, X is finite-dimensional, and we recover the usual PE notion.

Now, let us state the main theorem of this section that solves Problem 2.1.

Theorem 3.7 (Observer convergence). Suppose that Assumptions 2.4 and 3.1 are satisfied.

Define α * :=

2 12 w 12 2 (Xw 12 , • ) 2(1-2 22 w 22 2 (Xw 22 , • ) ) . Then, for all α > α * , for all u 1 , u 2 ∈ C 0 (R + × Ω, R n i ), any solution of (2)-(5) is such that lim t→+∞ z1 (t) Xz 1 = lim t→+∞ z2 (t) Xz 2 = 0
and w11 (t) (Xw 11 , • F ) and w12 (t) (Xw 12 , • F ) remain bounded for all t 0. Moreover, for any solution of (2), the corresponding error system (5) is uniformly Lyapunov stable at the origin, that is, for all ε > 0, there exists δ > 0 such that, if

z1 (t 0 ), z2t 0 , w11 (t 0 ), w12 (t 0 ) Xz 1 ×C 0 ([-d,0],Xz 2 )×Xw 11 ×Xw 12 δ
for some t 0 0, then z1 (t), z2 (t), w11 (t), w12 (t) Xz 1 ×Xz 2 ×Xw 11 ×Xw 12 ε for all t t 0 .

Furthermore, if S 11 and S 12 are differentiable, u 1 and u 2 are bounded 2 and if

ρ 1 ∈ L 2 (Ω 2 , R n 1 ×n 1 ) and ρ 2 ∈ L 2 (Ω 2 , R n 2 ×n 2 ) are self-adjoint positive-definite kernels such that the signal g : t → ((r, r ) → (S 11 (z 1 (t -d 11 (r, r ), r )), S 12 (z 2 (t -d 12 (r, r ), r )))) is PE with respect to P ∈ L(L 2 (Ω 2 , R n 1 +n 2 ), L 2 (Ω, R n 1 +n 2 )) defined by (P (x 1 , x 2 ))(r) := Ω 2 ρ 1 (r, r )x 1 (r , r )dr dr , Ω 2 ρ 2 (r, r )x 2 (r , r )dr dr for all (x 1 , x 2 ) ∈ L 2 (Ω 2 , R n 1 +n 2 ) and all r ∈ Ω, then lim t→+∞ w11 (t) • ρ 1 (Xw 11 , • F ) = lim t→+∞ w12 (t) • ρ 2 (Xw 12 , • F ) = 0. ( 8
)
The proof of Theorem 3.7 is postponed to Section 3.3.

Remark 3.8. In the case where all the state is measured, i.e., n 2 = 0, note that α * = 0. Hence, under the PE assumption on g, the convergence of ŵ11 • ρ 1 towards w 11 is guaranteed for any positive observer gain α. This means that the observer does not rely on any high-gain approach. This fact will be of importance in Section 4, to show that the controller answering Problem 2.2 is not high-gain when the full state is measured, contrary to the approach developed in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF].

Remark 3.9. The obtained estimations of the kernels w 11 and w 12 in (X w 11 , • F ) is blurred by the kernels ρ 1 and ρ 2 . The stronger is the semi-norm induced by ρ j (which is a norm if and only if ρ j is positive-definite), the stronger is the PE assumption, and the finer is the estimation of w 1j . In particular, if the counting measure replaces the Lebesgue measure over a finite set Ω and ρ j is a Dirac mass as suggested in Remark 3.6, then w1j • ρ j = w1j , hence the convergence of ŵ1j to w 1j obtained in Theorem 3.7 is in the topology of L 2 (Ω 2 , (R n i ×n j , • F )), i.e., coefficientwise.

Remark 3.10. The main requirement of Theorem 3.7 lies in the persistence of excitation requirement, which is a common hypothesis to ensure convergence of adaptive observers (see, for instance, [START_REF] Besançon | Remarks on nonlinear adaptive observer design[END_REF][START_REF] Farza | Adaptive observers for nonlinearly parameterized class of nonlinear systems[END_REF][START_REF] Sastry | Adaptive control: stability, convergence, and robustness[END_REF] in the finite-dimensional context and [START_REF] Curtain | Adaptive observers for slowly time varying infinite dimensional systems[END_REF][START_REF] Demetriou | Adaptive observers for a class of infinite dimensional systems[END_REF] in the infinite-dimensional case). Roughly speaking, it states that the parameters to be estimated are sufficiently "excited" by the system dynamics. However, this assumption is difficult to check in practice since it depends on the trajectories of the system itself. In Section 5, we choose in numerical simulations a persistently exciting input (u 1 , u 2 ) in order to generate persistence of excitation in the signal (S 11 (z 1 -d 11 ), S 12 (z 2 -d 12 )). This strategy seems to be numerically efficient, but the theoretical analysis of the link between the persistence of excitation of (u 1 , u 2 ) and that of (S 11 (z 1 -d 11 ), S 12 (z 2 -d 12 )

) remains an open question, not only in the present work but also for general classes of adaptive observers. This issue is further investigated in Section 4.2, where we look for a feedback law allowing simultaneous kernel estimation and practical stabilization. Another approach could be to design an observer not relying on PE, inspired by [START_REF] Pyrkin | An adaptive observer for uncertain linear time-varying systems with unknown additive perturbations[END_REF][START_REF] Wang | On robust parameter estimation in finite-time without persistence of excitation[END_REF] for example. These methods, however, do not readily extend to the infinite-dimensional delayed context that is considered in the present paper. They could be investigated in future works.

Remark 3.11. According to Definition 3.3, the PE assumption on g in Theorem 3.7 can be rewritten as follows: there exist positive constants T and κ such that, for all

(x 1 , x 2 ) ∈ L 2 (Ω 2 , R n 1 +n 2 ), T 0 2 j=1 Ω 2 g j (t + τ, r, r ) T x j (r, r )drdr 2 dτ κ 2 j=1 Ω Ω 2 ρ j (r, r )x j (r , r )dr dr 2 dr ( 9
)
This characterization will be used in the proof of Theorem 3.7.

Remark 3.12. The choice of the operator P is a crucial part of Theorem 3.7. First, its null space is given by Ker

P = {(x 1 , x 2 ) ∈ L 2 (Ω 2 , R n 1 +n 2 ) | Ω x j (r, •)dr ∈ Ker P j , ∀j ∈ {1, 2}} = {(x 1 , x 2 ) ∈ L 2 (Ω 2 , R n 1 +n 2 ) | Ω x j (r, •)dr L 2 (Ω,R n j ) = 0, ∀j ∈ {1, 2}} (
where P j denotes the HS integral operator of kernel ρ j ), since ρ 1 and ρ 2 are positive-definite. Secondly, remark that P can be written as a block-diagonal operator, with two blocks in L(L 2 (Ω, R n 1 )) and L(L 2 (Ω 2 , R n 2 )), respectively. Roughly speaking, this means that the PE signal g must excite "separately" on its two components so that we are able to distinguish them and to reconstruct separately w 11 and w 12 . Finally, note that if d 1j does not depend on r (i.e., d 1j (•, r ) is constant for all r ∈ Ω), then neither does g j (we write g j (t, r ) := g j (t, r, r ) by abuse of notations), and the kernel of P simply means that we do not require to excite the system along r. In other words, in that case, (9) can be rewritten as

T 0 2 j=1 Ω g j (t + τ, r ) T X j (r )dr 2 dτ κ 2 j=1 Ω Ω ρ j (r, r )X j (r )dr 2 dr
where X j defined by X j (r ) = Ω x j (r, r )dr spans L 2 (Ω, R n j ) as x j spans L 2 (Ω 2 , R n j ), which means that g is PE with respect to the HS integral operator having kernel diag(ρ 1 , ρ 2 ), which is a self-adjoint positive-definite endomorphism of L 2 (Ω, R n 1 +n 2 ).

Remark 3.13. One of the drawbacks of the observer (4) is that Theorem 3.7 does not guarantee input-to-state stability (ISS; see, e.g., [START_REF] Mironchenko | Input-to-state stability[END_REF][START_REF] Sontag | Input to State Stability: Basic Concepts and Results[END_REF]) of the error system with respect to perturbations of the measured output z 1 or model errors. This is an important issue in the context of neurosciences since model parameters are often uncertain. In particular, the assumption that S ij and d ij are known is based on models [] that may vary with time and with individuals.

From a mathematical point of view, it is due to the fact that the Lyapunov function V (see [START_REF] Brivadis | Existence of an equilibrium for delayed neural fields under output proportional feedback[END_REF]) used to investigate the system's stability cannot easily be shaped into a control Lyapunov function. Numerical experiments are performed in Section 5 to investigate robustness to measurement noise. From a theoretical viewpoint, in order to obtain additional robustness properties, new observers should be investigated in order to obtain global exponential contraction of the error system, for example, inspired by [START_REF] Burghi | Online estimation of biophysical neural networks[END_REF]. In any case, the Lyapunov analysis performed in Section 3.3.1 is still a bottleneck for proving the convergence of observers of this kind.

3.3 Proof of Theorem 3.7 (observer convergence)

3.3.1 Step 1: Proof that lim t→+∞ z1 (t) Xz 1 = lim t→+∞ z2 (t) Xz 2 = 0
In order to obtain the first part of the result, we seek a Lyapunov functional V for the estimation error dynamics. Inspired by the analysis performed in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF], let us consider the following candidate Lyapunov function:

V (z 1 , z2t , w11 , w12 ) := V z 1 (z 1 ) + V z 2 (z 2 ) + V w 1 ( w11 ) + V w 2 ( w12 ) + W 1 (z 2t ) + W 2 (z 2t ) (10)
where, for all (z 1 , z2t , w11 , w12 )

∈ X z 1 × C 0 ([-d, 0], X z 2 ) × X w 11 × X w 12 , V z i (z i ) := 1 2 Ω zi (r) τ i (r)z i (r)dr, (11)
V w j ( w1j ) := 1 2 Ω 2 Tr( w1j (r, r ) τ j (r) w1j (r, r ))dr dr,

W i (z 2t ) := Ω 2 γ i (r) 0 -d i2 (r,r ) |z 2t (s, r )| 2 dsdr dr, (12)
and γ i ∈ L 2 (Ω, R), for i ∈ {1, 2}, are to be chosen later.

Computing the time derivative of these functions along solutions of (5), we get:

d dt V z 1 (z 1 (t)) = -α Ω |z 1 (t,
i (z 2t ) = Ω 2 γ i (r)(|z 2 (t, r )| 2 -|z 2 (t -d i2 (r, r ), r )| 2 )dr dr.
Combining the previous computations, we obtain that

d dt V (z 1 (t), z2t , w11 (t), w12 (t)) = -α z1 (t) 2 Xz 1 -1 - 2 i=1 Ω γ i (r)dr z2 (t) 2 Xz 2 + 2 i=1 N i (z i (t), ẑ2t , z 2t ) - 2 i=1 Ω 2 γ i (r)|z 2 (t -d i2 (r, r ), r )| 2 dr dr
where

N i (z i (t), ẑ2t , z 2t ) := Ω 2 zi (t, r) w i2 (r, r )(S i2 (ẑ 2 (t-d i2 (r, r ), r ))-S i2 (z 2 (t-d i2 (r, r ), r ))) dr dr.
Let us provide a bound of N i by applying Cauchy-Schwartz and Young's inequalities.

|N i (z i (t), ẑ2t , z 2t )| = Ω 2 zi (t, r) w i2 (r, r ) S i2 (ẑ 2 (t -d i2 (r, r ), r )) -S i2 (z 2 (t -d i2 (r, r ), r )) dr dr Ω zi (t, r) Ω w i2 (r, r ) S i2 (ẑ 2 (t -d i2 (r, r ), r ) -S i2 (z 2 (t -d i2 (r, r ), r )) dr dr Ω zi (t, r) 2 dr Ω Ω w i2 (r, r ) S i2 (ẑ 2 (t -d i2 (r, r ), r ) -S i2 (z 2 (t -d i2 (r, r ), r )) dr 2 dr 1 2ε i zi (t) 2 Xz i + ε i 2 Ω Ω w i2 (r, r ) S i2 (ẑ 2 (t -d i2 (r, r ), r ) -S i2 (z 2 (t -d i2 (r, r ), r )) dr 2 dr 1 2ε i zi (t) 2 Xz i + ε i 2 Ω Ω w i2 (r, r ) 2 dr Ω S i2 (ẑ 2 (t -d i2 (r, r ), r ) -S i2 (z 2 (t -d i2 (r, r ), r )) 2 dr dr 1 2ε i zi (t) 2 Xz i + ε i 2 i2 2 Ω Ω w i2 (r, r ) 2 dr Ω |z 2 (t -d i2 (r, r ), r )| 2 dr dr for all ε i > 0. Now, set γ i (r) := ε i 2 i2
2 Ω w i2 (r, r ) 2 dr for all i ∈ {1, 2}. We finally get that

d dt V (z 1 (t), z2t , w11 (t), w12 (t)) -α - 1 2ε 1 z1 (t) 2 Xz 1 -1 - 2 i=1 Ω γ i (r)dr - 1 2ε 2 z2 (t) 2 Xz 2 . ( 14
)
In order to make V a Lyapunov function, it remains to choose the constants ε i . By definition of γ i , we have .

2 i=1 Ω γ i (r)dr = 1 2 2 i=1 ε i 2 i2 w i2 2 (Xw i2
Then, there exist two positive constants c 1 and c 2 , given by c 1 = α -α * and c 2 = 1 4 (1 -

2 22 w 22 2 (Xw 22 , • )
), such that for any solution of (5),

d dt V (z 1 (t), z2t , w11 (t), w12 (t)) -c 1 z1 (t) 2 Xz 1 -c 2 z2 (t) 2 Xz 2 . ( 15
)
Since τ i ∈ C 0 (Ω, D n i ++ ), V z i and V w j define norms that are equivalent to the norms of X z i and (X w 1j , • F ), respectively. Hence, the error system is uniformly Lyapunov stable, since V is non-increasing. Moreover, zi , and w1j are bounded for i, j ∈ {1, 2}. Moreover, we have for all .

t 0, d dt V z 1 (z 1 (t)) -α z1 (t) 2 Xz 1 + Ω |z 1 (t,
Hence d dt V z 1 (z 1 ) and d dt V z 2 (z 2
) are bounded since zi , and w1j are bounded for all i, j ∈ {1, 2}. Thus, according to Barbalat's lemma applied to

V z i (z i ), V z i (z i (t)) → 0 as t → +∞, hence zi Xz i → 0. Remark 3.14. At this stage, note that d dt (V z 2 (z 2 ) + W 2 (z 2t )) -c 2 z2 2 Xz 2
whenever γ 2 and ε 2 are chosen as above. Hence, Assumption 3.1 implies that z2 is converging towards 0. This justifies that Assumption 3.1 is indeed a dissipativity assumption of the z 2 sub-system, i.e. a detectability assumption, since ẑ2 has the same dynamics than z 2 .

3.3.2

Step 2: Proof that lim t→+∞ w11 (t)P X Xw 11 = lim t→+∞ w12 (t)P Y Xw 12 = 0 Now, assume that t → (S 11 (z 1 (t -d 11 )), S 12 (z 2 (t -d 12 ))) is PE with respect to ρ. The error dynamics (5) can be rewritten as

τ 1 (r) ∂ z1 ∂t (t, r) = f 0 (t, r) + Ω w11 (t, r, r )g 1 (t, r, r )dr + Ω w12 (t, r, r )g 2 (t, r, r )dr τ 1 (r) ∂ w11 ∂t (t, r, r ) = f 1 (t, r, r ) τ 2 (r) ∂ w12 ∂t (t, r, r ) = f 2 (t, r, r )
where g j (t, r, r

) := S 1j (z j (t -d 1j (r, r ), r )), ∀j ∈ {1, 2}, f 0 (t) 2 Xz 1 := Ω -αz 1 (t, r) + Ω ŵ12 (t, r, r )(S 12 (ẑ 2 (t -d 12 (r, r ), r )) -S 12 (z 2 (t -d 12 (r, r ), r )))dr 2 dr 2α 2 z1 (t) 2 Xz 1 + 2 2 12 Ω Ω ŵ12 (t, r, r ) 2 dr Ω z2 (t -d 12 (r, r ), r ) 2 dr dr 2α 2 z1 (t) 2 Xz 1 + 2 2 12 ŵ12 (t) 2 (Xw 12 , • ) sup s∈[-d,0] z2t (s) 2 Xz 2 , f 1 (t) (Xw 11 , • ) := Ω 2 |z 1 (t, r)S 11 (z 1 (t -d 11 (r, r ), r ))| 2 S2 11 µ(Ω) z1 (t) 2 Xz 1 , and
f 2 (t) (Xw 12 , • ) := Ω 2 |z 1 (t, r)S 12 (ẑ 2 (t -d 12 (r, r ), r ))| 2 S2 12 µ(Ω) z1 (t) 2 Xz 1 .
Recall that, according to Step 1, zi (t) Xz i → 0 for all i ∈ {1, 2}, and ŵ12 (t

) (Xw 12 , • ) w12 (t) (Xw 12 , • ) + w 12 (t) (Xw 12 , • ) remains bounded as t → +∞. Hence f 0 (t) (Xz 1 , • ) → 0 and f j (t) (Xw 1j , • ) → 0 as t → +∞ for all j ∈ {1, 2}. Moreover, t → (g 1 (t), g 2 (t)) is PE with respect to ρ by assumption.
Applying twice Duhamel's formula (once on z1 , then once on w1j ) , we get that for all t, τ 0,

τ 1 (r)z 1 (t + τ, r) = τ 1 (r)z 1 (t, r) + τ 0 f 0 (t + s, r)ds + 2 j=1 Ω τ 0 w1j (t + s, r, r )g j (t + s, r, r )dsdr = τ 1 (r)z 1 (t, r) + τ 0 f 0 (t + s, r)ds + 2 j=1 Ω w1j (t, r, r ) τ 0 g j (t + s, r, r )dsdr + 2 j=1 Ω τ 0 s 0 τ j (r) -1 f j (t + σ, r, r )g j (t + s, r, r )dσdsdr . For any t, T 0, define O(t, T ) := T 0 Ω z1 (t+τ, r) τ 1 (r) 2 z1 (t+τ, r)drdτ . Since z1 (t) Xz 1 → 0, O(t, T ) → 0 as t → +∞ for all T 0. Moreover, O(t, T ) = T 0 Ω τ 1 (r)z 1 (t, r) + τ 0 f 0 (t + s, r)ds + 2 j=1 Ω τ 0 s 0 τ j (r) -1 f j (t + σ, r, r )g j (t + s, r, r )dσdsdr 2 drdτ + T 0 Ω 2 j=1 Ω w1j (t, r, r ) τ 0 g j (t + s, r, r )dsdr 2 drdτ + 2 T 0 Ω 2 j=1 Ω w1j (t, r, r ) τ 0 g j (t + s, r, r )dsdr τ 1 (r)z 1 (t, r) + τ 0 f 0 (t + s, r)ds + 2 j=1 Ω τ 0 s 0 τ j (r) -1 f j (t + σ, r, r )g j (t + s, r, r )dσdsdr drdτ Since z1 (t) Xz 1 → 0, f 0 (t) (Xz 1 , • ) → 0 and f j (t) (Xw 1j , • ) → 0 as t → +∞, |g j (t, r, r )|
S1j for all t 0 and all r, r ∈ Ω, and t → w1j (t) Xw 1j is bounded, we get that for any T > 0,

lim t→+∞ T 0 Ω 2 j=1 Ω w1j (t, r, r ) τ 0 g j (t + s, r, r )dsdr 2 drdτ = 0. ( 16
)
For all t, τ 0 and all r ∈ Ω, define h(t, τ, r) := 2 j=1 Ω w1j (t, r, r ) τ 0 g j (t+s, r, r )dsdr . By [START_REF] Curtain | Adaptive observers for slowly time varying infinite dimensional systems[END_REF], h(t, •) L 2 ((0,T ),Xz 1 ) → 0 as t → +∞. Note that ∂h ∂τ (t, τ, r) = 2 j=1 Ω w1j (t, r, r )g j (t + τ, r, r )dr hence h(t, •) ∈ W 1,2 ((0, T ), X z 1 ) and is bounded since g j 's are bounded. Moreover, since u j is supposed to be bounded for all j ∈ {1, 2}, dz i dt is also bounded according to Proposition 2.5. Hence, since S 1j 's are differentiable with bounded derivative, ∂ 2 h ∂τ 2 (t, τ, r) = j=1 Ω w1j (t, r, r ) ∂g j ∂τ (t + τ, r, r )dr is well-defined and bounded. Therefore, for all t 0, h(t, •) ∈ W 2,2 ((0, T ), X z 1 ) and h(t, •) W 2,2 ((0,T ),Xz 1 ) c 3 for some positive constant c 3 independent of t. According to the interpolation inequality (see, e.g., [37, Section II.2.1]),

h(t, •) 2 W 1,2 ((0,T ),Xz 1 ) c 3 h(t, •) L 2 ((0,T ),Xz 1 ) . Thus ∂h ∂τ (t, τ ) L 2 ((0,T ),Xz 1 ) → 0, meaning that lim t→+∞ T 0 Ω 2 j=1 Ω w1j (t, r, r )g j (t + τ, r, r )dr 2 drdτ = 0. ( 17
)
Let (e k ) k∈N be a Hilbert basis of X z 1 . We have, for all t 0, Now, since g is PE with respect to ρ over L 2 (Ω 2 , R n 1 +n 2 ), there exist positive constants T and κ such that (9) holds. for all x j ∈ L 2 (Ω 2 , R n j ), j ∈ {1, 2}, and all t 0. Choosing x j (r, r ) = w1j (t, r, r ) e k (r), we get

T 0 Ω 2 j=1 Ω w1j (t, r, r )g j (t + τ, r, r )dr 2 drdτ k∈N κ 2 j=1 Ω Ω 2
ρ j (r, r ) w1j (t, r , r ) e k (r )dr dr 2 dr.

On the other hand,

2 j=1 w1j (t) • ρ j 2 (Xw 1j , • F ) = 2 j=1 w(t) • ρ j 2 L 2 (Ω 2 ,(R n 1 ×n j , • F )) = 2 j=1 ρ j • w(t) * 2 L 2 (Ω 2 ,(R n j ×n 1 , • F )) = 2 j=1 k∈N Ω Ω (ρ j • w(t) * )(r, r )e k (r )dr 2 dr = 2 j=1 k∈N Ω Ω 2 ρ j (r, r ) w1j (t, r , r ) dr e k (r )dr 2 dr = 2 j=1 k∈N Ω Ω 2 ρ j (r, r ) w1j (t, r , r ) e k (r )dr dr 2 dr 1 κ T 0 Ω 2 j=1 Ω w1j (t, r, r )g j (t + τ, r, r )dr 2 drdτ.
Thus, by [START_REF] Demetriou | Adaptive observers for a class of infinite dimensional systems[END_REF],

lim t→+∞ w11 (t) • ρ 1 (Xw 11 , • F ) = lim t→+∞ w12 (t) • ρ 2 (Xw 12 , • F ) = 0,
which concludes the proof of Theorem 3.7.

Adaptive control

In order to tackle Problem 2.2, we now introduce an adaptive controller based on the previous observer design.

Exact stabilization

Let z 1,ref ∈ X z 1 be a constant reference signal at which we aim to stabilize z 1 . The z 2 dynamics of (2) can be written as When z 1 is constantly equal to z 1,ref and u 2 = 0, we have v 2 = 0. Hence, according to [22, Proposition 3.6], system (18) admits, in that case, a stationary solution, i.e., there exists z

τ 2 (r) ∂z 2 ∂t (t, r) = -z 2 (t, r) + Ω w 22 (r, r )S 22 (z 2 (t -d 22 (r, r ), r ))dr + v 2,ref (r) + v 2 (t, r) (18)
2,ref ∈ X z 2 such that z 2,ref (r) = Ω w 22 (r, r )S 22 (z 2,ref (r ))dr + v 2,ref (r), ∀r ∈ Ω.
Moreover, we have the following stability result. Lemma 4.1 ([14, Proposition 1]). Under Assumption 3.1, system (18) is input-to-state stable ( ISS) at z 2,ref with respect to v 2 , that is, there exist functions β of class KL and ν of class K ∞ such that, for all z 2,0 ∈ C 0 ([-d, 0], X z 2 ) and all v 2 ∈ C 0 (R + , X z 2 ), the corresponding solution z 2 of (18) satisfies, for all t 0,

z 2 (t) -z 2,ref Xz 2 β( z 2 (0) -z 2,ref Xz 2 , t) + ν sup τ ∈[0,t] v 2 (τ ) Xz 2 . ( 19
)
Remark 4.2. The result given in [14, Proposition 1] holds for systems whose dynamics is of the form (after a change of variable)

τ 2 (r) ∂z 2 ∂t (t, r) = -z 2 (t, r) + S 22 Ω w 22 (r, r )z 2 (t -d 22 (r, r ), r )dr + v 2,ref (r) + v 2 (t, r),
which is slightly different from [START_REF] Detorakis | Incremental stability of spatiotemporal delayed dynamics and application to neural fields[END_REF]. However, one can easily check that this modification does not impact at all the proof given in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]. In particular, the control Lyapunov function given in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF] remains a Lyapunov function for (18) (where v 2 is the input).

In particular, Assumption 3.1 implies that z 2,ref is unique. In the case where z 1,ref = 0 and S 2j (0) = 0 for all j ∈ {1, 2}, we also have z 2,ref = 0.

We aim to define a dynamic output feedback law that stabilizes (z 1 , z 2 ) at the reference (z 1,ref , z 2,ref ). We propose the following feedback strategy: for all t 0 and all r ∈ Ω, set

u 1 (t, r) = -α(z 1 (t, r) -z 1,ref (r)) + z 1 (t, r) - r ∈Ω ŵ11 (t, r, r )S 11 (z 1 (t -d 11 (r, r ), r ))dr - r ∈Ω ŵ12 (t, r, r )S 12 (ẑ 2 (t -d 12 (r, r ), r ))dr , u 2 (t, r) = 0, (20)
where α > 0 is a tunable controller gain. The motivation of this controller is that, when w 1j = ŵ1j and ẑ2 = z 2 , the resulting dynamics of z 1 is τ 1

∂z 1 ∂t = -α(z 1 -z 1,ref ).
Hence, since z 2 has a contracting dynamics by Assumption 3.1, this would lead (z 1 , z 2 ) towards (z 1,ref , z 2,ref ). Note that for the controller [START_REF] Detorakis | Structure of receptive fields in a computational model of area 3b of primary sensory cortex[END_REF], the resulting dynamics of ẑ1 in ( 4) is τ 1

∂ ẑ1 ∂t = -α(ẑ 1 -z 1,ref ).
Hence, since the choice of the initial condition of the observer is free, one particular instance of the closed-loop observer is given by ẑ1 (t, r) = z 1,ref (r) for all t 0 and all r ∈ Ω. For this reason, in the stabilization strategy, we can reduce the dimension of the observer by setting ẑ1 = z 1,ref , i.e., z1 = z 1,ref -z 1 . Finally, the closed-loop system that we investigate can be rewritten as system (2) coupled with the controller (20) and the observer system given by

τ 2 (r) ∂ ẑ2 ∂t (t, r) = -ẑ 2 (t, r) + Ω w 21 (r, r )S 21 (z 1 (t -d 21 (r, r ), r ))dr + Ω w 22 (r, r )S 22 (ẑ 2 (t -d 22 (r, r ), r ))dr τ 1 (r) ∂ ŵ11 ∂t (t, r, r ) = (z 1 (t, r) -z 1,ref (r))S 11 (z 1 (t -d 11 (r, r ), r )) τ 1 (r) ∂ ŵ12 ∂t (t, r, r ) = (z 1 (t, r) -z 1,ref (r))S 12 (ẑ 2 (t -d 12 (r, r ), r )) . (21)
First, let us ensure the well-posedness of the resulting closed-loop system.

Proposition 4.3 (Closed-loop well-posedness). Suppose that Assumption 2.4 is satisfied.

For any initial condition (z 1,0 , z 2,0 , ẑ2,0 , ŵ11,0 , ŵ12,0 )

∈ C 0 ([-d, 0], X z 1 ) × C 0 ([-d, 0], X z 2 ) 2 × X w 11 × X w 12
, the closed-loop system (2)-( 20)-( 21) admits a unique corresponding solution

(z 1 , z 2 , ẑ2 , ŵ11 , ŵ12 ) ∈ C 1 ([0, +∞), X z 1 × X 2 z 2 × X w 11 × X w 12 ) ∩ C 0 ([-d, +∞), X z 1 × X 2 z 2 × X w 11 × X w 12 ).
Proof. We adapt the proof of Proposition 3.2. The main difference is that, since the system is in closed loop, the state variables (z 1 , z 2 ) cannot be taken as external inputs in the observer dynamics. Therefore, we consider the map

F : R + × C 0 ([-d, 0], X z 1 ) × C 0 ([-d, 0], X z 2 ) 2 × X w 11 × X w 12 → X z 1 × X 2 z 2 × X w 11 × X w 12 such that (4) can be rewritten as d dt (z 1 , z 2 , ẑ2 , ŵ11 , ŵ12 )(t) = F (t, z 1t , z 2t , ẑ2t , ŵ11 , ŵ12
). Since τ i are continuous and positive, S ij are bounded and w ij are square-integrable over Ω 2 and d i are continuous, the map F is well-defined by the same arguments than [START_REF] Faye | Some theoretical and numerical results for delayed neural field equations[END_REF]Lemma 3.1.1]. Let us show that F is continuous, and globally Lipschitz with respect to (ẑ 1 , ẑ1t , ŵ11 , ŵ12 ), so that we can conclude with [26, Lemma 2.1 and Theorem 2.3]. Define F 1 taking values in X z 1 , F 2 and F 3 taking values in X z 2 , F 4 taking values in X w 11 and F 5 taking values in X w 12 so that F = (F i ) i∈{1,2,3,4,5} From the proof of [24, Lemma 3.1.1], F 1 , F 2 , and F 3 are continuous and globally Lipschitz with respect to the last variables. From the boundedness and global Lipschitz continuity of S, F 4 , and F 5 are also continuous and globally Lipschitz with respect to the last variables. This concludes the proof of Proposition 4.3. Now, the main theorem of this section can be stated. (Xw 22 , • ) ) . Then, for all α > α * , any solution of (2)-( 20)-( 21) is such that ), that is, for all ε > 0, there exists δ > 0 such that, if

lim t→+∞ z 1 (t) -z 1,ref Xz 1 = lim t→+∞ z 2 (t) -z 2,
z 1 (t 0 )-z 1,ref , z 2t 0 -z 2,ref , ẑ2t 0 -z 2,ref , ŵ11 (t 0 )-w 11 , ŵ12 (t 0 )-w 12 Xz 1 ×C 0 ([-d,0],X 2 z 2 )×Xw 11 ×Xw 12 δ
for some t 0 0, then Proof. As in Section 3, define the observer error z2 = ẑ2 -z 2 and w1j = ŵ1j -w 1j . Define also z1 := z 1,ref -z 1 . Then (z 1 , z2 , w11 , w12 ) satisfies [START_REF] Besançon | On adaptive observers for systems with state and parameter nonlinearities[END_REF]. Thus, Theorem 3.7 can be applied. It implies that (z 1 , z2 ) → 0 in X z 1 × X z 1 , ( w11 , w12 ) remains bounded, and this autonomous system is uniformly Lyapunov stable at the origin. Moreover, z 2 satisfies the dynamics [START_REF] Detorakis | Incremental stability of spatiotemporal delayed dynamics and application to neural fields[END_REF] with v 2 (t, r) = Ω w 21 (r, r )(S 21 (z 1,ref (r ) -z1 (t, r )) -S 21 (z 1,ref (r )))dr . In other words, (z 1 , z2 , w11 , w12 ) coupled with z 2 is a cascade system where z 2 is driven by the other variables. According to Lemma 4.1, ( 18) is ISS with respect to v 2 . Thus, z 2 → z 2,ref in X z 2 and is uniformly Lyapunov stable at z 2,ref .

z 1 (t) -z 1,ref , z 2 (t) -z 2,ref , ẑ2 (t) -z 2,
Remark 4.5. Let us consider the case where the full state is measured, i.e., n 2 = 0. Then α * = 0, which means that the controller does not rely on any high-gain approach (see Remark 3.8). This is the main difference between the present work and the approach developed in [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF], where the gain has to be chosen sufficiently large even in the fully measured context. Remark 4.6. Convergence of the kernel estimation ŵ1j towards w 1j is not investigated in Theorem 4.4. Actually, one can apply the last part of Theorem 3.7 to show that, under a PE assumption on the signal g = (S 11 (z 1 (t -d 11 )), S 12 (z 2 (t -d 12 ))), the kernel estimation is guaranteed in the sense of (8). However, since the feedback law is made to stabilize the system, z 1 and z 2 are converging towards (z 1,ref , z 2,ref ). Hence, for g to be PE with respect to P , one must have that for some positive constant T and κ,

T |(S 11 (z 1,ref ), S 12 (z 2,ref )) * v| 2 κ P x 2 Y , ∀x ∈ X , ∀t 0.
i.e. that (S 11 (z 1,ref ), S 12 (z 2,ref )) * induces a semi-norm on X stronger than the one induced by P . This operator being a linear form, it implies that the kernel of P must contain a hyperplane, which makes the convergence of ŵ1j towards w 1j too much blurred by P to claim that any interesting information on w 1j can be reconstructed by this method, except when dim X = 1.

In particular, if (S 11 (z 1,ref ), S 12 (z 2,ref )) = (0, 0), then P = 0, hence no convergence of ŵ1j towards w 1j is guaranteed. The objective of simultaneously stabilizing the neuronal activity while estimating the kernels is investigated in the next section.

Simultaneous kernel estimation and practical stabilization

The goal of this section is to propose an observer and a controller that allow to simultaneously answer Problems 2.1 and 2.2. In order to do so, we make the following set of restrictions (in this subsection only) for reasons that will be pointed in Remark 4.11:

(i) Exact stabilization is now replaced by practical stabilization, that is, for any arbitrary small neighborhood of the reference, a controller that stabilizes the system within this neighborhood has to be designed.

(ii) All the neuronal activity is measured, i.e., n 2 = 0. Hence, we set z = z 1 , w 11 = w, u = u 1 , S = S 11 , τ = τ 1 and d = d 11 to ease the notations.

(iii) The state space X z is finite-dimensional. To do so, as suggested in Remark 2.3, we replace the Lebesgue measure with the counting measure and take Ω as a finite collection of N points in R q , so that X z R (N ×n) .

(iv) The delay d is constant. We write d := d(r, r ) for all r, r ∈ Ω by abuse of notations.

(v) The decay rate τ (r) is constant and all its components are equal, that is, τ (r) = τ Id R n for all r ∈ Ω for some positive real constant τ .

(vi) The reference signal is z ref is constant and S(z ref ) = 0.

(vii) S is locally linear near z ref , that is, there exists ε > 0 such that S(z)

= dS dz (z ref )(z -z ref ) for all z ∈ R n satisfying |z -z ref | ε. Moreover, dS dz (z ref ) is invertible.
Under these restrictions, we now consider the following problem. Problem 4.7. Consider the system (3). From the knowledge of S, τ and d and the online measurement of u(t) and z(t), for any arbitrary small neighborhood of 0 ∈ X z , find u in the form of a dynamic output feedback law that stabilizes z in this neighborhood, and estimate online w.

As explained in Remark 4.6, there is no hope of estimating w when employing the feedback law [START_REF] Detorakis | Structure of receptive fields in a computational model of area 3b of primary sensory cortex[END_REF], since stabilizing z at z ref prevents the persistency condition from being satisfied. For this reason, we suggest adding to the control law a small excitatory signal, whose role is to improve persistency without perturbing too much the z dynamics. Of course, this strategy prevents obtaining exact stabilization. This is why this condition has been relaxed in (i). More precisely, to answer Problem 4.7, we suggest to consider the feedback law

u(t, r) = v(t, r) -α(z(t, r) -z ref ) + z(t, r) - r ∈Ω ŵ(t, r, r )S(z(t -d, r ))dr , (22)
coupled with the observer

τ ∂ ẑ ∂t (t, r) = -αẑ(t, r) + v(t, r) τ ∂ ŵ ∂t (t, r, r ) = -(ẑ(t, r) -z(t, r))S(z(t -d, r )) (23)
where α > 0 is a tunable controller gain and v is a signal to be chosen both small enough in order to ensure practical convergence of z towards 0 and persistent in order to ensure convergence of ŵ towards w. In practice, v can also be seen as an external signal arising from interaction with other neurons whose dynamics are not modeled. Since v is supposed to be continuous, the proof of well-posedness is identical to the proof of Proposition 4.3 and we get the following. Proposition 4.8 (Well-posedness of ( 1)-( 22)-( 23)). Suppose that Assumption 2.4 is satisfied. For any v ∈ C 0 (R + , X z ) and any initial condition (z 0 , ẑ0 , ŵ0 ) ∈ C 0 ([-d, 0], X z ) 2 × X w , the closed-loop system (1)-( 22)-( 23) admits a unique corresponding solution (z, ẑ, ŵ) ∈

C 1 ([0, +∞), X 2 z × X w ) ∩ C 0 ([-d, +∞), X 2 z × X w ).
Theorem 4.9 (Simultaneous estimation and practical stabilization). Suppose that Assumption 2.4 is satisfied. Then, for all α > 0 and all input v ∈ C 0 (R + , X z ), any solution of (2)-( 22)-( 23 Moreover, there exists ε > 0 such that if v ∈ C 1 (R + , X ) is bounded by ε, has bounded derivative, and is PE with respect to Id X , then we also have

lim t→+∞ ŵ(t) -w (Xw, • F ) = 0.
The proof of Theorem 4.9 relies on the following lemma that states standard properties of the PE condition. Lemma 4.10 (Properties of PE). Let X be a Hilbert space and Y be a Banach space. Let g ∈ C 0 (R + , X ) and P ∈ L(X , Y).

(a) If (6) is satisfied only for t t 0 for some t 0 0, then g is PE with respect to P .

(b) If g is PE with respect to P , then t → g(t -d) is also PE with respect to P for any positive delay d 0.

(c) If g is PE with respect to P and W ∈ L(X ), then W g is PE with respect to P W * .

Moreover, if X = Y is finite-dimensional and P = Id X , we also have:

(d) If g is PE with respect to Id X with constants T and κ in (6), and g is bounded by some positive constant M , then M κ T . Conversely, for all positive constants T , κ and M such that M 2κ dim X T

, there exists g ∈ C 1 (R + , X ) that is bounded by M and PE with respect to Id X with constants (T , κ).

(e) If g is bounded and PE with respect to Id X and δ ∈ C 0 (R + , X ) is such that δ(t) → 0 as t → +∞, then g + δ is also PE with respect to Id X .

(f ) If µ is a positive constant and g ∈ C 1 (R + , X ) is bounded, with bounded derivative, and PE with respect to Id X , then any solution z ∈ C 1 (R + , X ) of dz dt = -µz(t) + g(t) is also PE with respect to Id X . Lemma 4.10 is proved in A. r, r ∈ Ω.

In order to test the observer (4), the inputs u i are chosen as spatiotemporal periodic signals with irrational frequency ratio, i.e., u i (t, r) = µ sin(λ i tr) with µ = 10 3 and λ 1 /λ 2 irrational. This choice is made to ensure persistency of excitation of the input (u 1 , u 2 ), which in practice seems to be sufficient to induce persistency of excitation of t → (S 11 (z 1 (t -d)), S 12 (z 2 (td))). Note that for u 1 = u 2 = 0, the persistency of excitation assumption seems to be not guaranteed. Hence the observer does not converge (the plot is not reported). For testing the controllers, the inputs are respectively chosen as [START_REF] Detorakis | Structure of receptive fields in a computational model of area 3b of primary sensory cortex[END_REF] for exact stabilization and ( 22) with v(t, r) = µ sin(λ 1 tr), µ ∈ R, for simultaneous practical stabilization and estimation. In the latter case, we must fix n 2 = 0 as imposed by Section 4.2 (ii).

The parameters of the system (2), the observer (4), and the controller [START_REF] Detorakis | Structure of receptive fields in a computational model of area 3b of primary sensory cortex[END_REF] are set as in Table 1, so that Assumptions 2.4 and 3.1 are fulfilled. The convergence of the observer error (5) towards zero is verified in Figure 2. In particular, the estimation of w 11 by the observer is shown at several time steps in Figure 1. The convergence of the state towards zero ensured by the controller ( 20) is shown in Figure 3. As explained in Remark 4.6, no convergence of the kernels estimation can be hoped for in Figure 3 since stabilizing the state prevents PE.

Figure 4 enlightens the compromise made by the feedback law ( 22) between observation and estimation: by choosing v(t, r) = µ sin(λ 1 tr) with µ = 100, the asymptotic regime of the state remains in a neighborhood of zero (practical stabilization), which allows the convergence of the kernel estimation ŵ towards w. When increasing µ, the estimation rate increases (one obtains a plot similar to Figure 2 for µ = 10 3 ) but the asymptotic regime of the state moves away from zero. On the contrary, when decreasing µ towards zero, one obtains an asymptotic regime of z closer to zero (one obtains a plot similar to Figure 3 for µ = 0.1), but the convergence of ŵ towards w is slower.

Conclusion

In this paper, a new adaptive observer has been proposed to estimate online the synaptic strength between neurons from partial measurement of the neuronal activity. We proved the convergence of the observer under a persistency of excitation condition by designing a Lyapunov functional taking into account the infinite-dimensional nature of the state due to the spatial distribution of the neuronal activity and to the time-delay. We have shown that this observer can be used to design dynamic feedback laws that stabilize the system to a target point, even without persistency of excitation. From the theoretical viewpoint, the main open question remains to extend our result on simultaneous estimation and stabilization. It currently relies on important limitations on the system, that cannot be lifted without a deeper analysis of the PE condition proposed in the paper. In particular, sufficient conditions ensuring that choosing a PE input signal guarantees PE of the state of the neural fields should be sought. Hence, g is PE with respect to Id X with constants T and κ.

(e) Denote by M a bound of g. Denote by T and κ the PE constants of g with respect to Id X . Let ε = κ 4M T . Let t 0 > 0 be such that δ(t) X ε for all t t 0 . Then for all t t 0 and all x ∈ X , Hence, g + δ is PE with respect to Id X by (a).

(f ) This proof follows the one given in [START_REF] Loria | A nested matrosov theorem and persistency of excitation for uniform convergence in stable nonautonomous systems[END_REF]Property 4]. We give it here for the sake of completeness. Denote by M a bound of g and ġ. Let z be a solution of dz dt = -µz(t) + g(t). By Duhamel's formula, z(t) X e -µt z(0) X + t 0 e -µ(t-τ ) g(τ ) X dτ e -µt z(0) X + M µ .

Hence there exists t 0 0 such that z(t) X

2M

µ for all t t 0 . For any x ∈ X , define φ x : R + → R by φ x (t) = -z(t), x X g(t), x X . Then φ x is continuously differentiable and φx = z, x X µg -ġ, x X -| g, x X | 2 . Hence, for all t t 0 and all T > 0, φ x (t + T ) -φ x (t) = t+T t z(τ ), x X µg(τ ) -ġ(τ ), x X dτ

- t+T t | g(τ ), x X | 2 dτ. ( 24
)
Since g is PE with respect to Id X , there exist T, κ > 0 such that, for any k ∈ N, x 2 . Choose k large enough that kκ > 4M 2 µ . Combining ( 24) and ( 25) yields, for all t t 0 and all x ∈ X , t+kT t z(τ ), x X µg(τ ) -ġ(τ ), x X dτ kκ -

4M 2 µ x 2 X . (26)
Finally, we get by Cauchy-Schwartz inequality that

Combining ( 26)-( 27)-( 28), we obtain that for all x ∈ X and all t t 0 , t+kT t | z(τ ), x X | 2 dτ (kκ -4M 2 µ ) 2 kT (µ + 1) 2 M 2 x 2 X which implies that z is persistently exciting with respect to Id X by Lemma 4.10 and (a).

Assumption 3 . 1 (

 31 Strong dissipativity). It holds that 22 w 22 (Xw 22 , • ) < 1. Assumption 3.1 yields that for any pair (z a 2 , z b 2 ) of solutions of (1) (replacing τ , w, S and d by τ 22 , w 22 , S 22 and d 22 ), the distance z a 2

2

 2

2 g

 2 , r, r )g j (t + τ, r, r )dr , r, r )g j (t + τ, r, r ) e k (r)dr dr j (t + τ, r, r ) w1j (t, r, r ) e k (r)dr dr 2 dτ.

where v 2 ,

 2 ref (r) := Ω w 21 (r, r )S 21 (z 1,ref (r ))dr and v 2 (t, r) := u 2 (t, r) + Ω w 21 (r, r )(S 21 (z 1 (t, r)) -S 21 (z 1,ref (r )))dr .

Theorem 4 . 4 ( 2 12 w 12 2 (

 4422 Exact stabilization). Suppose that Assumptions 2.4 and 3.1 are satisfied. Define α * := Xw 12 , • )

S ij = tanh τ i = 1 λ 1 = 100 λ 2

 12

= 10 Figure 1 :

 101 Figure 1: Evolution of the kernel estimation ŵ11 (t, r, r ) when running the observer (4).

Figure 2 :

 2 Figure 2: Evolution of the estimation errors w1i Xw 1i and zi Xz i for i ∈ {1, 2} of the observer (4).

Figure 3 :

 3 Figure3: Evolution of the norm of the state z i Xz i and of the estimation errors w1i Xw 1i and zi Xz i for i ∈ {1, 2} for the control law[START_REF] Detorakis | Structure of receptive fields in a computational model of area 3b of primary sensory cortex[END_REF].

Figure 4 :

 4 Figure 4: Evolution of the norm of the state z Xz and of the estimation errors w Xw and z Xz for the control law (22) with v(t, r) = 100 sin(λ 1 tr).

|

  g(τ ), x X | 2 dτ T M 2 x 2 X , hence M κ T . Conversely, set g(τ ) = 2κ T dim X =1 sin( 2 πτT )e , where e is a basis of X . Then g is bounded by2κ dim X Tand has a bounded derivative. Moreover, for all t 0 and allx = dim X =1 x e ∈ X , t+T t | g(τ ), x X | 2 dτ = 2κ T

  τ ) + δ(τ ), x X | 2 dτ t+T t | g(τ ), x X | 2 dτ -2 t+T t | g(τ ), x X δ(τ ), x X | dτ (κ -2M T ε) x 2

| 2 µ x 2 .

 22 g(τ ), x X | 2 dτ kκ x 2 X , ∀x ∈ X , ∀t 0. (25)Moreover, by Cauchy-Schwartz inequality, |φ x (t)| 2M Hence φ x (t+T )-φ x (t)-4M 2 µ

  ), x X µg(τ ) -ġ(τ ), x X dτ (µ + 1)M x X t+kT t z(τ ), x X dτ (27)andt+kT t z(τ ), x X dτ √ kT t+kT t | z(τ ), x X | 2 dτ .

  1,ref . Assuming adequate contraction properties of the z 2 dynamics, z 2 would tend towards some reference z 2,ref that depends on z 1,ref . In particular, pathological oscillations would vanish in state-state, without any large-gain assumption on the control policy. This motivates us to investigate Problem 2.1 and use the observer in dynamic output feedback to address Problem 2.2.

  we denote by x t the history of x over the latest time interval of length d, i.e., x t (s) = x(t + s) for all t 0 and all s ∈ [-d, 0]. For any n ∈ N, denote by D n ++ ⊂ R n×n the set of positive diagonal matrices. For the sake of reading, if d i2 does not depend on r, i.e. if d i2 (•, r ) is constant for all r ∈ Ω, we simply write d ij (r ) := d ij (r, r ). For any globally Lipschitz map S ij : R n j → R n i , denote by ij its Lipschitz constant and set Sij

  ref Xz 2 = lim t→+∞ ẑ2 (t) -z 2 (t) Xz 2 = 0 and ŵ11 (t) (Xw 11 , • F ) and ŵ12 (t) (Xw 12 , • F ) remain bounded for all t 0. Moreover, the system (2)-(20)-(21) is uniformly Lyapunov stable at (z 1,ref , z 2,ref , z 2,ref , w 11 , w 12

  ref , ŵ11 (t) -w 11 , ŵ12 (t) -w 12 Xz 1 ×X 2 z 2 ×Xw 11 ×Xw 12

	ε
	for all t t 0 .

Table 1 :

 1 System and observer parameters for the numerical simulation of Figures2-4

  | g(τ ), x X | 2 dτ κ P x 2 Y , ∀x ∈ X ,which shows that g is PE with respect to P . ), x X | 2 dτ κ P x 2 Y for all x ∈ X , then for all t d,| g(τ ), x X | 2 dτ κ P x 2 Y , ∀x ∈ X ,which shows that t → g(t -d) is PE with respect to P by (a).which shows that W g is PE with respect to P W * by (a).(d) If g is bounded by M and for all t 0, t+T t| g(τ ), x X | 2 dτ κ x 2 X for all x ∈ X , then Cauchy-Schwartz inequality yields

	t+t 0 +T	t+t 0 +T
		t+t 0
	(b) If, for all t 0, | g(τ t+T t+T t	t+T -d
	| g(τ -d), x X | 2 dτ =
	t	t-d

A Proof of Lemma 4.10

(a) If, for all t t 0 , t+T t | g(τ ), x X | 2 dτ κ P x 2 Y for all x ∈ X , then for all t 0, t | g(τ ), x X | 2 dτ (c) If, for all t 0, t+T t | g(τ ), x X | 2 dτ κ P x 2 Y for all x ∈ X , then for all t 0, t+T t | W g(τ ), x X | 2 dτ = t+T t | g(τ ), W * x X | 2 dτ κ P W * x 2 Y , ∀x ∈ X ,

Actually, there is a typo in the condition stated in [14,Proposition 3]. The mistake is corrected in the proof, making it equivalent to our Assumption 3.1. See also[START_REF] Detorakis | Incremental stability of spatiotemporal delayed dynamics and application to neural fields[END_REF] Theorem 3] for a corrected version of the hypothesis.

This assumption is missing in[START_REF] Brivadis | Online estimation of hilbert-schmidt operators and application to kernel reconstruction of neural fields[END_REF] while it is implicitly used in the proof.

Proof of Theorem 4.9. Clearly, from [START_REF] Faugeras | Persistent neural states: stationary localized activity patterns in nonlinear continuous n-population, q-dimensional neural networks[END_REF] and Grönwall's inequality, lim sup t→+∞ ẑ(t) Xz lim sup t→+∞ v(t) Xz . Moreover, setting z1 = ẑ -z and w11 = ŵ -w, we see that, due to the choice of u given by [START_REF] Faugeras | Absolute stability and complete synchronization in a class of neural fields models[END_REF], (z 1 , w11 ) satisfies [START_REF] Besançon | On adaptive observers for systems with state and parameter nonlinearities[END_REF]. Hence, according to Theorem 3.7, lim sup t→+∞ ẑ(t) -z(t) Xz = 0 and ŵ(t) remains bounded. This yields the first part of the result. To show the second part of the result, it is sufficient to show that t → (r → S(z(t-d, r)) is PE with respect to Id Xz in order to apply the second part of Theorem 3.7.

To end the proof of Theorem 4.9, we use Lemma 4.10 as follows. Using Assumption (vii), let ε > 0 be such that S(z)

) is bounded by ε, has bounded derivative, and is PE with respect to Id X . Such a signal v exists by (d). By (c), τ -1 v is also PE with respect to Id X . Since ẑ -z ref satisfies [START_REF] Faugeras | Persistent neural states: stationary localized activity patterns in nonlinear continuous n-population, q-dimensional neural networks[END_REF], (f ) with µ = τ -1 α shows that ẑ -z ref is also PE with respect to Id X . Moreover, by Grönwall's inequality, there exists

) is PE with respect to Id X . Since ẑ(t) -z(t) → 0 as t → +∞ and S has bounded derivative, S(ẑ(t)) -S(z(t)) → 0. Hence, according to (e), t → S(z(t)) is also PE with respect to Id X . Hence, by (b), t → S(z(t -d)) is PE with respect to Id Xz , which concludes the proof of Theorem 4.9.

Remark 4.11. All assumptions (i)-(vii) have been used in the second part of the proof of Theorem 4.9 at crucial points where, without them, one cannot conclude. In particular, these assumptions allow to use the properties stated in Lemma 4.10. Without them, stronger versions of these properties should be required to show that S(ẑ(• -d)) is PE with respect to Id Xz . For example, without (iii), the property (e) would be required in an infinite-dimensional context, which is impossible since counter-examples can easily be found. Without (v), (f ) would be required for filters of the form dz dt = -Σz(t) + g(t) where Σ is a positive definite matrix, which is also known to be false (see [START_REF] Narendra | Persistent excitation in adaptive systems[END_REF]Example 7]). Similarly, (ii) is required to have that ẑ is a filter of v in the form of [START_REF] Faugeras | Persistent neural states: stationary localized activity patterns in nonlinear continuous n-population, q-dimensional neural networks[END_REF], which is necessary to use (f ). Without (iv), (b) would be required for non-constants delays, which is not possible due to counter-examples such as R + t → (sin(t -d 1 ), cos(t -d 2 )) ∈ R 2 with d 1 = 0 and d 2 = π 2 . Properties (vi) and (vii) are used in the end of the proof of Theorem 4.9. Without them, passing from ẑ -z ref being PE to S(ẑ) being PE remains an open problem.

Numerical simulations

We provide numerical simulations of the observer and controllers proposed in Theorems 3.7, 4.4, and 4.9. The observer simulations are in line with those presented in [START_REF] Brivadis | Online estimation of hilbert-schmidt operators and application to kernel reconstruction of neural fields[END_REF], while the two controllers (exact stabilization and simultaneous practical stabilization and estimation) are new. We consider the case of a two-dimensional neural field (namely, n 1 = n 2 = 1) over the unit circle Ω = S 1 with constant delay d. The kernels are given by Gaussian functions depending on the distance between r and r , as it is frequently assumed in practice (see [START_REF] Chaillet | Robust stabilization of delayed neural fields with partial measurement and actuation[END_REF]): w ij (r, r ) = ω ij g(r, r )/ g L 2 (Ω 2 ;R) , with g(r, r ) = exp(-σ|r -r | 2 ) for constant parameters σ and ω ij given in Table 1. Simulations code can be found in repository [START_REF] Brivadis | KernelEstimation Project[END_REF]. The system is spatially discretized over Ω with a constant space step ∆r = 1/20, and the resulting delay differential equation is solved with an explicit Runge-Kutta (2, 3) method. Initial conditions are taken as z 1 (0, r) = z 2 (0, r) = ẑ1 (0, r) = 1, ẑ2 (0, r) = ŵ11 (0, r, r ) = ŵ12 (0, r, r ) = 0 for all


Sebastian Hörl

Jakob Puchinger

Modeling the ecological and economic footprint of last-mile parcel deliveries using open data: A case study for Lyon

Keywords: parcels, urban, last-mile, logistics, optimization, VRP

Modeling the ecological and economic footprint of last-mile parcel deliveries using open data: A case study for Lyon

Sebastian Hörl, Jakob Puchinger

INTRODUCTION

The amount of parcels delivered in the urban space is expected to increase strongly in the coming years. Today, cities already reflect upon strategies to regulate urban logistics, understanding the complex interplay between its economic, ecological and social impacts becomes ever more important. While ideas and research efforts on sustainable urban logistics policies are gaining traction [START_REF] Mucowska | Trends of Environmentally Sustainable Solutions of Urban Last-Mile Deliveries on the E-Commerce Market-A Literature Review[END_REF][START_REF] Neghabadi | Systematic literature review on city logistics: Overview, classification and analysis[END_REF][START_REF] Patella | The Adoption of Green Vehicles in Last Mile Logistics: A Systematic Review[END_REF]. Recent advances in transport simulation aim to model urban logistics on a systemic level [START_REF] De Bok | Application of an empirical multi-agent model for urban goods transport to analyze impacts of zero emission zones in The Netherlands[END_REF][START_REF] Sakai | SimMobility Freight: An agent-based urban freight simulator for evaluating logistics solutions[END_REF][START_REF] Toilier | Freight transport modelling in urban areas: The French case of the FRETURB model[END_REF], but reliable data remains scarce. The present short paper is an attempt to model one specific sector of urban logistics -home parcel deliveries -solely based on open data for a use case of Lyon.

METHODOLOGY

Our approach follows various steps from generating the parcel demand for a territory and defining the supply in terms of operators and distributions center. We then define cost structures to obtain the used vehicles and driven distances to deliver all parcels based on a cost-minimization and vehicle-routing approach. The individual steps are described below.

1

Demand data

Our demand generation process is based on a synthetic population for the Rhône-Alpes region around Lyon. Such a synthetic population, which is a digital representation of households and persons in a region, along with their socio-demographic attributes can be generated based on open data in France. We make use of a replicable data processing pipeline that can be applied anywhere in France [START_REF] Hörl | Synthetic population and travel demand for Paris and Île-de-France based on open and publicly available data[END_REF]. For the present study, we only consider households in our study area, which comprises the city of Lyon, the Grand Lyon metropolitan region and bordering municipalities with relevant logistics infrastructure (Figure 1). For this perimeter, the population synthesis pipeline generates 1.6 million persons in about 795,000 households.

We fuse the synthetic population data with surveys on the purchasing behavior of the local population. Specifically, [START_REF] Gardrat | Méthodologie d'enquête: Le découplage de l'achat et de la récupération des marchandises par les ménages[END_REF] provides statistics on the annual number of orders made per household based on various socio-demographic characteristics. In [START_REF] Hörl | From synthetic population to parcel demand: Modeling pipeline and case study for last-mile deliveries in Lyon[END_REF] we have proposed a method to make use of this information to generate the probable daily parcel demand for the synthetic households using Iterative Proportional Fitting. Applying the model to our study area yields 16,252 geolocated parcels to be delivered during an average day (Figure 1).

Operator model

The goal of our methodology is to let operators minimize their cost by choosing relevant vehicle types for delivering their assigned parcels and optimizing the vehicle routes. Unfortunately, information on the cost structures of parcel operators is scarce. However, we can assume that the main cost components for offering their service are salaries, vehicle maintenance and investment costs, and per-distance costs. A substantial part of our research was to collect information from gray literature on these cost components. While a detailed analysis of our sources and ag-gregation procedures exceed the scope of this paper, they will be detailed in an extended publication.

Figure 2: Monthly rent versus transport volume

For personnel costs, we assume an effective net salary of 1,300 EUR per month per driver, leading to an approximate gross salary of 1,700 EUR, and to monthly costs of about 3,400 EUR per full-time employee and month for the operator. Divided by 25 operating days, we arrive at daily salary costs of 136 EUR.

In terms of per-vehicle costs, we have examined the long-term rental offers of the major French vehicle manufacturers, along with the characteristics of the advertised vehicles. This analysis has yielded distinct vehicle classes (of about 3.3m 3 , 5m 3 , 10m 3 ) for which costs increase linearly with the transport volume. This is true for thermal and electric vehicles while the slope of the latter is higher (Figure 2). We document the daily unit costs per prototypical vehicle that are used in our model in Table 1 ranging from 210 EUR for a small thermal vehicle up to 800 EUR for large electric truck. In the final optimization we divide these cost by 25 active days per month.

The per-distance costs depend strongly on the consumption of the individual vehicle types.

Based on our analysis of manufacturer offers, we have attached representative values for thermal vehicles (in L/100km) and electric vehicles (in Wh/km) to our prototypical vehicle types in Table 1. The per-distance costs are calculated by multiplying the driven distance per vehicle type with the respective consumption factor and the price for fuel (in EUR/L) and electricity (in ct/kWh), respectively.

Additionally, we have noted down representative CO2 equivalent emissions rates (in gCO2eq/km) for each vehicle type. The rates for thermal vehicles are based on our manufacturer analysis, while the rates for electric vehicles are based on the French average of 90gCO2eq/kWh for electricity production1 .

Finally, Table 1 shows the values for a prototypical cargo-bike (Be) based on current rental offers in France and typical consumption rates.

Operator assignment

To link the demand and the operators, we need to assign an operator to each generated parcel in the synthetic population. For that, we perform weighted random draws from the set of operators based on their market shares. Those market shares have been elaborated from gray literature and a dedicated model. These steps cannot be covered in detail but will be explained in an extended publication. Table 2 shows the resulting parcels assigned to each operator.

Distribution centers

To know from where parcels are dispatched, we make use of the SIRENE database2 , which lists all enterprises and their facilities in France, along with their address and the number of employees. From this database we have extracted all facilities belonging to any of the parcel distributors listed in Table 2. The resulting distribution centers are shown in Figure 1 and Table 2 indicates the number of centers per operator.

For each parcel, we select the distribution center of the respective operator that is the closest in terms of road distance. These road distances are calculated using a network extracted from OpenStreetMap data3 and the osmnx library [START_REF] Boeing | OSMnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks[END_REF]. The process results in nine distribution centers out of 104 with more than 300 assigned parcels (Figure 1).

Heterogeneous Vehicle Routing Problem

Based on the inputs above, we define a Heterogeneous Vehicle Routing Problem (HVRP) per distribution center with the following characteristics:

• The goal is to minimize the overall cost which is the sum of the unit costs and a daily salary per chosen vehicle and the total distance-based cost of the vehicle trajectories. • The operator can vary the number of vehicles of each of the seven types (Table 1) and the individual vehicle routes.

• Vehicles start at the distribution center and must return before the end of the day. Their total active time cannot exceed a daily duration of 10h. It consists of the travel times between parcels and depot; service times of 120s at delivery; and service times of 60s per pick-up. • Vehicles cannot carry more parcels than their capacity allows. We assume 10 parcels per m 3 in Table 1. We allow multiple tours per vehicle during one day.

For each distribution center, we obtain a distance matrix and a travel time matrix from our extracted OpenStreetMap network using the osmnx library [START_REF] Boeing | OSMnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks[END_REF]. Since osmnx calculates travel times based on the speed limits of the road segments, we further inflate these values using averaged factors from the TomTom Congestion Index [START_REF] Cohn | The TomTom Congestion Index[END_REF] factors of Lyon to arrive at approximately congested travel times.

Finally, we solve the resulting Heterogeneous Vehicle Routing Problems using the open-source VRP solver VROOM4 .

RESULTS AND DISCUSSION

We define three individual scenarios:

• Baseline: The scenario is based on our synthetic population for 2022. The prices are chosen such that they reflect the long-term cost structures of the operators that have given rise to the distribution schemes that we see today.

• Today: The scenario considers recent increases in energy prices beginning of 2023 with fuel prices of about 1.90 EUR/L5 and 14 ct/kWh6 . It hence shows how the distribution system could develop in case prices stay at this level in the long term. • Future: The scenario is a future scenario in which we consider an updated synthetic population that considers population growth7 and a general increase of parcels per capita by a factor of two. We assume that prices have increased by +80% for fuel and +60% for electricity.

The results are shown in Table 3. For the Baseline case, we obtain 139 thermal vehicles being used, but only eight electric vehicles and 24 cargo bikes. This reflects today's reality where electric vehicles do not have a large share in the transport system. In terms of vehicle sizes, large vehicles (133) dominate, followed by smaller ones (24) of which the majority are cargobikes. Medium-sized vehicles are rarely used. The total distance driven for last-mile deliveries is almost 10,000 km per day for thermal vehicles, the distance for electric vehicles and cargobikes is ten times smaller. Only 7% of all parcels are delivered by electric vehicles or cargobikes. In terms of consumption, 780 liters of fuel are needed and 236 kWh of electricity. This consumption translates into about 2100 kg of CO2 equivalents emitted for the last-mile deliveries during one day which makes 131g per parcel. To calculate the total consumed energy we assume a conversion rate of 10 kWh/L and arrive at a total of 8000 kWh per day with 497 Wh per parcel.

For the Today scenario with adjusted prices (increase of 30% for fuel and 55% for electricity), we see a slight shift of electric vehicles from eight to 14. Still, this shift represents a doubling of the driven distance of electric vehicles and a doubling in parcels delivered by electric vehicles while their overall percentage remains low with 7% for electric transporters and 12% for cargobikes. Accordingly, electricity use doubles while fuel consumption drops by 10%. These shifts lead to a reduction in emissions by 8% in total to 120g per parcel. Total energy use is also reduced by 5% while no large shifts in the cost structures can be observed. Despite electricity prices having increased stronger than fuel prices, the observed shifts can be explained by the different ratios of capital expenses versus operational expenses between thermal and electric vehicles. The latter have higher vehicle prices with lower per-distance costs. In the Baseline case, the break-even daily distance at which a single electric vehicle becomes cheaper in total is at about 34 km, while the point shifts to 26 km in the Today scenario (see Figure 3). In the Future scenario, both the numbers of thermal and electric vehicles increase because of the higher demand. However, electric vehicles increase strongly from 8 to 101. In terms of vehicle size especially large vehicles double in count. At the same time the distance for thermal vehicles goes down and the distance for electric vehicles increases tenfold. Interestingly, the total distance is not doubled, but increases by only 50%, which shows that there are scale effects with respect to the transported parcel volumes. While parcels delivered by electric vehicles are rare Today they make up 50% of all flows in the Future scenario. In the latter, fuel consumption goes down by 37% while electricity use increases by a factor of 13 and total used energy by a factor of 10. On the contrary, total emissions decrease, but only by 24% despite a reduction of 63% per parcel. This effect is due to the generally increased demand. Total costs increase by 75% but not equally on all cost components (67% on salaries, 144% on vehicles, 109% on distance), which puts a higher influence on operational costs on the overall costs. Per parcel, there is a margin of 26ct per parcel between the Today and Future scenario.

In all scenarios, we see that cargo-bikes are rarely used because of their limited capacity.

CONCLUSION

In this paper we have documented a model on the economic and ecological characteristics of last-mile parcel deliveries in a city. While the model makes use of a multitude of assumptions, its main value lies in the comparison of scenarios. In the future, more detailed and distinct scenarios should be evaluated. In terms of validation, system-level reference data is not likely to emerge in the near future. Hence, we are engaging actively in discussing our operational assumptions with domain experts and practitioners to compile a comprehensive list of limitations and future improvements, which will be detailed in an extended publication on the model.

Figure 1 :

 1 Figure 1: Map of the study area, generated parcels, and distribution centers (Background: OpenStreetMap)

Figure 3 :

 3 Figure 3: Break-even points for a small electric vehicle

Table 1 : Operator model

 1

	Vehicle type	St	Mt	Lt	Se	Me	Le	Be
	Size	S	M	L	S	M	L	S
	Propulsion	T	T	T	E	E	E	E
	Capacity	33	50	100	33	50	100	14
	Fuel consumption	5	6	8	-	-	-	-
	[L/100km]
	Electricity consumption	-	-	-	160	200	300	42
	[Wh/km]
	Unit cost	210	260	370	260	400	800	160
	[EUR/month]
	Distance cost*	304.5.0	377.00	522.00	14.00	18.00	27.00	3.80
	[EUR/100km]	0
	Emissions**	130	160	215	14.4	18	27	3.8
	[gCO2eq/km]

Be: Cargo-bike; Size: S -Small, M -Medium, L -Large; Propulsion: T -Thermal, E -Electric; *Indicative distance costs based on 1.45 EUR/km and 9ct/kWh; **Electric vehicle emissions based on 90gCO2eq/kWh

Table 2 : Operator statistics

 2

	Operator	Distribution centers	Market share [%]	Parcels
	La Poste (Colissimo)	72	40.08	6,384
	Chronopost	6	14.98	2,430
	UPS	2	13.55	2,210
	DPD	3	9.94	1,632
	DHL	8	8.95	1,477
	GLS	2	6.93	1,169
	Colis privé	2	5.36	917
	Fedex	9	0.21	33
	Total	104	100	16,252

Table 3 : Optimization results

 3

		Baseline	Today	Future
	Scenario
	Population	2022	2022	2030
	Demand factor	1.0	1.0	2.0
	Fuel price [EUR/L]	1.45	1.90	3.40
	Electricity price [EUR/kWh]	0.09	0.14	0.23
	Vehicles by type
	Thermic	139	133	164
	Electric	8	14	101
	Cargo-bike	24	24	21
	Vehicles by size
	Small (S)	32	33	32
	Medium (M)	6	6	8
	Large (L)	133	132	246
	Distances [km]
	Thermic	9,835	8,916	6,185
	Electric	1,194	2,294	11,438
	Cargo-bike	750	785	861
	Total	11,778	11,995	18,484
	Parcels
	Thermic	15,001	14,493	21,926
	Electric	500	1,009	11,236
	Cargo-bike	751	750	761
	Total	16,252	16,252	33,923
	Consumption
	Fuel [L]	783	710	494
	Electricity [kWh]	236	546	3,261
	Environment
	Emissions [kgCO2eq]	2,127	1,956	1,622
	Per parcel [gCO2eq]	131	120	48
	Energy [kWh]	8,071	7,642	8,205
	Per parcel [Wh]	497	470	242
	Cost [EUR]
	Salaries	23,256	23,256	38,896
	Vehicles	2,217	2,302	5,416
	Distance	1,159	1,426	2,427
	Total	26,632	26,985	46,740
	Per Parcel	1.64	1.66	1.38

https://www.rte-france.com/eco2mix

https://www.sirene.fr

https://download.geofabrik.de

https://github.com/VROOM-Project/vroom

Diesel, France, 09/01/2023, https://www.tolls.eu/fuel-prices

EUROSTAT, Non-household, S2 2022, https://ec.europa.eu/eurostat/cache/infographs/energy_prices/enprices.html

Based on INSEE prediction scenarios https://www.insee.fr/fr/information/6536990

ACKNOWLEDGEMENTS

This paper presents work developed in the scope of the project LEAD, which has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement no. 861598. The content of this paper does not reflect the official opinion of the European Union. Responsibility for the information and views expressed in this paper lies entirely with the authors.