Agustín Somacal
email: a.somacal@aristas.com.ar

Yamila Barrera
email: y.barrera@aristas.com.ar

Leonardo Boechi
email: lboechi@ic.fcen.uba.ar

Matthieu Jonckheere

Vincent Lefieux
email: vincent.lefieux@rte-france.com

Dominique Picard
email: picard@math.univ-paris-diderot.fr

Ezequiel Smucler
email: e.smucler@aristas.com.ar

dpto. N Piso

Uncovering differential equations from data with hidden variables

Keywords: Differential equations, Dynamical systems, Lasso, Latent variables, Machine learning

come

Introduction

Many branches of science are based on the study of dynamical systems. Examples include meteorology, biology and physics. The usual way to model deterministic dynamical systems is by using (partial) differential equations. Typically, differential equations models for a given dynamical system are derived using apriori insights into the problem at hand; then the model is validated using empirical observations. In an era in which massive data-sets pertaining to different fields of science are widely available, an interesting problem is whether it is possible for a useful differential equations model to be learned directly from data, without any major modeling effort required by the researcher.

The SINDy method is an approach of sparse identification of nonlinear dynamical systems that consists of linking the dynamical system discovery problem to a statistical regression problem [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF], [START_REF] Rudy | Data-driven discovery of partial differential equations[END_REF] and [START_REF] Quade | Sparse identification of nonlinear dynamics for rapid model recovery[END_REF]). The main idea of the SINDy method is to consider a set of differential operators (possibly in both time and space if appropriate), discretize them, for example by using finite differences, and then regress the outcome of interest on the discretized differential operators. By solving the regression problem using an ad-hoc thresholded least-squares algorithm, they can build sparse, interpretable models, that use mostly low order derivatives. The authors explored the applicability of their method on simulated data, but only in situations in which all the variables of the simulated models are observed.

Our goal is to extend the SINDy model for the case in which not all relevant variables are observed, that is, in cases in which the main variable of interest depends on other variables of which no measurements are available. As an example of application we consider the climate time series of the Réseau de Transport d' Électricité (RTE). RTE is the main electricity network operational manager in France, who is interested in understanding the behavior of climate time series because of their impact on energy consumption. RTE uses high-level simulations of hourly temperature series to study the impact different climate scenarios have on electricity consumption, and hence on the French electrical power grid. The underlying simulations are based on the Navier-Stokes equations and include variables as wind velocity, density, pressure, etc. The resulting dynamic system is known to be chaotic, see [START_REF] Kida | Route to Chaos in a Navier-Stokes Flow[END_REF]. For that reason, our goal is to learn a system of differential equations that adequately models the dynamics of the temperature time series if the only observable variable is temperature, that is, if pressure, wind velocity, etc, are hidden variables.

To accommodate the possibility of hidden variables we note that, for a large class of dynamical systems, it is possible to reconstruct a trajectory (equivalent to the original one) given only one of the model variables, using its higher order derivatives ( [START_REF] Takens | Detecting strange attractors in turbulence[END_REF]).

Related to our approach is the Generalized Polynomial Modeling method (known as GPoMo), that addresses the recovery problem via a combinatorial search among a predefined set of polynomial functions of the observable variables [START_REF] Mangiarotti | Polynomial search and global modeling: Two algorithms for modeling chaos[END_REF] and [Bongard and Lipson, 2007a]. The GPoMo method proceeds by choosing iteratively a family of combination of terms that minimize the Akaike or the Bayesian information criterion. Finally, it returns the set of best models. The authors also discuss the ability of their algorithm to find equations able to capture the dynamics when only some variables are observed. However, we will show that unfortunately this approach does not scale to large problems. Other approaches for learning dynamical systems from data available in the literature, 2 such as those based on symbolic regression ( [Bongard and Lipson, 2007b]), also have the drawback of being too computationally expensive.

The article is organized as follows. In Section 2 we review the GPoMo and SINDy methods and afterwards, describe our methodology (named L-ODEfind) in detail. Section 3 presents the results of our experiments. In particular, in Section 3.2, we compare the performance of GPoMo and L-ODEfind in recovering differential equations using empirical data in the case in which all relevant variables are observed. In Section 3.3 we compare them in the harder case in which at least one relevant variable driving the dynamical system is latent. We apply our proposed method to real world temperature times series in Section 3.4. Finally, in Section 4 we discuss future work and possible extensions.

Methods

GPoMo: Differential equations recovery as a combinatorial search problem

GPoMo is a method proposed and implemented by [START_REF] Mangiarotti | Polynomial search and global modeling: Two algorithms for modeling chaos[END_REF] that addresses the differential equations recovery problem via a combinatorial search in the space of differential equations that can be expressed as polynomial functions of the observed variables. The method uses a genetic algorithm in which at each step new test models are generated by randomly choosing some polynomial terms to be included in the equations. This choice is made by making small variations (take or add a few terms) over the best previously seen models. Then, to select the winning models at each step, they are integrated and compared to the original data. This combination of combinatorial search and integration steps makes the method slow (as we will see in the experiments). Moreover, except for polynomial combinations of the variables, it does not allow other types of regressors to be included, such as functions of the time variable. To overcome the aforementioned limitations of GPoMo, we based our approach on SINDy, which uses sparse regression to discover governing physical equations from measurement data. We briefly review the SINDy method in the following section.

SINDy: Differential equations recovery as a linear regression problem

In our work we focus on ordinary differential equations even though SINDy can also tackle partial differential equations. In particular, consider a dynamical system represented by functions f 1 (t) . . . f H (t) satisfying a set of differential equations of the form

Df = U(f , Ef ), (1)
where f = (f 1 , . . . , f H ), D, E are differential operators in the temporal variables (t) and U : R J → R H is an unknown map. Suppose we have a series of T equally spaced in time measurements, that is, we observe f h (t i ) i ∈ {1, . . . , T }, h ∈ {1, . . . , H}. An example of a dynamical system we will study in this paper is the classical Rossler system (2) ( [Rossler, 1976]). This system was originally designed to have similar properties and be simpler than the Lorenz system ( [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]) which was, at the same time, a simplified model for atmospheric convection. The system is given by

df 1 (t) dt = -f 2 (t) -f 3 (t), df 2 (t) dt = f 1 (t) + αf 2 (t), ( 2
)
df 3 (t) dt = β + f 3 (f 1 (t) -γ),
for constants α, β, γ. This system can be written in the form (1) by taking

f = (f 1 , f 2 , f 3 ), D = (d/dt, d/dt, d/dt) and U = (U 1 , U 2 , U 3 ) where U 1 (v 1 , v 2 , v 3 ) = -v 2 -v 3 , U 2 (v 1 , v 2 , v 3 ) = v 1 + αv 2 and U 3 (v 1 , v 2 , v 3 ) = β + v 3 (v 1 -γ).
For certain values of the parameters α, β, ρ, the system is known to have chaotic solutions [START_REF] Amaral | Piecewise affine models of chaotic attractors: The rossler and lorenz systems[END_REF][START_REF] Ibrahim | Chaotic behaviour of the rossler model and its analysis by using bifurcations of limit cycles and chaotic attractors[END_REF].

Suppose now that we have access to a particular time series that was generated by this system.

Our objective is to find some system of differential equations that can explain the behaviour of the measurements. The SINDy method works by choosing a large dictionary of functions and regressing discretisations of ∂f1 ∂t , . . . , ∂f H ∂t on the dictionary. The dictionary in question can be formed, for example, by collecting polynomial powers of f h , h = 1, . . . , H, spatial derivatives of f 1 . . . , f H and trigonometric functions of t. A concrete simple example of such a dictionary in the case in which H = 1 is the following:

A = t, t 2 , sin(t), f 1 , f 2 1 , tf 1 .
Of course in practice all derivatives are replaced by the corresponding finite differences taken from the measurements represented in f . Having chosen a dictionary, we let A = (A 1 , . . . , A p ) be the vector collecting all members of the dictionary.

Using the observations of the dynamical system, a regression model can be fitted to find the combination of the elements of the dictionary of functions that adequately explains the behavior of ∂f1 ∂t , . . . , ∂f H ∂t . That is, we look for a vector of regression coefficients c = (c 1 , . . . , c p ) such that for all

h = 1, . . . , H ∂f h (t) ∂t ≈ p i=1 c i,h .A i (t). (3)
The regression model has to be learned using the available data. This regression problem could be solved in principle using least-squares. However, the ordinary least-squares regression estimator is ill-defined in cases in which the number of predictor variables p is larger than the number of observations. Since the analyst is usually uncertain about the number of elements in the dictionary needed to adequately model the system of interest, the method used to solve the regression problem at hand should allow for large number of predictor variables (possibly larger than the number of observations) and automatically estimate sparse models, that is, generate accurate models that only use a relatively small fraction of predictor variables. The Lasso regression technique is perfectly suited for this task. The Lasso is an 1 -regularised least-squares regression estimator, defined as follows. For

h = 1, . . . , H such that f h is observable we let c * h = arg min c h ∈R p l ∂ n f h ∂t n (t l ) - p i=1 c i,h A i (t l ) 2 + λ||c h || 1 , ( 4
)
where λ > 0 is a tuning constant, measuring the amount of regularization. It can be shown [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF]) that the 1 penalty encourages sparse solutions and that the larger λ > 0 the sparser the solution vector c * h will be. In practice, λ is usually chosen by cross-validation. Note that any other sparse regression technique could have been used to estimate the coefficients. We prefer the Lasso due to its simplicity and the wide availability of efficient algorithms to compute it. See for example [START_REF] Friedman | Pathwise coordinate optimization[END_REF].

The main assumption behind this methodology is that the dynamical system that generated the data at hand can, in reality, be at least approximated using a sparse model, that is, that the vectors c * h in (4) are either exactly or approximately sparse. This hypothesis is known to hold for several dynamical systems of interest in different fields of science. See [START_REF] Quade | Sparse identification of nonlinear dynamics for rapid model recovery[END_REF]. If the hypothesis holds, we can expect the Lasso estimates to select only a few elements of the dictionary, namely, those that do a good job at explaining variations in the response variable [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF]).

In [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF], [START_REF] Rudy | Data-driven discovery of partial differential equations[END_REF] and [START_REF] Quade | Sparse identification of nonlinear dynamics for rapid model recovery[END_REF] the authors propose to use an ad-hoc linear regression estimator based on iteratively thresholding the least-squares estimator and applied this method only to first order systems. Through extensive numerical experiments, they show that this methodology is able to learn systems of partial differential equations that adequately model the dynamical system that generated the data. Unfortunately, if some variables are latent, that is, if one is unable to measure at least one of f 1 , . . . , f H , the approach described above cannot be used directly. Next, we describe a way of extending this methodology to deal with the case in which some variables are latent.

Our proposal: L-ODEfind

To accommodate the possibility of latent variables we note that, for a large class of dynamical systems, it is possible to reconstruct a trajectory (equivalent to the original one) given only one of the model variables, using its higher order derivatives ( [START_REF] Takens | Detecting strange attractors in turbulence[END_REF]). Moreover, we recall that in the case of a linear system of n ordinary differential equations there is an equivalence between this multidimensional system and a single differential equation of order n, which we can interpret as latently including the information of the other n -1 unobserved variables. [START_REF] Mangiarotti | Polynomial search and global modeling: Two algorithms for modeling chaos[END_REF] also makes use of higher order time derivatives to deal with unobserved variables.

Based on this ideas, we propose to augment the methodology developed in [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF], [START_REF] Rudy | Data-driven discovery of partial differential equations[END_REF] and [START_REF] Quade | Sparse identification of nonlinear dynamics for rapid model recovery[END_REF] by choosing the target variable to be a higher-order time derivative, to tackle situations in which not all relevant variables are observed. We estimate the coefficients of the dynamical system using the Lasso estimator (4). As mentioned earlier, we chose to use the Lasso due to its simplicity, the abundance of theoretical guarantees on its performance [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF]) and the availability of efficient algorithms to solve the convex optimization problem that defines the estimator. The choice of the tuning constant λ in ( 4) is done by 10-fold cross validation, using the LassoCV method from sklearn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] with a maximum number of steps equal to 10000, and 100 candidate λs. After estimating the regression coefficients, we build a forecasting method by integrating the retrieved differential equation. We call this method Latent ODE find (L-ODEfind).

For instance, suppose we have observed a single time series f (H = 1). If we choose as a target variable the third time derivative and we use polynomial combinations up to degree 2 of the series and derivatives up to order 2 as regressors, equation ( 3) becomes

∂ 3 f (t) ∂t 3 ≈ p i=1 c i A i (t), (5)
and

A = 1, f, ∂f (t) ∂t , ∂ 2 f (t) ∂t 2 , f 2 , ∂f (t) ∂t 2 , ∂ 2 f (t) ∂t 2 2 , f ∂f (t) ∂t , f ∂ 2 f (t) ∂t 2 , ∂f (t) ∂t ∂ 2 f (t) ∂t 2 .
This regression model is then fitted using the Lasso (4), as described earlier.

Evaluation with predictions

Evaluating the equations found by these methods can be challenging. The evaluation of a method that aims at recovering the differential equation behind the observed data needs to be done in different ways depending on the information available.

In the case of a simulation where the coefficients of the differential equation are known and the variables are fully observed, the adjusted coefficients and the true ones can be compared using the mean squared error. If not all variables are observed, the adjusted coefficients refer to a different differential equation that in most cases cannot be obtained analytically (exceptions are, for instance, linear ODEs). In this case, the coefficients comparison cannot be done and another way of evaluating the method is needed. The same happens when working with real world data where the differential equation behind is not known. In order to evaluate a method in this context, we propose to integrate the fitted differential equation to make predictions for different time horizons and compare them to the observed data by using the symmetric mean absolute percentage error (SMAPE). Using a time horizon of n, the SMAPE is defined as:

SM AP E(A, F ) = 1 n n t=1 |F t -A t | (|A t | + |F t |)/2
, where A t is the real value and F t the forecasted value. If the predictions made by method A are better (lower values of SMAPE) than the ones made by method B, the differential equation found by the method A is deemed better than the one found by method B.

Implementation details

All the experiments in this paper were performed using Python 3.8, except for GPoMo , which was performed using R 3.6 [R Core Team, 2020] and the GPoM package [START_REF] Mangiarotti | Recherche pour le Développement[END_REF]. For L-ODEfind we use our own Python implementation (available in https://github.com/agussomacal/ L-ODEfind). The integration of the differential equations was done using odeint from python Scipy library [START_REF] Virtanen | SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF].

Results

In this section we first compare L-ODEfind with GPoMo for the problem of learning an ordinary differential equation with no hidden variables. Then, we compare the performance of these methods in simulated systems with hidden variables and finally in temperatures series provided by RTE.

Names abbreviation

The names abbreviation used in this section can be found in Table 1. Target time derivative is the degree of the fitted differential equation and poly degree is the maximum degree of the polynomial combinations of the derivatives used as regressors. An example of the regression problem with target derivative 3 and poly degree 2 can be found in equation ( 5). Notice that the number in the model name refers to the degree of the differential equation.

Method

Simulated data with fully observed variables

We compare the performance of L-ODEfind with that of GPoMo for the task of learning an ordinary differential equation with no hidden variables. Note that since in this case there are no hidden variables, our method coincides with SINDy. The goal of this comparison is to highlight the fact that, because L-ODEfind solves a continuous optimization problem and GPoMo approximately solves a combinatorial optimization problem, L-ODEfind can be orders of magnitude faster that GPoMo.

We generated data using the Lorenz attractor equations:

dx dt = σ(y -x), dy dt = x(ρ -z) -y, ( 6
)
dz dt = xy -βz,
using the coefficients σ = 10, ρ = 28 and β = 8 3 , for twenty different random starting conditions, with distribution N (0, 1). The differential equation was integrated using a discretization step of 0.01.

Then, we fitted the resulting datasets using GPoMo and L-ODEfind. For GPoMo we set a maximum number of steps of 10240. We compare the true coefficients of the differential equations with the ones found by each model using as a metric the mean square error (MSE). All the coefficients were taken into account for computing the MSE, including the ones that are supposed to be zero. We average the MSE corresponding to different random starting conditions and report this as our goodness of fit metric.

Figure 1 shows that L-ODEfind method is nearly two orders of magnitude more accurate than the best GPoMo case. On the other hand, L-ODEfind is also orders of magnitude faster, taking less than 10 seconds to compute accurate approximations of the true coefficients while GPoMo takes hours.

Simulated data with hidden variables

In this section we use three ODE systems, an oscillator, and the Rossler and Lorenz systems, as examples to evaluate and compare, using the methodology explained in 2.4, the accuracy of several L-ODEfind and GPoMo models. For a given differential equation, we generated time series of length 5000 points and integration step of dt = 0.01 for 20 different initial conditions. Following the methodology described in 2.4 we fit each series with GPoMo or L-ODEfind and then integrate the equation found in each case to predict the values of the series for several time horizons and compute the corresponding SMAPE. The number of maximum iterations for GPoMo is set to 5120.

Example 1: Oscillator

We start by considering an oscillator, which is a first order linear system of two variables. An equation of order two involving only one of the variables can be derived. So, in this case, the problem of having hidden variables can be effectively tackled by choosing a higher order time derivative (order two) as the target. We used the two variables oscillator equation which can be written in general form as:

d dt   x y   =   a b c d   •   x y   (7)
where x and y are the variables and a, b, c, d are the coefficients of the linear equation that links the variable with its derivatives. If we only had access to the variable x, this system could be rewritten in a second order differential equation taking the form:

d 2 x dt 2 = (a + d) dx dt -(ad -bc)x = β dx dt + αx. ( 8
)
In our experiments we set a = 0.1, b = -1, c = 1, d = 0 (α = -1, β = 0.1) and we only observed the variable x so the corresponding second order equation derived as in 8 is 0.1 dx dt -x. We can see in Figure 2 (a) that using the second derivative in time as target gives the lowest prediction error for both GPoMo and L-ODEfind, as expected from equation ( 8). We also see that the model found by L-ODEfind manages to approach the true model much better than GPoMo as it maintains an SMAPE below 0.02 in horizons were GPoMo has already arrived to 0.1. Moreover, when looking at the coefficients found by GPoMo (α = -0.99506 ± 0.00372 and β = 0 ± 0) and L-ODEfind (α = -0.99965 ± 0.00030 and β = 0.09963 ± 0.00013) are within a small tolerance the expected from equation ( 8). Finally, when looking at the fitting time (Fig. 2 (b)), we find that L-ODEfind is around 50 times faster than GPoMo. 10

Example 2: Rossler

Next, we consider a more complicated case, the Rossler system, which is a non-linear (quadratic) system:

dx dt = -x -z, dy dt = x + ay, ( 9
)
dz dt = b + z(x -c).
If only variable y is observed (x and z hidden) an equation of order 3 and polynomial degree 2 can be deduced for y [START_REF] Letellier | Relation between observability and differential embeddings for nonlinear dynamics[END_REF] .

In our experiments, we set a = 0. Therefore, when considering systems where an analytical solution can be deduced, such us the oscillator (x observed, y hidden) and Rossler system (y observed, x and z hidden), both methods perform very well in terms of SMAPE prediction.

Next, we consider the Rossler system in the case of variable x observed (y and z hidden). In this scenario, there is no analytical solution using only polynomials [START_REF] Letellier | Relation between observability and differential embeddings for nonlinear dynamics[END_REF]. Interestingly, we can see that using higher order time derivatives as target, both GPoMo and L-ODEfind find a differential equation that is able to provide predictions whose accuracy is comparable (for short-term horizons) to the previous case (y observed) where there was an analytical solution (Figure 4). The fitting times continue to show that L-ODEfind is consistently faster than GPoMo. x as observed variable (y and z hidden).

Example 3: Lorenz attractor

Here we consider the Lorenz system [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF], discussed in Section 3.2, where no differential equation using only polynomials can be derived for x as the only observed variable.

Following the same methodology as before, we tried different target derivatives for both GPoMo and L-ODEfind to fit observed variable x (y and z served). We found that using higher order time derivatives (in particular second order) helps to find models that can approximate better the observed time series when integrated, although the prediction accuracy suffers from the added complexity of the problem (Fig. 5). As a consequence, the accuracy degrades faster reaching an SM AP E = 0.5 as soon as 35 time steps while in the previous cases this was attained around 185 steps. For a narrow difference again L-ODEfind outperformed GPoMo while also keeping fitting times 4 to 40 times faster than GPoMo. The time series are not measured temperatures nor the output of a simulation, but rather the result of a reanalysis process. Also, some relevant variables for modeling the atmospheric system are not available to us, for instance, wind and pressure. Our temperature time-series are not historical measurements but can be thought as a possible realization of the temperature in Paris. These temperature time series have 365 × 24 = 8760 time points.

In order to evaluate the methods, we select 39 out of the 200 temperature time series (due to GPoMo's computation time), fit the different models and use them to predict for short time horizons.

The fit was done with the first 8560 time points and the prediction was evaluated with the following time points (up to a time horizon of 15 time points).

In tackling this complex problem we want to analyze the performance of L-ODEfind and GPoMo in the forecasting task and compare their behavior to classical forecasting methods naive predictor (for every time horizon, predicts the average of the last 24 time points) and exponential smoothing (ES) (triple exponential smoothing with an additive seasonal component6 ).

In Figure 6 the SMAPE of the prediction for different time horizons in hours is displayed. L-ODEfind2 gives the lowest SMAPE for time horizons lower or equal to 6 hours, whereas GPoMo3 is better for time horizons greater than 6 hours. Notice that for time horizons greater than 6 hours, the naive predictor performs better than exponential smoothing, giving a rough idea of the reasonable predictability horizon that forecasting methods can give. In any case, both L-ODEfind and GPoMo give better forecasts than both ES and naive when the target derivative is 3. As can be seen in

Discussion

In this paper we addressed the problem of recovering differential equations from data where not all variables are observed by enhancing the approach outlined in [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF]. We tested this approach in simple and complex ODE systems and consistently found that the proposed approach of using time derivatives of higher order as target regressing variables allows to find models whose future predictions are more reliable than only using first order derivatives. We also compared our method to GPoMo and found that: (i) our proposal is orders of magnitude faster that GPoMo, and (ii) our proposal learns models with comparable or even higher prediction accuracy for several dynamical systems. Finally, we faced the challenge of addressing a real world problem and found that both L-ODEfind and GPoMo used as forecasting methods gave comparable results while at the same time outperforming classical forecasting methods, exponential smoothing and naive. In summary, L-ODEfind proved to be an accurate and fast method for recovering ordinary differential equations from data with hidden variables.

  7

Figure 1 :

 1 Figure 1: Comparison between L-ODEfind and GPoMo when finding the coefficients of data generated by the Lorenz attractor differential equations. We plot the average MSE across different random initial conditions for the system, and the time (in seconds) needed by the algorithm to fit the model.

( a )

 a Mean SMAPE versus time horizon. (b) Fitting time.

Figure 2 :

 2 Figure 2: Prediction accuracy and fitting time for L-ODEfind and GPoMo models when data comes from the oscillator system with x as observed variable (y hidden).

  52, b = 2, c = 4 and only observe y (x and z latent). L-ODEfind is consistently faster than GPoMo although, in this example, GPoMo3 and L-ODEfind3 have almost the same prediction SMAPE. Notice that the SMAPE is less than 0.2 for all models with target time derivative 2 or 3 up to 125 time horizon steps.

( a )

 a Mean SMAPE versus time horizon. (b) Fitting time.

Figure 3 :

 3 Figure 3: Prediction accuracy and fitting time for L-ODEfind and GPoMo models when data comes from the Rossler system with y as observed variable (x and z hidden).

( a )

 a Mean SMAPE versus time horizon. (b) Fitting time.

Figure 4 :

 4 Figure 4: Prediction accuracy and fitting time for L-ODEfind and GPoMo models when data comes from the Rossler system with

( a )

 a Mean SMAPE versus time horizon. (b) Fitting time.

Figure 5 :

 5 Figure 5: Prediction accuracy and fitting time for L-ODEfind and GPoMo models when data comes from the Lorenz system with x as observed variable (y and z hidden).

3. 4 .

 4 Temperature series provided by RTE In this section, we fit different models to temperatures series provided by the Réseau de Transport d' Électricité (RTE). The data consist of 200 hourly measured temperatures time series. These time series correspond to temperatures in Paris along a year for 200 different possible years or scenarios.

  figure 6 (b), the fitting times for GPoMo are much greater than L-ODEfind, whereas L-ODEfind and exponential smoothing have similar fitting times.

( a )

 a SMAPE vs time horizons in hours for different predictive methods. (b) Time to fit the different models in log scale.

Figure 6 :

 6 Figure 6: SMAPE and fitting times for different models applied on a real world data set of temperatures.

Table 1 :

 1 Models names abbreviations used in the graphics above.

We used the implementation available in the Python library sktime[Löning et al.]

Acknowledgements

Funding: This work was supported by Aristas S. R. L. and RTE (Réseau de Transport d' Électricité).


Edi Kouassi

Hilaire
email: edi.hilaire@yahoo.fr

Duquenne Philippe
email: philippe.duquenne@ensiacet.fr

Un modèle à évènements pour étudier la flexibilité opérationnelle d'un flow shop flexible

Keywords: Ordonnancement, Flexibilité opérationnelle, MINLP, Convergence Scheduling, Operational Flexibility, MINLP, Convergence

La gestion des ateliers de production recherche le meilleur compromis entre les impératifs commerciaux, financiers et industriels. L'objectif de cet article est la caractérisation et la modélisation des leviers de flexibilité d'un atelier de production à cheminement unique (« flow shop »), afin d'assurer la flexibilité opérationnelle du système. Cette flexibilité porte sur le fait que chaque poste de la chaîne de production dipose d'une ou plusieurs machines identiques, sur des durées opératoires variables à chaque poste selon le nombre des opérateurs affectés, et sur les dates de livraison des produits. La convergence entre les ressources requises pour l'exécution d'un plan de production et celles mises en oeuvre pour sa réalisation est favorisée. La fonction objectif intègre l'évaluation des stocks et en-cours, le coût du travail et de l'inactivité, et les conséquences des perturbations. Le problème est résolu en traitant l'ensemble des variables simultanément, selon une formulation de programmation mathématique non-linéaire en variables mixtes (MINLP).

Introduction

En vue de relever le défi de la performance et de la compétitivité dans un environnement économique en perpétuel changement, les entreprises sont de plus en plus amenées à s'adapter, et accroître leur réactivité vis-à-vis du marché afin d'assurer leur pérennité. Un aspect de cette réactivité concerne la convergence entre la charge de travail et la capacité en ressources humaines d'une entreprise à repondre à une demande d'un plan de production. Cette capacité de s'adapter, reconnue sous le nom de flexibilité [START_REF] Edi | Affectation flexible des ressources dans la planification des activités industrielles : prise en compte de la modulation d'horaires et de la polyvalence[END_REF][START_REF] El-Awady | Flexible resources allocation techniques : characteristics and modelling[END_REF], est la traduction d'un ensemble de décisions de gestion prises aux niveaux stratégique, tactique et opérationnel. Ces décisions permettent d'obtenir une organisation efficace du système de production, d'assurer le respect des méthodes de fabrication et des contraintes sur l'utilisation des ressources. Elle portent également sur les moyens de production ; les moyens de gestion de la production ; les moyens d'achats et ventes ; du système d'information et de conception du produit. De ce point de vu, les chercheurs ont développé et proposé des travaux d'ordonnancement en vue de maximiser la production tout en minimisant les coûts de production. En plus des coûts liés à l'utilisation des moyennes de production (salaires, énergies et coûts d'organisation), les coûts de stockages et les pénalités de retards prennent parfois des proportions démesurées. La recherche d'un ordonnancement capable d'obtenir un niveau de stock zéro ou minimum et la production de biens dans un délai acceptable s'avèrent décisive pour réduire le coût de production et augmenter le gain de l'entreprise. Après une présentation de l'état des travaux d'ordonnacement d'un flowshop, nous exposerons notre modèle et enfin nous présenterons une resolution à parir d'une methode exacte : procedure par séparation et évaluation (« Branch and bound »).

Méthodologie

2-1. Situation du problème : ordonnancement d'un flow shop

L'ordonnancement est un processus très important utilisé dans la gestion et la production des biens et services. Il consiste en une allocation des ressources dans le but d'exécuter un ensemble de travaux sur une période de temps donnée. Depuis les travaux de Salveson [START_REF] Salveson | On a quantitative method in production planning and heduling[END_REF] et de Johnson [START_REF] Johnson | Optimal two-and three-stage production schedules with setup times included[END_REF], un grand nombre d'études ont été publiées sur l'ordonnancement. Dans un problème d'ordonnancement de flow shop simple, chaque poste de travail est constitué d'une seule machine [START_REF] Campbell | A heuristic algorithm for the n job, m machine sequencing problem[END_REF][START_REF] Nawaz | A heuristic algorithm for the mmachine, n-job flowshop sequencing problem[END_REF]. Si un poste au moins comporte plus d'une machine, on parle alors de flow shop flexible ou flowshop hybride [START_REF] Logendran | Group scheduling in flexible flow shops[END_REF][START_REF] Kyparisis | Flexible flow shop scheduling with uniform parallel machines[END_REF][START_REF] Linn | Hybrid flow shop scheduling : A survey[END_REF]. Ainsi, un flow shop flexible n'est qu'une généralisation du concept de flow shop simple [START_REF] Chung | Scheduling flexible flow shops with no setup effect[END_REF]. Un flow shop flexible est une configuration d'ateliers à cheminement unique dans lequel les travaux sont exécutés de façon linéaire à travers des postes de travail constitués de machines en parallèle [START_REF] Shaukat | Scheduling in a Flow Shop with Multiple Processors[END_REF]. Des travaux effectués dans [START_REF] Graham | Optimization and approximation in deterministic sequencing and scheduling theory: a survey[END_REF] ont permis d'élaborer une classification standard des problèmes de flow shop (notation standard α|β|γ). Cette classification a été étendue par Vignier [START_REF] Vignier | Les problèmes d'ordonnancement de type flow shop hybride : état de l'art[END_REF] dans le but de prendre en compte, les nouvelles architectures d'organisation des machines dans les ateliers (machines parallèles …). En réalité, les unités de production utilisent rarement une seule machine à chaque poste pour exécuter une opération. La chaîne de production est constituée d'une série de postes de travail constituée de machines identiques en parallèle. Cette duplication des machines à chaque poste a pour but d'accroître la flexibilité et l'efficacité de la chaîne de production. De plus, lorsque la durée opératoire d'un travail donné est la même sur toutes les machine (resp. tous les postes de travail), on parle alors de flow shop proportionnel (resp. de flowshop flexible proportionnel). Des travaux portent sur une configuration de flowshop flexible dans laquelle, l'exécution d'une opération est précédée (resp. suivie) d'un temps de montage (resp. démontage) dépendant de la séquence des travaux sur la machine. Ce temps mobilise uniquement la ressource [START_REF] Kurz | Scheduling flexible flow lines with sequence dependent set-up times[END_REF]. Il existe également des études qui portent sur l'ordonnancement avec des contraintes de blocage [START_REF] Wajdi | Ordonnancement des systemes de production flexibles soumis a differents types de contraines de blocage[END_REF]. Le concept de blocage vient du fait que dans les lignes flexibles de production, les postes de travail sont reliés par un convoyeur. Ainsi, un poste de travail ne peut-il libérer un produit que si le convoyeur d'accumulation aval a une position libre. A l'inverse, un poste de travail en aval peut être en attente par manque de produit provenant du poste amont.

Dans ce travail, nous supposerons une variation de la durée opératoire sur chaque poste. Cette variation est obtenue grâce à des modifications du nombre des opérateurs (ressources humaines disponibles).

2-2. Complexité

L'ordonnancement de travaux sur une chaîne de production de type flow shop flexible est considéré comme NP-difficile [START_REF] Sriskandarajah | Scheduling algorithms for flexible flow shops : worst and average case performance[END_REF]. Même pour le cas de flow shop simple avec une seule machine à chaque poste de travail (2, 3 postes de travail ou davantage). Le problème est toujours considéré comme NP-difficile au sens fort du terme comme le montre Garey [START_REF] Garey | Two-processor scheduling with start times and deadlines[END_REF]. De plus, selon Gupta, tout flowshop de deux postes de travail disposant plus d'une machine à un poste est NP-difficile au sens fort du terme [START_REF] Gupta | Schedules for the two-stage hybrid flowshop with parallel machines at the second stage[END_REF].

2-3. Méthodes de résolution

Plusieurs méthodes de résolution ont été élaborées pour aborder les différents problèmes de flow shop flexible. On distingue les méthodes exactes, les heuristiques et les métaheuristiques.

Parmi les méthode exactes, la méthode de séparation et évaluation est la plus utilisée. Brah et Hunsucker [START_REF] Brah | Branch and bound algorithm for the flow shop with multiple processors[END_REF] ont-ils présenté un algorithme basé sur la méthode de séparation et évaluation dans une configuration de flowshop hybride avec un nombre arbitraire de poste de travail et la prise en compte de stocks intermédiaires . Portmann et al. [START_REF] Portmann | Branch and bound crossed with GA to solve hybrid flowshops[END_REF] ont proposé une version améliorée de l'algorithme de séparation et évaluation de Brah et Hunsucker. Dans leur travail, ils ont amélioré leur borne inferieure et ont réduit le nombre de branches utilisé dans l'arbre de recherche. Moursli et Pochet [START_REF] Moursli | Q branch and bound algorithm for the hybrid flowshop[END_REF] ont également présenté un algorithme de séparation et évaluation pour la résolution d'un problème de flow shop hybride dans lequel chaque poste de travail est constitué de plusieurs machines identiques en parallèle. En outre, Sawik quant à lui présente une formulation par programmation en nombre entier pour la résolution d'un problème d'ordonnancement d'une ligne de production flexible avec une capacité de stockage limitée [START_REF] Sawik | Mixed integer programming for scheduling surface mount technology lines[END_REF].

Bien que les méthodes dites exactes fournissent des solutions optimales, elles se sont avérées de moins en mois efficaces pour la résolution des problèmes de grandes tailles et nécessitant un temps de calcul dépassant l'attente des responsables de la planification. Cette situation a conduit certains chercheurs à élaborer des méthodes dites approchées. Il existe deux types de méthode approchée dont les heuristiques et les méta-heuristiques. Les chercheurs [START_REF] Chen | A tabu search algorithm for the integrated scheduling problem of container handling systems in a maritime terminal[END_REF][START_REF] Wardono | A tabu search algorithm for the multistage parallel machines problem with limited buffer capacities[END_REF] ont proposé une méthode tabou pour aborder le même problème. Une étude réaliste a été effectuée par Low [START_REF] Low | Simulated annealing heuristic for flow shop scheduling problems with unrelated parallel machines[END_REF] sur le model du flow shop flexible. Dans cette étude, il utilise le recuit simulé pour l'optimisation du temps d'écoulement total dans le système. Houari et Mhallah ont étudié un problème de flow shop hybride à deux postes de travail dans lequel ils utilisent une procédure de recherche tabou et une procédure de recuit simulé pour sa résolution [START_REF] Haouri | Heuristic algorithms for the two-stage hybrid flowshop problem[END_REF]. Ishibuchi et Misaki [START_REF] Ishibuchi | Modified simulated annealing algorithms for the Flow-Shop sequencing problem[END_REF] se sont intéressés à la méthode de recuit simulé pour résoudre ces problèmes. D'autres auteurs se sont focaliser sur les algorithmes génétiques [START_REF] Imran | Minimizing makespan for a no-wait flowshop using genetic algorithm[END_REF][START_REF] Torabi | A hybrid genetic algorithm for the finite horizon economic lot and delivery scheduling in supply chains[END_REF] pour optimiser la durée d'exécution des travaux sur les postes de travail.

Dans la plupart des études menées, la durée de traitement d'un article sur chaque poste de travail ou machine (durée opératoire ou cycle) est fixe et connue. Aussi, elles ne font pas allusion de l'utilisation des effectifs des ressources humaines comme un aspect de la flexibilité d'un flow shop. De plus, nous considérons que la date de livraison d'un article est une variable qui est défini dans un intervalle prédeterminé.

2-4. Modélisation du problème 2-4-1. les variables et les données

Les indices i : indice pour l'article et m le nombre total d'articles sur la période, 1 ≤ i ≤ m j : indice pour le poste et n le nombre total de postes, 1 ≤ j ≤ n k,o : indices pour les événements et T le nombre total d'événements, 1

≤ k ≤ T et 1 ≤ o ≤ T avec m n T m n + ≤ ≤ 2. .
Les variables continues : t k : date d'occurrence de l'événement k, L i 0 : date de livraison théorique de l'article i.

Les variables discrètes sont les suivantes : N j,k : nombre de positions disponibles sur le poste j entre les événements k et k+1, as i,j,k : variable binaire désignant l'affectation du début d'opération à l'événement k (as i,j,k = 1 si l'article i arrive sur le poste j à la date de l'événement k et = 0 sinon), ae i,j,k : variable binaire désignant l'affectation de fin d'opération à l'événement k (ae i,j,k = 1 si l'article i quitte le poste j à la date de l'événement k et = 0 sinon), ω i,j : variable binaire matérialisant un changement de durée opératoire sur le poste j entre les articles i et i+1, E k : effectif présent dans l'atelier entre les événements k et k+1, E tot : effectif maximum présent dans l'atelier sur l'horizon de temps E tot = max(E k ), E i,j : effectif présent sur le poste j pour l'article i, y k : variable d'écart pour un événement k possible et non affecté, c i,j : cycle de l'article i sur le poste j, δ i,j : date de disponibilité des éléments de l'article i pour le poste j.

Les données : c j min : valeur minimale du cycle pour le poste j, c j max : valeur maximale du cycle pour le poste j, γ : nombre d'articles sur lesquels se répercutent les perturbations, N j min : nombre minimal de positions disponibles sur le poste j (en principe égal à 0), N j max : nombre maximal de positions disponibles sur le poste j, N j 0 : nombre de positions initialement disponibles sur le poste j (par défaut, N j 0 =N j max) , L i -: date de livraison au plus tôt pour l'article i (par défaut, L i -= L i 0 -2), L i + : date de livraison au plus tard pour l'article i (par défaut, L i + = L i 0 + 2), χ j : charge de travail sur le poste j, v i,j : valeur de l'article i à l'entrée du poste j (en coût matière uniquement), w i,j : valeur de la matière ajoutée à l'article i sur le poste j. il établit la relation suivante, pour assurer la cohérence des données, entre v i,j et w i,j :

j i j i j i w v v j n i , , , 1 , 1 + =

2-4-4. Les contraintes du problème

Dans ce paragraphe nous exprimons les contraintes du système de production par un ensemble d'équations comme présentées ci-dessous : -Chaque événement k (fin ou début d'opération pour un article i sur un poste j) doit se produire une et une seule fois sur l'ensemble de l'horizon.

1 , 1 , 1 , , 1 , , = ∀ = ∀ ∑ ∑ = = T k k j i T k k j i ae j i et as j i (1)
-Relation de précédence pour un article i entre deux postes (j et j+1) de la gamme : Pour qu'une opération j débute au k ème événement, l'événement déclenchant la fin de l'opération j-1 doit être affectée à un événement antérieur (donc compris entre 1 et k).

∑ - = - ≤ > ∀ > ∀ ∀ 1 1 1 1 , 1 , k o ,o i,j i,j,k ae as k j i ( 2
)
-Cohérence des événements de début et de fin de chaque opération : pour qu'une opération se termine au k ème événement, elle doit avoir débuté par une affectation à un événement antérieur (de 1 à k-1).

∑ - = ≤ > ∀ > ∀ ∀ 1 1 1 , 1 , k o i,j,o i,j,k as ae k j i (3)
-Respect de la date de livraison échue : la mise à disposition de l'article peut avoir lieu dès la fin de la dernière opération, cette date étant définie par le membre de gauche de l'inégalité ci-dessus.

+ - = ≤ ≤ ∀ ≤ ∀ ∑ i i i i T k k k n i L L L i et L t ae i . 1 , , (4)
-Affectation cohérente des postes : le poste dont le nombre de positions occupées initialement est connu, fait l'objet d'un bilan des actions qui débutent et qui s'achèvent lors de l'occurrence de chaque événement. La variable désignant le nombre de positions disponibles est contrainte par le nombre maximum de positions utilisables. max , min

j k j j m i k j i m i k j i k j k j N N N avec ae as N N k j ≤ ≤ + - = > ∀ ∀ ∑ ∑ = = - 1 , , 1 , , 1 , , 1 ,
-Définition de la durée opératoire et des dates de réalisation : la durée de chaque opération est calculée comme la différence entre les dates de fin et de début d'opération, exprimée en fonction des dates des événements et des variables binaires d'affectation. Chaque durée est contrainte entre deux bornes. De plus, les événements sont classés par ordre chronologique

k k j j i j j i T k k k j i T k k k j i t t k et c c c avec c t as t ae j i ≤ ≥ ∀ ≤ ≤ = - ∀ ∀ - = = ∑ ∑ 1 max , min , 1 , , 1 , , 1 . . , (6)
-Relation entre l'effectif et le cycle sur chaque poste et pour chaque article : l'effectif requis pour exécuter une opération sur le poste est égale au rapport de la charge de travail par le cycle.

∀ ∀ - = i j E c i j i j j , . , , χ 0 (7)
-Détermination de l'effectif dans l'atelier : l'effectif total de l'atelier à chaque état (entre deux événements consécutifs) est calculé comme la somme des effectifs sur les postes entre les événements qui délimitent l'état. La variable E tot interviendra dans la fonction coût, afin de chiffrer le coût total du travail.

tot k m i n j i,j i,j,k m i n j i,j i,j,k k- k E E k .E ae .E as E E k et E ≤ ≤ ∀ - + = ∀ = ∑∑ ∑∑ = = = = 0 0 1 1 1 1 1 0 (8)
-Assignation effective des événements possibles et classement par ordre chronologique : la variable d'écart y k est forcée à prendre la valeur 1 si aucune action n'est observée lors de l'événement k ; les variables d'écart sont classées par valeurs décroissantes afin de regrouper tous les événements non affectés (pour lesquels la variable d'écart est égale à 1) en début de liste. Si un événement est associé à une action, la variable d'écart prend une valeur nulle.

y k ≤ > ∀            ≤ + ≤ + ≥ + + ∀ ≤ ≤ ∀ - = = = = = = = = ∑∑ ∑∑ ∑∑ ∑∑ 1 1 1 , , 1 1 , , 1 1 , , 1 1 , , 1 1 1 1
et 1 0 (9)
-Ecart minimal entre deux événements affectés : afin d'éviter que deux événements distincts ne soient associés à deux dates égales, nous définissons un écart minimal ε entre deux dates consécutives (la valeur de ε doit être petite par rapport aux cycles de production, nous avons retenu une valeur de 0,5 unités de temps). Si les événements ne sont pas affectés, ils peuvent être associés à des dates égales.

( )

1 1 . 1 + ≤ - + > ∀ k k k t y t k ε (10)
-Respect des dates de disponibilité : une opération ne peut commencer que lorsque les éléments matériels nécessaires (fournitures, approvisionnements) à sa réalisation sont disponibles sur le site de l'atelier. De plus, l'hypothèse d'une alimentation à débit constant permet de calculer l'ensemble des dates d'alimentation pour tous les articles et pour tous les postes à partir de la première date (δ 1,1 ), en utilisant la durée entre deux approvisionnements dans l'atelier, noté α, et la durée entre l'approvisionnement de deux postes consécutifs, noté β. La date δ 1,1 intervient dans le terme de la fonction objectif ; elle prendra donc la valeur maximale permettant de respecter la contrainte, afin de minimiser le coût total.

( ) ( )

1 . 1 . . , 1 , 1 , , 1 , , - + - + = ≥ ∀ ∀ ∑ = j i avec t as j i j i j i T k k k j i β α δ δ δ (11)

2-4-5. La fonction objectif à minimiser

La fonction objectif à optimiser est la somme de plusieurs critères, exprimés en unités monétaires ; à savoir : -F 1 : coûts liés au stockage amont, -F 2 : coûts liés aux en-cours et stocks intermédiaires, -F 3 : coûts liés au stockage de produits finis, -F 4 : pénalités de retard, -F 5 : coût du travail effectif, -F 6 : coût de l'inactivité, -F 7 : coût des perturbations. -Les stocks de matière avant chaque poste (F1) : les coûts liés à l'immobilisation de la matière première et des articles stockés avant leur entrée dans l'atelier sont calculés en fonction de la valeur de chacun des éléments, sur la durée comprise entre leur date de disponibilité et leur date d'entrée sur le poste correspondant.

(

) ( )

( ) ( ) ( ) ( ) ( ) ∑ ∑ ∑ ∑ ∑ ∑ = = = = = =               - + - +       - + - + - =         - + - = m i T k k k i n j T k k k j i j i m i i i n j j i j i j i i t as v j i t as w d v d w F i i 1 1 , 1 1 , 1 , 1 1 , 1 1 , , , 1 1 , 1 , 1 , , , 1 1 .
. . 1 . 1 . . . . . . . 0 , 0 , α δ β α δ τ δ δ τ
avec d i,j : date de début de l'opération j pour l'article i, (12) v i,0 : valeur de l'article i avant le début de l'opération 1.

-Les articles en cours de production et les stocks intermédiaires (F2) : les encours sont calculés, article par article, entre la date d'entrée de l'article dans l'atelier (d 1,1 ) et la date de mise à disposition de celui-ci (f i,n ), en prenant en compte l'évolution de sa valeur. Soit f i,j , la date de fin de l'opération j pour le poste i, on a :

( ) ( ) ( ) ( ) ( ) F f d v d f v d d v f d v v as as t i j i j i j j n i m i j i j i j j n i m i j i j i j j n i n i n i n i m i j i j k i j k k k 2 1 1 1 1 1 1 1 1 1 1 1 = - + -       = - + -               = - = = + = - = + = - = + = ∑ ∑ ∑ ∑ ∑ ∑ τ τ τ . . . . . . . . , , , , , , , , , , , , , , , , , ( ) 1 1 1
1 1 1 T j n i n i n k i n k k k T i m i m v ae as t ∑ ∑ ∑ ∑ ∑ = - = = = + -               , , , , , . . (13)
-Stockage des produits finis avant la livraison (F3) : le stockage des produits finis s'étale sur la période comprise entre la date de mise à disposition de l'article et sa date de livraison (L i ), au cours de laquelle la valeur n'évolue plus, et le risque d'obsolescence sera considéré comme négligeable.

( ) ∑ ∑ ∑ = = =       - = - = m i n i T k k k n i i m i n i n i i v t ae L v f L F 1 , 1 , , 1 , , 3 . . . . . τ τ (14)
-Pénalités de retard (F4) : la fonction P i (t) est la loi de pénalité applicable pour un retard de livraison de l'article i ; elle peut varier selon le client. Elle est a priori croissante au sens strict.

( )

F t dt i f L P i i 4 0 = - ∫ (15)
-Coût du travail (F5) et de l'inactivité (F6) : les composantes permettant de chiffrer le coût du travail réellement effectué F 5 et l'inactivité F 6 , et peuvent être évaluées ensemble ou indépendamment l'une de l'autre. En effet, la somme de ces deux termes correspond à la rémunération du nombre maximum d'opérateurs entre le début et la fin du plan de production ; on peut ainsi écrire, en notant U le coût du travail par unité de temps :

( ) ( ) ∑ = - = - = + T k k k n m k tot tot n m as ae t U E U E d f F F 1 , 1 , 1 , , 1 , 1 , 6 5 .
. . . . (16)
Si l'on souhaite dissocier activité et inactivité, par exemple en leur associant des coûts différents, notés respectivement U et S (si l'on considère que l'inactivité n'est pas totalement perdue mais peut être utilisée à d'autres activités, telles que formation, action d'amélioration, etc.), on obtient alors :

( )

F U E t t k k k k T 5 1 1 1 = - + = - ∑ . . (17) ( ) ( )
F S E E t t tot k k k k T 6 1 1 1 = - - + = - ∑ . . (18)
-Coût des perturbations (F7) : chaque changement de cycle (durée opératoire) sur un poste entraîne une modification de la composition de l'équipe qui lui est affectée ; des perturbations sont ainsi générées selon le sens et l'amplitude de la variation. La diminution d'un cycle n'a pas la même influence sur le coût qu'une augmentation de même amplitude.

Pour un poste j donné, le coût d'un changement de cycle entre deux articles consécutifs aura la forme suivante (où γ est un coefficient pondérateur caractérisant l'importance donnée aux perturbations et f désigne le coût « local » de la perturbation entre les articles i et i+1) : ( ) La modélisation aboutissant à un problème de programmation non-linéaire en variables mixtes, sa résolution nécessite l'utilisation d'un code capable de traiter l'arborescence des solutions, par une méthode de séparation et évaluation. Le problème est tout d'abord traité comme un problème de programmation non-linéaire en variables réelles, qui permet d'obtenir la solution relaxée du problème ; la résolution du problème de programmation non-linéaire à chaque noeud de la procédure de « Branch & Bound » permet de passer progressivement du problème relaxé au problème en variables mixtes. La totalité des calculs ont été effectués avec OPTMIX (code en Fortran utilisable sur PC [START_REF] Reneaume | Formulation et résolution du problème d'optimisation nonlinéaire en variables mixtes dans un environnement modulaire : Application à la synthèse optimale des procédés[END_REF]. Compte tenu de l'organisation et de l'indexation des variables, l'algorithme va procéder comme suit :

si c i+1,j > c i,j , alors f U c c j i j i j 7² 1 1 = -         + . . .
F U c c c c j i j i j i j i j i j i j j n i m 7 1 1 1 1 1 1 1 = - --         + + = = - ∑ ∑ . . . . . , , , , , , γ χ ω ω (19)
-Résolution du problème relaxé (en considérant que toutes les variables sont continues) ; -Séparation et évaluation des sommets, en commençant par les variables de types ae i,k , puis as i,j,k , puis c i,j et enfin ω i,j ;

-Obtention d'une solution admissible, retenue comme optimale ou non, selon la valeur de la fonction objectif en ce point ; -Exploration de l'arborescence tant qu'il reste des sommets non développés et dont l'évaluation est meilleure que la solution optimale courante.

La détermination de l'ordre de séparation et d'évaluation des sommets répond à un impératif de réduction du temps de calcul ; en effet, le chemin total parcouru dans l'arborescence pour atteindre un sommet peut connaître d'importantes différences selon l'ordre retenu.

Dans le cas présent, la logique retenu répond aux considérations suivantes :

-Détermination de l'enchaînement des livraisons et de leur affectation aux événements ; -Détermination de l'ordre des débuts de poste, en fonction des livraisons prévues pour la branche de l'arborescence en cours d'exploration ; -Choix des cycles de production, compatibles avec les décisions d'affectation préalables ;

-Détermination des valeurs des variables ω i,j d'après les valeurs des c i,j ;

-Détermination de la valeur des variables continues. Il est cependant important de noter les remarques suivantes :

-Les perturbations engendrées par les variations de durées opératoires (inactivité F 6 , et coûts de changements de cycles, F 7 ) ne sont compensées par l'ensemble des gains sur les coûts d'immobilisation que dans le cas de produits dont la valeur est élevée par rapport aux coûts salariaux.

-De même, le taux d'actualisation des encours τ joue un rôle majeur ; si sa valeur est élevée, la réorganisation de la production est aisément justifiable par les gains potentiels sur les frais financiers, alors qu'une valeur faible de τ diminuera notablement l'intérêt de telles modifications.

-Enfin, la répartition temporelle des dates de livraison est primordiale ; si celles-ci sont réparties presque de manière homogène, il est préférable de stocker les produits finis durant de courtes périodes.

Conclusion

La description du temps selon un principe de discrétisation non-uniforme permet la modélisation du problème en prenant en compte simultanément l'ensemble des variables de décision. De plus, elle possède un caractère pédagogique sur la présentation du problème et les relations existant entre les variables. La résolution du problème de programmation non linéaire en variables mixtes (MINLP) conduit à la solution optimale ; celle-ci a ainsi pu être comparée à la solution à débit constant, qui ne fait appel à aucun des potentiels de flexibilité. Au regard des critères complémentaires (taux d'utilisation des ressources et taux de variation des effectifs), la solution optimale présente également des caractéristiques satisfaisantes. Cependant, la taille du problème exprimé en variables mixtes sur un exemple de dimensions très modestes, ainsi que le temps nécessaire à la résolution, constituent un obstacle sérieux quant à son utilisation pour un exemple de dimensions industrielles ; il est donc nécessaire d'étudier une autre méthode qui pourra être l'objet d'une autre publication.

2 .Figure 1 :

 21 Figure 1 : Effectifs et durée opératoire admissible

  Notre objectif, est l'obtention d'une solution optimale, non contrainte par des limites sociales, et de mesurer l'intensité des perturbations. En désignant par τ le taux d'actualisation par unité de temps, les coûts imputés à l'immobilisation d'un capital V durant une période d sont de la forme : V.[(1+τ) d -1]. Lorsque cette durée d est courte, nous pouvons nous livrer à l'approximation suivante : V.[(1+τ) d -1] ≈ V.τ'.d et Lorsque τ est petit, l'approximation τ = τ' est acceptable.

,

  et si c i,j > c i+1,j , alors f

2- 5 . Méthodologie de résolution du problème 2 - 5 - 1 .

 5251 Choix de la méthode et stratégie de resolution

3 .

 3 Resultats et DiscussionLa résolution du problème défini ci-dessus à l'aide de la formulation MINLP conduit à l'obtention de la solution optimale représenté sur le diagramme Figure5, obtenue en 183.10 3 secondes CPU, soit environ 50 heures (sur un PC équipé d'un microprocesseur P400 MHz) ; l'axe des abscisses représente l'échelle des temps, les deux postes considérés sont représentés en ordonnée : on obtient donc le diagramme de Gantt des différents postes.

Figure 1 :

 1 Figure 1 : Diagramme de Gantt des postes de la solution optimale

Tableau 5 :

 5 Solution optimale (Cycles et dates de livraison)La figure 6 constitue le plan de charge des opérateurs par poste associé à la solution optimale ; il représente le nombre d'opérateurs nécessaire sur chaque poste par unité de temps.

Figure 6 :

 6 Figure 6 : Plan de charge de la solution optimale

Tableau 6 :

 6 Décomposition de la fonction coût pour la solution optimale

2-5-2. Données de simulation du problème

  La formulation décrite est appliquée à la résolution d'un exemple ; l'atelier étudié est constitué de deux postes de travail, P1 et P2, sur lesquels doivent être traités dix articles, suivant une gamme identique. Chaque article doit d'abord passer sur le poste P1 avant de passer sur le poste P2 ; de plus, chaque poste ne compte qu'une seule position. L'ensemble des pièces devant être intégrées dans l'atelier arrivent dans le stock amont selon un principe « à débit constant» ; le rythme d'arrivée (intervalle entre deux arrivées) est connu a priori, mais la date d'arrivée du premier lot dans le stock (δ 1,1 ) est une variable, qui détermine l'arrivée de tous les autres lots. Les valeurs discrètes autorisées pour les cycles des postes P1 et P2 sont au nombre de deux, et sont présentées dans le tableau 1 ci-dessous, ainsi que les charges de travail sur chaque poste (en hommes.unités de temps) : Cycles et charge de travail par poste De plus, l'effectif total disponible sur l'atelier ne peut en aucun cas dépasser 15 opérateurs, qui est donc une borne supérieure pour E tot . Les intervalles autorisés pour la livraison de chacun des articles sont énumérés dans le tableau 2 ci-dessous. Date de livraison théorique et intervalles de livraison pour chaque article

	Cycle minimum (en u. t.) min j c Cycle maximum (en u. t.) max j c Charge de travail (en h. u. t.) 1 2 3 4 5 6 7 8 9 10 Tableau 1 : N° de l'article Tableau 2 :	L i 26 -35 38 55 65 73 78 87 93 101	Poste P1 6 8 48 L i 0 28 37 40 57 67 75 80 89 95 103	Poste P2 5 8 60 L i + 30 39 42 59 69 77 82 91 97 105


Segregation of the terminal region of the Escherichia coli chromosome

Cohesion period

Replication Segregation

Mathieu Stouf 1 , Jean-Christophe Meile 1 , Marie Deghorain 2 , Carine Pages 1 and François Cornet 1 . 1 Laboratoire de Microbiologie et de Génétique Moléculaire, CNRS, Université de Toulouse, France 2 Université Catholique de Louvain, Louvain-la-Neuve, Belgium

The FtsK DNA-translocase controls the last steps of chromosome segregation in E. coli. It is associated with the cell division septum and translocates sister chromosomes using the KOPS DNA motifs to orient its activity, thereby controlling both decatenation by TopoIV and chromosome dimer resolution by XerCD/dif recombination. Using XerCD/dif recombination to probe the interaction of chromosome loci with FtsK, we have shown that FtsK activity is restricted to a ~400 kb region around the dif position, the FtsK High Activity Region (FHAR). Measurement of the relative segregation times of different loci tagged with fluorescent proteins by microscopy shows that this region is subjected to an elevated level of post-replicative cohesion between sister-chromosomes. Using pairs of loci tagged with two different fluorescent proteins, we showed that loci of the FHAR segregate in an ordered manner, the dif site being the last point of contact between sister chromosomes. These data suggest that the oriented translocase activity of FtsK controls the ordered segregation of a specific terminal region of the chromosome.

A model for the segregation of the terminal region A model for the segregation of the terminal region --> Releasing of cohesion in the FHAR occurs from the outside towards the --> Releasing of cohesion in the FHAR occurs from the outside towards the dif dif site site --> The segregation of FHAR is ordered and oriented --> The segregation of FHAR is ordered and oriented --> The --> The dif dif site appears to be the last point of contact between the two sister chromosomes. site appears to be the last point of contact between the two sister chromosomes.

--> What is the last contact point between the two sister chromosomes ? --> What is the last contact point between the two sister chromosomes ? --> How is cohesion released during segregation of the --> How is cohesion released during segregation of the ter ter region ? region ?

--> FtsK activity and cohesion time of

I I N N T T R R O O D D U U C C T T I I O O N N R R E E S S U U L L T T S S C C O O N N C C L L U U S S I I O O N N
We used a recombination system with a cassette contaning lacI gene flanked by two dif sites. When FtsK translocates this sequence, it activates XerCD recombination between the two dif site. The loss of the lacI gene is measured as the apparition of dark blue colonies on Xgal-containing plates.

At the end of replication, the sister ter region stay together.

During the past few years, it has been shown that, between replication and segregation, there is a period where the two replicated loci stay together. This period of colocalisation called "Cohesion" is measured by fluorescence microscopy as the proportion of cells with a single focus of a given tagged locus. Cohesion was shown to be most important in the ori and ter region, suggesting these two regions host specific activities for cohesion and segregation.

  The mecanism of segregation leads to the repartition of the sister chromosomes in the two daughter cells. This process occurs progressively from the ori to the ter region and concomitantly with DNA replication. Chromosome dimers and FtsK Chromosome dimers and FtsK Segregation of the Echerichia coli chromosome Segregation of the Echerichia coli chromosome FtsK acts in a restricted region and maintains an high level of cohesion FtsK acts in a restricted region and maintains an high level of cohesion Replication of the Echerichia coli chromosome proceeds bidirectionally along the two chromosome arms from the oriC site to the ter region, where it is stopped by terminator site. During replication, recombinaision can occur between sister chromosomes leading to the formation of chromosome dimers. These structures are resolved by the activation of the tyrosine recombinase XerCD at the dif site by the DNA translocase FtsK. FtsK is associated to the division septum by its N-Ter domain, which is essential for cell division, and translocates DNA via its C-Ter domain which plays a central role in segregation of the ter region. FHAR segregation is oriented and ordered FHAR segregation is oriented and ordered --> Four loci were chosen in the FHAR --> Two loci visualized in a same cell using systems ParB/parS P1 and pMT1

  FHAR are not controlled by terminaison of replication and are rather --> FtsK activity and cohesion time of FHAR are not controlled by terminaison of replication and are rather linked to postreplicative processing of the terminal region linked to postreplicative processing of the terminal region XerCD/dif recombination system The FHAR : FtsK High Activity Region The FHAR is not where replication terminates --> Where does FtsK acts? --> Where does FtsK acts? --> Is there a relation between FtsK activity and the terminaison of the replication ? --> Is there a relation between FtsK activity and the terminaison of the replication ? FHAR cohesion level depends on the FtsK translocase activity

  dif dif region is the region with the highest level of region is the region with the highest level of cohesion of all the FHAR cohesion of all the FHAR -FHAR = High recombination frequence --> The FHAR shows higher level of cohesion than its adjacent regions.--> Whithout translocase activity, cohesion is lower in the FHAR Loci away from dif are segregated in first.The level of cohesion between loci closest from the dif site is the highest of all the FHAR. Cohesion is released progressively on the right arm Here, we have shown that FtsK activity is restricted to a region of about 400 kb in the ter high cohesion region. Whe also showed that in this region, cohesion is released progressively from the outside to the dif site which appears to be the last point of contact between the two sister chromosomes. Recently, an important actor for the organisation and condensation of the ter region was discovered. The MatP/matS organizes the ter region into a Macro Domain. The order of segregation and the time of cohesion of the FHAR is being studied in the lab in different mutant contexts : p-topoIV, FtsK ATP-, FtsK KOPS blind , matP-, mukB-and XerCD-in order to understand and characterize the mecanism of segregation of the ter region and precise the last steps of the cell cycle.When the septum is closing, the two dif sites are segregated by the final releasing of the FHAR cohesion.