Antoine Caillebotte
email: caillebotte.antoine@inrae.fr

Estelle Kuhn
email: estelle.kuhn@inrae.fr

Sarah Lemler
email: sarah.lemler@centralesupelec.fr

Estimation and variable selection in a joint model of survival times and longitudinal outcomes with random effects

Keywords: Joint model, non-linear mixed effects model, Cox model, high dimension, preconditioned stochastic gradient, proximal operator

This paper considers a joint survival and mixed-effects model to explain the survival time from longitudinal data and high-dimensional covariates. The longitudinal data is modeled using a nonlinear effects model, where the regression function serves as a link function incorporated into a Cox model as a covariate. In that way, the longitudinal data is related to the survival time at a given time. Additionally, the Cox model takes into account the inclusion of high-dimensional covariates. The main objectives of this research are two-fold: first, to identify the relevant covariates that contribute to explaining survival time, and second, to estimate all unknown parameters of the joint model. For that purpose, we consider the maximization of a Lasso penalized likelihood. To tackle the optimization problem, we implement a pre-conditioned stochastic gradient to handle the latent variables of the nonlinear mixed-effects model associated with a proximal operator to manage the non-differentiability of the penalty. We provide relevant simulations that showcase the performance of the proposed variable selection and parameters' estimation method in the joint modeling of a Cox and logistic model.

Introduction

A very current issue in many fields is better understanding the interactions between dependent dynamic phenomena. For example, in medicine, this may involve the dynamics of a patient's tumors in oncology and the effects of anti-cancer treatments administered to the patient. Another example in plant science is the dynamics of plant development in a plot and the spread of an epidemic disease or pests in that plot. The phenomena considered are often complex, both in terms of their modes of interaction and their temporal and spatial dynamics. Moreover, these phenomena are often observed in populations of heterogeneous or structured individuals, such as patients or plants.

Mathematical modeling has proven to be a powerful tool for understanding the interactions between multiple dynamic phenomena. It also allows for considering variabilities present in the observed population of individuals. Joint modeling of several phenomena has demonstrated its effectiveness in several fields, including medicine, pharmacology, and biology ( [START_REF] Kerioui | Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models[END_REF]). A particular case of joint models concerns the simultaneous modeling of longitudinal data and survival data observed on the same individual. In this type of joint model, longitudinal data are often modeled by a mixed-effects model ( [START_REF] José | Mixed-Effects Models in S and S-PLUS[END_REF][START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF]), and survival data by a survival model such as the Cox model ( [START_REF] Cox | Regression models and life-tables[END_REF]). The latter allows for modeling the instantaneous risk of the survival variable as a function of covariates. It is also possible to include longitudinal data modeling as a covariate in the Cox model via a linking function. The objective is then to estimate the model parameters from the observations and to select relevant covariates. Several authors have proposed such approaches ( [START_REF] Wulfsohn | A joint model for survival and longitudinal data measured with error[END_REF], [START_REF] Rizopoulos | Joint models for longitudinal and time-to-event data: With applications in R[END_REF], [START_REF] Król | Joint model for left-censored longitudinal data, recurrent events and terminal event: Predictive abilities of tumor burden for cancer evolution with application to the ffcd 2000-05 trial[END_REF]). Due to the presence of latent variables in the mixed-effects model, inference by maximum likelihood can be made via Expectation Maximization (EM) like algorithms ( [START_REF] Wulfsohn | A joint model for survival and longitudinal data measured with error[END_REF], [START_REF] Guedj | Joint modeling of the clinical progression and of the biomarkers' dynamics using a mechanistic model[END_REF], [START_REF] Rizopoulos | Jm: An r package for the joint modelling of longitudinal and time-to-event data[END_REF], [START_REF] Desmée | Using the saem algorithm for mechanistic joint models characterizing the relationship between nonlinear psa kinetics and survival in prostate cancer patients[END_REF]). The EM-type algorithms, such as the classical Stochastic Approximation Expectation Maximization (SAEM), are the most classical approaches for inferring parameters in the presence of latent variables. They have been developed for estimation in general latent variable models. They are particularly easy to implement in the context of a curved exponential family based on sufficient statistics of the model. Moreover, theoretical convergence results have been established in this context. However, when the model does not belong to the exponential family, which is the case in our context, the methodology is not generic in practice, and the theoretical results fail.

Some exponentialization trick has then been proposed to face this restrictive assumption of the curved exponential family. It consists in considering some unknown parameters as random population variables. However, [START_REF] Debavelaere | On the curved exponential family in the stochastic approximation expectation maximization algorithm[END_REF] have shown that, in general, the parameter returned by the SAEM on the modified model is not a maximum likelihood of the initial model, and they have suggested the use of this exponentialization trick with variances of the new random population variables that decrease as the iterations of the algorithm progress. This approach also has limitations in practice due to complex algorithmic settings and tuning. The gradient-based methods are another type of approach, often omitted for estimating parameters in latent models. Recently, [START_REF] Baey | Efficient precondtionned stochastic gradient descent fro estimation in latent variables models[END_REF] suggested using a preconditioned stochastic gradient algorithm to deal with parameter estimation in the presence of latent variables. This approach is particularly interesting when considering a model that does not belong to the exponential family, as is the case for the joint model. [START_REF] Baey | Efficient precondtionned stochastic gradient descent fro estimation in latent variables models[END_REF] showed that this algorithm performs well for the nonlinear logistic growth mixed-effects model, which can be used to represent some longitudinal data. Note that Bayesian numerical methods have also been proposed in parallel ( [START_REF] Brown | A bayesian semiparametric joint hierarchical model for longitudinal and survival data[END_REF], [START_REF] Rizopoulos | The r package jmbayes for fitting joint models for longitudinal and time-to-event data using mcmc[END_REF], [START_REF] Kerioui | Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models[END_REF]).

Besides, in many applications, current technological means allow for collecting highdimensional explanatory covariates. These may include, for example, genetic markers or omics data. In addition to the wealth of information provided by these covariates, they also generate difficulties in the statistical analysis of models as it is necessary to adapt statistical and numerical approaches to their high dimensionality. One possible approach is to consider a penalized estimator, such as the Lasso ( [START_REF] He | Simultaneous variable selection for joint models of longitudinal and survival outcomes: Variable selection in joint models[END_REF], [START_REF] Yi | Simultaneous variable selection and estimation for joint models of longitudinal and failure time data with interval censoring[END_REF]), and adapted numerical methods, such as stochastic proximal gradient ( [START_REF] Achab | Learning from sequences with point processes[END_REF], [START_REF] Fort | Stochastic proximal gradient algorithms for penalized mixed models[END_REF]).

In this paper, we consider a joint model which combines, through a link function, a nonlinear mixed effect model for longitudinal data and a Cox model for the survival times, including covariates of high dimension. Our work aims to select the relevant variables among the high-dimensional covariates in the Cox model part of the joint model based on the whole dataset and then to estimate the model's unknown parameters. For that purpose, we propose an estimate for model parameters, which include a Lasso penalization for the regression parameter of the Cox model. To calculate this estimate in practice, we develop an algorithm combining a preconditioned stochastic gradient to deal with the latent variables in the joint model out of the exponential family and a proximal gradient to handle the non-differentiability of the Lasso penalty used for variable selection in the Cox model part. The proposed algorithm is easy to implement in general joint models without assuming that model density belongs to the curved exponential family.

The paper is organized as follows. In Section 2, we detail the joint model constructed from a nonlinear mixed effects model for longitudinal data and a Cox model for survival data, with high-dimensional covariates and a link function. In Section 3, we present the proposed inference method based on a Lasso penalized estimator and numerical procedure based on a stochastic proximal gradient algorithm. Finally, we illustrate the methodology in Section 4 through a simulation study. The paper ends with a conclusion.

Joint survival and mixed-effects model

We consider N individuals and study, for each individual i, the survival time T i , corresponding to the duration until the occurrence of an event of interest, and longitudinal data, more precisely repeated observations J times denoted by Y i,j with i P t1, . . . , N u and j P t1, . . . , Ju. Note that our work can easily be generalized to the case where there are different number of longitudinal observations for each individual of the population. The following describes the joint model we considered.

Survival model

The survival time T i of individual i is the time between a fixed initial moment and the occurrence of an event of interest. It is a positive random variable. To characterize the distribution of T i , we use the hazard function defined by:

h i ptq :" lim dtÑ0 Ppt ď T i ă t `dt|T i ě tq dt ; @t ě 0. (1)
The Cox model ( [START_REF] Cox | Regression models and life-tables[END_REF]) is one of the most classical models in survival analysis. It allows us to relate the hazard function of the survival time T i to covariates U i P R p with p being the number of covariates. In our approach, we will consider the high-dimensional setting with a large number of covariates, such as p is very large with respect to N . The Cox model for individual i is written as follows:

hpt|U i q " h 0 ptq exppβ T U i q, (2)
with β P R p a regression parameter and h 0 the baseline hazard function that characterizes a common behavior in the observed population. In the sequel, we will consider a parametric baseline function denoted by h θ base where θ base P R b are its parameters. Therefore, the Cox model's unknown parameters are β and θ base .

In addition to the covariates, we consider explaining some of the survival time variability using the longitudinal data dynamic, which will be modeled using a nonlinear mixed effects model. Let us present the mixed-effects model before explaining the integration of this new component into the Cox model.

Nonlinear mixed-effects model

The longitudinal data are observed J times for each individual i P t1, . . . , N u. Let us denoted by Y i,j the j-th observation of the i-th individual for j P t1, . . . , Ju and i P t1, . . . , N u. We model this longitudinal observation using a nonlinear function m that depends on individual parameters represented by the latent variable Z i as follows:

# Y i,j " mpt j ; Z i q `εi,j ; ε i,j " i.i.d. N p0, σ 2 q, Z i " i.i.d. N pµ, Γq. @1 ď i ď N, 1 ď j ď J (3)
where, t j is the j-th observation time, and ε i,j is an additive noise assumed centered Gaussian with unknown variance σ 2 . The latent variable Z i describes the inter-individual variability of the population. It is assumed that Z i follows a Gaussian distribution with unknown expectation µ and variance Γ. The unknown parameters of the nonlinear mixedeffects model are therefore µ, Γ, and σ 2 . Let us introduce in the following the link function, which will combine the two previous models by modeling the influence of the dynamic of the longitudinal observation of the hazard function.

Joint survival and mixed-effects model

We assume that the hazard of the survival time is related to the longitudinal data dynamic through the link function m as follows: hpt|Mpt, Z i q, U i q " h θ base ptq exppβ T U i `αmpt, Z i qq, @t ě 0,

where Mpt; Ziq " tmps; Ziq|@s, 0 ď s ă tu describes the past values of the longitudinal dynamic up to time t. The parameter α represents the influence of the longitudinal dynamic on the survival data. The joint model can be written as follows:

$ ' & ' % hpt|Mpt, Z i q, U i q " h θ base ptq exppβ T U i `αmpt, Z i qq Y i,j " mpt j ; Z i q `εi,j Z i " i.i.d. N pµ, Γq ; ε i,j " i.i.d. N p0, σ 2 q. @1 ď i ď N, 1 ď j ď J (5)
The unknown parameters for the joint model include the parameters of the Cox model and those of the nonlinear mixed effects model, as well as the link function parameter of the joint model. We note θ " pθ base , β, µ, Γ, σ 2 , αq P Θ the vector of unknown parameters with Θ Ă R d being the parameter space. In the following section, we propose an estimation method for these parameters.

Definition of the marginal likelihood

We consider the maximum likelihood estimator to infer the joint model parameters. In the context of latent variable models, the marginal likelihood, denoted by L marg , is obtained by integrating the complete likelihood over the latent variables, which are not observed. Let D " pY, Tq be the observed variables:

L marg pθ; Dq " ż f θ pD, ZqdZ " ż g θ pD|Zqp θ pZqdZ (6)
where f θ , g θ , p θ are respectively the density of the pair pD, Zq, the density of Z conditionally to D, the density of Z. Due to the integral, it is difficult to directly compute the maximum of the marginal likelihood, which does not have an analytical form in this latent variable model. Therefore, we will use numerical methods to solve this maximization problem.

Definition of the penalized estimator for variable selection

We introduce a penalty and consider a penalized maximum likelihood estimator to deal with the high dimension of the covariates. We aim to select relevant variables among the covariates of the survival model. We use the Lasso (Least Absolute Shrinkage and Selection Operator) procedure which was initially proposed for linear regression models ( [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]) and the Cox model ( [START_REF] Tibshirani | The lasso method for variable selection in the cox model[END_REF]). This method enables us to handle high-dimensional data and select a subset of explanatory covariates from a large collection. We consider a Lasso penalty which only depends on the parameter β:

penpθq " }β} 1 " p ÿ k"1 |β k | ,
Our goal is then to maximize the logarithm of the marginal likelihood where the penalty is integrated as follows. Let us define the penalized maximum likelihood estimator by: θ " arg max θPΘ plog L marg pθ; Dq ´λ penpθqq ,

where Θ denotes the parameter space and where λ is a positive parameter called the regularization parameter. The larger the value of λ, the more β will be constrained to have zero components. Conversely, the smaller the value of λ, the more free the components of β will be. It is customary to determine the value of λ using cross-validation ( [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]).

Usually, when we deal with latent variables, since the marginal likelihood is nonanalytic, classical methods used to infer the unknown parameters are Expectation Maximization like algorithme ( [START_REF] Shu | Handbook of computational statistics: concepts and methods[END_REF]). The inconvenience of these procedures is that it is welladapted to models belonging to the curved exponential family, which is not the case for the joint model we consider. Recently [START_REF] Baey | Efficient precondtionned stochastic gradient descent fro estimation in latent variables models[END_REF] have proposed a preconditioned stochastic gradient descent for estimation in a latent variable model adapted to general latent variables models. Moreover due to the non-differentiability of the considered penalty, we will use a proximal algorithm as proposed by [START_REF] Achab | Learning from sequences with point processes[END_REF] and [START_REF] Fort | Stochastic proximal gradient algorithms for penalized mixed models[END_REF]. Thus, we add a proximal gradient in the procedure proposed in [START_REF] Baey | Efficient precondtionned stochastic gradient descent fro estimation in latent variables models[END_REF] and implement a preconditioned stochastic proximal gradient algorithm to calculate the estimator.

Implementation of the inference procedure

We deal simultaneously both with unobserved random effects of the mixed-effects model and the penalty term by implementing a preconditioned stochastic proximal gradient, called SPG-FIM in the sequel. A forward-backward splitting algorithm such that a Stochastic Proximal Gradient can compute the estimate 7. The latter has a stochastic approach to replace missing data with simulations. The algorithm is divided into three steps; a realization of the latent variables is sampled with a first step called Simulation, which uses a Metropolis-Hastings sampler ( [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF]). The second step is the classical gradient descent on the approximate complete likelihood, the Forward step. Following the procedure proposed in [START_REF] Baey | Efficient precondtionned stochastic gradient descent fro estimation in latent variables models[END_REF], we have chosen to use a preconditioning of the gradient with an estimate of the Fisher information matrix. The latter is updated during the iterations using the estimate proposed by [START_REF] Delattre | Estimating fisher information matrix in latent variable models based on the score function[END_REF]. The last step, called Backward, deals with the penalty term. We apply the classical proximal operator ( [START_REF] Jacques | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF][START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF], defined below Prox pen pβq " arg min

β 1 PR p ˆpenpβ 1 q `1 2 }β ´β1 } 2 2 ˙. ( 8
)
With the Lasso penalty, the proximal operator has an explicit form:

pProx Lasso pβqq i " $ & % 0 si |β i | ă λ β i ´λ si β i ě λ β i `λ si β i ď ´λ ; @i P t1, ..., pu. (9)
The Backward step corresponds to the application of the proximal operator on the result of the Forward step.

As the penalty only depends on β, the proximal operator selects the β components that seem to be the most explanatory of the data. It computes a sparse solution for β but also applies shrinkage on the non-zero components so that the Lasso estimator is biased. Therefore, we detail a method to obtain an unbiased estimator in what follows.

Algorithm 1 provides the steps of the stochastic proximal gradient, where pγ k q kě1 is a step size such that @k P N, γ k P r0, 1s,

Compute v k " 1 N N ÿ i"1 ∇ log f θ k pD i Z pkq i q 7 ' FIM computation : 8 ' Compute the stochastic approximation 9 @i P t1, ..., N u, ∆ pkq i " p1 ´γk q∆ pk´1q i `γk ∇ log f θ k pD i Z pkq i q 10 ' Compute the FIM : 11 F IM k " 1 N N ÿ i"1 ∆ pkq i `∆pkq i ˘T 12 ' Gradient descent : 13 ' Forward step : ω k`1 " θ k ´γk F IM ´1 k v k 14 ' Backward step : 15 θ k`1 " Prox γ k pen pω k`1 q " arg min θ 1 PΘ " γ k penpθ 1 q `1 2 }ω k`1 ´θ1 } 2 2
* 16 end 17 return θ " θ K

Simulation study

We generated data according to the joint model presented previously in Equation 5. We consider N " 100 individuals, each individual being observed J " 20 times. We use the classical logistic function in the nonlinear mixed effect model detailed in Equation 3:

m : t Þ Ñ Z 1 1 `exp ˆZ2 ´t Z 3 ˙, (10)
where Z 1 represents the asymptotical maximum value of the curve, Z 2 represents the value of the sigmoid's midpoint, and Z 3 represents the logistic growth rate. We model for each individual i the corresponding individual parameter Z i P R 3 through a Gaussian random variable with expectation µ P R 3 and a diagonal variance Γ " diagpγ 2 1 , γ 2 2 , 0q meaning that the third parameter µ 3 is modeled as a fixed effect. We consider a fixed Weibull baseline defined as h a,b ptq " ba ´bt b´1 , where a " 80 and b " 35 are fixed (i.e. not estimated) in the simulation study.

We fix p " 100 and choose the vector β such that the first four components are equal to p´2, ´1, 1, 2q and the rest are equal to zero. Additionally, we generate the matrix of covariates U with N rows and p columns, following a uniform distribution U i,l " Upr´1, 1sq; @i P 1, ..., N , l P 1, ..., p. We choose finally the link function parameter α " 11.11.

All the parameter values are detailed in Table 4. We focus in this simulation study on the selection of variables and on the inference of the parameters of the mixed-effect model as well as α the multiplicative parameter of the Cox model. The proximal operator ( 9) has a shrinking effect on the estimator after its application, meaning that the values found for β are smaller than expected. Therefore as it is usually the case the estimator of β is biased. We thus divide the inference into two steps, an exploratory one that allows us to select the support of the vector β through a Lasso penalization estimation, and a second step of inference without penalization, where we have restricted the number of covariates with respect to the selected support. We conduct the following inference methodology: ' Run the SPG-FIM algorithm in order to compute θLasso pλq " arg max θPΘ plog L marg pθ; Dq ´λ pen Lasso pθqu , for different values of λ on a fixed grid.

Parameter µ 1 µ 2 µ 3 γ 2 1 γ 2
' Select the best regularization parameter such that λ m " arg min λ BICpλq according to the BIC criterion (see [START_REF] Schwarz | Estimating the dimension of a model[END_REF]):

BICpλq " ´2 logpL marg p θLasso pλq; Dqq `k logpN q

where k is the number of non-zeros components in β. Note that the quantity L marg pθ; Dq " ż Z f θ pD, ZqdZ is computed by approximating the integral using a Monte-Carlo procedure.

' Choose the reduced support of β according to the estimate θLasso pλ m q obtained previously from a run of the SPG-FIM (Algorithm 1).

' Compute θMLE the maximum likelihood estimate in the reduced model without the penalty term and therefore without bias with the SG-FIM (SPG-FIM without the Backward step).

Following those steps, we consider 50 runs. Each time, the initial parameters θ 0 were randomly drawn. We initialize each component of β 0 according to a uniform distribution Upr´1, 1sq. We choose the following step sequence:

γ k " " 1 si k ď 1100 pk ´1100q ´1 si k ą 1100 , (11)
The step sequence starts to decrease after a warm-up period of 1100 iterations.

The components of β selected by the Lasso in the full model are presented in Figure 1. We can see that out of the 100 variables considered, only 8 were selected by the Lasso procedure. Note that the four variables that are additionally selected are only output once out of the 50 runs, with a very small estimated value of 10 ´4. Figures 4 and5 display the estimated parameter and gradient values as a function of iterations during the first execution of the SPG-FIM algorithm. Figure 2 shows their estimation in the reduced model without penalization. We can observe from Figure 2 that the estimates are much better in particuliar have less bias, once the support has been properly selected. Table 4 summarizes the results of the 50 algorithm runs. The mean of the estimates over the 50 runs is near to the true value, and the relative RMSE is small. We can see that the estimate of the regression parameter of the link function α is more difficult to be estimated than the other parameters. Figure 6 illustrates the evolution of the 100 components of the vector β throughout the iterations of the SPG-FIM algorithm. For each iteration, each component is represented by a colored square, where the square's color corresponds to the component's magnitude. Shades of blue represent small values, and shades of red larger values. Overall, this graphical representation provides a concise and visual depiction of how the individual components of the vector evolve over time in the SPG-FIM algorithm. We can see that a large set of components converge to 0, while the first four components converge to the expected values. So, at the final iteration, we obtain the support of β.

Conclusion and perspectives

In this work, we jointly addressed variable selection and parameter estimation in a joint model for survival data and longitudinal outcomes. We connected a survival model with a nonlinear mixed-effects model through a linking function and covariates eventually in high dimension. To estimate the unknown parameters of this joint model, we conduct a preconditioned proximal stochastic gradient to deal with the latent variables and the Lasso penalty. Our methodology has been thoroughly evaluated on simulated data to demonstrate its performance.

One interesting perspective of this work consists in addressing the prediction task carefully. It would be interesting to set up a method for predicting survival time from some start of longitudinal data observation. Previous studies, such as [START_REF] Desmée | Nonlinear joint models for individual dynamic prediction of risk of death using hamiltonian monte carlo: application to metastatic prostate cancer[END_REF] and [START_REF] Kerioui | Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models[END_REF]), have explored this task within the context of joint models.

Then, we assumed no censorship in the observed data to simplify the inference procedure. However, it could be essential to consider the presence of censorship in real-world scenarios. Therefore another important perspective is to adapt our procedure to survival times that can be censored.

Furthermore, considering the partial likelihood for the survival part rather than the complete likelihood would be an intriguing avenue to explore. This approach would eliminate the need to make assumptions about the form of the baseline risk, which is especially relevant when dealing with the Cox model. Indeed, when considering the Cox model, the partial likelihood allows for estimating the regression parameter without knowing the baseline function. It would be interesting to see how our algorithm performs with this partial likelihood instead of the complete one.

Algorithm 1 : 3 ' 4 Draw 5 '

 1345 Stochastic proximal gradient with FIM preconditionning (SPG-FIM)Require : Number of iterations K ě 1 ; sequence of step-size pγ k q kě1 1 Initialize Starting point θ 0 P R d , ∆ 0 2 for k " 1 to K do Simulation step : Z pkq using a single step of a Hastings Metropolis procedure Gradient computation : 6

2 σ 2 αTable 1 :

 221 True value 0.3 90 7.5 2.5 10 ´3 20 10 ´3 11.11 Parameter β 1 β 2 β 3 β 4 β 5 . . . β p True value -2 -1 1 2 0 . . . 0 True parameter values used for the simulation

Figure 1 :

 1 Figure 1: Box plots of the selected variables during the 50 runs of the SPG-FIM in the full model. The colored stars indicate the true value used for the data simulation

Figure 2 :

 2 Figure 2: Box plots of the estimates for the selected variables during the 50 runs of the SG-FIM in the reduced model. The colored stars indicate the true value used for the data simulation

Figure 3 :

 3 Figure 3: Distribution of estimates for each of the parameters during 50 runs of the SG-FIM. Vertical black lines represent the true values of each parameter

Figure 4 :

 4 Figure 4: Parameter estimates during SPG-FIM iterations Figure 5: Gradient values during SPG-FIM iterations

Figure 6 :

 6 Figure 6: Tile representation of the selected variables during iterations of SPG-FIM

Table 2 :

 2 Mean and relative RMSE for 50 runs of the algorithm.

		value mean estimate relative RMSE
	α	11.1	11.9	0.209
	γ 2 1 γ 2 2	0.0025 20	0.00254 20.1	0.182 0.160
	µ 1	0.3	0.300	0.0196
	µ 2	90	90.1	0.00577
	µ 3	7.5	7.52	0.0232
	σ 2	0.001	0.001	0.0366

Inference methodNote that there is often censoring in survival analysis, which leads to partially observed data: survival times are not directly observed. Available information is censored times and indicators, making the estimation task more complex. For the sake of simplicity, since we focus on the high-dimensional covariates selection task, we will not consider censoring in our approach for the moment. However, it will be part of further work.

Funding and Acknowledgements

This work was funded by the Stat4Plant project.


Lucien Arnaud

Ndjie Ngale

Eddy Caron
email: eddy.caron@ens-lyon.fr

Yulin Zhang

Mélanie Fontaine
email: melanie.fontaine@crispi-upjv.fr

Simulation-based Estimation of Resource needs in Fog Robotics Infrastructures

Keywords: Robotics, Fog Computing, IoT, Simulation, Machine Learning

Embedded devices are increasingly connected to the Internet to provide new and innovative applications in many areas. These devices (Edge devices or "Things" in IoT) are heterogeneous sensors, cameras and even robots performing sometimes certain tasks locally. Fog Computing (or Fog robotics) optimizes the management of these tasks, offering data management mechanisms (computation and storage) closer to the data source. Nevertheless, many aspects remain closed to Fog computing environments like resource needs estimation in such environments. Indeed, such a topic remains a critical challenge, as it falls under either solving very complex optimization problems or comparing hypothetical scenarios very time consuming and/or expensive for deployment in a real environment. To help on this challenge we built SERFRI, an approach to estimate the resource needs in Fog robotics environments based on simulation. This approach optimizes simultaneously the duration and the Fog resources utilization cost in order to determine the minimum resource requirements compromising both metrics. We validated this approach on an existing robotics use case. This one aims at deploying a human face detection service on streaming images.

INTRODUCTION

The Internet of Things (IoT) has transformed the Internet, enabling communication between every kind of objects (Things). The growing number of sensors and smart devices increased the possibilities of data generation and collection (R. [START_REF] De Oliveira | A real-time and energy-aware framework for data stream processing in the internet of things[END_REF]. These smart devices are robots in our case, making the previous observation more important in the field of robotics [START_REF] Turnbull | Cloud robotics: Formation control of a multi robot system utilizing cloud infrastructure[END_REF][START_REF] Norman | On data-limited and resource-limited processes[END_REF]. Indeed robotics is becoming more and more topical, as it facilitates the automation of tasks and guarantees better accuracy and speed of execution compared to human processing [START_REF] Kwon | When humans aren't optimal: Robots that collaborate with risk-aware humans[END_REF][START_REF] Moniz | Robots working with humans or humans working with robots? searching for social dimensions in new human-robot interaction in industry[END_REF].

Basically, the high level architecture of a robot (Figure 1) consists of: sensors which get information about the robot environment, actuators which allow the robot to move in its environment, micro-controller which manages the actuators' motion, a microprocessor which computes tasks and the battery alimenting the whole system. In the basic running process of a robot, there are two main problems related to the computational capabilities of the robot: (i) the first one is about energy consumption, since this process runs frequently and asynchronously the battery will be used by the assembly of robot components simultaneously, which makes the activity not feasible in the long term; (ii) the second one is about the fact that in most cases, robots do not have enough computing power or storage capacity to process (for a current service) and post-process (e.g: for the updating of a service) data generated.

In some cases it has been shown that it is more advantageous relative to time constraints to deploy robotic applications (which must be real-time) on remote resources. Figure 2 presents for example a comparison between the ETL (Extract -Transform -Load) process duration of deploying a robotic application locally and remotely (36 times more powerful than the robot). This application was made using a Reachy humanoid robot [START_REF] Mick | Reachy, a 3d-printed humanlike robotic arm as a testbed for human-robot control strategies[END_REF] that captured 100 images (at the rate of 7 images per second) that were processed locally or on remote resources, to determine whether there was a human face or not. Then by the fiftieth request, the local deployment becomes in average 50s slower than the remote one. Based on the fact that industrial robotics applications often consist of millions of requests and based on the benefits of the Fog Computing (in terms of latency and power computing), it is then really more advantageous to deploy such applications on Fog-based environments.

Fog Computing offers a large computing and storage capacity near data sources (robots in this case, where we talk about Fog Robotics Platform (FRPs) [START_REF] Gudi | Fog robotics: An introduction[END_REF]). This is a recent approach that solves the latency problem present in native Cloud environments, for long-term and real-time applications (as in the case of industrial robotics). The architecture of Fog computing consists of three main layers: an IoT layer which contains connected data sources (such as smartphones, sensors, robots, etc.) sending requests (periodic, asynchronous and heterogeneous in the case of robotics applications), a Fog layer which contains resources (static/dynamic and heterogeneous/homogeneous) near the data sources and a Cloud layer which contains static resources far the data sources.

We are interested in two main problems related to such platforms: (i) indeed, in general, setting up certain hypothetical scenarios of Fog platforms is not possible at the real scale (this is sometimes due to high price of industrial robots and the financial cost of reserving or using a certain number of resources); (ii) moreover, the global makespan 1 of the platform often decreases with the number of resources increasing while the financial cost of using these resources does not have a known and specific evolution a pri-1 global queries' completion time ori. Two questions then arise, respectively to the two problems mentioned above: (a) How to simulate the deployment of FRPs in a realistic way? (b) how to estimate, based on a realistic simulation, the sufficient (minimum) number of Fog resources needed to guarantee a fast execution and at a lower financial cost of a robotic application?

This paper presents an approach for Simulationbased Estimation of Resource needs in Fog Robotics Infrastructures (SERFRI). It aims at providing a FRP configuration reducing the activity duration as well as the resources usage cost. More precisely, it determines the minimum number of Fog resources reducing significantly simultaneously: the makespan and the financial utilization cost of Fog Computing nodes. As a proof of concept, this approach has been used on an existing robotics use case. For this use case the IT architecture specification has been estimated to support its activity.

In the following sections of the paper, we present: the case study in Section 2, the related works in Section 3, the core principle of our contribution in Section 4, the experimental settings, results and comments in Section 5 and the conclusion in Section 6.

THE CASE STUDY

This section presents the study and the specifications of an existing robotics use case which is about deploying a human face detection service on a stream of images by a Reachy humanoid robot [START_REF] Mick | Reachy, a 3d-printed humanlike robotic arm as a testbed for human-robot control strategies[END_REF]. Reachy (Figure 3) is a human-scale robotic arm with seven joints from shoulder to wrist, which can help researchers explore, develop and test innovative control strategies (such as tele-operation or gaze control) and interfaces on a human-like robot. Its 3D printed structure and ready-to-use actuators make it inexpensive compared to the price of an industrial-grade robot. Using an open-source architecture, its design makes it widely pluggable and customizable, so it can be integrated into many applications. In Figure 3 we can see that Reachy has two high definition cameras and several articulations. This allows realizing through its use several use cases. The present use case employs one of Reachy's cameras to capture a stream of images to determine if there is at least one human face in each image. In this case, Reachy will raise the left hand and if not, the right hand. In Section 1 we have shown that for a number of resources that are globally 36 times more powerful than Reachy's NUC, the ETL process (photo capture, image computation, response to Reachy as shown in Figure 5) was faster. Two main challenges appear: (i) Since the latency be- tween Reachy and the remote resources represents a major parameter (more than 200 ms) in the duration of the ETL process, what type of distributed system would be best suited to deploy such a use case? (ii) To perform the experiment highlighted in Figure 2, 3 Grid'5000 nodes of the chetemi cluster 2 (Intel Xeon E5-2630 v4) located in Lille were used (Figure 4), however with 2 nodes of the same type after a certain number of requests or with 1 node of another type of resources, the same type of behavior could be observed; what would be the sufficient (taking less time and less financial cost to run) number of remote computing nodes for this use case? Part of Section 3 gives some clarification on the type of distributed system 2 https://www.grid5000.fr/w/Lille:Hardware#chetemi adequate for the application underlying this use case and Section 4 will provide an element of answer regarding the estimation of the resource requirements on the adequate platform.

RELATED WORK

Modeling and simulation technology has become a useful and powerful tool in Large Scale Distributed Systems (LSDS) deployments (particularly in the Cloud computing research community) to deal with performance evaluation and security problems in such environments [START_REF] Zhao | Modeling and simulation of cloud computing: A review[END_REF]. It is necessary to model both the platform and the application concerned. In this Section, we present: (a) A brief background on the deployment of robotic ecosystems on LSDS in terms of brief comparative study to justify the use of Fog Computing Platforms, (b) some simulators of application deployments on Fog-based environments and (c) some Fog resource needs estimation approaches.

Deployment of Robotic Ecosystems on Large Scale Distributed Systems

Many works present the deployment of robotic ecosystems on large-scale distributed systems in various paradigms as Service Oriented Architecture (SOA), Cloud Computing, Edge Computing, Fog Computing and peer-to-peer. Cloud and SOA [START_REF] Papazoglou | Service oriented architectures: approaches, technologies and research issues[END_REF]) offer a high level of security, large computing capacities and interoperability. An initial type of deployment to virtualize robotic hardware and software resources and offer them as services on the Cloud via the web is more and more topical. For example, RobotWeb [START_REF] Koubaa | A service-oriented architecture for virtualizing robots in robot-as-a-service clouds[END_REF]) is a SOAP-based web service middleware-based system that binds robot computing resources as services and publishes them to end users. In [START_REF] Koubaa | A service-oriented architecture for virtualizing robots in robot-as-a-service clouds[END_REF], the robots run on the Robotic Operating System (ROS), as it provides a hardware abstraction that facilitates application development.

Edge Computing enables the decentralization of the management of robot functionalities, thus increasing the reliability and reducing the latency of the centralized system [START_REF] Chen | Internet of intelligent things and robot as a service[END_REF]. The work presented in [START_REF] Ray | Internet of robotic things: Concept, technologies, and challenges[END_REF] proposes a new concept that addresses the problems of supporting control and monitoring activities at deployment sites and industrial automation. In that work, smart objects can monitor peripheral events, induce sensor data acquired from various sources, use ad hoc, local and distributed Artificial Intelligence (AI) to determine the appropriate course of action. The goal is to control or disseminate static or dynamic robotic objects that are aware of their position in the physical world in a seamless manner by providing a way to operate them as suggested by the Internet of Robotic Things (IoRT) concept [START_REF] Ray | Internet of robotic things: Concept, technologies, and challenges[END_REF].

The use of pairing technologies [START_REF] Kapitonov | Robot-as-a-service: From cloud to peering technologies[END_REF] presents a wide decentralization of robotic devices, giving them economic independence by using Peer-to-peer technologies, distributed registry and cryptography. Table 1 [START_REF] Sakovich | Fog computing vs. cloud computing for iot projects[END_REF] gives a more explicit comparison of these different approaches. The concept of Robot as a Service (RaaS) has made it easier to manage robotic functionalities. Everything is done so that the end user does not have to worry about updating the functionality. The RaaS concept directs the work related to Fog robotics towards the ROS improvement in order to favor in a first step, a large decentralization of the management of robotic functionalities. ROS is a widely used software development framework with packaging tools and communication infrastructure for robotic scenarios, often easily deployed on virtual machines running on Android devices, smart mobile gateways and especially robotic hardware [START_REF] Quigley | Ros: an open-source robot operating system[END_REF].

The work done by [START_REF] Mushunuri | Resource optimization in fog enabled iot deployments[END_REF], consists in linking the NLopt optimization library (Johnson, 2014) (which implements several types of optimization algorithms) to ROS. Through this library, they optimize the communication infrastructure of robotic scenarios, thus promoting the management of limited resources in various other mobile IoT environments.

Simulation in Fog-Based Environments

Simulating application execution in LSDS is very common these days generally speaking in Cloudbased environments. In this Section we will discuss some of the most relevant and existing simulators in Fog environments.

iFogSim [START_REF] Gupta | ifogsim: A toolkit for modeling and simulation of resource management techniques in the internet of things, edge and fog computing environments[END_REF] is a simulation framework allowing to model IoT, Edge and Fog environments and to measure the impact of resource management techniques on latency, energy consumption, network congestion and the cost of using resources. it also demonstrates strong scalability in terms of memory consumption and execution time. However, it is limited with regard to the network communication model, since it does not deal with device-to-device communications because it considers a hierarchical organization of Fog equipment.

MyiFogSim [START_REF] Lopes | Myifogsim: A simulator for virtual machine migration in fog computing[END_REF]) is a simulation framework that is based on iFogSim, so it inherits all its assets. In addition to these strengths, it supports the analysis of the impact of virtual machine migration policies on the quality of service of applications. It therefore also inherits the network communication modeling limitations present in iFogSim.

YAFS [START_REF] Lera | Yafs: A simulator for iot scenarios in fog computing[END_REF]) is a simulation framework for IoT applications in Fog environments. It has broader functionality than iFogSim and Myi-FogSim, in that it deals with both Fog resource management and IoT application placement on Fog devices. In addition, its modeling of Fog environments is much more realistic and strong. However, its cost model is not available, its network model is limited and customization of the mobility model is not supported [START_REF] Kunde | Comparison of simulators for fog computing[END_REF].

Fog Resource Needs Estimation for IoT Applications

Fog Computing brings a recent approach to data management, resource management is very little discussed in Fog-based environments. Previous studies are limited to Cloud and Edge scenarios only.

Reference [START_REF] Aazam | Fog computing micro datacenter based dynamic resource estimation and pricing model for iot[END_REF] presents a service oriented resource management model for IoT devices, through Fog, which can help in efficient, effective, and fair management of resources. Their work mainly focused on customer type and device based resource estimation and pricing, but they do not take into account IoT device capabilities, like CPU and memory in their model.

Reference [START_REF] Garcia | Do we really need cloud? estimating the fog computing capacities in the city of barcelona[END_REF] analyzes and estimates the real Fog computational capacity in the city of Barcelona, which is one of the world's most popular smart cities. Thus, they do not focus on specific applications or specific types of applications like robotics applications.

Reference [START_REF] Kattepur | A-priori estimation of computation times in fog networked robotics[END_REF] mainly deals with the estimation of the execution times of robotics applications in heterogeneous Fog environments. They implement their approach on applications pro-cessing images and videos. However, they do not address the estimation of Fog resource requirements, moreover the model does not take into account the CPU degradation of Fog resources.

This work primarily discusses the Fog resource requirements of robotic devices of the Edge layer, since this kind of deployment is a compromise for the majority of the comparison criteria (good latency, better security, good data management) in Table 1. It differs from existing works about multi-objectives optimization in that it does not focus on resource utilization optimization but on resource needs estimation (whose use can be further optimized). This estimation targets the reduction of two metrics simultaneously: the makespan and financial Fog resource utilization cost; Through robust and reliable simulation based on the well-known simulator, SimGrid3 [START_REF] Casanova | Versatile, scalable, and accurate simulation of distributed applications and platforms[END_REF], which has one of the best network communication models to date.

SERFRI'S CORE PRINCIPLE

To estimate the resource requirements of an application based on the simulation of its deployment in a FRP, it would be necessary to: (i) model this FRP, (ii) model the robotic application, (iii) set up the simulation and (iv) simulate application deployments on multiple platform scenarios (since simulation is less time-consuming) and compare the metrics (makespan, cost) to determine the best type of resources to use (among numerous types). These steps represent the phases of SERFRI and are the subject of presentation in this section.

Fog Robotics Platform Modeling

Figure 6 takes up the high-level architecture of a robot already presented in Figure 1 without on-board computing and storage capabilities. In this work, we deal only with homogeneous Fog environments containing static resources.

The first objective of this approach is to determine the global resource needs (Workload4 (W ), and Storage (S) capacity) of the application, then to estimate the configuration of resources (storage capacity, and computing power) which satisfy at best the global needs and simultaneously compromises the makespan and the cost of using Fog resources.

To parameterize this situation, we represent a Fog resource of order j by two types of information: • the fixed parameters (without dependency on its activity): the cost of using the resource per unit of time (c j ); resource memory (r j ); and its storage capacity (s j ).

• the calculated parameters (depend on its activity): the calculation speed of the node (v j ), which is the product of the number of cores, the clock frequency and the number of operations per instructions cycle of the Fog node; the start computing time (b j ); the completion time of all the requests entrusted to this resource during the activity (m j ). Finally, a Fog resource of order j (FN j ) can be written according to the expression FN j = (r j , s j , c j , v j , b j , m j ).

Considering n Fog resources, C the resource utilization cost and M the makespan of the Fog layer, then appears the following constraints and relations:

S < n ∑ j=1 s j
(1)

C = n ∑ j=1 c j × m j (2) M = Max({b j + m j } n j=1 ) -Min({b j } n j=1 ) (3)

Robotic Application Modeling

Figure 6 suggests that the application be divided into two sub-applications: a Client application to issue requests and a Server application to process these requests. In the context of robotics, these queries are:

• periodicals: for each pair of consecutive requests, there is a constant duration between their emission dates (Equation ( 5));

• asynchronous: sending and processing a request does not prevent the same operations for other requests;

• heterogeneous: the quantity of data communicated varies according to the requests.

The i-order request (R i ), is characterized by: a date of issue (e i ), a quantity of data that it communicates (d i ). Finally we write: R i = (e i , d i ). And for a set of q periodical (t-periodical (Equation ( 5))) and asynchronous (Equation ( 6)) requests we define a function T giving the completion time of a i-order request processed by the j-order Fog resource by:

T : {R i } q i=1 × {FN j } n j=1 → R + (R i , FN j ) → t i j (4)
So that :

t i j = value > 0 if R i is processed by FN j 0 else. e i+1 = e i + t, ∀ 1 ≤ i < q (5) e i+1 ≤ e i + t i , ∀ 1 ≤ i < q (6)
The workload of the i-order request will then be defined by Equation ( 7). This workload being the same for the same application even if the resources specification change, helps in compute the corresponding completion time of a request by the new resource, having the completion time of the previous one.

W i j = t i j × v j (7)
Figure 6 suggests also three kind of events5 : Sending event (which sends the requests), a Computing event (which processes the data communicated by the requests) and the Receiving event (Which load to the data source/an IoT device the results of requests). These events happen in this order: Sending, computing and Receiving; The robot runs the Sending and the Receiving events when the resources handle the Computing event.

Simulation Settings

Even though there are a very large number of simulators in Fog based environments, there are some limitations on the most known of them, as they are based on other simulators that are limited in terms of network communication models for example. To circumvent this limitation, we used SimGrid which is very powerful for simulations in many types of LSDS and in Cloud and Grid computing in particular. SimGrid has an accurate network communication model, which helps make the simulation realistic [START_REF] Velho | Accuracy study and improvement of network simulation in the simgrid framework[END_REF]. It has been used in hundreds of projects and helps in both distributed applications and LSDS evaluation. Moreover, there is a very large community around him since it represents more than twenty years of code improvement.

A SimGrid study consists of:

• The description of the platform concerned: defining the zones of the platform studied, the entities (hosts) which composed them as well as a description of the specifications of these hosts and the connections between them. This can be done through an XML file or an assembly of code instructions, built from the following parameters: resource specifications {FN j } n j=1 , latency (lt) and bandwidth (bw) of the connection between the IoT layer and the Fog layer;

• The description of the application: a set of functions or events that must occur during the execution of a deployment scenario. In SimGrid, they are represented by Classes or Event Functions;

• The deployment: which specifies which host executes which event. It can be expressed through an XML file or an assembly of code instructions, built taking into account the workflow described above (Sending -Computing -Receiving) executed by the corresponding hosts (robot -Fog resources -robot).

Our simulation is therefore likened to a simulate() function which finally takes as parameters: {FN j } n j=1 , lt, bw, the requests sending period t and the function T in Equation ( 4)). This function (simulate()) returns again {FN j } n j=1 , with suitable values of calculated parameters of each resource related to the execution of the application, since before the execution those parameters were all set to 0. However the workload generated by a request is not known, it is then necessary to find what the function T corresponds to in our simulation. T can be assimilated (depending on the use case) to a mathematical formula, machine learning model or mapping object and so on, that give exactly or estimate the completion time of a request by a Fog resource according to the characteristics of the request and those of the Fog resource. Then an estimation of the query workload is given by Equation ( 7)).

Algorithm of the Approach

The objective of SERFRI is to propose a resource type and the optimal number of resources of this type among a set of resource types, in order to simultaneously reduce as far as possible the makespan and the global cost of resources use. The goal is to apply it on each type of resource studied and select the one that optimizes the makespan and the cost. Indeed, according to a specific resource type, for each resource of order j, FN j = FN 1 and FN 1 = (r, s, c, v j , b j , m j ) where r, s, c are constant. We do not deal with resource management and assume that the scheduling used for the distribution of requests to resources is the Round-Robin algorithm, and that the scheduling of requests on a node follows the logic of the FIFO algorithm.

The problem that this approach tries to solve can be reformulated as follows: Assuming we have l types of resources ({hn j } l j=1 ) and p configurations 6 ({con f ig k } p k=1 ) possible for each type of resource, what would be the pair (hn j , con f ig k ) that would best reduce the makespan and financial cost of using the resources (by serving requests in order of arrival, without attempting on line scheduling on them)? The results of the simulations of each pair in {hn j } l j=1 × {con f ig k } p k=1 can be represented by two matrices of order p × l: the makespan matrix A (A k j is the makespan when con f ig k resources of type hn j are used) and the financial cost matrix B (B k j is the financial cost of using con f ig k resources of type hn j ).

The problem as formulated can be solved by a Pareto optimization [START_REF] Ngatchou | Pareto multi objective optimization[END_REF], whose solution (hn j 0 , con f ig k 0 ) verifies that there is no pair (A k j , B k j ) so that A k j < A k 0 j 0 and B k j < B k 0 j 0 . This can lead mathematically speaking to several solutions, some of which may not be acceptable in practice.

Algorithm 1 of the SERFRI approach aims at selecting the most practically acceptable solution in the theoretical Pareto front (of couples verifying the Equation ( 1)). It proceeds by calculating the ratios (Equation ( 8)) relative to a basic resource type and a basic configuration for each metric and by performing optimization operations on them. This algorithm highlights some functions that we define below:

φ : M p×l → M p×l M → { M k j -min cell val(M) min cell val(M) } (k, j) ∈ 1,p × 1,l (8)
In the φ function definition, min cell val(M) returns the minimum value of the elements of M. This function (φ) returns a matrix where each cell has for value, the relative error of the value of the cell of the same position as this one in matrix M, compared to the cell of matrix M whose value is minimal.

σ : M p×l 2 → M p×l (M, N) → {max(M k j , N k j )} (k, j) ∈ 1,p × 1,l
(9) The σ function returns a matrix where each cell has, 6 the possible numbers of resources for value, the maximum value of two cells of the same position as this one in the matrices M and N.

δ : M p×l 2 → M p×l (M, N) → {min(M k j , N k j )} (k, j) ∈ 1,p × 1,l
(10) In contrast to the function σ, the function δ returns a matrix where each cell has a value equal to the minimum value between two cells at the same position in matrices M and N.

Algorithm 1 Determining sufficient specifications for Fog resources.

Require: A ∨ B Ensure: SP 0 λ ← σ(φ(A), φ(B)) µ ← δ(φ(A), φ(B)) SP ← argmin k, j ({λ k j } (k, j) ∈ 1,p × 1,l ) SP 0 ← argmin k, j ({µ k j } (k, j) ∈ SP ) return SP 0
This algorithm is indeed of complexity O(p × l), which is faster than first applying a Pareto optimization of exemplary complexity (O(p × l × log(p × l)) [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-ii[END_REF]) before implementing a selection algorithm to find the solution of the front that is in practice acceptable. From the pair (k 0 , j 0 ), we then have the type of resources (hn j 0 ) and the number of resources of this type (con f ig k 0 ) which compromises the simultaneous reduction of the makespan and financial cost of using the resources. Mathematically it can be demonstrate that SERFRI algorithm solutions are always Pareto front solutions.

EXPERIMENTS

In this section, we first validate the simulation of the deployment of the use case expressed in Section 2. And then, we present the estimated resource type specifications (among many types of resources) to optimally set up the infrastructures of this case. The simulator7 is publicly available on GitLab and contains dataset, source code and instructions to reproduce our results.

Simulation Validation

To make a simulation reliable, it should be as realistic as possible. One way to validate it is to evaluate its accuracy and its scalability compared to the same scenario in a real environment. We evaluate the accuracy of our simulator by comparing the time-based traces obtained in real situations and in simulation, about the completion time and we evaluate the scalability by evaluating the machine learning models predicting the services' completion times when the number of cores is increasing.

Simulation Scalability Evaluation

Equation (4) (Section 4.2) presents a generic function for estimating the completion time of a i-order request on a j-order Fog resource. We said in Section 4.3 that this function can be similar to a mathematical formula, machine learning model or mapping object and so on in our simulation (For each platform setting).

In our use case, we built machine learning models trying to find the optimal correspondence between the data (e i , d i ) and t i j according to the platform setting. We thus first collected execution traces of the human face detection service on 400 images for each configuration of the Fog layer (1, 2, 3, 4 and 5 chetemi nodes). Then a clustering algorithm (KMeans in our case) was applied to 280 random occurrences of the couple (e i , d i ) in the dataset containing the 400 starting couples. The value of the predicted completion time for each triplet of a cluster was then the mean value of the completion times of the couples of this cluster. We evaluate the accuracy of the model like the percentage of predictions whose absolute errors relative to the actual values are less than a given tolerance. This tolerance should be as small as possible. The CDF on Fig. 7 gives the completion time estimation model accuracy according to the tolerance (and according to the number of nodes), and Table 3 gives the accuracy corresponding to 15 s according to the number of cores. This accuracy increases with the number of nodes, meaning that the simulation reaches the scale-out property.

Simulation Accuracy Evaluation

Figure 8 gives the evolution of the life-cycle duration (in real and in simulation) according to the order of the request sent for each configuration of the Fog layer. This life-cycle consists of three steps (extract, transform and load) and is traced for each request in simulation and in real deployment to compare both deployments. According to the number of cores, the difference between both life-cycles duration in average is short (22.97 s by computing mean value of CT MAE column of Table 2 giving an accuracy of 95% in average for the completion time estimation model) . This tolerance can be considered as short because it does not really have an impact on the global makespan, completion time and duration of the activity, since requests are asynchronous. Despite the fact that the data communicated by the queries are of the same order (because standard deviation of image size distribution (Table 2) is very lower than the minimum value of this distribution (Table 4)), the completion times in real life vary greatly (Table 2). This is because at certain times there are many more requests sharing the CPU than at other times. The simulation follows this same observation well with respect to the mean absolute errors of the completion time for each configuration of the Fog layer highlighted in Table 2.

Fog Resource Needs Estimation

To implement SERFRI on the use case presented in Section 2, we consider that Fog instances have similar characteristics to OVHCloud instances 8 designed for data analysis and data-science uses: r2-15, r2-30, r2-60 , r2-120 and r2-240. We simulate the platform scenario made-up of : 2, 4, 8, 16, 32 and 64 resources for each of OVHCloud instance types considered. We present first the simulation settings, then the pareto optimization and finally the SERFRI steps: (i) phitransformation 9 , (ii) sigma-transformation 10 and result expression.

Settings

To perform the simulations with these new instances, we do not rebuild an estimation model but we rely on the models built with the chetemi nodes. Indeed, with regard to Equation ( 7)), the workload of a request i on a resource j with speed v j is written W i j = t i j × v j , in the same way, this request processed by another type of resource of order j with speed v ′ j will have a workload that will be written W i j = t ′ i j × v ′ j . This allows us to deduce the desired completion time: -15, r2-30, r2-60, r2-120, r2-240 Table 4 gives the values of the simulator inputs. Regarding the characteristics of the instances, the level of granularity is very low, because we consider: the memory, the storage, the number of cores, the resources cost per unit of time, the clock frequency, the maximum number of operations per instruction cycle 8 Cloud pricing: Comparison of public Cloud offers-OVH, https://www.ovhCloud.com/fr/public-Cloud/prices/ 9 The application of φ function 10 The application of σ function and bandwidth in a public network. Requests sizes are randomly generated between 420000 B and 480000 B, and we keep this same distribution for all deployment scenarios of the use case. Table 5 and Table 6 present the makespan and financial costs of using Fog resources, respectively. They are therefore the matrices we referred to in Section 4.4. The graphical transcription of Table 5 (Fig. 9) shows that the makespan decreases as the number of resources increases and more so when the type of resources is more powerful. The graphical transcription of Table 6 (Fig. 10) , on the other hand, grows as the number of resources increases and more so when the type of resources is more powerful, that makes the optimization very difficult. By naively ap- plying Pareto optimization to reduce the makespan and the cost of using resources as much as possible, we obtain the following Pareto front: {(r2-15, 8), (r2-15, 64), (r2-60, 64), (r2-240, 8), (r2-240, 32), (r2-240, 64)}. 8)) tables from Table 5 and Table 6 respectively. In other words, we look for the smallest element in Table 5 and we calculate for each element of this table the relative error with respect to this smallest element, and the same thing is done for the Table 6. This allows to bring the makespan and the financial cost into the same dimensional space. The makespan decreases according to the number of resources and the power of resources, but it is an inverted observation for the financial cost ratio. The Pareto front is the same front presented in Section 5.2.2.

t ′ i j = t i j ×v j v ′ j .

The Pareto Optimization

Phi-transformation

Sigma and delta transformations

The sigma-transformation applies the σ function (Equation 9) on normalized matrices. Maximizing these matrices (Table 7 andTable 8) leads to the creation of the Table 9, according to the Equation ( 9), from which we can determine the minimum element representing the couple (r2-240, 64). Then among results cells of the sigma-transformation we choose the one minimizing the equivalent cells in the deltatransformation (Table 10), therefore we extract the solution of the Pareto front that in practice is the most efficient in terms of time and financial cost: (r2-240, 64). It is indeed, because compared to the other Pareto solutions, it saves much more time than it costs.

CONCLUSION

In this article we address the estimation of Fog resource requirements (number of resources, memory, computing and storage capacity of each resource) of robotics applications. We then proposed an approach (SERFRI) for estimating static and homogeneous Fog resource needs based on the simulation of the deployment of Fog Robotic infrastructures. This approach aims at providing the best resource type and the best number of them, compromising the reduction at best of the makespan and the one of the financial cost of using these Fog resources. We implement the approach on an existing use case (Reachy) which could also have industrial vocation. The validation of the simulator on which the approach is based consisted mainly on the validation of the machine learning models used (accurate on average at 95 % with a prediction margin of 23 s for each request, not really impacting the global life-cycle duration of the application) to estimate the request completion time. This approach can therefore, also greatly help industrial companies (and/or those who base their activities on the use of robots) to have an approximate idea of the budget required for the activities of their robots on Fog resources. It utilizes reliable and realistic simulation results to determine a subset of the Pareto front with linear complexity relative to the number of simulations required. This approach also saves the financial budget and the use of the robot's battery, allowing it to perform more tasks. However, a similar approach should be implemented in the generic case of platform, where Fog resources are heterogeneous and robots are mobile. The study indeed deserves reflection because it is not just an implementation of the knapsack problem as in the case of static heterogeneous resources.

Figure 1 :

 1 Figure 1: Five main components of a robotic architecture.

Figure 2 :

 2 Figure 2: Comparing the ETL process of deploying a human face recognition service locally and remotely.

Figure 3 :

 3 Figure 3: Reachy's Full kit picture: 3D printed MJF Painted -Flexible PU molded -Aluminium Construction; 670x1800x200mm with arms outstretched; 180 W of Max Power consumption and 6.2 kg of weight.

Figure 4 :

 4 Figure 4: The Reachy use case architecture.

Figure 5 :

 5 Figure 5: The Reachy use case workflow modeling.

Figure 6 :

 6 Figure 6: Fog Robotic Platform modeling.

Figure 7 :

 7 Figure 7: CDF according to the mean absolute error between real completion and prediction. Each chetemi node has 20 cores. Reachy sends 400 requests to Fog nodes.

  5 chetemi nodes: 100 processing cores.Figure 8: Requests life-cycle durations. Comparison between Real and Simulation deployments.

Figure 9 :

 9 Figure 9: Makespan by the number of resources for each type of resource.

Figure 10 :

 10 Figure 10: Financial cost by the number of resources for each type of resource.

Table 1 :

 1 Comparison of distributed system paradigms.

	Paradigm	Latency	Data management	Scalability	Computing Capabilities
	SOA and Cloud	High	Centralized	High	High
	Edge
	computing	High	Decentralized	Low	Low
	(IoT)
	Fog computing	Low	Decentralized	Low	Low
	Peer-to-peer	Very Low	Widely decentralized	Very Low	Very High

Table 2 :

 2 Mean Absolute Errors (MAE) and Standard deviations (STD) according to the number of cores. CT stands for Completion Time.

	Number of cores	CT MAE (s) CT STD (s)	Image size STD (KB)
	20	22.83	26.07	23.76
	40	21.01	17.81	25.79
	60	24.12	15.95	16.8
	80	20.25	13.4	23.7
	100	26.62	9.39	21.8

Table 3 :

 3 completion time accuracy for a tolerance of 15 s, according to the number of nodes.

	Number of cores Accuracy (%)
	20	80.83
	40	80.92
	60	88.3
	80	90.0
	100	96.5

Table 4 :

 4 Simulation settings.

	Number of requests	5000
	Frequency
	(number of	7
	requests/second)
	Maximum number of instances	64
	Bandwidth between
	Edge and	20
	Fog Layers (Mbps)
	Latency between
	Edge and	0.005
	Fog Layers (s)
	Data sizes range (B)	420000 -480000
	Response json size (B) 62
	Instances	r2
	(ovh Cloud instances)

Table 5 :

 5 Makespan (in second) by resource type and configuration.

	Configuration	R2-15	R2-30	R2-60	R2-120	R2-240
	2	82182.11	82182.11	41138.88	20610.64	10253.34
	4	41138.88	41138.88	20610.64	10592.45	5591.26
	8	20610.64	20610.64	10592.45	5591.26	2690.48
	16	10592.45	10592.45	5591.26	2772.42	1530
	32	5591.26	5591.26	2772.42	1530	869.8
	64	2772.42	2772.42	1530	869.8	494.5

Table 6 :

 6 Financial cost (in euro) by resource type and configuration.

	Configuration	R2-15	R2-30	R2-60	R2-120	R2-240
	2	4.43	5.12	5.01	5.04	4.95
	4	4.45	5.14	4.98	5.09	5.13
	8	4.41	5.1	5.01	5.14	5.05
	16	4.46	5.15	5.03	5.12	5.22
	32	4.45	5.14	5.04	5.28	5.44
	64	4.44	5.13	5.16	5.46	6.4

Table 7 :

 7 Makespan ratios by resource type and configuration.

	Configuration	R2-15	R2-30	R2-60	R2-120	R2-240
	2	165.19	165.19	82.19	40.68	19.73
	4	82.19	82.19	40.68	20.42	10.31
	8	40.68	40.68	20.42	10.31	4.44
	16	20.42	20.42	10.31	4.61	2.09
	32	10.31	10.31	4.61	2.09	0.76
	64	4.61	4.61	2.09	0.76	0

Table 8 :

 8 Financial cost ratios by resource type and configuration.

	Configuration	R2-15	R2-30	R2-60	R2-120	R2-240
	2	0	0.16	0.14	0.14	0.12
	4	0.01	0.17	0.13	0.15	0.16
	8	0	0.16	0.14	0.17	0.15
	16	0.01	0.17	0.14	0.16	0.18
	32	0.01	0.17	0.14	0.2	0.23
	64	0.01	0.16	0.17	0.24	0.45

Table 7 and

 7 Table 8 represent the normalized (phitransformed tables obtained by applying φ function of Equation (

Table 9 :

 9 Maximized ratios by resource type and configuration.

	Configuration	R2-15	R2-30	R2-60	R2-120	R2-240
	2	165.19	165.19	82.19	40.68	19.73
	4	82.19	82.19	40.68	20.42	10.31
	8	40.68	40.68	20.42	10.31	4.44
	16	20.42	20.42	10.31	4.61	2.09
	32	10.31	10.31	4.61	2.09	0.76
	64	4.61	4.61	2.09	0.76	0.45

Table 10 :

 10 Minimized ratios by resource type and configuration.

	Configuration	R2-15	R2-30	R2-60	R2-120	R2-240
	2	0	0.16	0.14	0.14	0.12
	4	0.01	0.17	0.13	0.15	0.16
	8	0	0.16	0.14	0.17	0.15
	16	0.01	0.17	0.14	0.16	0.18
	32	0.01	0.17	0.14	0.2	0.23
	64	0.01	0.16	0.17	0.24	0

https://simgrid.org/usages.html

in FLOPs(FLoating Point OPerations)

An event means an action (Sending, Computing or Receiving) in this case

Our simulator is available at: https://gitlab.inria.fr/lndjieng/fridsim/-/tree/master

ACKNOWLEDGEMENTS

We are deeply grateful to all those who played a role in the success of this project. We would like to thank Interreg AiBLE and REACT-EU UV-Bot for their support throughout the research process.