Quentin Clairon

Adeline Samson

Optimal control for estimation in partially observed elliptic and hypoelliptic stochastic differential equations

Keywords: Stochastic Differential Equations, Ellipticity, Hypoellipticity, Estimation, Optimal Control Theory, Linear-Quadratic Theory, Pontryagin maximum principle

Multi-dimensional Stochastic Differential Equations (SDEs) are a powerful tool to describe dynamics of several fields (pharmacokinetic, neurosciences, ecology, etc). The estimation of the parameters of these systems has been widely studied. We focus in this paper in the case of partial observations, only a one-dimensional observation is available. We consider two families of SDE, the elliptic family with a full-rank diffusion coefficient and the hypoelliptic family with a degenerate diffusion coefficient. The estimation for the second class is much more difficult and only few references have proposed estimation strategies in that case. Here, we adopt the framework of the optimal control theory to derive an estimation contrast (or cost function) based on the best control sequence mimicking the (unobserved) Brownian motion. We propose a full data-driven approach to estimate the parameters of the drift and of the diffusion coefficient. The estimation reveals to be very stable in a simulation study conducted on different examples (Harmonic Oscillator, FitzHugh-Nagumo, Lotka-Volterra).

Introduction

We focus on the statistical inference for d-dimensional stochastic dynamical systems modeled by a stochastic differential equation (SDE). We are interested in the case of partial observations: only the first one-dimensional coordinate, denoted V t , of the system is observed while the other (d -1)-dimensional coordinates, denoted U t , are unobserved. The system is written as follows:

dV t = a 1 (V t , U t , t; θ)dt + σ 1 dW 1t
(1)

dU t = a 2 (V t , U t , t; θ)dt + σ 2 dW 2t
where a 1 and a 2 are the two drift functions that depend both on V t and U t , θ are the drift parameters, (W 1t ) t and (W 2t ) t are two independent Brownian motions and σ 1 , σ 2 are the two diffusion coefficients.

We consider two classes of models (1). The first class, called elliptic, corresponds to an SDE with a full non degenerate diffusion coefficient. This means that denoting B σ = (σ 1 , σ 2 ) the diffusion coefficient, the matrix B σ B T σ is full rank, X T being the transposed matrix of X. The second class, called hypoelliptic, corresponds to an SDE with a degenerate stochastic noise: the diffusion coefficient B σ B T σ is not invertible. For example when σ 1 = 0. These two specificities, partial observations and hypoelliptic/elliptic properties, are of increasing importance in many applications, and we give some examples below. But before describing the examples, let us remark that these two specificities are not of the same nature. The first one is linked to the type of observations. In many examples, the system is complex and is modeled by a multi-dimensional system, while the experimentalists are only able to measure, often at discrete times, a one-dimensional signal. As will be recalled later when reviewing the literature, this increases the difficulty of estimating the parameters of model (1). The second specificity is not linked to any experimental constraint, but is a mathematical way of describing the intrinsic noise of the process (V t , U t ). The two systems, elliptic and hypoelliptic, may not have the same interpretation depending on the applications (see below examples in neurosciences). It might nevertheless be difficult for the modeler to know in advance if the system is more likely elliptic rather than hypoelliptic. Unfortunately, estimation methods are often strongly different depending on the nature of the noise (see more details below), and may fail down when applied on the 'wrong' class of models. This is an advantage of our method which is the same for elliptic and hypoelliptic SDE.

Let us now give some examples of applications. Partially observed SDEs have been used in pharmacokinetics for modeling the concentration of a drug in the body, either in a elliptic or a hypoelliptic version [START_REF] Ditlevsen | Parameter estimation in a stochastic model of the tubuloglomerular feedback mechanism in a rat nephron[END_REF][START_REF] Cuenod | Parameter estimation and change-point detection from dynamic contrast enhanced MRI data using stochastic differential equations[END_REF][START_REF] Donnet | A review on estimation of stochastic differential equations for pharmacokinetic -pharmacodynamic models[END_REF]. In system biology, the famous stochastic Lotka-Volterra model [START_REF] Lotka | Elements of Mathematical Biology[END_REF][START_REF] Meeds | Optimization monte carlo: Efficient and embarrassingly parallel likelihood-free inference[END_REF][START_REF] Graham | Asymptotically exact inference in differenciable generative models[END_REF][START_REF] Mao | Environmental brownian noise suppresses explosions in populations dynamics[END_REF] describes the interaction between two species, predator and prey, through a two-dimensional elliptic system. It is often possible to observe only one of the two species, leading to partial observations. In neurosciences, several stochastic systems have been proposed to model the dynamic of one single neuron. The first equation V t corresponds to the dynamics of the membrane potential of the neuron and U t to a recovery variable, or a synaptic conductance, that can not be measured. We can cite the synaptic-conductance based models [START_REF] Pospischil | Extracting synaptic conductances from single membrane potential traces[END_REF][START_REF] Paninski | A new look at state-space models for neural data[END_REF], 2012[START_REF] Ditlevsen | The morris-lecar neuron model embeds a leaky integrate-andfire model[END_REF][START_REF] Ditlevsen | Hypoelliptic diffusions: discretization, filtering and inference from complete and partial information[END_REF] or the FitzHughNagumo model [START_REF] Gerstner | Spiking Neuron Models[END_REF]. These models have been proposed with stochastic noise on the synaptic conductance dynamic (U t ) only, leading to hypoelliptic SDEs [see e.g. [START_REF] Paninski | Inferring synaptic inputs given a noisy voltage trace via sequential Monte Carlo methods[END_REF] references therein], or on both coordinates leading to elliptic SDE [START_REF] Ditlevsen | The morris-lecar neuron model embeds a leaky integrate-andfire model[END_REF]. A last class of models is the stochastic hypoelliptic Damping Hamiltonian system where the first coordinate represents the position of a particle and the second its velocity (d = 2). It is natural that the (Brownian) noise appears only in the velocity coordinate, the position being defined as the deterministic infinitesimal integral of the speed. The position can be measured with precision, but the speed is not directly available. In these models, the first equation of the dynamical system reduces to dV t = U t dt.

Let us now review the main estimation methods that have been proposed in the literature to estimate the drift parameter θ and the diffusion coefficient B σ .

Estimation of elliptic SDE has been widely studied. In the complete observations cases (both V t and U t observed), we can cite among others [START_REF] Bibby | Martingale estimating functions for discretely observed diffusion processes[END_REF], [START_REF] Pedersen | A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations[END_REF], [START_REF] Kessler | Estimation of an ergodic diffusion from discrete observations[END_REF], [START_REF] Aït-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF], [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF], [START_REF] Sørensen | Parametric inference for diffusion processes observed at discrete points in time: a survey[END_REF], [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes[END_REF], [START_REF] Jensen | Markov chain Monte Carlo approach to parameter estimation in the FitzHugh-Nagumo model[END_REF], [START_REF] Ditlevsen | Estimation in the partially observed stochastic morris-lecar neuronal model with particle filter and stochastic approximation methods[END_REF], van der Meulen and Schauer [2016a]. The case of partial observations has also been considered with several approaches. The unobserved coordinates are treated as missing data and are imputed, see for examples [START_REF] Elerian | Likelihood inference for discretely observed nonlinear diffusions[END_REF], [START_REF] Bjørn | MCMC analysis of diffusion models with application to finance[END_REF], Golightly andWilkinson [2006, 2008], [START_REF] Ditlevsen | Estimation in the partially observed stochastic morris-lecar neuronal model with particle filter and stochastic approximation methods[END_REF], van der Meulen and Schauer [2016b]. Most methods propose to approximate the transition density by the Euler-Maruyama scheme and consider a Monte-Carlo approximation to impute and filter the unobserved coordinates. Therefore, they are computationally intensive. We will show that the methodology we develop is less demanding in terms of time of computation.

Let us now explain why the estimation of hypoelliptic systems is more difficult. Let us imagine that the complete observations of (V t , U t ) are available in continuous time. Estimating θ would be naturally performed through the Girsanov formula, that gives directly the likelihood [START_REF] Lipster | Statistics of random processes I : general theory[END_REF]. However, the Girsanov formula requires the matrix B σ B T σ to be invertible. Because of the singularity of B σ in the hypoelliptic case, this inverse does not exist and the likelihood is not properly defined. The same problem occurs for the estimation methods developed for elliptic systems that have been cited above, as they also generally require B σ B T σ to be invertible. They can thus not be applied to hypoelliptic systems. There are thus only few references for hypoelliptic SDEs. The stochastic Damping Hamiltonian system (with dV t = U t dt) has been the most studied. In the parametric framework, [START_REF] Gloter | Parameter estimation for a discretely observed integrated diffusion process[END_REF], [START_REF] Samson | A contrast estimator for completely or partially observed hypoelliptic diffusion[END_REF] propose Euler contrasts with a correction of the bias due to partial observations. [START_REF] Pokern | Parameter estimation for partially observed hypoelliptic diffusions[END_REF] propose a Gibbs sampling in a bayesian approach, but do not correct the bias. In the non-parametric framework, Cattiaux et al. [2014a,b], [START_REF] Comte | Adaptive estimation for stochastic damping hamiltonian systems under partial observation[END_REF] consider the estimation of the drift, the diffusion coefficient and the invariant density using kernel estimators. For hypoelliptic SDEs that are more general than the stochastic Damping Hamiltonian system, we are only aware of the work of [START_REF] Ditlevsen | Hypoelliptic diffusions: discretization, filtering and inference from complete and partial information[END_REF]. Their approach is based on a discretization scheme of order 1.5, a particle filter to approximate the unobserved coordinate and a maximization of a statistical contrast by stochastic approximation. The main drawbacks are the computational time induced by the particle filter and the fastidious calculations to exhibit the sufficient statistics of the likelihood that are then stochastically approximated. There is thus a need to develop new approaches.

In this paper, we propose a strategy based on the optimal control theory. Let us recall the main question the optimal control theory aims to address. For a given dynamical system in a given initial state, which input, or control, do we have to apply to it in order to steer it to a desired behavior in an optimal way ? This problem is formulated as an optimization problem under constraint where the optimality is defined through the introduction of a cost function and the proposed dynamical model belongs to the contraints. Under this theory a large variety of theoretical and numerical tools have been developed to solve this kind of optimization problems, that is to find the so-called optimal control which minimizes the proposed cost function. Recently, this theory has been advantageously used for statistical purpose. We can cite among others the pioneer work of [START_REF] Martin | Optimal control, statistics and path planning[END_REF] for non-parametric estimation of B-splines. Parametric approaches have been proposed more recently by Brunel and [START_REF] Brunel | A tracking approach to parameter estimation in linear ordinary differential equations[END_REF], Clairon andBrunel [2016, 2017], [START_REF] Iolov | Optimal design for estimation in diffusion processes from first hitting times[END_REF], [START_REF] Zhang | Quantization and stochastic ccontrol of trajectories of underwater vehicle in bearings-only tracking[END_REF]. One way of using the optimal control theory is to rewrite an estimation criteria, typically a likelihood or a posterior distribution as a tracking problem. Tracking problems are a specific class of optimal control theory problems: the aim is to find the optimal control which leads a coordinate of the system the closest possible to a target trajectory, here the observations, on a given observation interval. A way of solving an optimal control problem is the Pontryagin maximum principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF][START_REF] Trelat | Controle optimal : Theories et applications[END_REF][START_REF] Sontag | Mathematical Control Theory: Deterministic finite-dimensional systems[END_REF]] which allows to calculate this optimal control as well as the corresponding system response.

This idea has already been successfully developed by Brunel and [START_REF] Brunel | A tracking approach to parameter estimation in linear ordinary differential equations[END_REF] to estimate the parameters of ordinary differential equations. It proves to be numerically efficient and stable, especially when the problem is ill-conditioned. Moreover, the consistency and the rate of convergence of the corresponding estimator have been proven [Clairon andBrunel, 2016, 2017]. However, their method, based on the deterministic theory of optimal control, is restricted to the case of deterministic system. In this paper, we adapt this idea to the problem of SDE estimation. To do so, we resort to the framework of the discrete optimal control theory. Theoretical and numerical results have been fully developed for linear model, this is the discrete linear-quadratic (LQ) theory, a particular case of the Pontryagin maximum principle. Indeed, this theory ensures the existence, uniqueness and gives the closed form of the solution of the control problem defining our estimation criteria. The main advantage of this theory is that it applies without any hypothesis on the diffusion coefficient B σ . Especially, it applies also to the hypoelliptic case. This is the approach that we use for linear SDE. When estimating parameters of a nonlinear SDE, we propose to rewrite the estimation criteria to enter this theory as well. Unfortunately, the consistency and convergence results presented for ODE models in Brunel and [START_REF] Brunel | A tracking approach to parameter estimation in linear ordinary differential equations[END_REF] do not apply here and we were not able to prove a theoretical result for our estimator. Nevertheless, it reveals to be efficient in practice.

Our estimation procedure follows three steps. First we define a criteria to estimate the drift parameter θ alone, σ being considered known. This criteria is minimized with the help of the linear-quadratic theory. This method introduces a weighting parameter w which needs to be selected. The second step consists in constructing an external criteria based on moments of the Brownian process allowing to data-select the weight w. The third step is the estimation of σ by profiling the functional used in step 2.

The paper is organized as follows. Section 2 presents the elliptic and hypoelliptic models. Section 3 introduces the estimation of θ when σ is known. In section 4, the method used to select the weighting parameter w is presented. Section 5 explains the estimation of σ. Section 6 illustrates the procedure on four elliptic or hypoelliptic models: an elliptic Lotka-Volterra and FitzHugh-Nagumo model, and then an hypoelliptic Harmonic Oscillator and FitzHugh-Nagumo model.

Models and objectives

As explained in the introduction, the estimation procedure that we propose, based on the discrete optimal control theory, is fully developed for linear models. We thus introduce a linear SDE and the corresponding estimation criteria. In the simulation section, we will show an adaptation of the method for a non-linear SDE.

Elliptic and hypoelliptic stochastic differential equations

We consider a d-dimensional state variable Z t ∈ R d , d ≥ 2, defined for t in a time interval [0, T ]. We distinguish in the following the first observed state variable denoted V t ∈ R from the last d -1 other unobserved variables denoted U t ∈ R d-1 . The dynamic of (Z t = (V t , U t )) t≥0 is described by the following stochastic dynamical system:

dV t dU t = A θ (t) V t U t dt + B σ dW t (2)
with a known initial condition (V 0 , U 0 ) and where W t is a m dimensional Brownian motion.

The drift is assumed linear with respect to Z t = (V t , U t ). The d × d-matrix A θ (t) depends on the unknown parameter vector θ and may be a function of time t. The d × m-matrix B σ is called the diffusion coefficient and depends on an unknown parameter vector σ. We consider two cases:

• Elliptic SDE: m = d and the d × d matrix B σ is not singular: det(B σ B T σ ) > 0 (3)
• Hypoelliptic SDE: m < d and the matrix B σ is singular:

B σ = 0 m b σ , (4)
with 0 m the m-dimensional row vector of zeros and b σ a (d -1) × m-matrix such that det(b σ b T σ ) > 0 A noticeable feature in the hypoelliptic SDE is that the equation ruling V t does not contain a stochastic part. The matrix b σ models the way the stochastic disturbance acts on the unobserved variables U t of the system and indirectly on V t through U t .

In the hypoelliptic case, we consider the following assumption:

(H1) Let A 1j (t) denote the jth element of the first column of A θ (t) and B ,j the jth column of B σ . For any t, there exists at least one j = 2, . . . , d such that

A 1,j (t)B ,j = 0
Under assumption (H1), the noise is propagated also to the first coordinate V t . This property that the noise generates the entire space R d is a characteristic of the hypoelliptic property [see e.g. [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. Note that the standard assumptions to ensure the hypoellipticity are different from (H1) [see e.g. [START_REF] Samson | A contrast estimator for completely or partially observed hypoelliptic diffusion[END_REF], but it reduces to (H1) for a linear system.

Objectives and issues

The state variables Z t = (V t , U t ) are split in two because we observe only the first onedimensional state variable V t . We denote t -→ Y (t), the realization of V t . We assume that Y is discretely observed on the interval [0, T ] at times 0 = t 0 < . . . < t n = T without measurement error and denote (Y 0 , . . . , Y n ) these observations. We thus have

Y i = CZ ti , with C a 1 × d- matrix C = (1, 0 d-1
) T and 0 d-1 is the row vector of zeros of size d -1. The aim of the paper is to estimate the unknown parameters of model ( 2) in the elliptic and hypoelliptic cases using the discrete observations (Y 0 , . . . , Y n ). We will distinguish two cases: 1/ σ is known and only θ is estimated; 2/ both θ and σ are estimated.

Before introducing our approach based on the optimal control theory, let us explain why the estimation problem is difficult. When one wants to estimate θ and σ of an elliptic, partially observed SDE, the standard approach starts with the discretization of the diffusion (2). The Euler-Maruyama discretization scheme at time (t 1, , . . . t n ) is defined as follows for i = 0, . . . , n-1 and Z i = (V i , U i ):

Z i+1 = Z i + ∆ i A θ (t i )Z i + B σ η i =: A θ (t i )Z i + B σ η i (5) where Z i , U i , V i stand for Z ti , U ti , V ti , ∆ i = t i+1 -t i , A θ (t i ) = I d + ∆ i A θ (t i
) and the η i are independent variables distributed as N (0, ∆ i ). When the system is elliptic, that is when B σ B T σ is invertible, minus twice the log-likelihood of this discretized process, assuming both V t and U t discretely observed, is then

L Euler ((V i , U i ) i=0,...,n , θ, σ) = n-1 i=0 (Z i+1 -A θ (t i )Z i ) T (B σ B T σ ) -1 (Z i+1 -A θ (t i )Z i ) +n log(det(B σ B T σ )). ( 6
)
For partial observations, the log-likelihood (6) has to be integrated with respect to (U 0 , . . . , U n ).

The estimator is thus defined as arg min

θ,σ L Euler ((V i , U i ) i=0,...,n , θ, σ)d(U 0 , . . . , U n )
This integral can be viewed as a filtering problem: the unobserved trajectory (U 0 , . . . , U n ) is filtered with respect to the observations (V 0 , . . . , V n ). The Kalman filter can be applied for linear SDEs. In the non-linear case, the extended Kalman filter or a particle filter has to be used, at the price of an important computing time. For hypoelliptic SDE, the criteria ( 6) can not even be computed, B σ B T σ being not invertible. Therefore, there is a need to propose alternatives. In this paper, we take advantage of the optimal control theory to filter the unobserved coordinate or more precisely to find a surrogate value for the (unobserved) realization of the Brownian motion, under the form of a control sequence that drives the trajectory. This sequence allows to define an estimation criteria, that involves only the non-degenerate part of the stochastic noise. This criteria is thus suitable for elliptic as well as for hypoelliptic SDEs. The estimation procedure is presented in Sections 3-5.

Estimation of θ via optimal control theory

In this section, we assume σ known and fixed at its true value. We expose our estimation procedure: its main principle is to turn the statistical problem into an optimal control one. The procedure is the same whether B σ is singular or not. We start this section by deriving the general optimal control problem from which our estimator is defined. The optimal control problem is then solved by using a result coming from the linear quadratic theory (Section 3.2).

Optimal control problem and definition of θ

The principle of the control problem is to introduce a control sequence such that the data (Y 0 , . . . , Y n ) are close to a solution of the discretized model ( 5), close meaning of the order of the Euler scheme error [START_REF] Bally | The law of the euler scheme for stochastic differential equations i. convergence rate of the distribution function[END_REF]. In the context of SDEs, the natural control sequence is the increment of the Brownian motion. Let us introduce u i = W i+1 -W i , i = 0, . . . , n -1, which will play the role of the control value at time t i . Note that u i ∈ R m , i = 0, . . . , n -1 with a dimension m that can be different from d, thus allowing the noise to be degenerate.

The discretized model ( 5) can be reformulated under the form of a discrete controlled system:

V i+1 U i+1 = A θ (t i ) V i U i + B σ u i (V t0 , U t0 ) = (V 0 , U 0 ) (7)
We denote:

• u the vector of discrete values taken by the control: u = (u 0 , . . . u n-1 ).

• Z i,θ,σ,u = (V i,θ,σ,u , U i,θ,σ,u ), the solution of (7) corresponding to the given θ, σ and u.

When model ( 2) is true, there exists one realization of the Brownian motion such that the data (Y 0 , . . . , Y n ) are a sample of model (2). The control sequence u can be viewed as this specific realization of the Brownian motion. The objective is to infer this sequence u by extracting knowledge from the data. For that purpose, a cost function C(u, Y ; θ, σ) is defined such that the optimum in u corresponds to this realization of the Brownian motion. We say that the sequence u is designed through the cost function C(u, Y ; θ, σ).

A natural cost function is the conditional posterior probability P (u|Y ; θ, σ) of u knowing the data Y . Filtering consists in computing the conditional expectation of this posterior distribution. An alternative is to compute the maximum a posteriori (MAP) by maximizing P (u|Y ; θ, σ) with respect to u. That is the way we consider to compute the control in this paper. A direct computation of the MAP reveals to be numerically unstable or too strict for SDEs. Therefore, we introduce a weight to stabilize the optimization problem. Let us introduce the basic MAP and then the weighted cost function.

For a fixed value of the parameters (θ, σ), the MAP would be defined as follows: where P (Y |u; θ, σ) is the density of the data given the u, P (u; θ, σ) is the density of the Brownian motion and P (Y ; θ, σ), the likelihood, does not depend on u. Let us now detail the two terms that depend on u. The conditional density P (Y |u; θ, σ) is described by the error induced by the Euler discretisation of the SDE. It is well known that this error has a variance of order the time step [START_REF] Bally | The law of the euler scheme for stochastic differential equations i. convergence rate of the distribution function[END_REF]. The second term P (u; θ, σ) is the density of a discretized Brownian motion, which has a variance . Therefore, one can formulate the MAP as

ûMAP = arg max u - n i=0 (V i,θ,σ,u -Y i ) 2 i - n-1 i=0 u T i u i i = arg min u C(u, Y ; θ, σ) (8) with C(u, Y ; θ, σ) = n i=0 (V i,θ,σ,u -Yi) 2 i + n-1 i=0 u T i ui
i . The optimisation of the MAP (8) reveals to be too strict and intractable in practice. Therefore, we relax the constraints by introducing a weight between the two terms. To ease the writing of the cost function, we replace the scaling 1 i of the first term by a weight w > 0. The cost function of the control is finally defined as

C w (u, Y ; θ, σ) = w n i=0 (V i,θ,σ,u -Y i ) 2 + n-1 i=0 u T i u i i .
The weight w has to be chosen by the user. In Section 4 a procedure is proposed to select w adaptively from the data. To enter the optimal control theory, we exhibit the last observation Y n which plays a specific role in a control problem. The cost function can thus be rewritten:

C w (u, Y ; θ, σ) = w (V n,θ,σ,u -Y n ) 2 + n-1 i=0 w (V i,θ,σ,u -Y i ) 2 + u T i u i i .
For a given weight w, the best control is the sequence u that minimizes C w (u, Y ; θ, σ) under the constraint of model ( 7). It is defined as the solution of the following optimal control problem:

Minimize in u: C w (u, Y ; θ, σ) = w (V n,θ,σ,u -Y n ) 2 + n-1 i=0 w (V i,θ,σ,u -Y i ) 2 + 1 i u T i u i Subject to:    V i+1,θ,σ,u U i+1,θ,σ,u = A θ (t i ) V i,θ,σ,u U i,θ,σ,u + B σ u i (U 0,θ,σ,u V 0,θ,σ,u ) = (V 0 , U 0 ) . (9
) It is called the optimal control and denoted u θ,σ . For a fixed value of σ and given the discrete observations (Y 0 , . . . , Y n ), we then define the estimator of θ as:

θ w (σ) = arg min θ S w (Y ; θ, σ) (10)
where

S w (Y ; θ, σ) := min u C w (u, Y ; θ, σ) = C w (u θ,σ , Y ; θ, σ)
is the profiled cost C w over the set of possible sequences of the Brownian motion increments.

The computations of u θ,σ and S w (Y ; θ, σ) require to solve the optimal control problem (9). We thus have to prove the existence and uniqueness of the solution of problem (9), and that this unique solution is numerically computable to obtain u θ,σ and S w (Y ; θ, σ). These two results are given by the linear quadratic theory [START_REF] Trelat | Controle optimal : Theories et applications[END_REF][START_REF] Sontag | Mathematical Control Theory: Deterministic finite-dimensional systems[END_REF]] and are exposed below.

Linear quadratic theory

The linear quadratic theory is derived from the Pontryagin maximum principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF] in the restricted framework of linear models. For a given (θ, σ), the linear quadratic theory ensures, for the problem (9):

• the existence and uniqueness of u θ,σ , the solution of (9),

• that u θ,σ and S w (Y ; θ, σ) can be computed by the introduction of a backward finite difference equation, called the Riccati equation.

Let us be more precise. For a given (θ, σ), and a given weight w, let us denote the following matrices, for i = 0, . . . , n -1:

A i,θ = A θ (t i ) 0 T d 0 d 1 ∈ R d+1 × R d+1 , Q i = C T C -C T Y i -CY i Y 2 i w ∈ R d+1 × R d+1 , and
B 0 σ = B σ 0 m ∈ R d+1 × R m .
Set R n,θ,σ = Q n and let R i,θ,σ be the positive solution of the discrete backward Riccati equation, for i = n -1, . . . , 0:

R i,θ,σ = Q i + A T i,θ R i+1,θ,σ -R i+1,θ,σ B 0 σ 1 i I m + B 0 σ T R i+1,θ,σ B 0 σ -1 B 0 σ T R i+1,θ,σ A i,θ (11)
We can establish the following theorem. Theorem 1. The optimal control problem (9) has a unique solution u θ,σ = (u 0,θ,σ , . . . , u n-1,θ,σ ) with the value u i,θ,σ ∈ R m at time t i given by:

u i,θ,σ = - 1 i I m + B 0 σ T R i+1,θ,σ B 0 σ -1 B 0 σ T R i+1,θ,σ A i,θ V i,θ,σ,u θ,σ U i,θ,σ,u θ,σ . (12)
Moreover, the minimum value of C w is equal to:

S w (Y ; θ, σ) = min u∈R m×(n-1) C w (u, Y ; θ, σ) = C w (u θ,σ , Y ; θ, σ) (13) = (V 0 , U T 0 , 1) R 0,θ,σ (V 0 , U T 0 , 1) T The proof is given in Appendix.
Theorem 1 provides a closed expression for S w that depends on the matrix R 0,θ,σ which is easily obtained by solving the Riccati equation ( 11) and on (V 0 , U 0 ) the initial condition which is assumed to be known. The estimator of θ can thus be rewritten as

θ w (σ) = arg min θ (V 0 , U T 0 , 1) R 0,θ,σ (V 0 , U T 0 , 1) T
Note that using the expression of u θ,σ given by ( 12) in the discretized model (7), the system dynamics under the optimal control is given by:

Z i+1,θ,σ,u θ,σ = A θ (t i ) -B σ ( 1 i I m + B 0 σ T R i+1 B 0 σ ) -1 B 0 σ T R i+1 A i Z i,θ,σ,u θ,σ Z 0,θ,σ,u θ,σ = (V 0 , U 0 ) . (14)
Let us now enlighten the noticeable advantage of our approach. Theorem 1 holds for elliptic (m = d) SDEs as well as for hypoelliptic (m < d) SDEs. The estimation criteria does not involve the inverse of the matrix B σ B T σ , contrary to the standard estimation approaches for SDEs. Moreover, the linear quadratic theory itself does not require the matrix B σ B T σ to be full-conditioned. This can be explained by our choice of the cost function C ω (u, Y ; θ, σ), which is strictly convex with respect to the sequence u due to the term u T i u i . The cost function is thus automatically adapted to the dimension of the noise. Once again, it is noticeable that the solution to the optimal control problem is the same for the two classes of SDEs and that no specific calculations are required for the hypoelliptic case.

The relaxation made on the original cost C coupled with the Euler-Maruyama discretization scheme allows us to make our estimation problem fits in the framework of the discrete LQ theory. This ensures in turn the well posedness nature of the optimization problem (uniqueness, existence and continuity of the solution w.r.t problem parameters) defining the cost S w as well as an efficient method to compute it no matter the elliptic nature of the sde model. However, it comes with the price that the theoretical methods used for deriving consistency and convergence rate in Brunel and [START_REF] Brunel | A tracking approach to parameter estimation in linear ordinary differential equations[END_REF] cannot be applied anymore. In Brunel and [START_REF] Brunel | A tracking approach to parameter estimation in linear ordinary differential equations[END_REF], the control problem was a continuous one, the original cost C was not approximated and the original ODE model was not discretized, only perturbed by the addition of a control function u. Thus, the source of error was due to the noisy nature of observation and thus theoretical results have been derived by classic methods for M-estimator (see [START_REF] Van Der | Asymptotic Statistics. Cambridge Series in Statistical and Probabilities Mathematics[END_REF] for example). Nevertheless, our procedure reveals to be efficient in practice, see Section 6.

The next step is the selection of the weight w involved in the definitions of the optimal control problem (9) and the estimator θ w (σ). This is explained in the next section.

Adaptive w selection

The choice of the weight w is critical to ensure the stability of the procedure. We propose to select it adaptively by minimizing an external functional criteria G(Y ; σ, w, θ) that depends on w, the data and the parameters θ and σ. We propose hereafter two possible choices for G: G (1) uses a property of the average quadratic growth of a Brownian process and gives importance to the term of adequacy of the control sequence u; G (2) relies on a new contrast that gives more importance to the data fidelity term.

For the two functionals G (1) and G (2) , the selected weight w and the final estimator of θ are computed through the following steps. Definition 1 (Estimator θ(σ), σ known). Assume σ known. The estimation of θ is defined as follows:

1. For any weight parameter w, compute the estimator θ w (σ) = arg min θ S w (Y ; θ, σ). The two functionals are presented in the next sections.

Minimize the functional

G = G (1) or G = G (2)

w selection via a quadratic growth moment condition

We now present the first functional G (1) . It comes from a constraint linked to the distribution of the control. The control u θ,σ mimics the increments of the Brownian motion. Thus, it is natural to impose that the values u i,θ,σ √ i are independent and distributed with N (0, I m ). We can then derive constraints related to this distribution. For example, let us denote u i,j,θ,σ the j-th component of the m dimensional Brownian motion at time t i , for j ∈ 1, m . The law of large number implies that, almost surely,

1 n n-1 i=0 u 2 i,j,θ,σ i -→ 1.
We thus propose the following functional to select w:

G (1) (Y ; σ, w, θ) = m j=1 n-1 i=0 u 2 i,j,θ,σ i -n 2 . (17)
Here, the dependency of G (1) (Y ; σ, w, θ) w.r.t w is made through the optimal control sequence u θ,σ . This functional gives importance to the control sequence property.

w selection via a data fidelity criteria

The second functional is driven by the idea of giving importance to the data fidelity term. We can show by recurrence that for i = 0, . . . , n -1 and any integer 0 ≤ l ≤ i -1:

Z i+1,θ,σ,u = l j=0 A θ (t i-j ) Z i-l,θ,σ,u + l k=0 k-1 j=0 A θ (t i-j ) B σ u i-k
where k-1 j=0 A θ (t i-j ) = 1 when k = 0. Thus, by multiplying the left and right hand side by C, we can explicitly link the observations Y and the control:

Y i+1 -C l j=0 A θ (t i-j ) Z i-l,θ,σ,u = C l k=0 i-k k-1 j=0 A θ (t i-j ) B σ u i-k √ i-k . The functional G (2) is based on the distribution of Y i+1 -C l j=0 A θ (t i-j )Z i-l,θ,σ,u .
To ease the reading, let us introduce the m-vector

Γ θ,σ (k, t i ) = C i-k k-1 s=0 A θ (t i-s ) B σ
and the sequence of scalars, for i = 0, . . . , n -1:

X i,θ,σ,u = Y i+1 -C l j=0 A θ (t i-j )Z i-l,θ,σ,u . (18)
We can derive the law followed by

X i,θ,σ,u = l k=0 Γ θ,σ (k, t i ) u i-k √ i-k : X i,θ,σ,u ∼ N (0, γ 2 l,θ,σ (t i )) with γ 2 l,θ,σ (t i ) = l k=0 Γ θ,σ (k, t i ) Γ θ,σ (k, t i ) T .
Let us choose l such that γ 2 l,θ,σ (t i ) is non 0. Let us explain the intuition why choosing a sufficiently large index l yields γ 2 l,θ,σ (t i ) > 0. The idea is to give time to the stochastic elements u i to diffuse. In other words, after "enough time" (characterized by l), the elements u i are able to perturb the observations Y i+l+1 and σ becomes univocally identified. Now set X l i,θ,σ,u := ti) . Then E X l i,θ,σ,u 2 = 1. Note that X l i,θ,σ,u depends on (u i-l , . . . , u i ) and X l i+l+1,θ,σ,u depends on (u i+1 , . . . , u i+l+1 ) thus they are independent. For the i.i.d sequence of X l (l+1)i+1,θ,σ,u θ,σ 0≤i≤ n-1 l+1 of size L := n+l l+1 , the law of large number implies that, almost surely,

X i,θ,σ,u γ l,θ (
1 L L-1 i=0 X l (l+1)i+1,θ,σ,u 2 -→ 1.
We can now define the functional G (2) as:

G (2) (Y ; σ, w, θ) = L-1 i=0 X l (l+1)i+1,θ,σ,u θ,σ 2 -L 2 . ( 19
)
In practice, we choose the smallest value l which ensures γ 2 l,θ,σ (t i ) > 0 in order to use the largest sequence X l (l+1)i+1,θ,σ,u θ,σ 0≤i≤ n-1 l+1 .

• For elliptic SDEs, it is sufficient to take l = 0. Indeed Γ θ,σ (0,

t i ) = C √ i B σ = 0 m since B σ is of full rank.
• For hypoelliptic SDEs, we have Γ θ,σ (0,

t i ) = C √ i B σ = 0 m . With l = 1, we have γ 2 1,θ,σ (t i ) = Γ θ,σ (1, t i ) Γ θ,σ (1, t i ) T = i-1 CA θ (t i )B σ B T σ A θ (t i ) T C T and assumption (H1) implies γ 2 1,θ,σ (t i ) > 0
The computations of γ 2 0,θ,σ (t i ) and γ 2 1,θ,σ (t i ) are detailed in Section 6 for some examples of SDEs.

Estimation of σ

So far, the vector σ of the parameters involved in the diffusion coefficient B σ has been considered known. We now present its estimation. The standard estimator of a diffusion coefficient is the quadratic variation of the stochastic process. It can however not be applied in the context of partial observations because only the first coordinate is observed. We therefore propose to minimize an estimation criteria. It turns out that the two functionals G (1) and G (2) presented in Section 4 can be used to estimate σ. It might seem counter-intuitive that the same criteria could be used to estimate a weight parameter and a diffusion coefficient. However, let us recall that the two functionals have been constructed as constraints of the optimal control problem. These constraints are not specifically linked to the weight parameter but to the model itself. They can thus be used judiciously to estimate also the diffusion coefficient B σ .

Let us now present the procedure. It is a nested procedure that provides the final estimation of θ, σ and the data-driven selection of the weight w. Definition 2 (Estimators σ and θ). The estimation of (θ, σ) is defined as follows:

1. For any weight parameter w and diffusion coefficient σ, compute the estimator θ w (σ) = arg min θ S w (Y ; θ, σ). 2) with respect to σ, using the plugged value θ w (σ), and define σ(w) 2) with respect to w, using the plugged value θ w (σ), and define

Minimize the functional

G = G (1) or G = G (
:= arg min σ G(Y ; σ, w, θ w (σ)) (20) 3. Minimize the functional G = G (1) or G = G (
w := arg min w G(Y ; σ(w), w, θ w ( σ(w))) (21)
4. Define the final estimators of θ and σ as

σ := σ( w) (22)
θ := θ w ( σ)
The two functions G (1) and G (2) can be used in this nested procedure.

Simulation study

We want to evaluate the two estimation procedures given in definition 1 when σ is known and definition 2 when σ is also estimated. A simulation study is conducted on several SDEs, elliptic (Section 6.1) and hypoelliptic (Section 6.2). For each model, a hundred trajectories with n = 1 000 equidistant points are simulated. Then the mean and standard error (SE) of the estimators are computed. The mean required CPU time (computed with the MATLAB function cputime) is also reported. When σ is known, two estimators are computed using the nested procedure given in definition 1. We denote θ (1) (σ) ( θ (2) (σ), respectively) the estimator defined by ( 16) with the weight w adaptively selected by choosing the functional G = G (1) (G = G (2) , respectively). When σ is also estimated, two estimators are also computed. We denote θ (1) and σ(1) ( θ (2) and σ(2) , respectively) the estimators given in definition 2 using G = G (1) (G = G (2) , respectively).

Three models are used: the Lokta-Volterra model, the FitzHugh-Nagumo model and the Harmonic Oscillator (HO). The two firsts are non linear and HO is linear. As explained in the introduction, the procedures that we propose apply to linear SDE only. This allows to take advantage of the explicit solution to the optimal control problem thanks to the linear-quadratic theory. For the non-linear SDE, we thus propose strategies to bypass this limit and adapt our estimators to these non-linear models.

Simulation study for elliptic systems

Two models are studied in their elliptic version, the Lotka-Volterra system and the FitzHugh Nagumo model.

Elliptic Lotka-Volterra process

Stochastic versions of the Lotka-Volterra model have been proposed, for examples by [START_REF] Lotka | Elements of Mathematical Biology[END_REF], [START_REF] Meeds | Optimization monte carlo: Efficient and embarrassingly parallel likelihood-free inference[END_REF], [START_REF] Graham | Asymptotically exact inference in differenciable generative models[END_REF], [START_REF] Mao | Environmental brownian noise suppresses explosions in populations dynamics[END_REF][START_REF] Mao | Asymptotic behaviour of the stochastic lotka-volterra model[END_REF]. Originally, it is a predator-prey model that describes the dynamics of biological systems with two populations, named predator and prey, which interact together. The preys are assumed to have an unlimited food supply. The perturbed model is defined through a system of two differential equations disturbed by a stochastic noise. Let V t and U t denote the number of predators and preys at time t, respectively. The dynamics of V t and U t is described as:

dV t = (-θ 1 V t + θ 2 V t U t )dt + σdW 1,t dU t = (θ 3 U t -θ 4 V t U t ) dt + σdW 2,t (V 0 , U 0 ) = (0.5, 1) (23)
where θ 1 is the death rate of the predator, θ 2 is the growth rate of the predator population, θ 3 is the exponential growth of the prey, θ 4 is the rate of predation upon the prey which reflects the interaction between the two species, (W 1,t ) t≥0 and (W 2,t ) t≥0 are two independent Brownian motions, and σ is the diffusion coefficient assumed to be the same for both coordinates. The system is thus elliptic.

To simulate the trajectories, a Euler-Maruyama discretization scheme is used with = 0.01. The parameters are set to (θ 1 , θ 2 , θ 3 , θ 4 ) = (0.5, 0.5, 0.5, 0.5) and σ = 0.2. An example of trajectory for the model ( 23) is presented in Figure 1.

The objective is to estimate these parameters from discrete observations of the predator population only. For the sake of parametric identifiability, θ 2 is considered as known.

The system ( 23) is non linear due to the interaction term between predators and preys. The estimation procedure proposed in Section 3 can thus not be applied directly. To deal with the nonlinear interaction term V t U t , we take advantage of the discrete observations Y of the first coordinate. The estimation procedure only uses the discretized version of the continuous model. Instead of a discretized version of model ( 23), we consider a discrete model where the interaction term

V i U i is replaced by Y i × U i . V ti+1 = V ti + ∆ i (-θ 1 V ti + θ 2 Y i U ti ) + σ W 1,ti+1 -W 1,ti U ti+1 = U ti + ∆ i (θ 3 U ti -θ 4 Y i U ti ) + σ W 2,ti+1 -W 2,ti (V 0 , U 0 ) = (0.5, 1) . ( 24
)
The estimation procedures are then applied to the model ( 24). We have B σ = σ 0 0 σ . Note that to ensure that G (2) is well defined, it is sufficient to take l = 0, as γ 2 0,θ,σ (t i ) = i σ 2 > 0. Results are presented in Table 1. The four procedures give overall good results. The two procedures with σ known are very similar, in terms of mean, SE of the estimators and computation time. When σ is also estimated, the second method performs better, especially to estimate the diffusion parameter σ. This could be explained by the fact that the G (2) -based procedure relies more on the observations {Y i } than G (1) . As the diffusion coefficient is the same for both coordinates, the estimation of σ is more accurate with G (2) . Note that the computation time of the procedures estimating the σ's is about hundred times larger for G (1) and twenty times larger for G (2) .

θ 1 θ 3 θ 4 σ CPU True values
0.5 0.5 0.5 0.2 time (s) σ known, estimation procedure given in Definition 1 θ(σ) (1) 0.64 (0.34) 0.58 (0.15) 0.62 (0.40) -158.3 θ(σ) (2) 0.65 (0.35) 0.59 (0.17) 0.64 (0.42) -158.5 σ unknown, estimation procedure given in Definition 2 θ (1) (2) . Mean individual CPU times are given for each estimation procedure.

Elliptic FitzHugh-Nagumo

The FitzHugh-Nagumo model (FHN) describes the dynamic of an excitable neuron, modeling the characteristic spikes of the neuron [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF]. It is defined by a system of two differential equations, that we consider perturbed by a random noise. The two coordinates model the membrane potential of the neuron and a recovery channel mimicking the opening/closing of ion channels. More formally, let V t and U t denote the membrane potential at time t and the value of the recovery channel at time t, respectively. The stochastic FHN model is defined as:

dV t = 1 V t -V 3 t -U t dt + σ 1 dW 1,t dU t = (γV t -U t + β) dt + σ 2 dW 2,t (V 0 , U 0 ) = (0, 0) (25)
where is a time scale parameter, γ and β are kinetic parameters, σ 1 and σ 2 are the two diffusion coefficients, (W 1,t ) t and (W 2,t ) t are two independent Brownian motions. Only the membrane potential V t can be measured experimentally.

The trajectories are simulated with a Euler-Maruyama scheme with = 0.02. The parameters are set to = 0.1, γ = 1.5, β = 0.8, σ 1 = 0.1, and σ 2 = 0.3. An example of trajectory for the model ( 25) is presented on Figure 2. The time scale is difficult to estimate and is considered as known (and set to 0.1) in the following. The objective is to estimate θ = (γ, β, σ 1 , σ 2 ) from the discrete observations of the first coordinate (Y 1 , . . . , Y n ).

The estimation procedures require a homogeneous linear model. As done for the Lotka-Volterra model, we will propose a linear model close to (25) and apply the estimation procedures to that linearized model. In FHN model, there is a non-linear term (V 3 ) in the first equation and a non-homogeneous term (β) in the second equation. Both have to be replaced. The trajectory V can be replaced by Y , the observed realization of V . As we will see, several choices can be done to replace the cubic term. For the non-homogeneous term, a strategy is to introduce a new variable R t assumed to be constant equal to 1. This can be done by increasing the number of coordinates by one and stating dR t = 0.

A first attempt to propose a linear model consists in replacing V 3 by Y 3 in a discretized version of (28). But this reveals to be very unstable, due to the dramatic propagation of the discretization error committed by using Y 3 instead of V 3 . An alternative is replacing V -V 3 by (1 -Y 2 )V . The non-homogeneous term is replaced by a new variable Z t . We then consider the following discretized linear system:

V ti+1 = V ti + i 1 1 -Y 2 i V ti -U ti + σ 1 (W 1,ti+1 -W 1,ti ) U ti+1 = U ti + i (γV ti -U ti + βR i ) + σ 2 (W 2,ti+1 -W 2,ti ) R ti+1 = 1 (V 0 , U 0 , R 0 ) = (0, 0, 1) (26)
The estimation procedures are then applied to the model ( 26) with d = 3, C = ( 1 0 0 ) and

B σ =   σ 1 0 0 σ 2 0 0   .
To estimate σ(2) defined by ( 19), note that again the index l = 0 is sufficient as γ 2 0,θ,σ (t i ) = i σ 2 1 > 0. Results are presented in Table 2. When σ = (σ 1 , σ 2 ) is known, the two estimators θ(σ) (1) and θ(σ) (2) are very close. This is not the case when σ is also estimated. This could be explained by the fact that the G (1) criteria uses the 2 dimensional control sequence in a direct way to construct an estimator of (σ 1 , σ 2 ). At the contrary the G (2) criteria only uses the observation to estimate the diffusion, the control sequence being used to estimate the parameters of the unobserved coordinate. As expected, the computational time is much larger when σ is unknown.

γ β σ 1 σ 2 CPU true value 1.5 0.8 0.1 0.3 time (s) σ known, estimation procedure given in Definition 1 θ (1)
1.52 (0.16) 0.82 (0.12) --84.3 θ (2) 1.53 (0.17) 0.82 (0.12) --86.3 σ unknown, estimation procedure given in Definition 2 θ (1) and σ(1) 1.51 (0.15) 0.81 (0.10) 0.06 (0.02) 0.18 (0.10) 12746 θ (2) and σ(2) 1.57 (0.13) 0.85 (0.11) 0.04 (0.01) 0.12 (0.12) 2800

Table 2: Estimation results for the elliptic FitzHugh-Nagumo model. From a hundred simulated trajectories, mean and standard error (in brackets) of the different estimators of the estimators given in Definition 1 (σ known) and Definition 2 (σ unknown) obtained with the two functionals G (1) and G (2) . Mean individual CPU times are given for each estimation procedure.

Simulation study for hypoelliptic systems

We now compare the estimation procedures on two hypoelliptic systems, the Harmonic Oscillator and a hypoelliptic FitzHugh-Nagumo neuronal model. As explained before, the estimation procedures we propose are robust to the hypoellipticity and we thus apply the same four procedures.

Harmonic oscillator

The Harmonic Oscillator is a mechanistic model describing oscillations governed by a white noise. It is described by a system of two equations, denoted V t and U t , the noise entering only the second equation [START_REF] Pokern | Parameter estimation for partially observed hypoelliptic diffusions[END_REF]. The model is defined as follows:

   dV t = U t dt dU t = (-DV t -δU t )dt + σdW t (V 0 , U 0 ) = (0, 0) (27)
with D, δ, σ > 0 and (W t ) t a Brownian motion.

The trajectories are simulated with a Euler-Maruyama scheme with = 0.02. The parameters are set to D = 4, δ = 0.5 and σ = 0.5. The objective is to estimate θ = (D, δ, σ) from the discrete observations of the first coordinate (Y 1 , . . . , Y n ).

The four estimation procedures are applied. Again, let us give some details on the procedure for obtaining σ(2) . We have

B σ = 0 σ and A θ (t i ) = 1 + i i -i D 1 -i δ . So Γ θ,σ (0, t i ) = 0, Γ θ,σ (1, t i ) = 3 2 i σ.
We can deduce that γ 2 1,θ (t i ) = 3 i σ 2 = 0. This order of variance for V t is the true one, as explained in [START_REF] Samson | A contrast estimator for completely or partially observed hypoelliptic diffusion[END_REF]. Thus our procedure automatically propagates the noise from the second coordinate to the first one, yielding to a invertible covariance matrix of the process. In that sense, it is close to what has recently been proposed by [START_REF] Ditlevsen | Hypoelliptic diffusions: discretization, filtering and inference from complete and partial information[END_REF].

Results are presented in Table 3. When σ is known, the first approach gives better results than the second. Especially δ is estimated with bias by θ (2) and not by θ (1) . This is the opposite when σ is also estimated. The second procedure performs better for δ and more convincingly for σ which is accurately estimated. As noted for the elliptic case, the computation time is larger when σ is estimated than when it is fixed at its true value. (1)

4.12 (0.44) 0.52 (0.13) -119.3 θ(σ)

(2) 4.06 (0.51) 0.16 (0.07) -117.5 σ unknown, estimation procedure given in Definition 2 θ (1) and σ(1) 4.11 (0.51) 0.10 (0.12) 0.15 (0.09) 2869 θ (2) and σ(2) 4.03 (0.50) 0.19 (0.15) 0.50 (0.35) 3067

Table 3: Estimation results for the hypoelliptic Harmonic Oscillator model. From a hundred simulated trajectories, the mean and standard error (in brackets) of the estimators given in Definition 1 (σ known) and Definition 2 (σ unknown) obtained with the two functionals G (1) and G (2) . Mean individual CPU times are given for each estimation procedure.

Hypoelliptic FitzHugh-Nagumo

We now consider the FHN model without noise on the first coordinate, as studied by [START_REF] Ditlevsen | Hypoelliptic diffusions: discretization, filtering and inference from complete and partial information[END_REF]. The model is thus defined by: and we refer to Section 6.1.2 for the description of the variables and parameters.

dV t = 1 V t -V 3 t -U t dt dU t = (γV t -U t + β) dt + σdW t (V 0 , U 0 ) = (0, 0) (28)
The trajectories are simulated with a Euler-Maruyama scheme with = 0.02. The parameters are set to = 0.1, γ = 1.5, β = 0.8, and σ = 0.3. An example of trajectory for the model ( 28) is presented in Figure 3. The time scale is difficult to estimate and as in the elliptic case, we consider it known (and set to 0.1) in the following. The objective is to estimate θ = (γ, β, σ) from the discrete observations of the first coordinate (Y 1 , . . . , Y n ).

Model ( 28) is non-linear and time-inhomogenous. Thus, we have to consider a linear version, as done in the elliptic case. For the instability reason evoked in Section 6.1.2, we directly consider the linear discretized model:

V ti+1 = V ti + i 1 1 -Y 2 i V ti -U ti U ti+1 = U ti + i (γV ti -U ti + βR ti ) + σ(W ti+1 -W ti ) R ti+1 = 1 (V 0 , U 0 , R 0 ) = (0, 0, 1) (29)
The four estimation procedures are applied on model (29). Let us give some details on the procedure σ(2) defined by ( 19). We have C = ( 1 0 0 ),

B σ =   0 σ 0   and A θ (t i ) =   1 + i 1 -Y (t) 2 -i 0 i γ 0 β i 0 0 1   .
So Γ θ,σ (0, t i ) = 0, Γ θ,σ (1, t i ) = -3 2 i σ . We can deduce from that γ 2 1,θ (t i ) = 3 i σ 2

ε 2 = 0. This corresponds to the first term of the exact variance of V t , as proved by [START_REF] Ditlevsen | Hypoelliptic diffusions: discretization, filtering and inference from complete and partial information[END_REF].

Results are presented in Table 4. The most striking and somewhat counter-intuitive result in the hypoelliptic case is the required computational time, in particular when σ needs to be estimated. This computational time is smaller than in the elliptic case. Indeed, for a fixed σ, the complexity of the optimization problem linked to the drift parameter estimation θ does not depend on the singularity of B σ . It only relies on the dimension of θ and the dimension of the control space. So the smaller the dimension of σ, the quicker the minimizer of S is found. Since the dimension of σ is smaller in the hypoelliptic example than in the elliptic one, it explains the observed difference in the computational time. (1)

1.49 (0.10) 0.81 (0.08) -147.9 θ(σ)

(2)

1.49 (0.10) 0.81 (0.08) -152.2 σ unknown, estimation procedure given in Definition 2 θ (1) and σ(1) 1.46 (0.34) 0.81 (0.28) 0.38 (0.13) 3151 θ (2) and σ(2) 1.48 (0.12) 0.79 (0.10) 0.29 (0.20) 3822

Table 4: Estimation results for the hypoelliptic FitzHugh-Nagumo model. From a hundred simulated trajectories, mean and standard error (in brackets) of the different estimators given in Definition 1 (σ known) and Definition 2 (σ unknown) obtained with the two functionals G (1) and G (2) . Mean individual CPU times are given for each estimation procedure.

Conclusion and discussion

In this work, we propose a new method based on control theory to estimate parameters in SDEs. Its main feature is to propose a unified framework for both elliptic and hypoelliptic models by using a criteria focusing on estimating the Brownian motion realization given the observation rather than solely on the observation. By doing so, we manage to construct a criteria S w well defined no matter the structure of B σ . Another interest of the method is the computational time. The use of the discrete LQ theory allows to avoid the use of MCMC and/or stochastic approximation computationally costful. These two points enlighten the benefit of using control theory for dealing with statistical problems. However, because of the nested structure of our estimation procedure, we can see that the computational time drastically increases when σ has to be estimated. In particular it is very sensitive to the dimension of σ. This can be a limitation for high dimensional systems. Interestingly, at the contrary to classic statistical approaches, the hypoelliptic nature of a system is an advantage for us because the dimension of σ is smaller than the one in its elliptic counterpart. A perspective could be to investigate a criteria which allows us to simultaneously estimate θ and σ instead of using our current nested procedure. However, so far we did not manage to find a criteria which fits in the framework of the discrete LQ theory. That leads us to our second limitation, our method is currently restricted to linear models or non-linear ones with a specific structure (i.e models only nonlinear w.r.t the observed state variable). Our main challenge would be to extend our method to general nonlinear SDEs.

  |u; θ, σ)P (u; θ, σ) P (Y ; θ, σ)

  with respect to w, using the plugged value θ w (σ), and define w :the final estimator of θ as θ(σ) = θ w (σ) (16)

Figure 1 :

 1 Figure 1: Lotka-Volterra model: example of a simulation with the dynamics along time of the number of predators (green plain line) and the number of preys (blue dotted line).

Figure 2 :

 2 Figure 2: Elliptic FitzHugh-Nagumo model: example of a simulation with the dynamics along time of the membrane potential (left) and the recovery variable (right)

Figure 3 :

 3 Figure 3: Hypoelliptic FitzHugh-Nagumo model: example of a simulation with the dynamics along time of the membrane potential (left) and the recovery variable (right)

Table 1 :

 1 Estimation results for the Lotka Volterra model. From a hundred simulated trajectories, mean and standard error (in brackets) of the estimators given in Definition 1 (σ known) and Definition 2 (σ unknown) obtained with the two functionals G (1) and G

	and σ(1)	0.68 (0.35) 0.57 (0.15) 0.65 (0.45) 0.16 (0.01)	11405
	θ (2) and σ(2)	0.68 (0.35) 0.59 (0.16) 0.65 (0.44) 0.19 (0.01)	2438

Acknowledgements

Adeline Samson has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissements davenir.

A Discrete LQ theory: Main theorem

In section 3, the theorem which gives us the profiled cost value, as well as the corresponding control sequence, is a particular case of a more general theorem. This theorem, fundamental in LQ theory, ensures the existence, uniqueness and gives a closed form for the global minimizer of a cost under the form:

where u = {u 0 , . . . u N -1 } and the state variable sequence x = {x 0 , . . . , x N } are linked by the finite difference equation:

The derivation of the next theorem by the optimality principle can be found in [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF].

Theorem 2. Let us assume that S is positive semi-definite, for all i ∈ 0, N -1 , Q i is positive semi-definite and R i is positive definite. Then the cost (30) reaches its global minimum for the control sequence u * given by:

x k and the minimal cost value is equal to:

where P k is given by the discrete time Riccati difference equation:


Veronica Karina Quintuna Rodriguez
email: veronicakarina.quintunarodriguez@orange.com

Fabrice Guillemin
email: fabrice.guillemin@orange.com

Performance Analysis of Resource Pooling for Network Function Virtualization

Keywords: Scheduling, virtualization, VNF, cloud computing, fog computing. BBU EPC sub-function computing resources Scheduler

In the framework of network function virtualization, we consider in this paper the execution of Virtualized Network Functions (VNFs) in data centers whose computing capacities are limited. We assume that each VNF is composed of subfunctions to be executed on general purpose hardware, each subfunction requiring a random amount of processing time. Because of limited processing capacity, we investigate the relevance of resource pooling where available cores in a data center are shared by active VNFs. We study by simulation various algorithms for scheduling sub-functions composing active VNFs (namely Greedy, Round Robin and Dedicated Core algorithms). We additionally introduce an execution deadline criterion, which means that VNFs can renege if their sojourn time in the system exceeds by a certain factor their service time. This feature is especially relevant when considering the processing of real-time VNF. Simulations show that sub-functions chaining is critical with regard to performance. When sub-functions have to be executed in series, the simple Dedicated Core algorithm is the most efficient. When sub-functions can be executed in parallel, Greedy or Round Robin algorithms offer similar performance and outperform the Dedicated Core algorithm. Enabling as much as possible parallelism and avoiding chaining when designing a VNF are fundamental principles to gain from the available computing resources.

I. INTRODUCTION

The emergence of virtualization technology plays a crucial role in the evolution of telecommunications network architectures, notably by enabling the virtualization of network functions. While such functions were so far implemented on proprietary and closed platforms offered by equipment providers, virtualization allows network functions to be implemented on Commercial off-the-shelf (COTS) servers. This is clearly a great revolution in the design of networks and IT infrastructures, which can eventually be completely merged.

For today's network operators, virtualization promises more flexibility and scalability as well as resource efficiency in the deployment of network services [START_REF] Aleksic | Network virtualization: Paving the way to carrier clouds[END_REF]. In addition, the use of virtualization techniques allows the network operator to deal with traffic growth without constantly investing in costly infrastructure, since Virtualized Network Functions (VNFs) can be hosted in traditional data centers and customized according to specific needs.

In this context, a VNFs is a piece of software which can even be available as an open source program and executed on general purpose processors (e.g., ARM processors, x86 servers) in multi-platform environments (i.e., operating systems); this guarantees portability and compatibility. Furthermore, virtualization techniques and data center processing enable on-demand provisioning, elasticity and resource pooling. In terms of flexibility and cost efficiency, network functions should not be simply virtualized but also redesigned to gain from the available computing resources.

An end-to-end network service (e.g., video streaming, mobile voice) can be represented as a forwarding graph of network functions, which is commonly referred to as a service chaining [START_REF]Network Functions Virtualisation ETSI. Architectural framework. Technical report[END_REF]. Today's hardware-based approach makes services implementation extremely complex and time-consuming [START_REF]Openflow-enabled sdn and network functions virtualization[END_REF]. In a virtualized environment, a VNF Forwarding Graph (VNF FG) runs on the top of the virtualization layer and can be managed more efficiently [START_REF]Network Functions Virtualisation ETSI. Architectural framework. Technical report[END_REF]. For instance, to update a network service, it is enough to add a new VNF on a virtual machine; in the same way, to scale the network, it is only necessary to instantiate the underlying VNF on the computing facility.

Furthermore, each VNF consists of a number of software components, referred to in the following as sub-functions. The flexibility brought by virtualization should be exploited to realize modular and parallel sub-functions. We assume that cores of a computing facility are assigned for executing sub-functions according to a specific scheduling algorithm. This algorithm plays a key role in the system performance, especially when there are real-time constraints or deadlines for the execution of a particular VNF.

In this paper, we shall consider three scheduling algorithms, referred to as Round Robin, Greedy and Dedicated Core. The goal of this contribution is to identify which algorithm is the most efficient in terms of execution time of a VNF and reneging rate (when a deadline for the execution of a VNF applies). In like manner, we aim to determine the impact of VNF design in the global performance of a network function.

The organization of this paper is as follows: Potential use cases are presented in Section II. In Section III, we introduce the model and performance measures to characterize the efficiency of the system. Scheduling algorithms are described in Section IV. Simulation results of scheduling algorithms with or without reneging are discussed in Section V. Finally, concluding remarks are presented in Section VI.

II. PROBLEM STATEMENT AND USE CASES

Network Function Virtualization (NFV) approach enables the deployment of a large number of VNFs, which can be executed on cloud-computing as well as fog computing environments. While large centralized data centers are commonly operated by the cloud framework, fog computing enables data center dissemination on the edge of the network infrastructure. The advantage of small data centers disseminated at the network edge is that they are close to end-users and ensure lower latency and better customer experience. The key difference between fog and cloud computing is in the limitation of computing capacity [START_REF] Stojmenovic | The fog computing paradigm: Scenarios and security issues[END_REF].

This work pays special attention to the execution of VNFs in data centers with limited capacity in terms of computing. We more precisely investigate the feasibility of resource pooling. We consider a computing facility composed of a number of cores, which can be dynamically allocated for the execution of sub-functions of an active VNF.

Functions of mobile core network (e.g., those of the Evolved Packet Core (EPC) of 4G mobile or future 5G networks) are good candidates for virtualization in the form of virtualized EPC (vEPC), Cloud-EPC (C-EPC) or EPC as a Service (EPCaaS) [START_REF] Taleb | Ease: Epc as a service to ease mobile core network deployment over cloud[END_REF], [START_REF] Pentikousis | Mobileflow: Toward software-defined mobile networks[END_REF]. As shown in Fig. 1 several functions can be handled by software, such as, session setup, user authentication and access authorization, currently handled by the Mobility Management Entity (MME). This is also the case of packet filters, policy control and charging, supported until now by P-Gateway as well as mobility functions (e.g., in case of handover between eNodeBs or between LTE and other 3GPP access) held by S-Gateway. Each VNF can be executed on Virtual Machines (VM) or containers. In this way, operators can instantiate each EPC component in order to cope with the incoming traffic demands. On the other hand, it is also possible to optimize resources during non-peak hours.

Moreover, a virtual EPC could be instantiated for a specific need, e.g., a company which is willing to operate its own private mobile network for its employees. The same applies for other traditional network functions, e.g., IP Multimedia Subsystem (IMS) for private Voice over IP (VoIP) networks.

In spite of the promises of virtualization, several challenges should be dealt with in future implementations, particularly because some network functions have strict requirements in terms of real-time processing. This is notably the case of Radio Access Network virtualization for short vRAN. For instance, in 4G mobile networks, the Time Division Duplex (TDD) process requires a response ACK/NACK (i.e., ACK to proceed to the transmission of a new frame, or NACK to attempt a retransmission) to the User Equipment (UE) or eNodeB within 3ms after a frame is arrived, where 2ms are available for reception process and 1ms for transmission process [START_REF] Chih-Lin | Recent progress on c-ran centralization and cloudification[END_REF]. This operation named Hybrid Automatic Repeat-Request (HARQ) is performed at the MAC level, after demodulation and decoding are done. If successful decoding is not possible the packet re-transmission is scheduled. Nevertheless, this fact degrades the customer experience, because it reduces the achievable UE peak rate [START_REF] Wubben | Benefits and impact of cloud computing on 5g signal processing: Flexible centralization through cloud-ran[END_REF].

Considering the realization of vRAN, traditional eNodeBs with integrated radio and base-band processing are replaced with shared processing and distributed radio elements, i.e., Remote Radio Head (RRH). So, virtual Base Band Unit (BBU) sub-functions are located at a centralized site, named below vBBU, where a pool of computing resources could be dynamically allocated based on traffic conditions. In addition, BBU virtualization can dramatically improve the global network performance; since several radio elements are handled by a single BBU pool, inter-cell coordination facilitates cooperative multi-point processing and massive MIMO implementation while avoiding interference [START_REF] William Dawson | On the benefits of ran virtualisation in c-ran based mobile networks[END_REF].

In this way, the entire BBU functionality [START_REF] Nikaein | Processing radio access network functions in the cloud: Critical issues and modeling[END_REF], [START_REF] Wang | Potentials and challenges of c-ran supporting multi-rats toward 5g mobile networks[END_REF] including RAN L1, L2, and L3 protocol layers could be represented as a suite of virtual sub-functions available for treating Physical Resource Blocks (PRB), as illustrated in Fig. 1. Modulation and encoding process as well as data convergence process (e.g., packet compression and encryption) and radio bands allocation might be considered as virtual sub-functions of the base-band processing module. Thus, the processing load mainly depends of the Modulation and Coding Scheme (MCS), the execution time of these sub-functions is determined by the channel conditions between the eNodeB and the current UE. Depending on which value is reported by the UE, network transmits data with different MCS. If network gets high Channel Quality Indicator (CQI) from UE, (i.e., good radio quality) it transmits the data with high-order modulation and a simple channel coding scheme. Likewise, for a low CQI, the network constructs the transmission block with a robust MCS. All sub-functions except modulation/demodulation and encoding/decoding might have small run time enabling the possibility of sharing computing resources [START_REF] Wubben | Benefits and impact of cloud computing on 5g signal processing: Flexible centralization through cloud-ran[END_REF].

Decoding function (e.g., LDPC, turbo codes) is the most complex in terms of processing load, since its number of iterations depends on the signal quality. Therefore, decoders are promising candidates for parallelization to meet real-time constraints; for instance, multiple code words may be decoded in parallel or even the decoder itself might be decomposed into multiple threads that run in parallel [START_REF] Wubben | Benefits and impact of cloud computing on 5g signal processing: Flexible centralization through cloud-ran[END_REF].

In view of the two above use cases (i.e., vBBU, vEPC), we see that virtualization requires an appropriate decomposition of virtual network functions in order to guarantee an efficient resource utilization in terms of computing by adequate scheduling algorithms. This fact enables resource pooling and statistical multiplexing in the execution of sub-functions on multi-core platforms, performing the same tasks with less hardware or capacity. It is set out in more detail in the next section.

III. MODEL DESCRIPTION A. Model setting

At present resource allocation in virtual platforms have near static behavior; virtual machines or containers are reserved for a specific VNF (e.g., vEPC, vBBU) even though a VNF is sporadically invoked. As a consequence, efficiency in resource utilization is not achieved, since computing resources are frozen but no used. It would thus be better to perform statistical multiplexing on computing resources. More precisely, in the following, we assume that a set of cores is available to execute VNFs with dynamic resource orchestration. Cores are allocated for the execution of a VNF when that function is invoked.

As discussed in Section II, in most cases, the runtime of a virtual network sub-function is deterministic. Nevertheless, other aspects discussed in [START_REF] Haugh | Examining factors of the nfv-i impacting performance and portability[END_REF] are intrinsic to the execution of a sub-function such as memory access time, caching policies, disk access time, inter-processor communication, system buses and I/O buses behavior, among others. In fact, the multiprocessing architecture plays an important role in the whole network virtualization performance. In the following, we include these issues in the holding time of cores, which becomes a random variable. This is a classical assumption in network performance modeling when various factors influence the execution of a job.

Taking into account a pool of cores which are statistically shared by several active VNFs, a significant gain in resource savings is expected. The counterpart is that the execution of a VNF might be delayed until some core is available in comparison with the case when computing resources are dedicated to a single VNF. In this paper, we assume that a VNF is composed of sub-functions, each of them being executed on the multi-core platform. The goal of this work is to investigate how the underlying sub-functions should be scheduled or even conceived to improve the VNF performance.

In this context, the functional disaggregation in the virtualization process should take into account the correspondence between sub-functions, because it determines the behavior in the execution process. Then, a particular VNF could be viewed as a process flow with sub-functions either running in sequence or else being executed in simultaneous threads. The present work analyzes the behavior of scheduling algorithms for both configurations, i.e., when VNFs are executed as a chain of sub-functions in contrast with the case when sub-functions are allowed to be executed in parallel.

B. Queuing system formulation

Let us consider a computing facility (i.e., a data center) equipped with a pool of c cores capable of executing elementary sub-functions. These sub-functions are part of a VNF, also referred to macro-function. Such a VNF is composed of k sub-functions; when a VNF is invoked the entire batch of k sub-functions has to be handled by the computing facility. The model proposed below takes k as a constant, however, it could be easily extended to the case when k is a random variable (e.g., with a geometric distribution).

VNF requests occur at the computing facility according to a Poisson process with rate λ. The execution time of a subfunction is random with an exponential distribution, with mean 1/µ j for the jth sub-function of a VNF, where j = 1, . . . , k. We assume that execution times of various sub-functions are independent, hence, the execution time of a macro-function is originally the sum of k exponential random variables. The mean execution time is

E[s] = K j=1 1 µ j .
When all µ j are equal to some µ > 0, the execution time is simply an Erlang random variable with mean k/µ and variance k/µ 2 , denoted, for short, by Er(k, µ).

In general, this model involves a pool of cores which execute only one sub-function at a time. The holding time of a core by the jth sub-function of a VNF is exponential, with mean 1/µ j . Hence, when a request occurs while all cores are busy, the function is queued, see Fig. 2 for an illustration. The total amount of time r to treat a VNF by the computing facility is composed by the queuing delay q and the execution time of sub-functions s, so that, the system response time for the execution of a VNF is given by r = s + q. The queuing delay is the amount of time that the VNF spends in the system while none of its sub-functions is executed. We assume that the waiting room of the computing facility is infinite. It is then clear that the system is stable if and only if the load ρ defined by

VNF i-1 VNF i 1 2 j-1 j k i,1 i,2 i,j i,k … … i,j-1 i-1 i
ρ = λE[s] c
is less than 1. In the following, we assume that this stability condition is always satisfied. That condition is not necessary in the case of reneging, since this phenomenon makes the system stable at the price of rejecting requests.

C. Performance measures

In view of the above queuing system formulation, we characterize its performance by the distribution of the sojourn time r of VNFs in the system. If r takes small values, then macrofunctions (VNFs) are rapidly executed and consequently, the overall efficiency is achieved.

On the other hand, if the sojourn time is too large, then some VNFs are slowly executed and in some cases (e.g., vRAN) completing the execution of a VNF becomes useless. As presented in Section II, a time budget of 3ms for BBU processing is available to allow continuous transmission per UE. This fact is equivalent to introduce the concept of reneging. This budget represents the time between the end of the up-link transmission and the start of the down-link subframe carrying the corresponding ACK/NACK [START_REF] Dötsch | Quantitative analysis of split base station processing and determination of advantageous architectures for lte[END_REF]. More precisely, a customer reneges if its sojourn time in the system exceeds too much its service time. Reneging in multipleservers queuing systems is discussed in [START_REF] Liu | Balking and reneging in m/g/s systems exact analysis and approximations[END_REF]. The adopted reneging criterion is given below:

Definition 1 (Reneging macro-function): A VNF execution request reneges if its sojourn time r > (1 + θ)s for some θ > 0. This is equivalent to the condition q > θs, where q is queuing delay of the macro-function.

In case of reneging, the VNF leaves the system and all resources used prior to its departure are wasted. Hence, the performance of the system is characterized beyond the sojourn time r by the reneging rate η(θ) as a function of the parameter θ.

Whether with the distribution of sojourn time or reneging, the system performance depends on the queuing delay of VNF requests, i.e., the time that a VNF spends in the system without receiving service. This queuing delay depends on how subfunctions are sorted for their execution. This is precisely the role of the scheduling algorithms described in the following section.

IV. SCHEDULING ALGORITHMS

The scheduling algorithm determines which VNF and, more precisely, which sub-function gets access to a particular core of the computing facility. In other words, the scheduling algorithm selects which sub-function is executed. The scheduling strategy has a direct impact on the system response time as well as on the resources utilization efficiency. Thus, the scheduling analysis aims at identifying conditions and constraints to improve the system performance.

The scheduler selects the next sub-function to be processed among VNFs that are ready to be executed, and allocates to it a processing unit (core). Let us consider three scheduling algorithms: Dedicated Core (DC); Round Robin (RR) and Greedy (G) presented in the following subsections.

1) Allocating the entire macro-function to a Dedicated Core: This method processes all sub-functions forming a particular VNF upon the same core. The scheduler selects the VNF keeping the arrival order. In this case, the computing facility can be described by an M/G/c queue system, where the service time is the sum of independent exponential random variables. If all µ j are the same, then an M/Er(k, µ)/c queue system is obtained.

2) Allocating sub-functions by Round Robin criterion: This algorithm handles the VNFs under the RR criterion. Incoming VNFs integrate the circular cycle upon arrival, while the scheduler selects them in first-in-first-out order. Hence, subfunctions are assigned to the pool of cores in a circular order. As detailed in [START_REF] Singh | A priority based round robin cpu scheduling algorithm for real time systems[END_REF], Round Robin generally employs timesharing, assigning to each task a time slot, and interrupting the task if it is not completed. In this paper, the proposed algorithm does not force the sub-function going out of the core if it is not finished, since we assume that sub-functions are not divisible.

3) Greedy allocation of sub-functions: Unlike the precedent algorithms, the arrival order of VNFs is prioritized at the moment of allocating sub-functions, in the sense that, if a VNF starts its service, then the system tries to finish it as soon as possible. Thus, when the first sub-function of a VNF starts its execution, the next core which becomes available is used to serve the second sub-function and so on.

As a consequence, those VNFs which have started their service have priority over those which have not begun their execution. Incoming VNFs are obliged to wait to be processed until precedent VNFs have completely finished their execution. In other words, the system does not try to share the computing resources among the VNFs which are in the system as in the case of RR algorithm. Greedy algorithm aims at completing the service of VNFs at the earliest possible way without fairness.

V. PERFORMANCE ANALYSIS A. Simulation setting

The fundamental issue is to analyze the performance of the scheduling algorithms introduced in the previous section. We are interested in the sojourn time distribution of VNFs as well as in the reneging rate in case of deadline for the execution of a VNF. Keeping in mind that a VNF is formed by subfunctions, which can be executed in parallel or in series, the behavior of scheduling algorithms is analyzed for both cases.

Different scenarios have been considered in relation with the computing pool size and the number of sub-functions per VNF, taking values where c < k, c > k and c k. Furthermore, we have considered different load conditions, namely ρ = 0.6 (moderate load) and ρ = 0.9 (heavy load). Finally, we analyze the behavior when all sub-functions have the same mean execution time. This means that for all j = 1, . . . , k, µ j = µ for some constant µ > 0. The parameter 1/λ equal to the mean inter-arrival time of VNFs is taken as time unit. The parameter µ is adjusted so that the load ρ = kλ/(µc) is equal to a prescribed value.

The performance results analyzed below correspond to those scenarios with heavy load and a number of cores greater than the number of VNF's sub-functions. The same configuration was applied when sub-functions belonging to a particular VNF have different mean execution times. We have also evaluated the performance when different classes of VNFs share the same multi-core platform. This latter scenario offers the possibility to support new services with the same infrastructure, having for instance the co-execution of RAN functionality with other network functions. The results obtained for these scenarios are qualitatively the same as those presented below and are not reported in this paper.

B. Scheduling performance without reneging

In the framework of VNFs where their components represent a forwarding graph of sub-functions, Greedy algorithm is not applicable because the chaining constraint can not be respected. For instance, BBU is a chained process whose functional blocks roughly are modulation/demodulation, encoding/decoding, radio scheduling, and RLC/MAC PDU generation/decomposition. Hence, we analyze the performance in terms of sojourn time of VNFs when these are scheduled following the RR criterion as well as the DC algorithm. Fig. 3 illustrates the behavior of both algorithms considering a highly loaded system with parameter ρ = 0.9, which means that the arrival rate of VNFs is important for the computing capacity. From Fig. 3, it turns out that the simple DC algorithm offers slightly better performance than the RR algorithm, for the tail the of the sojourn time distribution. This will be confirmed by considering the case when VNFs can renege. The additional delay introduced by a higher sojourn time represents in most cases the degradation of the customer experience. For instance, when executing a vBBU, all information sent from the physical layer to the MAC layer and vice-verse has to deal with this extra time. In this case, a critical effect of a high sojourn time is the caducity of channel measurements sent from the UEs to the eNodeB, which are used for radio scheduling and other fundamental issues in the characterization of radio signals. As a consequence, the mobile network losses both energy and spectral efficiency.

Let us consider now, VNFs whose sub-functions are allowed to be executed in parallel. This means that the execution of a particular sub-function is independent of the results of the previous one. In this context, Greedy (G) algorithm can be applied. Fig. 4 presents the behavior of three algorithms where the performance of both G and RR is notably better than the behavior obtained with the DC algorithm. Even more, G algorithm shows a slightly better performance than that of RR.

It is evident that the chaining constraint considered in the first simulation balks the advantages of RR algorithm in its attempt of fairness. In this way the simplicity embedded in DC criterion is the most appropriate for a computing pool environment, all the more as the complexity of RR does not improve the sojourn time of VNFs. In this subsection, we analyze the scheduling performance considering reneging with parameter θ = 1. As explained in Section III, this factor represents the deadline present in some network functions which require real-time execution as in the case of base band processing of mobile networks. Considering the scenario where c > k and applying the chaining constraint, Fig. 5 illustrates the behavior of RR and DC algorithms. Again their performances in terms sojourn time of a VNF are similar with an advantage for DC when considering the tail of sojourn time distributions. Table I shows reneging rate which represents the number of VNFs that have not finished their execution.

The utilization rate and waste rate show respectively the occupation of the computing platform by VNFs which have been completely processed and the occupation of cores by subfunctions belonging to VNFs which have reneged. DC offers a slightly better performance, in line with the observation made for the tail of the sojourn time distributions when there is no reneging. RR algorithm yields a higher reneging rate, and consequently worse utilization factor than DC algorithm. Hence, more computing resources are wasted.

Analyzing the case when sub-functions can be executed in parallel, the behavior of scheduling algorithms turns out again more favorable for G and RR criteria. It is shown in Fig. 6 where reneging rate was established with θ = 1, it means that a VNF interrupts its service and leaves the system if its sojourn time is greater than the double of the required execution time. The same kind of behavior has been observed for smaller and greater values of θ. As much as in terms of execution delay as in terms of reneging rate presented in Table II, G has the best performance. It is evident that the behavior of DC is relegated by RR and G performance. Nevertheless, the utilization rate of DC remains interesting when compared with that resulting of RR execution. Results show that the worst algorithm in terms of resources saving is RR although it keeps a better sojourn time. Greedy becomes the most efficient and notably the most suitable when the execution of sub-functions is not limited by the chaining constraint.

VI. CONCLUSION

In this work, we have studied a system executing VNFs on a computing platform composed by a pool of cores; each VNF is composed of several sub-functions. We have analyzed by simulation the performance of three algorithms for scheduling the execution of sub-functions of active VNFs (namely, Round Robin, Dedicated Core and Greedy). It turns out that when sub-functions can be executed in parallel, the Greedy algorithm ensures the best performance in terms of execution delay. We have also considered the case when VNFs may renege because the sojourn time in the system exceeds some threshold related to the required amount of service. Still in this case, the Greedy algorithm offers the best performance.

In the case of chained sub-functions Greedy algorithm can not be applied, and the performances observed with Dedicated Core and Round Robin are similar, so the complexity added by this latter is not justified.

This phenomenon has to be taken into account when designing VNFs executed on a pool of cores. In particular when decomposing a network service into components or microservices, the best choice is to decompose a function into subfunctions which can be executed in parallel and independently of each other.

Fig. 1 .

 1 Fig. 1. Network Function Virtualization.

Fig. 2 .

 2 Fig. 2. Elements of the queuing system.

Fig. 3 .

 3 Fig. 3. Scheduling performance considering chained sub-functions

Fig. 4 .

 4 Fig. 4. Scheduling performance considering no chained sub-functions

Fig. 5 .

 5 Fig. 5.Scheduling performance considering chained sub-functions and reneging

1 PFig. 6 .

 16 Fig. 6. Scheduling performance considering no chained sub-functions and reneging