Danilo Dursoniah
email: danilo.dursoniah@univ-lille.fr

Maxime Folschette
email: maxime.folschette@centralelille.fr

Cédric Lhoussaine
email: cedric.lhoussaine@univ-lille.fr

Franc ¸ois Pattou Rebecca Goutchtat
email: rebecca.goutchtat@vet-alfort.fr

Violeta Raverdy

Limits of a Glucose-Insulin Model to Investigate Intestinal Absorption in Type 2 Diabetes Work in progress

Keywords: Ordinary Differential Equation, Systems Biology, Type 2 Diabetes, Intestinal Glucose Absorption

Abnormal regulation of glucose absorption in the small intestine is an important cause of Type 2 Diabetes (T2D). Even if this hypothesis is clinically well-known, it has not been fundamentally validated yet, mainly due to a lack of reliable metabolic knowledge on the glucose regulation. The main objective of this paper is to test this hypothesis on a highly referenced model composed of ordinary differential equations. This model is tested on an original dataset featuring the observations of obese diabetic patients. It shows its limits to predict our post-prandial glycemia and insulinemia time series especially with regard to the crucial complexity of gastro-intestinal regulation.

I. INTRODUCTION

Diabetes is a chronic metabolic disease characterized by a lack of insulin secretion and a decreased peripheral insulin response. Insulin is a hormone that down-regulates blood sugar concentration. Consequently, the balance of glucose and insulin concentrations in different tissues, called homeostasis, is pathologically perturbed: hyperglycemia is observed both during fasting and post-prandial periods. It gradually leads to tissues damages and subsequent diseases, i.e., a high rate of comorbidity [START_REF] Ghasemi | Type 2 diabetes: An updated overview[END_REF].

More precisely, Type 2 Diabetes (T2D) results from the body's ineffective use of insulin. Most patients (∼ 90%) with diabetes have T2D. Around 400 millions of people are affected worldwide by the disease representing a major public health issue in most developed countries [START_REF] Goedeke | Emerging pharmacological targets for the treatment of nonalcoholic fatty liver disease, insulin resistance, and type 2 diabetes[END_REF]. It is commonly accepted that this type of diabetes is largely caused by physical inactivity combined with an high-carbohydrate diet. However, through bariatric surgery, obese patients with T2D have seen their physiological glycemia immediately restored, independently to their weight loss [START_REF] Baud | Sodium glucose transport modulation in type 2 diabetes and gastric bypass surgery[END_REF]. This observation leads us to consider Intestinal Glucose Absorption (IGA) as a critical cause of T2D, among others. Bariatric surgery, and more precisely Roux-en-Y Bariatric Surgery (RYGB), anatomically leads to the decrease of the glucose absorption surface, which would explain, at least partially, this unexpected clinical benefit. Furthermore, the gastro-intestinal tract includes:

• enzymatic and mechanical transformation of starch (amylopectine and amylose) into absorbable glucose, • incretin secretion and effects on the blood sugar,

• and the small intestine microbiota, which may modulate dietary responses. This landscape of hypothetical causal factors shows that fundamental research effort on T2D must continue despite precise clinical understanding of the disease. However, all representations of glucose-insulin homeostasis largely underestimate the importance of the gastro-intestinal tract into the blood sugar consequences. Instead, they tend to model with increasing details the interaction of insulin with its related tissues (pancreas, liver and insulino-dependent tissues). We want to investigate the contribution of IGA to glucose homeostasis and its potential role in diabetes. To this aim, and as a preliminary work, we consider a typical and state of the art homeostasis glucose-insulin model [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] formalized as a system of Ordinary Differential Equations (ODEs). Our objective is twofold:

• test if this model can predict a significant improvement of glucose homeostasis by simulating RYGB as is observed experimentally, • test if this model can predict the time-course data of an original dataset of diabetic patients. In Section 2, we briefly describe the model. In Section 3, we present our parameter fitting results both from the original parameters of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] and for our own dataset We discuss the partial results in Section 4 and present the on-going and future work in Section 5.

II. MODEL

Many simulation models of the glucose-insulin system for the postprandial period have been developed [START_REF] Mari | Mathematical modeling for the physiological and clinical investigation of glucose homeostasis and diabetes[END_REF]. In this work, we consider a highly cited model, proposed in [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF], to simulate the postprandial physiological events of their own cohorts of normal subjects and T2D patients. This model is made of 12 ODEs and 36 parameters describing fluxes of glucose and insulin between physiological compartments: gastro-intestinal tract, plasma, liver, pancreas, muscle and adipose tissues (Figure 1). We recall in the following, in informal terms, how the physiological modules interact.

The Gastro-Intestinal Tract module describes the digestion process, from the stomach to the gut, and can be considered as the input of the whole system. It includes the complexity of the gastric emptying depending on the proportion of the solid and liquid phases of the alimentary bolus [START_REF] Dalla Man | A system model of oral glucose absorption: validation on gold standard data[END_REF]. Only the liquid phase that ends up in the gut is absorbed by the intestine and discharged to the peripheral blood. Incretins are not modeled. The Liver module describes the hepatic activity responsible for the endogenous glucose production down-regulated by the insulin. The glucagon hormone is not modeled. The Pancreas secretes insulin, which is up-regulated with the amount of glucose in the blood. The uptake of glucose by the (muscle and adipose) Tissues is regulated by the insulin. The Kidney is responsible for the glucose excretion-reabsorption. The insulin degradation is due to its lifespan and liver clearance. Even though the insulin independent glucose uptake by the brain is modeled (not in shown in Figure 1) the regulation of the glucose by the brain is neglected in [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF].

Using tracer-to-tracee ratio clamp technics [START_REF] Cobelli | Tracer-to-tracee ratio for analysis of stable isotope tracer data: link with radioactive kinetic formalism[END_REF], the authors of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] measured the time course of the glucose concentration in various compartments. This was done following a mixedmeal received by several normal and diabetic subjects. The parameters were estimated to fit this experimental data, which resulted in two sets of parameter values modeling respectively normal and T2D behaviors. For practical reasons, we fully reimplemented the model in the Julia programming language (version 1.5.3) with the DifferentialEquations package (version 6.17.0).

III. RESULTS

A. From T2D to normal model

In the following, we call normal model, resp. T2D model, the model instantiated with the normal, resp. diabetic, parameter values given in [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] (Table 1). We consider several parameter subsets corresponding to the previous modules: gastro-intestinal tract (also denoted as Ra in [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF]), liver (EGP), tissues (GK+U), pancreas (IK+S) and kidney (RE).

Starting from the T2D model, we re-estimated, in turn, each of these subsets of parameters, while leaving the other parameters fixed, in order to fit the plasma glucose dynamics of the normal model. Based on the data of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF], our first objective is to evaluate the capability of the model to predict, for each module alone, its capability to restore a normal glycemia. We estimate the parameters twice: with and without estimation of the basal values for insulin secretion by the pancreas, glucose production by the liver and utilization by the tissues. From the resulting 10 inferred models, we plotted the time course of the observed variables.

In this short paper, we only report (Figures 4a to 4d) the most relevant plots for our purpose: plasma glucose (G), plasma insulin (I) and the rate of intestinal absorption (Ra) for the two models obtained from the estimations of the gastrointestinal tract and pancreas compartments, with and without basal estimation. In order to compare the models' performance in fitting the normal model, we collect in a bar plot (Figure 3) the residual sum of squares for each model.

B. Parameter estimation of obese diabetics and RYGB

Our second objective is to test whether the model of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] can predict the time course concentrations of glucose and insulin obtained from our own dataset of diabetic patients who underwent RYGB surgery. For each patient, we use data before (hereafter referred to as visit A) and 3 months after (visit B) surgery. We first estimate all the parameters in order to fit the time course data of glucose and insulin from the visit A dataset. The model that we obtain is called the visit A model. Figure 2 shows the glucose and insulin plasma concentrations predicted by this model as well as the fitted data points. Then, as previously, we estimate each subset of parameters in order to fit the visit B dataset. Here, we only consider the case where we also estimate the basal concentrations, which indeed changed 3 months after surgery. We report in Figure 4 the time course of G, I and Ra after estimation of the parameters of the gastro-intestinal tract and pancreas.

IV. DISCUSSION

Estimating the basal values for insulin secretion, endogenous glucose production and insulino-dependent glucose utilization, can be interpreted as a prediction of the "longterm" effect of the parameter changes in the best case (since the model does not incorporate any long-term recovery mechanisms). Thus, not estimating these basal values can be interpreted as a "short-term" (or worst case) prediction.

A. From T2D to normal model

Our re-estimations of the parameters based on the data of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] predict (Figure 3) that the best performing compartments to restore a normal glycemia are pancreas and tissues, and then the intestinal tract. As expected, estimating the basal concentrations (i.e., long-term effect) improves the performance especially for the intestinal tract (see also G curve in Figures 4b and4c), which is consistent with experimental observation. However, the performance of the pancreas should be modulated. Indeed, Figures 4b and4a show that, in order to improve the glycemia, a very high plasma insulin concentration is necessary if only the pancreas parameters are modified. This seems physiologically unrealistic, meaning that the good performance of this compartment is over-estimated. Similarly, the estimation of the gastro-intestinal tract parameters on the short term (Figure 4d) indicates an unrealistic decrease of plasma insulin concentration. Finally, the "long-term" estimation of the gastro-intestinal tract parameters (Figure 4c) allows a good improvement of the glycemia with a realistic concentration of insulin and a decrease of intestinal absorption, which is consistent with observation.

B. Parameter estimation of obese diabetics and RYGB

From Figure 2, it is primarily interesting to notice that the T2D data of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] is significantly different from ours. Indeed, glucose in visit A is higher than in T2D model whereas insulin in visit A is lower than in T2D model. Despite this difference, the fitting is satisfying with all parameters set as free. This fitting is sensitive to the parameters estimation methods.

In Figures 4e and4f, the parameters are set free for the pancreas (IK+S) and the gastro-intestinal tract (Ra), respectively. Such process can be interpreted as surgery simulations targeting respectively the pancreas and the gastro-intestinal tract. No fitting attempt seems satisfying. On the one hand, freeing IK+S parameters seems to be satisfying for fitting the glycemia but clearly overestimates the insulinemia (cf. Figure 4e). On the other hand, by freeing Ra parameters, the rate of appearance is decreased as observed in experimental data (based on our D-xylose data, an alimentary glucose marker). Still, the parameter estimation fails completely to fit visit B glycemia and insulinemia (cf. Figure 4f).

V. CONCLUSION AND FUTURE WORKS

The parameter estimation performed by the authors of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] is based on training data generated by tracer-to-tracee clamp technique [START_REF] Cobelli | Tracer-to-tracee ratio for analysis of stable isotope tracer data: link with radioactive kinetic formalism[END_REF], which, despite its efficiency, remains an uncommon and complex method to monitor exogenous solutes. In practice, plasma glycemia and insulinemia are usually the only accessible clinical data. However, and as our model assessment suggests, this may raise parameter identifiability issues. To overcome this problem, we first plan to use the available additional D-xylose data, a marker that can be used to fit the rate of appearance (Ra) [START_REF] Bjarnason | Clinical biochemistry of the gastrointestinal tract[END_REF]. We also plan, by exploiting profile likelihood and sensitivity analysis, to study model reduction in order to eliminate the potential sources of non-identifiability. Other original datasets are currently used for the parameter estimation, generated from experiments on minipigs. Such biological models allow for more experiments and reproducibility, and decreased individual variabilities thus improving the reliability of parameter estimation. In this direction, another possibility, is to use publicly available datasets.

It should be noted that the failure of parameter estimation may be due to structural problems inherent to the model of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] which sub-model of IGA is largely simplified. For instance, it ignores the spatial none uniform glucose absorption rate along the intestin and the secretion of incretins. We plan to extend the model of [START_REF] Dalla Man | Meal simulation model of the glucose-insulin system[END_REF] with these aspects while simplifying the others to overcome identifiability issues that could emerge from additional parameters related to the gastro-intestinal tract.

Figure 1 :

 1 Figure 1: Simplified interaction graph of physiological compartments describing the model of [5].

Figure 2 :

 2 Figure 2: Parameter estimation of all parameters, including basal concentrations, to fit "Visit A" data.

Figure 3 :

 3 Figure 3: Loss from the parameters estimation applied on the different compartments.

  (a) Parameter estimation of IK+S compartment and basal concentration to fit the normal model for the T2D model. (b) Parameter estimation of IK+S compartment only to fit the normal model for the T2D model. (c) Parameter estimation of Ra compartment and basal concentration to fit the normal model for the T2D model. (d) Parameter estimation of Ra compartment only to fit the normal model for the T2D model. (e) Parameter estimation of IK+S compartment and basal concentrations, to fit visit B model from visit A model. (f) Parameter estimation of Ra compartment and basal concentrations, to fit visit B model from visit A model.

Figure 4 :

 4 Figure 4: Glucose (G), insulin (I), and rate of appearance (Ra) after parameter estimation, with and without basal concentrations, of pancreas (a, b), gastrointestinal tract (c, d) compartment for fitting normal model from TD2 model and for fitting visit B from visit A model (e, f).


J Hem

C Crauste-Thibierge

F Clément

D R Long

S Ciliberto

Simultaneous memory eects in the stress and in the dielectric susceptibility of a stretched polymer glass

We report experimental evidence that a polymer stretched at constant strain rate λ presents complex memory eects after λ is set to zero at a specic strain λw for a duration tw, ranging from 100 s to 2.2 × 10 5 s. When the strain rate is resumed, both the stress and the dielectric constant relax to the unperturbed state nonmonotonically. The relaxations depend on the observable, on λw and on tw. Relaxation master curves are obtained by scaling the time and the amplitudes by ln(tw). The dielectric evolution also captures the distribution of the relaxation times, so the results impose strong constraints on the relaxation models of polymers under stress and they can be useful for a better understanding of memory eects in other disorder materials.

Aging of amorphous and heterogeneous materials is characterized by a slow relaxation toward a new state after a perturbation has been applied [1]. This relaxation usually spans several orders of magnitude in time and may present memory eects, which depend on the sample history (see Ref. [START_REF] Keim | [END_REF] for a recent review).

Memory eects have been observed in glassy and disordered materials [318], and also in many other disordered systems [1926]. Memory eects may take the form of a Kovacs-like nonmonotonic relaxation, in which the response of a state variable, to an external perturbation, evolves in one direction before turning around at a timescale and an amplitude which depend on the sample history. These Kovacs-like relaxations have been observed in the mechanical [612] and electrical [13] responses of polymer glasses and melts [START_REF] Zapas | Deformation and Fracture of High Polymers[END_REF], in the static coecient of friction [START_REF] Rubinstein | [END_REF]20], in the load on crumpled thin sheets [21,22], and in the mechanical response of granular systems [2326].

Several theoretical models [2732] have been developed to explain, with common underlying physical principles, the similar features of these complex relaxations observed in very dierent systems. One of the key assumptions of these models is the existence in these disordered systems of a distribution of relaxation times (DRT) which is used to explain the origin of aging and of nonmonotonic relaxations. However, it is hard to discern between the different descriptions using experiments in which only one observable is measured. Thus to give more insight into this problem, we performed experiments on a stretched polymer glass where one gets quantitative information on the memory of several observables measured simultaneously and on the evolution of the DRT.

In this Letter we describe these experiments and report the observation that a stretched polymer glass exhibits complex memory eects of both the stress and of the dielectric susceptibility. Specically, we show that, after a perturbation of the strain rate, the relaxation dynamics of the stress and the dielectric susceptibility present nonmonotonic responses which depend not only on the sample history, as in the above-mentioned Kovacs-like memory eect, but also on the observable. Furthermore, from the evolution of the dielectric spectra, we get insights on the evolution of the polymer DRT. Thus our results impose constraints for modeling not only the universal features of memory eects but also the polymer relaxation dynamics under deformation, which is widely studied, both experimentally [3340] and theoretically [START_REF] Mckenna | third Workshop on Non-Equilibrium Phenomena in Supercooled Fluids, Glasses, and Amorphous Materials[END_REF]4145].

In our experiment the investigated polymer is an extruded lm MAKROFOL DE ® 1-1 000000 from BAYER based on Makrolon ® polycarbonate (PC) of thickness 125 µm. The lms are cut into rectangular pieces of size 42 × 21 cm with a dog bone shape before being tied to the tensile machine. During the experiment, we measure the instantaneous stress σ(t) and deformation λ(t) = L(t)/L 0 of the lms. The experiments are performed at imposed λ 10 -4 s and at T = 25 • C much below the glass transition temperature T g = 150 • C. In parallel to the stress measurement we use the dielectric spectroscopy to investigate the dynamics of the relaxation processes at work in the stretched PC lm [START_REF]Polycarbonate is a B-type dipole polymer[END_REF][START_REF] Kremer | Broadband Dielectric Spectroscopy[END_REF].

Dielectric spectroscopy is indeed a convenient way to probe in situ and without destruction the relaxation times of stretched polymers [3840]. It is worth to mention that other techniques, such as NMR measurements [33] or probe molecule experiments [34,35], can be used but dielectric spectroscopy is very versatile. In the rest of the Letter, we will focus only on the loss tangent tanδ(f) = -(f)/ (f), which does not depend on the lm thickness, being the ratio of the imaginary ( ) and real ( ) parts of the dielectric susceptibility, at frequency f . For the dielectric measure of our samples, the lm is sandwiched between two electrodes and conductive gel layers which ensure a uniform electrical contact. The electrodes are connected to a homemade dielectric spectrometer, which measures the dielectric spectrum every 32 s in the frequency range [10 -1 -10 3 ] Hz, simultaneously on all of the frequencies and with an accuracy of a few percent [START_REF] Pérez-Aparicio | [END_REF]. Thus it allows us to follow the evolution of the dielectric spectra during the stretching protocols which lasts more than 2 h. In order to get the required repeatability, we average the results over at least three tests performed in dierent samples submitted to the same stretching protocol. More details about the experimental apparatus can be found in Refs. [38,39,[START_REF] Pérez-Aparicio | [END_REF].

Hz and σ(t).

We rst recall the behavior of the stress and the dielectric susceptibility as a function of λ measured at constant strain rate λ, that we have discussed in Ref. [38,39]. Figure 1(a) depicts the stress evolution for polycarbonate at λ = 2.5 × -4 s -1 (black curve). The stress rst increases until the yield stress at about λ = 1.08, then, successively decreases, a regime known as strain softening, reaches a plateau at about λ = 1.20, and nally increases again, which is known as strain hardening. The breaking occurs at about λ = 1.8 -2 [START_REF]The sample deformation is homogeneous and there is no necking until the sample failure which occurs for λ ≥ 1[END_REF]. The dielectric spectra in the range [10 -1 -10 3 ] Hz have been measured every 32 s, during which λ has a negligible increase of less than 0.8% at this λ. As described in Refs. [3840], dielectric measurements give access to the evolution of the dominant relaxation time τ α of the polymer [50] by identifying the lower part of the dielectric spectrum as the high frequency tail of the α relaxation peak. In particular, we have shown in Ref. [39] that the increase (decrease) of the dielectric susceptibility in the range [10 -1 -10 3 ] Hz reects a decrease (increase) of τ α , which depends on λ and the strain rate λ. Therefore, the evolution of tanδ at f = 0.125 Hz depicted in Fig. 1(a) shows an acceleration of the relaxation times until the softening or onset of plastic ow (maximum around λ = 1.15 -1.17) and a slowdown during the strain hardening regime that we have studied in Ref. [39]. The results of the measurements performed at constant λ (black curves in Fig. 1(a)) will be called reference curves.

To investigate the memory eect, we have modied the protocol at constant λ in order to include an aging step at a specic λ in the post yield regime of the lm. More precisely, the lm is rst stretched at λ = 2.5 × 10 -4 s -1 then held at the waiting strain λ w during the waiting time t w , and nally further stretched with the same initial λ. As an example for the case λ w = 1.2 (plastic ow) and t w = 1 h, the evolutions of σ and of tanδ(0.125) are illustrated as a function of λ in Fig. 1(a). The corresponding time evolution is illustrated in Fig. 1(c) and 1(d) where the instantaneous strain rate λ is indicated, too. Three dierent intervals can be distinguished: During the initial one (λ < λ w and t < t wi ), σ and tanδ are identical to the reference curves [Fig. 1(a)]. During the waiting interval (i.e, λ w = 1.2 and t wi ≤ t ≤ t wf ), they decrease logarithmically with time, which is the signature of an aging process [see Fig. 1(b)]. Finally as soon as λ is resumed for t ≥ t wf , they increase again. In Fig. 1(a), the comparison with the reference curves obtained with the constant λ protocol shows that σ and tanδ(0.125)nally recover their reference curves. However, notice that in spite of the fact that σ and tanδ(0.125) have a similar dependence on λ before the aging step, the recoveries of the two observables are rather dierent [Fig. 1(a)]. In addition, we shall see that these recovery processes also depend on t w and on λ w .

Before presenting these results, we point out that Fig. 1(d) shows that the response of tanδ to λ increases by reducing the frequency f . From this frequency dependence we get useful information on the DRT evolution by studying the whole spectrum of tanδ. Thus we plot the tanδ spectra in Fig. 2(a) during the waiting interval at λ w = 1.2 and in Fig. 2(b) during the recovery for t > t wf . We rst notice that the spectra evolve in a signicant way only for f < 100 Hz. Furthermore, we see in Fig. 2(a) that the spectrum at λ = 1.2, at the beginning of the aging (t = t wi ), is identical to the reference one obtained during the constant λ protocol. Then, during the waiting interval for increasing ∆t wi = t -t wi , the spectrum progressively decreases towards the equilibrium spectrum (i.e, the spectrum at λ = 1). Inversely, in Fig. 2(b), the spectra during the recovery interval increase again until they reach the reference spectrum around λ = 1.3, i.e, ∆t = t -t wf = 400 s.

The main relaxation time can be estimated from these dielectric spectra following the method presented in Ref. [39] in which we have shown that at low frequencies tanδ(f) (sin βπ/2)/(2π f τ eff ) β where the characteristic time τ ef f and the exponent β have to be determined. Within the Cole-Cole [START_REF] Havriliak | [END_REF] model, τ ef f is proportional to the main relaxation time τ α of the polymer [52] and the exponent β is related to the width of the distribution of the relaxation times, specically the smallest is β < 1, the larger is the distribution [37,39]. Thus the t of the low frequencies part of the spectra in Figs. 2( To study the evolution of tanδ as a function of t w and λ w , we now select tanδ measured at f = 0.125 Hz as a typical amplitude of the low frequency spectra, which depend on τ α and β. Moreover, we nd that the waiting time has a strong impact on the relaxation. In the case of σ, Figs. 3(a) and 3(b) show that, as long as t w increases, the overshoot peak shifts progressively to longer time and presents a higher and sharper amplitude. In the case of tanδ, Figs. 3(c) and 3(d) show that an increase of t w also shifts the peak time and slows down the relaxation to the reference curve. Also in the dielectric case, we observe that a longer t w aects the peak amplitude.

Besides the waiting time, the choice of the waiting strain also modies the relaxation.

Figures 3(a all the relaxation curves, measured after dierent waiting times t w , can be superimposed on a master curve by scaling ∆σ and ∆t by ln(t w ), for t w > 100 s. Moreover, we observe that this relation holds regardless of the waiting strain value (λ w = 1.2 or λ w = 1.5). As a consequence, this means that the peak time and the peak amplitude depend on ln(t w ). Regarding the dependence of the peak properties on t w and λ w , we note that the peak amplitude ∆σ p is dierent for the two waiting strains, being higher at λ w = 1.2 than at λ w = 1.5 for all t w . In contrast the peak time ∆t p is independent on the waiting strain [Figs.

3(e) and 3(f )].

Finally, it is interesting to see that this way of scaling the data also works for the dielectric susceptibility. In a similar way as before, we have studied the evolution of the dierence between the tanδ curve at f = 0.125 Hz and the reference one ∆tanδ = tanδ -tanδ ref as function of ∆t. Figures 3(e) and 3(f ) show that all the relaxation curves merge on a master curve by rescaling ∆tanδ and ∆t by ln(t w ). This is observed again for the two waiting strains investigated, either at λ w = 1.2 [Fig. 3(c)] or at λ w = 1.5 [Fig. 3(d)]. Note that for the dielectric case, the peak time at λ w = 1.5 is delayed compared to λ w = 1.2, which was not the case for the stress relaxation. The logarithmic dependence on t w shows, in particular, that the memory of the aging history can also occur in tanδ, so that we provide evidence that the memory appears simultaneously in dierent observables but in dierent ways.

Concerning the scaling ∆t/ ln(t w ), we point out that in several systems with nonmonotonous relaxations, the memory eects present a ∆t/t w scaling [START_REF] Keim | [END_REF]22,30], which however is not generic. In several cases [4], a heuristic form of scaling ∆t/t µ w is used, which has been theoreti- cally justied [29]. If µ << 1 then ln t w cannot be distinguished from t µ w , and our scaling ∆t/ ln(t w ) can be understood within this theoretical framework. In these models the DRT does not depend on time whereas our measurements show that it does evolve, thus our measure of the DRT evolution will help in constructing more precise models. It remains to understand why the time scale of the relaxations are dierent for σ and tanδ. It is interesting to notice that for polymers the amplitude of the main yield stress peak also undergoes a similar logarithmic dependence on the aging time tw , which in this case is the time spent in the glassy state after a temperature quench from the melt [5356]. However, in contrast to our experiment, the positions of the yield peak do not depend on tw . This dierence could be related to the way the waiting step is realized. Indeed, when the polymer ages after a temperature quench, no stress is applied to the sample before the measure of the loading curve, whereas in our experiment the polymer ages from the plastic regime at constant strain. From this perspective, it must be concluded that stretching substantially aects the aging and the recovery processes, which is conrmed by the dependence of the memory eects on λ observed in our experiment. The evolution of the dielectric susceptibility and of the DRT when a stress is applied after the time tw has never been measured, and no compari- son with our results can be done. This comparison would be very useful for giving new insights into the polymer relaxations under stress.

In conclusion, we have shown that a stretched polymer presents complex memory eects in the nonmonotonous relaxations of σ and tanδ. Although nonmonotonous responses of σ have been already observed in stretched polymers [712], it is interesting to notice that tanδ behaves in a similar way as a function of λ w and t w . Nonmonotonous relaxation of the dielectric susceptibility has been observed [13] after an electric eld perturbation but never after a mechanical perturbation. This new observation is indeed very useful because the measure of tanδ as a function of frequency allows us to have information on the evolution of the DRT during the stretching, the waiting, and the recovery intervals. The result of this measure strengthens previous observations on DRT [36] and brings new insight to the long-lasting debate on mechanical rejuvenation in polymers (see, e.g, [START_REF] Mckenna | third Workshop on Non-Equilibrium Phenomena in Supercooled Fluids, Glasses, and Amorphous Materials[END_REF]39]).

The observed timescaling and DRT evolution is also use-

FIG. 1 :

 1 FIG. 1: (a) Dielectric susceptibility tanδ at f = 0.125 Hz and nominal stress σ = F/So with So = 210 × 0.125 mm 2 as a function of the strain λ. Black lines: constant strain rate protocol λ = 2.5 × 10 -4 s -1 , called reference curves. Red lines: relaxation protocol λw = 1.2 and tw = 1 h. The red line thickness corresponds to one standard deviation. (b) Temporal evolution tanδ at f = 0.125, 1, 7 Hz and σ during the waiting step over tw = 1 h. (c) λ(t) and σ(t) during the relaxation protocol with λw = 1.2, tw = t wf -twi = 1 h. Times twi and t wf delimit the beginning and the end of the waiting interval, respectively. (d) tanδ(t) at f = 0.125, 1, 7 Hz and σ(t).

  FIG. 2: Dielectric spectra tanδ recorded at dierent times for the relaxation protocol at λw = 1.2 and tw = 1 h.(a) Evolution during the aging interval at λ = 1.2 with ∆twi = t -twi and (b) during the recovery interval for λ > 1.2. Dashed lines correspond to the reference spectra recorded at constant λ. The evolutions of τ ef f and β are plotted in (c) and (d) respectively. Black line: constant λ protocol. Red line: relaxation protocol.

Figs. 3

 3 (a)-3(f ) bring more details on the recovery processes by showing dierent results of tests performed with dierent waiting times and wait-ing strains λ w = 1.2 (plastic ow) in Figs. 3(a) and 3(c) and λ w = 1.5 (strain hardening) in Figs. 3(b) and 3(d). When λ is resumed at t wf , we observe that the relaxations of σ and tanδ are both nonmonotonic and exhibit an overshoot above the reference curve before relaxing to it, with dierent shapes of the evolutions of the two observables. On one side, the stress overshoot has a distinct and large amplitude while occurring rapidly after the resumption [Figs. 3(a) and 3(b)]. In contrast, the overshoot of tanδ has a weaker amplitude and occurs much later in time [Figs. 3(c) and 3(d)].

  ) and3(b) show qualitatively that the stress overshoot peak is higher and narrower at λ w = 1.2 than at λ w = 1.5 and at the same time Figs.3(c) and 3(d) show that tanδ recovers much faster the reference curve at λ w = 1.2 than at λ w = 1.5. The dierence is signicant. For example at t w = 1 h and λ w = 1.2, tanδ recovers the reference curve at about λ = 1.45, equivalent to ∆t = t -t wf = 1024 s after the resumption of the stretching. This is much faster than for the case at t w = 1 h and λ w = 1.5 for which the recovery is even too long to occur during our experiment (λ > 1.85 equivalent to ∆t > 1240 s). In particular, this enhancement of the mobility recovery delay with λ w fully supports the long one observed in multicreep experiments[36], in which the aging relaxation is performed in the strain hardening regime.To be more precise regarding the dependence of the relaxations on those aging parameters we analyze the difference between the stress curve and the reference one ∆σ = σ -σ ref as a function of the time ∆t = t -t wf elapsed after resumption. Figures 3(e) and 3(f ) show that

FIG. 3 :

 3 FIG. 3: (a)-(d) Stress σ and dielectric susceptibility tanδ at f = 0.125 Hz as a function of λ for dierent λw and tw. Curves obtained by averaging three tests with tw = 62 h and four tests with tw = 100 s and 1 h. (a),(c) σ and tanδ at λw = 1.2 (in the plastic ow regime). (b),(d) σ and tanδ at λw = 1.5 (in the strain hardening). The reference curves are indicated in black. Note the stress axis has been expanded in (a) and (b) with respect to Fig. 1(a). (e),(f) Evolution of ∆σ/ ln(tw) (continuous lines) and ∆tanδ/ ln(tw) (symbols) as a function of ∆t/ ln(tw) for t ≥ t wf , at (e) λw = 1.2 and at (f) λw = 1.5 (see text for details).

  ful when one tries to t this observation into the more general models of memory eects which have been developed to explain the experimental observation in glassy and disordered systems. Thus our results, by showing important features of the dynamical processes of two interrelated quantities and of the polymer relaxation dy-

namics, challenge the theoretical models of the memory eect in glassy materials and disordered systems.

J.H. acknowledges the funding by Solvay. We thank O. Razebassia, A. Petrosyan and the Bartolo team for technical assistance on the tensile machine.


Sylvain Tourancheau
email: sylvain.tourancheau@univ-nantes.fr

Florent Autrusseau
email: florent.autrusseau@univ-nantes.fr

Z M Parvez Sazzad

Yuukou Horita
email: horita@eng.u-toyama.ac.jp

IMPACT OF SUBJECTIVE DATASET ON THE PERFORMANCE OF IMAGE QUALITY METRICS

Keywords: Image quality, Quality assessment, Subjective database

The interest in objective quality assessment have significantly increased over the past decades. Several objective quality metrics have been proposed and made publicly available, moreover, several subjective quality assessment databases are distributed in order to evaluate and compare the metrics. However, several question arises: are the objective metrics behaviours constant across databases, contents and distortions? how significantly the subjective scores might fluctuate on different displays (i.e. CRT or LCD)? which objective quality metric might best evaluate a given distortion? In this article, we analyse the behaviour of four objective quality metrics (including PSNR) tested on three image databases. We demonstrate that the performances of the quality metrics can strongly fluctuate depending on the database used for testing. We also show the consistency of all metrics for two distinct displays.

INTRODUCTION

As subjective experiments are extremely tedious and time consuming, objective quality metrics have lastly been extensively studied, and their performances have significantly improved over the decades. Whereas the aim of an objective metric is to substitute to tedious subjective experiments, while designing a metric, it is mandatory to compare the objective scores with the subjective ones, and to use statistical tools to accurately evaluate the objective quality assessment. Objective quality metrics may be used for various purposes, but they are usually designed and utilized within image and video compression context. However, the interest in objective quality metrics for other image and video processing applications (such as digital watermarking) has also recently increased. Evidently, the performances of objective metrics stongly rely on the tested distortion. For instance, the quality range of distortions induced

The authors would like to thank Romuald Pépion for his assistance to obtain the results described in the paper.

by watermarking techniques will be very narrow compared to the quality range of compressed data. This study is motivated by the need to faithfully analyse the metrics performances, and to demonstrate that, to be relevant, such analysis has to be conducted on a wide set of distortions and contents. It is thus strongly recommended to compare the metrics to several of the publicly available databases. One of the question we hereby address here, is the monotonicity of the metrics performances for various subjective databases. Four metrics were tested, a statistical metric (PSNR), one advanced metric exploiting Human Visual System features (C4), one based on an information-theoretic framework (VIF) and one based on structural similarities (SSIM). Three subjective databases are used, linear correlation and root mean squared error are used to assess the metrics. Strong variations among the four tested metrics within, and across the databases are observed. This may be explain by a somewhat important disparity among the databases. Besides, we also propose here a subjective study pointing out the important similarities between subjective quality assessment on both CRT and LCD monitors, since two subjective experiments, on the same dataset, were conducted independently on both displays. The analysis of both experiments will show the very weak variations of subjective data, which moreover can easily be induced by cultural factors and lab effects, rather than by the displays themselves. Similarly, a comparison is performed between the release 2 of the LIVE database and an update after a raw scores realignment processing. and shows minor difference behaviour for certain metrics. This paper is organized as follows: Section 2 describes the experiments conducted to construct our datasets, the performances of the quality metrics on each dataset are computed in Section 3 and results are discussed in Section 4.

DESCRIPTION OF THE SUBJECTIVE DATASETS

IRCCyN/IVC database

The IRCCyN/IVC subjective database [START_REF] Horita | Subjective quality assessment IRCCyN/IVC database[END_REF] consists in 10 original colour images with a resolution of 512×512 pixels from which 235 distorted images have been generated, using 4 dif-ferent process (JPEG, JPEG2000, LAR coding, Blurring). These algorithms have the advantage to generate very different type of distortions. Each distortion type have been optimized in order to uniformly cover the whole range of quality.

Subjective evaluations have been performed in a normalized room with lighting conditions and display settings adjusted according to ITU recommendation BT.500-11. The viewing distance was set to six times the picture's height. Fifteen observers participated to the experiments, they have been checked for visual acuity and color blindness.

A double stimulus impairment scale (DSIS) method have been used. Both distorted and original pictures were displayed sequentially. At the end of the presentation, the observer was asked to assess the annoyance he/she felt on the distorted image with respect to the original one. The impairment scale contained five categories marked with adjectives and numbers as follows: 5 -Imperceptible, 4 -Perceptible but not annoying, 3 -Slightly annoying, 2 -Annoying and 1 -Very annoying. Mean opinion score was then computed over the observers after the potentiel rejection of observers according to recommendations.

Toyama database

The Toyama subjective database [START_REF] Horita | Image quality evaluation database[END_REF] contains 182 images of 768×512 pixels. Out of all, 14 were original images (24 bit/pixel RGB) in each group. The rest of the images were JPEG and JPEG2000 coded images (i.e. 84 compressed images for each type of distortion). Six quality scales and six compression ratios were respectively selected for the JPEG and JPEG2000 encoders. The following codec softwares were used to generate the compressed images: JPEG using cjpeg software , and JPEG2000 with JasPer software .

Subjective experiments were conducted in a normalized room with low lighting conditions and display settings adjusted according to ITU-R BT.500.11. The viewing distance was set to four times the picture's height. Prior to participating the session all subjects were screened for normal visual acuity with or without glasses, normal color vision and familiarity with language. Sixteen non-expert subjects were shown the database; most of them were college students.

Single stimulus absolute category rating (SSACR) method was used in these subjective experiments. The subjects were asked to provide their perception of quality on a discrete quality score that was divided into five and marked with the numerical value of adjectives: Bad (1), Poor (2), Fair (3), Good (4) and Excellent (5). Note that the numerical values attached to each category were only used for data analysis and were not shown to the subjects. At the end of each test presentation, observers provide a quality rating using the adjective scale. The test presentation order was randomized according to standard procedure and the raw scores were collected in a data file by the computer. As the original images has been assessed as well, scores was converted in difference scores (DMOS) for each observer by computing the difference between the score obtained by the original image and the one obtained by the distorted image. Difference mean opinion scores (DMOS) were then computed for each image, after the screening of post-experiment results (most subjects had no outliers) according to ITU-R Rec. 500-10.

The Toyama database has been assessed on a CRT display in the University of Toyama in Japan. In order to check if the display is of central importance in such an experiment, we decided to conduct the same subjective quality assessment in IRCCyN laboratory in France, using a liquid crystal display (LCD). However, design the same experiment in two different labs it's a real challenge and the so-called "lab effect" can occur. Actually, even if set as similar as possible, the viewing conditions can differ from on testing room to another. Furthermore, using two different pools of observers can also lead to slight differences. Also, some cultural differences can appear between France and Japan, in the way to assess quality. For example, the way that observers consider the adjectives on the quality scale can be different. By the way, these two distinct experiments permits to obtain two different subjective datasets from the same image database.

LIVE database

The LIVE database release 2 [START_REF] Sheikh | LIVE image quality assessment database release 2[END_REF] contains 779 distorted images with five distortions type: JPEG, JPEG2000, white noise, gaussian blur and bit errors in JPEG2000 bit stream. Subjective quality scores have been published in form of difference scores in a quality range from 0 to 100. For more precisions concerning the experiments, please see Ref. [START_REF] Sheikh | LIVE image quality assessment database release 2[END_REF]. Recent work [START_REF] Sheikh | A statistical evaluation of recent full reference image quality assessment algorithms[END_REF] from LIVE laboratory have presented a new method to realigned subjective quality scores obtained from each session with each other. This realignment process used Z-scores transform in order to attenuate the inter-observers differences, following by a inverse transform to go back to the 0-100 quality range. This inverse transform has been adapted for each session. Following this work, an update of the subjective quality scores has been published online. In this paper, we will work on both datasets (release 2 and update).

QUALITY METRICS

First of all, objective image and video quality metrics can be classified according to the availability of the distortion free image/video signal, which may be used as a reference to compare an original image or video signal against its distorted counterpart. Specifically, such metrics are usually of three kinds. Full Reference (FR) quality assessment metrics for which the exact original image is needed in order to assess the visual quality of any distorted version. Reduced Reference (RR) quality metrics, for which a reduced form of the original images is used. No reference (NR) metrics, where only the distorted image is needed. Several quality metrics (QM) have been used on the five datasets described previously:

• Structural SIMilarity (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] is an objective metric for assessing perceptual image quality, working under the assumption that human visual perception is highly adapted for extracting structural information from a scene. Quality evaluation is thus based on the degradation of this structural information assuming that error visibility should not be equated with loss of quality as some distortions may be clearly visible but not so annoying.

• C4 [START_REF] Carnec | Objective quality assessment of color images based on a generic perceptual reduced reference[END_REF] is a metric based on the comparison between the structural information extracted from the distorted and the original images. What makes this metric interesting is that it uses reduced references containing perceptual structural information and exploiting an implementation of a rather elaborated model of the HVS.

• Visual Image Fidelity (VIF) [START_REF] Sheikh | Image information and visual quality[END_REF] is a full reference quality metric using the concept of information fidelity measurement. The VIF metric groups Natural Scene Statistics as well as HVS features and a model of the considered distortion to quantify the loss of information.

Two classic criteria have been chosen to evaluate the performance of QM on each database: the linear correlation coefficient (Pearson's correlation) and the root mean squared error (RMSE). These two values have been computed after a nonlinear regression on the results of the QM. This regression is performed to map the output of each QM to the quality range of the DMOS. This regression has been done for the whole databases, and also separately for the JPEG and JPEG2000 subsets in each database. The non-linear function used to compute the regression was a logistic function with three parameters, as described in Eq. 1. The values of the three parameters are optimized in order to minimize the RMSE.

Q mapped = a 1 + e -b(Q-c) (1)
Figure 1 presents the linear correlation coefficient (top line) and the RMSE (bottom line) for the three QM as well as for PSNR, on the five datasets. Vertical bars indicate the 95% confidence interval of each correlation values (computed after a Fisher transformation in order to transpose into a normal distribution). This information permits to determine if two correlations are statistically distinguishable or not. In this figure, results are doubled for the LIVE database (corresponding to release 2 and updated version) and for the Toyama database (corresponding to CRT and LCD dataset). Linear correlation between PSNR and JPEG subset of both Toyama subjective datasets was under 0.6.

DISCUSSION

As expected, the performances of quality metrics show more or less the same tendency for the five subjective datasets. VIF and C4 appears to be the best QM, with some very close results as well for correlation (where 95% confidence intervals are highly overlapped) as for RMSE. SSIM obtains some lower results, particularly for IVC and Toyama databases and finally PSNR gets the lower results for all datasets. The differences over databases, particularly the decrease of performance on IVC and Toyama databases for all QM might be explained by the difference of quality range in the databases. As illutrated in Figure 2, low quality anchors in the LIVE database are indeed strongly distorted pictures with extremely low quality. The corresponding low anchors in Toyama and IVC databases appears to have a sensibly better quality. It is observed that quality metrics based on the statistics of images such as PSNR provide better results on the LIVE database than on databases with a narrower range of quality such as IVC and Toyama. Actually, statistical metrics may be less accurate in the high quality range since quality perception in these area is mostly due to perceptual HVS features, rather than to the statistics of the image. The good performance of advanced quality metrics (VIF and C4), on all databases, seems to confirm this assumption.

The results showed that the type of display does not significantly interfere on the subjective scores. Actually, this was an expected behaviour since both subjective datasets on CRT (performed in the University of Toyama) and on LCD (performed in the IRCCyN lab.) were highly correlated (Pearson's correlation of 0.957) and quality metrics results on both datasets were, of course, the same. However, it is somewhat remarkable to have such close results on two different subjective datasets, despite the different displays, different labs and cultural effect. This is an evidence that subjective quality assessment experiments are quite reliable and that differences between databases are mostly due to difference of contents, distortions type and quality range.

We analyzed in this paper various quality metrics across subjective databases. We showed that, to be relevant, the performances evaluation of a quality metric has to be experimented on several databases. Evaluating a metric with one single subjective database might not be sufficient as the quality range of the database seems to be of central importance. As a result, it is important to be aware of the differences between subjective datasets, as well as to know what is the profile of the tested quality metric.

Fig. 1 :

 1 Fig. 1: Top line: Linear correlation coefficient between QM and subjective scores for the three databases. Bottom line: Root Mean Square Error between QM and subjective scores for the three databases.Results are computed on the whole databases (white bars) and independently on the JPEG and JPEG2000 subsets, gray and black bars respectively.

Fig. 2 :

 2 Fig. 2: Illustration of the difference of low anchor between two databases. (b-c): lowest quality pictures for JPEG distortion respectively in Toyama and LIVE databases.


Anas Makdesi
email: anas.makdesi@l2s.centralesupelec.fr

Antoine Girard
email: antoine.girard@l2s.centralesupelec.fr

Laurent Fribourg
email: fribourg@lsv.fr

Efficient Data-Driven Abstraction of Monotone Systems with Disturbances

Keywords: Monotone transition systems, data-driven control, symbolic control, abstraction

In this paper, we present a novel approach for the abstraction of monotone systems with bounded disturbances. The approach is data-driven and uses a given set of samples of the (unknown) dynamics of the system to compute an abstraction defined on partitions of the state and input spaces. The proposed method is efficient as its computational complexity is linear in the number of samples and in the size of the partitions. Moreover, the abstraction is shown to be minimally conservative in the absence of disturbances. We show that the resulting symbolic model is itself a monotone transition system and is related to the original system by an alternating simulation relation. We present some numerical experiments to show the effectiveness of the approach and to show how the choice of the partitions or the number of samples affects the quality of the approximation.

INTRODUCTION

Symbolic control is a computational approach to controller synthesis, based on the use of discrete abstractions (a.k.a. symbolic models) of continuous dynamical systems and of discrete controller synthesis techniques (see e.g. [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]; [START_REF] Belta | Formal methods for discrete-time dynamical systems[END_REF]). The main advantages of symbolic control are as follows. Firstly, it can be applied to general classes of nonlinear systems with input and state constraints, as long as one is able to compute an overapproximation of the system's reachable sets, using e.g. approaches based on mixed-monotonicity [START_REF] Coogan | Finite abstraction of mixed monotone systems with discrete and continuous inputs[END_REF]) or on growth bounds [START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF]). Secondly, since controller synthesis is handled using algorithmic approaches developed in the fields of supervisory control [START_REF] Ramadge | Supervisory control of a class of discrete event processes[END_REF]) or formal methods [START_REF] Bloem | Synthesis of reactive (1) designs[END_REF]), the type of specifications that can be addressed is broad: safety and reachability (Girard (2012)), behaviors described by automata [START_REF] Pola | Control of cyber-physical-systems with logic specifications: a formal methods approach[END_REF]), or temporal logic formulas [START_REF] Belta | Formal methods for discrete-time dynamical systems[END_REF]).

While symbolic control is usually presented as a modeldriven technique, some approaches such as those for monotone systems [START_REF] Meyer | Safety control with performance guarantees of cooperative systems using compositional abstractions[END_REF]) only require us to be able to sample the system dynamics on a given grid of states and inputs, and it seems natural to extend these techniques towards a purely data-driven approach. In this paper, we present a data-driven approach to abstraction of monotone systems with disturbances. We assume that we are given a set of samples of the dynamics of the system. Contrarily to model-driven approaches, the samples are assumed to be taken at random and given a priori. In addition to the data, we assume that we know bounds on the amplitude of disturbances affecting the system. Using only this information, we present an approach to compute a symbolic model of the system defined on partitions of the state and input spaces. The symbolic model and the original system are shown to be related by an alternating simulation relation [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]) and therefore, the symbolic model can be used to synthesize controllers that are known to be correct by design. Moreover, we show that in the absence of disturbances, the abstraction is minimally conservative in the sense that it is alternatingly simulated by all monotone systems that are consistent with the data, but not more. In addition, we show that the symbolic model is itself a monotone transition system, which makes it amenable to efficient symbolic control synthesis techniques presented in [START_REF] Saoud | Efficient synthesis for monotone transition systems and directed safety specifications[END_REF]).

There are few approaches in the literature that deal with data-driven symbolic control. In [START_REF] Hashimoto | Learning-based safe symbolic abstractions for nonlinear control systems[END_REF]), an approach is presented for Lipschitz dynamical systems with partially unknown dynamics, the unknown dynamics being assumed to belong to some regression kernel. The approaches in [START_REF] Milanese | Set membership identification of nonlinear systems[END_REF]; [START_REF] Sadraddini | Formal guarantees in data-driven model identification and control synthesis[END_REF]) also apply to Lipschitz dynamical systems with known bounds on the Lipschitz constants but use piecewise affine abstractions. In comparison, our approach does not assume that the system is Lipschitz but that it is monotone. While monotonicity is a strong property, monotone systems can be found in many practical applications such as adaptive cruise control [START_REF] Saoud | Efficient synthesis for monotone transition systems and directed safety specifications[END_REF]), temperature regulation in buildings [START_REF] Meyer | Safety control with performance guarantees of cooperative systems using compositional abstractions[END_REF]) or power networks [START_REF] Zonetti | Decentralized monotonicity-based voltage control of dc microgrids with zip loads[END_REF]), to name a few.

In our previous work [START_REF] Makdesi | Datadriven abstraction of monotone systems[END_REF]), only the case without disturbance was considered, meanwhile, the approach presented in this paper can be applied to monotone systems with bounded disturbances. Moreover, the complexity of the approach presented in this paper is linear with respect to the number of samples of the dynamics and to the size of the partitions, when the complexity of the previous approach was polynomial of order n + m where n and m are dimensions of state and input spaces.

The paper is organized as follows. In Section 2, we provide some preliminaries on monotone transition systems and their abstractions. In Section 3, we give a formal statement of the problem under consideration in the paper. Section 4 presents the main contributions of the paper by describing an efficient data-driven approach to abstraction of monotone systems with disturbances. Finally, Section 5 presents some numerical experiments to assess the effectiveness of our approach and to show how the amount of data and the choice of state and input space partitions affect the quality of the approximation.

PRELIMINARIES

In this section, after defining some notations, we introduce the necessary background on monotone transition systems.

Notations

R, R + 0 and N denote the sets of real, non-negative real, and natural numbers, respectively. For a vector x ∈ R n , we denote by x i its i-th coordinate, i = 1, . . . , n. We define a partial order on R n as follows: for x, x ∈ R n x x if and only if x i ≤ x i , for all i = 1, . . . , n. Similarly, x ≺ x if and only if x i < x i , for all i = 1, . . . , n. We denote

[x, x ] = {y ∈ R n | x y x } and [x, x ) = {y ∈ R n | x y ≺ x }.
The empty set is denoted by ∅. Given a set X we denote 2 X to the set of subsets of X. A relation R ⊆ X ×Y is identified with the set-valued map R : X → 2 Y where R(x) = {y ∈ Y |(x, y) ∈ R}.

Transition systems

Transition systems provide a unifying framework for considering continuous and discrete systems:

Definition 1. A transition system Σ is a tuple Σ = (X, U, F, Y, H), where X is a set of states, U is a set of inputs, F : X × U → 2 X is a transition relation, Y is a set of outputs, and H : X → Y is an output map.
Σ is finite if X and U are finite and it is said to be deterministic if for all x ∈ X and all u ∈ U , F (x, u) is singleton. We call any x ∈ F (x, u) a u-successor of the state x. In this work, we are assuming, for simplicity, that for all x ∈ X, for all u ∈ U , F (x, u) = ∅.

Alternating simulation relations make it possible to relate the behaviors of transition systems (see e.g. [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]). Let us consider two transition systems,

Σ i = (X i , U i , F i , Y i , H i ), i = 1, 2 sharing the same sets of outputs (Y 1 = Y 2 = Y ). Definition 2. R ⊆ X 1 × X 2 is an alternating simulation relation from Σ 1 to Σ 2 if the following conditions are satisfied: • for all x 1 ∈ X 1 , there exists x 2 ∈ X 2 such that (x 1 , x 2 ) ∈ R; • for all (x 1 , x 2 ) ∈ R, H 1 (x 1 ) = H 2 (x 2 ); • for all (x 1 , x 2 ) ∈ R, for all u 1 ∈ U 1 , there exists u 2 ∈ U 2 such that for all x 2 ∈ F 2 (x 2 , u 2 ), there exists x 1 ∈ F 1 (x 1 , u 1 ) satisfying (x 1 , x 2 ) ∈ R.
We say that Σ 1 is alternatingly simulated by Σ 2 if there exists an alternating simulation relation from Σ 1 to Σ 2 . This is denoted by Σ 1 AS Σ 2 .

It can be shown [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]) that if Σ 1 is alternatingly simulated by Σ 2 , then any controller for Σ 1 can be refined to a controller for Σ 2 so that the closed loop output behavior of Σ 2 is included in that of Σ 1 . This means that a controller synthesis problem with a specification on the output space for Σ 2 can be solved by considering the same problem for Σ 1 .

Monotone systems

We now introduce the class of monotone transition systems. Let us consider the transition system Σ = (X, U, F, Y, H), where X ⊆ R n , U ⊆ R m . Considering the partial order on R n and R m , we define monotone systems as follow:

Definition 3. Σ = (X, U, F, Y, H) is a monotone transition system if for all x, x ∈ X, u, u ∈ U with x x , u u , ∀y ∈ F (x, u), ∃y ∈ F (x , u
), y y , and ∀y ∈ F (x , u ), ∃y ∈ F (x, u), y y . Figure 1 shows the transition relations of two pairs of states and inputs. We can see that under our definition for the monotonicity, we do not require all the u -successor of x to be bigger than u-successor of x. In the case of deterministic systems, the equation in Definition 3 becomes:

def2.pdf x 1 x 2 x 1 x 2 x x u u F (x , u ) F (x, u)
x x , u u =⇒ F (x, u) F (x , u )
(1) which coincides with the classical definition of monotone systems (see e.g. [START_REF] Coogan | Finite abstraction of mixed monotone systems with discrete and continuous inputs[END_REF]).

PROBLEM FORMULATION

Let us consider a transition system Σ = (X, U, F, Y, H), defined on sets of states X = R n and inputs U = R m and whose dynamics is given by the set-valued map

F : X × U → 2 X . F (x, u) = F(x, u) + W, x ∈ X, u ∈ U
(2) where F : X × U → X is an unknown monotone deterministic map (i.e. satisfying (1)), and W ⊆ R n is a bounded set of disturbances with known lower and upper bounds w, w ∈ R n , i.e. W = [w, w]. It is straightforward to show that Σ is a monotone transition system. The output set Y and the output map H will be specified later.

In the following, we consider the system Σ as a black box, and we only have access to data generated by the system. Let us assume that we are given data consisting of a set of transitions

D = {(x k , u k , y k ) | y k ∈ F (x k , u k ), k = 1 . . . K}.
The main problem considered in this paper is that of computing from the data D, a symbolic model Σ D that is alternatingly simulated by Σ. Indeed, as mentioned previously, that would make it possible to use Σ D to synthesize a controller for Σ.

For i = 1, . . . , n, j = 1, . . . , m, let be given finite partitions (χ i q i ) q i ∈Q i and (ν j p j ) p j ∈P j of R where

Q i = {0, . . . , N i }, P j = {0, . . . , M j } and    χ i 0 = (-∞, a i 1 ) χ i q i = [a i q i , a i q i +1 ), q i = 1, . . . , N i -1 χ i N i = [a i N i , +∞)    ν j 0 = (-∞, b j 1 ) ν j p j = [b j p j , b j p j +1 ), p j = 1, . . . , M j -1 ν j M j = [b j M j , +∞) where a i 1 < • • • < a i N i and b j 1 < • • • < b j M j . Let us define Q = Q 1 × • • • × Q n and P = P 1 × • • • × P m
and consider the finite rectangular partitions (X q ) q∈Q , (U p ) p∈P of X and U , given for q = (q 1 , . . . , q n ) and p = (p 1 , . . . , p m ) by

X q = χ 1 q 1 × • • • × χ n q n , U p = ν 1 p 1 × • • • × ν m p m .
We denote by x q , x q , u p , u p the infimum and supremum, with respect to partial order , of X q and U p . Note that some components may be infinite if the intervals are unbounded.

We now define the set of outputs of Σ as Y = Q and its output map as H given by H(x) = q if and only if x ∈ X q .

We now provide a formal statement of our problem: Problem 1. Given disturbance lower and upper bounds W , W , the data set D and partitions (X q ) q∈Q , (U p ) p∈P , compute a symbolic model Σ D = (Q, P, F , Y, H ) such that Σ D AS Σ.

DATA DRIVEN SYMBOLIC MODEL

In this section, we provide a solution to Problem 1. Based on the data set and the given partitions, we are looking for a minimal over-approximation of the unknown system dynamics.

Over-approximation of the transition relation

We define a box-shaped map G : Q × P → 2 X with an upper bound G and a lower bound G: i.e. G(q, p) = [G(q, p), G(q, p)], for all q ∈ Q, p ∈ P . G(q, p) is defined component-wise as follows for i = 1, . . . , n:

G i (q, p) =      min (x k ,u k ) (xq,up) y i k + w i -w i if ∃(x k , u k ) (x q , u p ), +∞ otherwise.
(3) Similarly, G(q, p) is defined component-wise as follows for i = 1, . . . , n:

G i (q, p) =        max (x k ,u k ) (x q ,u p ) y i k + w i -w i if ∃(x k , u k ) (x q , u p ), -∞ otherwise. ( 4
)
The following result shows that G provides an overapproximation of the unknown map F : Proposition 4. For all q ∈ Q, p ∈ P , for all x ∈ X q , u ∈ U p , F (x, u) ⊆ G(q, p).

Proof. Let us consider q ∈ Q, p ∈ P , x ∈ X q , u ∈ U p . If there exists (x k , u k ) (x q , u p ), we have for all i ∈ {1, . . . , n}:

F i (x q , u p ) + w i ≤ F i (x k , u k ) + w i for all k such that (x k , u k ) (x q , u p ). Therefore, F i (x q , u p ) + w i ≤ min (x k ,u k ) (xq,up) F i (x k , u k ) + w i .
We also have, for all k = 1, . . . , K,

F i (x k , u k ) + w i ≤ y i k . Hence, F i (x q , u p ) + w i ≤ min (x k ,u k ) (xq,up) y i k + w i -w i ≤ G i (q, p).
We have also F i (x q , u p ) + w i ≤ G i (q, p) when G i (q, p) = +∞. In a similar way we can show that F i (x q , u p ) + w i ≥ G i (q, p). Hence, F (x, u) ⊆ G(q, p). 2

Actually, we can claim a slightly stronger result than that of Proposition 4:

Claim 5. Let F : X × U → 2 X be of the form (2) and such that y k ∈ F (x k , u k ), for all k = 1, . . . , K, we have F (x, u) ⊆ G(q, p), for all q ∈ Q, p ∈ P , for all x ∈ X q , u ∈ U p .

The proof of the previous claim is identical to that of Proposition 4 and is therefore omitted. The previous claim essentially states that G is an over-approximation of all maps of the form (2) that are consistent with the data. The following result states that in the absence of disturbances, it is actually the minimal over-approximation: Proposition 6. Let w = w = 0, let G : Q × P → 2 X be a box-shaped map such that the statement of Claim 5 holds. Then, we have for all q ∈ Q, p ∈ P , G(q, p) ⊆ G(q, p).

Proof. Let us consider q ∈ Q, p ∈ P . From (3), we get that for all (x k , u k ) (x q , u p ), F(x k , u k ) G(q, p). By

(3) and monotonicity of F, we also have for all (x k , u k ) (x q , u p ), F(x k , u k ) G(q, p). Then, we can choose a monotone map F : X × U × X such that F(x k , u k ) = F(x k , u k ), for all k = 1, . . . , K and F(x q , u p ) = G(q, p). Since Claim 5 holds, it follows that F(x q , u p ) ∈ G(q, p). Therefore, G(q, p) ∈ G(q, p). Using a similar approach, we can show that G(q, p) ∈ G(q, p). Since G is a boxed-shaped map, we get [G(q, p), G(q, p)] ⊆ G(q, p). 2

The previous proposition shows that the conservatism of our approach is minimal in the absence of disturbances. In the presence of disturbances, the absence of information on the actual value of the disturbance at samples prevents us from achieving minimal conservatism.

Resulting symbolic model

We can now define the following transition system Σ D = (Q, P, F , Y, H ) where for all q ∈ Q, p ∈ P , F (q, p) = {q ∈ Q| X q ∩ G(q, p) = ∅} and H (q) = q.

Let us remark that, formally Q and P are finite subsets of R n and R m so they are equipped with partial orders. Let us state the main result of the paper: Theorem 7. Σ D AS Σ and the alternating simulation relation is given by R = {(x, q) ∈ X × Q| x ∈ X q }. Moreover, the symbolic model Σ D is a monotone transition system.

Proof. Let us start by the alternating simulation relation. It is obvious from the definitions of R, H, H that the first two conditions of the alternating simulation relation are satisfied. For the last condition, we have from Proposition 4, F (x, u) ⊆ G(q, p) for all q ∈ Q, p ∈ P , for all x ∈ X q , u ∈ U p . Therefore, for all (x, q) ∈ R, for all p ∈ P , there exists u ∈ U p such that for all x ∈ F (x, u), where x ∈ X q , x ∈ G(q, p) implies that we have q ∈ F (q, p). From the definition of x , q we have (x , q ) ∈ R.

We now prove the monotonicity of the symbolic model. Let q 1 , q 2 ∈ Q, p 1 , p 2 ∈ P with q 1 q 2 and p 1 p 2 . Let us consider the set of points

K + (q i , p i ) = {(x k , u k ) | (x k , u k ) (x qi , u pi )}, i = 1, 2.
Let us remark that (q 1 , p 1 ) (q 2 , p 2 ) implies that (x q1 , u p1 ) (x q2 , u p2 ), which in turn implies that we have K + (q 2 , p 2 ) ⊆ K + (q 1 , p 1 ). It follows from (3) that G(q 1 , p 1 ) G(q 2 , p 2 ). ( 5) Equation ( 4), with the same reasoning, leads to G(q 1 , p 1 ) G(q 2 , p 2 ). ( 6) Then, let q 1 ∈ F (q 1 , p 1 ), we have that q 1 q 1 , where q 1 is such that G(q 1 , p 1 ) ∈ X q 1 . Let q 2 be such that G(q 2 , p 2 ) ∈ X q 2 , then from (5), it follows that q 1 q 2 , and therefore that q 1 ∈ q 2 . Also, by definition of F , we get that q 2 ∈ F (q 2 , p 2 ). Along the same lines and using (6), we can prove that for all q 2 ∈ F (q 2 , p 2 ), there exists q 1 ∈ F (q 1 , p 1 ) such that q 1 ∈ q 2 . Hence, the transition system Σ D is monotone. 2 Knowing that Σ D is a monotone system, makes it possible to use efficient discrete controller synthesis techniques [START_REF] Saoud | Efficient synthesis for monotone transition systems and directed safety specifications[END_REF]).

We would like to end the section with a discussion regarding the computational complexity of the approach. The computational cost of computing the symbolic model is dominated by that of computing the functions G and G. Using equations ( 3) and ( 4), it appears that this can be done with a complexity in O(N × |Q| × |P |), i.e. linear with respect to the number of data samples and also with respect to the size of the partition.

NUMERICAL EXAMPLE

In this section, we report some numerical experiments to show the performance of our approach.

Let us consider a model with two vehicles moving in one lane on an infinite straight road. The leader is uncontrollable (vehicle 2) while the follower is controllable (vehicle 1). A discrete-time approximation of this model is given by equations:

d(k + 1) = d(k) + (v 1 (k) -v 2 (k))T 0 + w 1 (k) v 1 (k + 1) = v 1 (k) + α(u(k), v 1 (k))T 0 + w 2 (k) v 2 (k + 1) = v 2 (k) + w 3 (k) (7)
Here u(k) is the control input which is the torque applied to the wheels, and d(k) is the signed distance between the vehicles. w(k) = [w 1 (k), w 2 (k), w 3 (k)] are disturbances which belong to a rectangular set W with known lower and upper bounds. α is a nonlinear function given by α

(u, v) = u -M -1 (f 0 + f 1 v + f 2 v 2 ). The vector of parameters f = (f 0 , f 1 , f 2 ) ∈ R 3
+ describes road friction and vehicle aerodynamics whose numerical values are taken from [START_REF] Nilsson | Correct-by-construction adaptive cruise control: Two approaches[END_REF]): f 0 = 51 N, f 1 = 1.2567 N s/m, f 2 = 0.4342N s 2 /m 2 . For the rest of parameters, we chose M = 1370 kg, T 0 = 0.5 s.

The system can be made monotone by making the change of variable v 2 (k) = -v 2 (k). We only use this model to generate the random set of data points which then is used by the algorithm to calculate the symbolic abstraction.

We sample the component of the data points, input and disturbance from uniform distributions in the sets

X D = [-100, 0] × [10, 30] × [-22, -18] U D = [-3, 3] W D = [-0.1, 0.1] × [-0.1, 0.1] × [-0.1, 0.1]
We define the partitions for each coordinate of the state space as follow:

(-∞, a i 1 ), [a i 1 , 1 N i (a i N i -a i 1 ) + a i 1 ), . . . , [ N i -1 N i (a i N i -a i 1 ) + a i 1 , a i N i ), [a i N i , +∞) with [a 1 1 , a 1 N 1 ] × [a 2 1 , a 2 N 2 ] × [a 3 1 , a 3 N 3 ] ⊂ X D .
In a similar way we define the partitions for u:

(-∞, b 1 ), [b 1 , 1 M (b M -b 1 ) + b 1 ), . . . , [ M -1 M (b M -b 1 ) + b 1 , b M ), [b M , +∞) with [b 1 , b M ] ⊂ U D .
In order to measure the conservatism in our abstraction we introduce the performance criterion µ(D, Q, P ) = q,p vol(X q × U p ) × q ∈F (q,p) vol(X q )

q,p (vol

(X q × U p )) vol(W D ) ∀q ∈ Q, ∀q ∈ Q, ∀p ∈ P such that X q ⊆ X D , X q ⊆ X D , U p ⊆ U D
Intuitively, µ(D, Q, P ) measures the degree of conservatism of the abstraction. The denominator represents the volume of the unknown map F for the part of space where we can find an approximation, whereas the numerator represents the volume of our over-approximation. µ can take its value in the interval [1, ∞) and the the smaller its value is, the more accurate the abstraction.

Abstraction

For our first experiment, we fixed the partitions and computed the abstraction using different numbers of data points. For a given number of data points, we sampled five different random sets of data, and calculated the associated abstraction. For the partitions, we chose

N 1 = 30, N 2 = 30, N 3 = 20, M = 12, a 1 1 = -80, a 1 N 1 = -20, a 2 1 = 15, a 2 N 2 = 25, a 3 1 = -21.5, a 3 N 3 = -18.5, b 1 = 1 2 3 4 5 6 7 8
Number of data points -2.5, b M = 2.5. In Figure 2, the average times needed to compute the abstraction for the five sets of samples are plotted with respect to the number of data points. We can see that the time increases linearly with the number of data points. In Figure 3, the average values of µ are plotted with respect to the number of data points. We can see how increasing the number of data points can reduce the conservatism of our abstraction.

For the second experiment, we fixed the data set D, but changed the size of the partitions. For each instance, we increased all the values of N 1 , N 2 , N 3 , M and calculated the abstraction. We chose a number of data points N = 10 4 , and we fixed the values a

1 1 = -80, a 1 N 1 = -20, a 2 1 = 15, a 2 N 2 = 25, a 3 1 = -21.5, a 3 N 3 = -18.5, b 1 = -2.5, b M = 2.5.
In Figure 4, we can see how the time of execution increases linearly with the size of partitions. Figure 5 shows how increasing the number of partition elements also improves the abstraction.

Controller synthesis

We briefly show that our abstractions are suitable for controller synthesis purposes. Since the leader vehicle is The synthesis of symbolic controllers enforcing assumeguarantee contracts such as (8) has been considered in [START_REF] Saoud | Contractbased design of symbolic controllers for safety in distributed multiperiodic sampled-data systems[END_REF]) and we use similar techniques in the present case. The abstraction was computed based on a data set of 10 6 samples with the same partition as that reported in the first experiment above. The controllable set of the resulting symbolic controller is shown on Figure 6. In this paper, we presented a data-driven approach for computing abstractions of monotone systems. The approach is computationally tractable and the abstractions can be used to design controllers that are correct by design. In the future, we would like to develop data-driven approaches to estimate the bounds on the disturbances, which currently are assumed to be known. We also aim at developing similar approaches for other classes of systems, such as e.g. Lispchitz continuous systems. We would also like to extend our approach to refine symbolic models and controllers online while collecting samples of closed loop behaviors.

Fig. 1 .

 1 Fig. 1. Example of transitions of a monotone transition system.

Fig. 3 .

 3 Fig. 2. Average time of execution with respect to the number of data points

Fig. 4 .Fig. 5 .

 45 Fig. 4. Time of execution with respect to the number of elements of the partitions

Fig. 6 .

 6 Fig. 6. Controllable set of the symbolic controller enforcing the assume-guarantee contract (8)

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 725144).