G L Gattobigio

A Couvert

B Georgeot

D Guéry-Odelin

Exploring classically chaotic potentials with a matter wave quantum probe

Keywords: numbers: 03.75.Pp, 05.45.Mt, 37.25.+k, 41.85.Ew

We study an experimental setup in which a quantum probe, provided by a quasi-monomode guided atom laser, interacts with a static localized attractive potential whose characteristic parameters are tunable. In this system, classical mechanics predicts a transition from regular to chaotic behavior as a result of the coupling between the different degrees of freedom. Our experimental results display a clear signature of this transition. On the basis of extensive numerical simulations, we discuss the quantum versus classical physics predictions in this context. This system opens new possibilities for investigating quantum scattering, provides a new testing ground for classical and quantum chaos and enables to revisit the quantum-classical correspondence.

The realization and the manipulation of Bose-Einstein condensates (BECs) have led to many experiments on atom-atom interactions in the low energy regime: low partial waves elastic scattering has been analyzed to deduce accurate determination of collisional properties [1,2], the possibility of tailoring the scattering length using either magnetic [3] or optical Feshbach resonances [4] has been demonstrated, and the role of the dimensionality on atom-atom scattering processes through very anisotropic confinement has been studied [5,6].

In all these examples, the BEC is used to reveal the microscopic atom-atom interactions properties. BECs have also been placed into disordered [7,8] or noncommensurable bi-chromatic potentials [9]. Such experiments do not study a single matter wave scattering event. This is to be contrasted with experiments where a BEC bounces on a mirror [10]. However, in this case, the potential with which the BEC interacts is extremely simple.

In this article, we demonstrate the use of a guided matter-wave [11][12][13][14] to study, in a confined environment, the interaction of a matter wave with a static tunable localized attractive potential (LAP) associated with a complex dynamics. The transverse confinement of the guide can be well approximated by an harmonic confinement. The LAP superimposed to the guide breaks the harmonic character of the potential experienced by the atoms and couples the longitudinal and transverse degrees of freedom (see Fig. 1). We study this scattering problem theoretically from the weak to the strong coupling regimes for which the system exhibits classical chaos. This can be done by tuning the LAP power (strength) and position d. Finally, we show that these different regimes can be clearly identified in our experimental results.

Our experiment provides a direct access to the modulus squared of the wave function through absorption imaging. This is to be contrasted with experiments carried out with microwave in chaotic structure in which the measurement is performed at a well-defined position [15]. In mesoscopic physics, experiments exploit the conductance measurements to infer properties of the chaotic structure [16][17][18][19]. However, in these systems, the chaotic dynamics effects cannot be singled out from disorder and interactions effects. In contrast with pioneering experiments using cold atoms to investigate classical chaotic behavior in static optical billiards [20,21], in our experiment we prepared atoms in a nearly pure quantum state and we analyzed a single event rather than relying on average observables. Other studies on quantum chaos have been performed with cold atoms in the different context of 1D time-dependent lattice potentials [22][23][24].

The experimental setup that delivers guided matter waves has been described in detail elsewhere [13]. A rubidium-87 BEC is produced in a crossed dipole trap using two focused laser beams with a wavelength of 1070 nm and waists w G ≈40 µm. Atoms are prepared in the |F = 1, m F = 0 internal state through a spin distillation process implemented during the evaporative cooling stage [12]. The quantum probe is realized by magnetic or all-optical outcoupling from the trap (located at A in Fig. 1) [13]. Atom-atom interactions turn out to be negligible in the guided atom lasers produced in this manner [13]. Experimentally, the propagating matter wave is characterized by three quantities: the atom flux (a few 10 5 atoms/s), the mean velocity v (tunable from 10 to 30 mm.s -1 ) and the population of the transverse modes. Experimental parameters are chosen to have the mean excitation number n below one, which corresponds to a ground state population larger than 50% [13]. Such a probe therefore explores the phase space with a transverse resolution close to 2 . The horizontal beam used for the dipole trap provides the guiding potential which can be approximated by a harmonic transverse confinement with a measured frequency ω ⊥ /2π ranging from 170 Hz to 200 Hz (guide power, P G ≃ 100-150 mW). The LAP originates from the same laser source as for the crossed dipole trap and is slightly detuned to wash out interferences. It has a waist of w LAP ≈ 60 µm and an adjustable power, P LAP , from 0 to 2W. It intersects the guide beam at a distance AB = 700 µm from the trap (A) along the propagation direction (see Fig. 1).

The interaction of the quantum probe with a wellcentered LAP (d = 0) has been investigated theoretiin [25]. To analyze the experimental results of the scattering problem for d = 0 it is instructive, as a first guideline, to use classical mechanics.

Consider that atoms are classical point-like particles delivered at the trap location (A) and at the bottom of the guide (see Fig. 1). As a result of the longitudinal curvature of the guide beam, atoms are accelerated towards point B. Close to B, atoms feel the combination of the guide and the attractive potential. When the LAP is well centered (d = 0), it does not generate a transverse oscillation downstream. However, when the transverse symmetry is broken (d = 0), the LAP has a dramatic effect. In the weak coupling regime, the beam goes downstream and exhibits a transverse dipole oscillation (Fig. 2a). A concomitant reduction of the longitudinal velocity occurs due to energy conservation. This output channel is stable against small changes of d/w LAP and P LAP /P G . In the strong coupling regime, the output channel can be either upstream (Fig. 2b) or downstream depending on the different parameters of the LAP (strength and displacement) and turns out to be extremely sensitive to small parameters changes (Fig. 2d).

The complex trajectories that result from the strong coupling between longitudinal and transverse degrees of freedom illustrate the chaotic nature of the potential in this regime (see Fig. 2c) where the Lyapunov exponents can be shown to be positive and the distribution of trapping times is found to exhibit an exponential decay (data not shown). Figure 2d summarizes the variation of the asymptotic transverse amplitude as a function of the ratio P LAP /P G and the displacement d of the LAP. It exhibits a region with strong variations of output channels with parameters, that we identify as the chaotic zone. To further characterize this zone, we have verified that the set of boundaries between left and right output channels as a function of d displays a fractal structure with a fractal dimension ξ. For a wide variety of fixed ratio P LAP /P G in the chaotic zone, we find an exponent ξ ≈ 0.9 revealing a strong fractality of the basin boundaries, usually associated with chaotic scattering [START_REF] Ott | Chaos in dynamical systems[END_REF].

However, the incident guided atom laser is not a classical probe since it is in the ground state of the transverse confinement [13]. In addition, the potential experienced by the atoms is not harmonic and the dynamics cannot thus be exactly mapped onto the classical one.

To explore the quantum-classical correspondence, we have developed a full 2D numerical simulation of the wave packet that represents the incident elongated matter wave. The extremely large anisotropy of the system with a longitudinal length of a few mm and a transverse size initially of a few hundreds of nm and the long interaction time inherent to the complex dynamics makes this kind of calculation memory size and time consuming, and even unrealistic in 3D with the currently available supercomputers. We have thus compared the dynamics of 2D quantum wavepackets with simulations of classical probability densities of same initial position and momentum distributions. Parameters are chosen in close similarity with the experimental conditions. In the weak coupling regime, the scattering is regular, the quantum wave packet keeps its structure and performs coherent oscilla- tions downward the LAP (Fig. 3a). In this regime, classical and quantum predictions are indistinguishable. In the strong coupling regime for which the classical counterpart is chaotic, the wave packet loses its filamentary structure. In this regime, the numerically simulated wavepacket follows the central trajectory for short time (see Fig. 3b).

For larger time, it loses its structure and classical and quantum dynamics can become different (see Fig. 3c and Fig. 3d). The quantum and classical results differ significantly when the wave function overlaps with itself in the LAP region and therefore exhibits some interference patterns (see Fig. 3d).

The scattering problem depends on many parameters: those that characterize the incident wave (velocity, number of occupied transverse modes, longitudinal coherence length, ...) and those that concern the confinement (powers, waists and distance d). Experimentally, we have investigated this large parameter space by varying the ratio P LAP /P G which in turn slightly modifies the displacement d (see below). The experiment duration, τ ∼ 300 ms, is chosen so that typically half the guided atom laser has interacted with the LAP in the weak coupling regime. This duration was chosen in order to observe the asymptotic regime in the weak coupling regime and the complex dynamics associated with the strong coupling regime. Absorption images are taken after time-of-flight once the confinement potentials are switched off. In this way, the ratio between the final, ∆, and initial, δ, transverse sizes can readily be extracted. Examples of experimental results are shown in Fig. 4. The two regimes (regular and chaotic behaviors) are clearly observed when the power is increased (see Fig. 4a and, Fig. 4b andc respectively). The transverse oscillation observed after the LAP in the weak interacting regime corresponds to a coherent transverse dipole oscillation (Fig. 4a) in close similarity with the corresponding quantum simulation (Fig. 3a). The apparent decrease of the transverse size of the helix motion is due to the guide anisotropy and the projection of the trajectory along the imaging axis.

For given parameters, the displacement d of the LAP is reproducible from shot to shot. However, d varies slightly with the LAP power. This latter effect is due to thermal effect generated by the RF power in the acousto-optic modulator that control the LAP power. In the parameter space of Fig. 2d, the experimental results thus do not follow a vertical line. This is clearly illustrated in Fig 4d where we have represented ∆/δ, in the weak coupling regime, for both the experimental points, and 3D (classical) simulations that assume a fixed displacement d = 4 µm. Here, we use the fact that classical and quantum predictions coincide in this regime. The observed saturation for P LAP > 10P G suggests that in this parameter regime d is approximately constant d ≃ 4 µm ± 0.3 µm. For P LAP < 10P G , the initial curvature is opposite to the simulated one revealing a change of d as P LAP is increased. Using a linear increase of d up to P LAP = 10P G towards the value d = 4 µm and a constant value afterwards gives a very good agreement with the experimental data (see Fig. 4d).

Interestingly enough, the experiments turn out to be robust with respect to the velocity dispersion of the incoming matter wave. This effect, confirmed by numerical simulations, originates from the fact that the larger available energy for larger incident velocity is compensated for by the smaller interaction time with the LAP.

Experimentally, we observe two types of images in the chaotic regime: those for which the trajectory can be well-identified (Fig. 4b) and those with a more diffuse cloud (Fig. 4c). Numerical simulations (see Fig. 3b) suggest that the former can be attributed to situations for which the chaotic dynamics did not have enough time to fully develop. The latter may result from either the ergodic mixture as in Fig. 3c that occurs in the chaotic regime for sufficiently long time or situations in which the relatively high local atomic density favors elastic collisions and thus the possibility of exploring a larger phase space zone with various output channels.

In Fig. 5, we assess the reproducibility of the experiment by using different runs with same initial conditions, where small initial differences are inevitably present. Examples of the different experimental images are shown in the top inset. We define the quantity σ i,j (z) as the difference of the mean transverse size of the wave packet at coordinate z downward the LAP (located at z B ) for two different runs i and j taken in the same conditions. We note σ(z) = N -1 r σ (z) where N r the number of runs (N r = 5 in our experiment) and σ is σ(z) averaged over z > z B . In Fig. 5 (top), this latter quantity is represented vs P LAP /P G . When the system enters the chaotic (strong coupling) regime a large dispersion in the transverse size is observed as a result of the high sensitivity of the system to parameter changes and initial conditions. The transition to chaos becomes visible for P LAP /P G ∼ 8, compatible with values of d ≈ 7.5 µm (d/w LAP = 0.125) (see Fig. 2d).

To make more visible the chaotic regime in the experiments, in Fig. 5 (bottom) we show the average distance σ(z) . As long as the wavepacket moves along the z coordinate, this coordinate gives the time evolution of the quantum wave function. Unavoidably small differences were present in the initial conditions. In the regular regime, these differences remain approximately of the same order of magnitude, while in the chaotic regime they display a sharp exponential amplification when the wavepacket enters the chaotic potential. This is reminiscent of what happens for classical trajectories (see Fig. 2c). Strictly speaking, σ(z) does not give a direct access to the Lyapunov exponent since more runs would be required and we have only a projection of the density probability in position space (instead of the phase space trajectory in classical physics). However, according to our numerical studies, σ(z) in many cases can be used to give a lower bound to the Lyapunov exponent.

In conclusion, we have explored the effects of a LAP on a guided atom laser. A classical analysis enables to identify different regimes, from regular to chaotic and un- stable dynamics. These different regimes are also visible in the quantum simulations as well as in the experiments. Differences between classical and quantum dynamics are especially noticeable in the non-regular regime. This system enables to investigate both classical and quantum chaos in a well-controlled environment. Two situations should be especially adapted to magnified quantum features: (i) the 2D geometry which favors interferences of the wavefunction with itself and (ii) situations for which LAP size is on the order of the transverse size of the initial wave function. In this latter regime quantum diffraction and/or tunneling effects should play an important role. The experimental exploration of these two regimes provides a natural continuation of our work. Another direction for future work includes the study of the nonlinear influence of atom-atom interaction and in particular matter wave solitons [START_REF] Khaykovich | [END_REF]28]. We think our approach paves the way for the use of guided matter waves to study highly non-equilibrium dynamics and to use such systems as quantum simulators for complex quantum scattering problems.

(OMASYC project) and the Institut Universitaire de France. We thank CalMip for the use of their super-computers.

FIG. 1 :

 1 FIG. 1: (color online). Sketch of the experimental setup. Atoms outcoupled from a BEC at point A are guided by the horizontal beam towards point B. In the vicinity of B, atoms also experience the attractive force exerted by the LAP (see text) whose center is at a distance d from the guide axis.

6 FIG. 2 :

 62 FIG. 2: (color online). (a) Typical trajectory in the weak coupling regime (see text). The interaction with the LAP yields a downstream propagation with a transverse oscillation. (b) Example of trajectory in the strong coupling regime. The complex dynamics may lead to the reflexion of the particle. In both (a) and (b) the LAP potential is represented by the blue shaded area. (c) Example of exponential separation Σ of nearby trajectories in phase space revealing the chaotic nature of the potential. The red dotted curve is the best exponential fit. Average over 100 couples of trajectories close to the central one is taken. (d) Amplitude of the transverse excitation as a function of the ratio PLAP/PG and the displacement d of the LAP: transmitted (red colors) and reflected (blue colors) particles. A chaotic zone with a complex dynamics can be clearly identified.

FIG. 3 :

 3 FIG. 3: (Color online). (a) (b) (c) Color density plots of the wave function from 2D numerical simulations. Color scale is not linear. Point B is the same as in Fig. 1. Wave packets are prepared in the ground state of the initial trap, and released in the guide at time t = 0 (wG = 45 µm and ω ⊥ = 2π × 203 Hz) at 700 µm from the LAP (wLAP = 41 µm). The incident mean velocity of the wave packet is v ∼ 10 mm.s -1 close to B. (a) d = 2 µm (d/wLAP = 0.05), PLAP = 0.5 W (PLAP/PG = 2.78), propagation time t = 120 ms, (b): d = 15 µm (d/wLAP = 0.37) and PLAP = 1W (PLAP/PG = 5.56), t = 120 ms, (c) d = 20 µm (d/wLAP = 0.5), PLAP = 1 W (PLAP/PG = 5.56), t = 160 ms, (d) projected distribution (c) on the z axis of the wave function dx|ψ(x, z)| 2 (black) and the initially identical classical distribution (red/gray); classical and quantum curves are identical for cases (a) and (b) (data not shown).

FIG. 4 :

 4 FIG. 4: (color online). Interaction between a guided atom laser and a LAP results). Absorption image (a) is the weak coupling regime (b) (c) in the strong coupling regime after a 15 ms time-of-flight. (d) Measurement of the displacement d in the weak coupling regime: the ratio ∆/δ of the transverse size before and after the interaction is plotted vs PLAP/PG (v ≃ 25 mm.s -1 , and ω ⊥ /2π ≃ 200 Hz). The dot-dashed curve corresponds to 3D (classical) simulation with d = 4 µm (d/wLAP ≈ 0.07). Most experimental points (disks) are inside the pink area delimited by the simulated curves (gray lines) for 3.7 and 4.3 µm.

FIG. 5 :

 5 FIG.5:(color online) Experimental difference between quantum evolutions for different runs with same initial conditions (same parameters as in Fig.4). (Top) σ (see text) vs PLAP/PG. The onset of the chaotic regime is clearly visible on both the large value of σ and its error bar. Inset: Examples of these different runs in the weak (left) and strong (right) coupling regimes. (Bottom) σ(z) (see text) vs z. Solid line is PLAP = 0.06 W (PLAP/PG = 0.6), dotted line is PLAP = 0.6 W (PLAP/PG = 6), dashed line is PLAP = 1 W (PLAP/PG = 10). We observe the amplification of small initial deviations when the beam enters the chaotic zone (z ∼ zB ∼ 0) as in Fig.2c.

We thank R. Mathevet and T. Lahaye for useful comments. We acknowledge financial support from the Agence Nationale de la Recherche (GALOP project), the Région Midi-Pyrénées, the university Paul Sabatier


Sebastien Destercke
email: sebastien.destercke@cirad.fr

Valerie Guillard
email: guillard@univ-montp2.fr

Interval analysis on non-linear monotonic systems as an efficient tool to optimise fresh food packaging

Keywords: Interval analysis, robust design, fuzzy sets, sensitivity analysis, uncertainty analysis

When few data or information are available, the validity of studies performing uncertainty analysis or robust design optimisation (i.e., parameter optimisation under uncertainty) with a probabilistic approach is questionable. This is particularly true in some agronomical fields, where parameter and variable uncertainties are often quantified by a handful of measurements or by expert opinions. In this paper, we propose a simple alternative approach based on interval analysis, which avoids the pitfalls of a classical probabilistic approach. We propose simple methods to achieve uncertainty propagation, parameter optimisation and sensitivity analysis in cases where the model satisfies some monotonic properties. As a real-world case study, we interest ourselves to the application developed in our laboratory that has motivated the present work, that is the design of sustainable food packaging preserving fresh fruits and vegetables as long as possible.

Introduction

There are many sources of uncertainties in life science and in agronomy, the main reasons for it being the high variability of living organism and the error of measurement devices. Also, the number of available samples for a given experiment may be limited (sometimes even reduced to one sample), due to cost or practical limitations. In such situations, it can be hard to determine a meaningful probabilistic model of the parameters, let alone a joint probabilistic model over all parameters. In other situations, uncertainty around a parameter or a constant can be described by expert opinions, and whether these opinions can be faithfully translated by single probabilities is questionable (see, e.g., [START_REF] Walley | Statistical reasoning with imprecise Probabilities[END_REF]Sec. 4] or [START_REF] Sandri | Elicitation, assessment and pooling of expert judgments using possibility theory[END_REF]).

In such situations, it may be better to use interval modelling and interval analysis to perform uncertainty studies, simply because determining intervals requires less data and knowledge. Also, using interval analysis amounts to make no assumptions about parameter dependencies.

In this paper, we consider dynamical non-linear models describing the evolution of variables, with the aim to optimise some of the parameter values w.r.t. some given objective. In such systems, the values of initial conditions, non-modifiable parameter values or even of the objectives may be ill-known. It is then desirable to perform some uncertainty analysis to achieve robust design. Performing such analysis with classical probabilistic methods [START_REF] Hertog | The impact of biological variation on postharvest behaviour: The case of dynamic temperature conditions[END_REF][START_REF] Bedford | Probabilistic Risk Analysis. Foundations and Methods[END_REF] usually requires to:

• specify the distribution of each input variable,

• specify the dependency structures between input variables,

• perform (costly) numerical analysis to evaluate the output uncertainty.

Meeting such requirements necessitates an important amount of information and data. It also involves the use of techniques having a high computational cost. In practice, when not enough information is available, distribution shapes (e.g., normality) and dependence assumptions (e.g., independence between all variables) are often chosen accordingly to some practical criterion rather than to available information. However, the validity of such choices, not confirmed by experiments or available knowledge, may be questioned, as well as the validity of subsequent analysis results [START_REF] Ferson | Different methods are needed to propagate ignorance and variability[END_REF]. They then provide overly precise and misleading conclusions, which may in turn lead to unwarranted and non-robust design choices.

When few data are available, an alternative is to use interval analysis [START_REF] Jaulin | Applied Interval Analysis[END_REF] to perform the uncertainty analysis, that is to consider that only the bounds in which each parameter may vary are known (an information that is often available). Such an analysis comes down to consider that:

• variable and parameter distributions are unknown (up to their bounds),

• dependence structure between variables and parameters is unknown.

Compared to probabilistic analysis, interval analysis can therefore be seen as a conservative analysis, in the sense that it does not make any additional hypothesis with regard to the available information, and possibly ignores some of the available information. However, in scientific analysis as well as in robust design, it is safer to use such conservative methods than to make unsupported assumptions. Also note that, when function f has some monotonic properties [START_REF] Singer | Bounding the solutions of parameter dependent nonlinear ordinary differential equations[END_REF][START_REF] Fortin | Gradual numbers and their application to fuzzy interval analysis[END_REF] (a common case in life sciences and other domains [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear monotone systems[END_REF][START_REF] Baudrit | Joint treatment of imprecision and variability in food engineering: Application to cheese mass loss during ripening[END_REF], where simple models are often encountered), performing interval analysis may require very few computations compared to, say, probabilistic Monte-Carlo analysis. When these monotonic properties are not satisfied, performing interval analysis requires more complex techniques [START_REF] Ramdani | A hybrid bounding method for computing an over-approximation for the reachable set of uncertain nonlinear systems[END_REF] with increased computational complexity (the computational cost then becoming comparable to the one of probabilistic methods). However, whether the monotonic properties are satisfied have no impact on the results conservativeness of interval analysis.

In this paper, we introduce a set of methods to achieve uncertainty propagation, parameter optimisation and sensitivity analysis on monotonic dynamical systems when uncertainty is described by intervals. Notations and general problem formulations are introduced in Section 2. Section 3 then provides details about the method themselves.

Finally, we illustrate the method on the real-world case study that have motivated the present work and that is currently treated in our laboratory [START_REF] Charles | Modeling of active modified atmosphere packaging of endives exposed to several postharvest temperatures[END_REF] (see website http://www.tailorpack.com/). It concerns the design of sustainable fresh food packaging, with the objective to preserve food from decay as long as possible. We use our method to, first, perform uncertainty analysis of a model describing gas exchanges between the packaging atmosphere and exterior atmosphere and, second, optimise oxygen and carbon dioxide permeances of packaging materials for a given fruit or vegetable, here chicory. The whole case study is described in Section 4.

Problem setting and notations

Vectors of model variables and parameters will be denoted by bold letters (x, p, . . .), while specific values of these vectors will de denoted by non-bold letters (x, p, . . .). As the models considered in this paper are multivariate, indexed letters x j will denote the jth element of vector x. Sets will be denoted by calligraphic letters (X , P, . . .). The real line will be denoted by R.

We consider a dynamical non-linear model ẋ = f (x, p E , p D ,t) with f : R m+n+1 → R m a time-dependent function that describes the evolution of m state variables x ∈ R m , i.e. x = (x 1 , . . . , x m ). The evolution of these variables depends on n E environmental parameters p E (whose values cannot be controlled) such as temperature or external pressure and n D design parameters p D (which values can be modified) such as mechanical or chemical properties of synthetic compounds. They form a vector p = (p D , p E ) of n parameters p = {p 1 , . . . ,

p n } ∈ R n .
Given some values p ∈ R n of the parameters and some initial conditions x(0) ∈ R m of the state variables, the solution of the system represented by f at time t is x(t) = (x 1 (t), . . . , x m (t)) where x i (t) describes the state of the ith state variable x i at time t.

When initial conditions x(0) and environmental parameter values p E are known, a classical design problem [START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF] consists in identifying the values p D ∈ R n D of the design parameters p D so that the solutions x(t) are as close as possible to a given objective x(t) on the state variables. However, both the exact values of initial conditions x(0), of parameters p or of the objective x(t) to reach are seldom known with certainty. In the next section, we detail how the problem can be treated when these values become interval and when the model satisfies some monotonic properties.

Interval analysis and design optimisation

In this section, we start by giving some refreshers on classical interval analysis and interval analysis on dynamical monotonic models, illustrating them on a simple example. We then detail our proposed optimisation and sensitivity analysis methods for such models.

Basics of interval analysis

In the computational literature, interval analysis was first developed to take account of numerical errors [START_REF] Moore | Methods and applications of Interval Analysis[END_REF]. However, interval analysis is now mostly used to perform robustness analysis in applications (robotics [START_REF] Jaulin | Guaranteed robust nonlinear estimation with application to robot localization[END_REF], chemical, biological, . . . ) where variable values are imprecisely known [START_REF] Jaulin | Applied Interval Analysis[END_REF].

A real interval [x] := [x -, x + ] of a variable X is a connected and closed subset of R. The set of real intervals is usually denoted IR. An interval vector [x] over R n (also called box) is the Cartesian product of n intervals. The classical problem of interval analysis consists in replacing, in a given function

Y = f (X 1 , . . . , X n ) from R n to R m , the point values x = (x 1 , . . . , x n ) of variables X = {X 1 , . . . , X n } by intervals [x] = ([x 1 ], . . . , [x n ]
) and to compute the range

f ([x]) = { f (x)|x ∈ [x]}. (1)
[x] Usually, f ([x]) is not a box over R m , but a complicated subset of it (i.e., it cannot be expressed as a Cartesian product of intervals). Rather than computing the exact propagation f ([x]), one can use guaranteed approximation techniques and compute a box

f ([x]) [ f ]([x])
[ f ]([x]
) over R m that will be an inclusion function, i.e.,

f ([x]) ⊆ [ f ]([x]). (2)
This comes down to compute an outer approximating interval for every dimension of f . The notion is illustrated in Figure 1. Among possible techniques to compute such inclusion functions is interval arithmetics [START_REF] Moore | Methods and applications of Interval Analysis[END_REF], where classical arithmetic operations {+, ×, /, -} are replaced by their interval equivalent. Depending on the technique and on the characteristic of the model f , [ f ]([x]) will be more or less close to f ([x]). In the next section, we recall such techniques that can be applied to the particular dynamical models we are interested in.

Interval analysis to propagate uncertainties in dynamical systems

Let us consider the more complex problem of evaluating the solution of the dynamical system Ẋ = f (X, P E , P D ,t) with f : R m → R m . In this model, the evolution of each state variable X i , i = 1, . . . , m is described by an ordinary differential equation (ODE) such that Ẋi = f i (X, P E , P D ,t). The solutions of such systems are m functions x i (t), i = 1, . . . , m. Here values of function x i (t) are discretised and computed for a finite number of time steps. We denote by T = [0,t] the time domain of the model and assume that each x i (t) is computed for T different times values t k , k = 1, . . . , T (with t kt k-1 being constant for every k).

For such models, the classical problem of interval analysis is formulated as follows: given m initial conditions

x i (0) ∈ [x i ](0) := [x i -, x i + ](0), i = 1, . . . ,
:= × m i=1 [x i ](0)× n j=1 [p j ].
Interval analysis on Ẋ = f (X, P E , P D ,t) then consists in finding, among all configurations in X , those reaching the bounds of x(t k ) for each time step t k , k = 1, . . . , T . As in the non-dynamical case (see Equation ( 2) and Figure 1), [x](t) can have a very complex shape, hence a slightly easier problem is to find, for each variable

X i , i = 1, . . . , m, some bounds [x i ](t) such that [x](t) ⊆ × m i=1 [x i ](t)
and then to consider the resulting boxes × m i=1 [x i ](t) as the solutions. This is the problem we will consider here.

Among the set X of configurations, denote by

H = × n i=1 {x i -, x i + } × m j=1
{p j -, p j + } the set of extreme ones, that is the set of configurations whose elements are bounds of the intervals. H therefore contains only combinations of such bounds, i.e., 2 (n+m) elements, while X contains an infinity of them. To illustrate the approach used in this paper, we will use the following simple example:

Example 1. Consider the very simple system from R → R Ẋ1 = f 1 (P 1 ) = P 1 ,
where p 1 ∈ [p 1 ] and x 1 (0) ∈ [x 1 ](0) are positive numbers, and the (analytical) solution of the system is x 1 (t) = P 1 • t + x 1 (0). Here, both X = X 1 and P = P 1 are reduced to a single element each. Our knowledge about p 1 and initial condition is given by the intervals p 1 ∈ [0.5, 1] and x 1 (0) ∈ [START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF][START_REF] Baudrit | Joint treatment of imprecision and variability in food engineering: Application to cheese mass loss during ripening[END_REF]. The sets of all possible combinations of [p 1 ] and [x 1 ](0), together with their extreme combinations are illustrated in Figure 2. In this example, the set X = [1, 3] × [0.5, 1], while the set H = {(1, 0.5), (1, 1), (3, 0.5), (3, 1)} is reduced to four points.

As for interval analysis on classical functions [START_REF] Fortin | Gradual numbers and their application to fuzzy interval analysis[END_REF], there exist specific cases [START_REF] Delanoue | A new method for integrating ode based on monotonicity[END_REF][START_REF] Ramdani | Computing reachable sets for uncertain nonlinear monotone systems[END_REF][START_REF] Singer | Bounding the solutions of parameter dependent nonlinear ordinary differential equations[END_REF] where the bounds of [x i ](t) on each dimension correspond to solutions of the system Ẋ = f (X, P E , P D ,t) corresponding to specific points in the set H of extreme configurations, thus reducing the number of necessary computations to find them. Before recalling what are these cases, let us first define the notion of monotonic for dynamical system Ẋ = f (X, P E , P D ,t). Definition 1. f i is said to be dynamically increasing w.r.t. variable X j , j = i or parameter and it said to be dynamically decreasing w.r.t. variable X j , j = i or parameter

P k if ∂ f i ∂ X j ≥ 0 or ∂ f i ∂ P k ≥ 0 (3) X 1 P 1 1 0.5 3 1 X : All possible couples {p 1 , x 1 (0)} H : Extreme configurations
P k if ∂ f i ∂ X j ≤ 0 or ∂ f i ∂ P k ≤ 0 (4)
For a given function f i (or, equivalently, for a given variable X i ), let us denote by P ,i , X ,i and P ,i , X ,i the set of parameters for which f i is dynamically increasing and decreasing, respectively. Note that they are disjoint, i.e. P ,i ∩ P ,i = / 0 and X ,i ∩ X ,i = / 0. We also assume from now on that these sets form partitions of the parameters and variables, that is for a given i, P ,i ∪ P ,i = P and X ,i ∪ X ,i = X. When a function f i is dynamically increasing or decreasing in each variable and parameter, the following proposition [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear monotone systems[END_REF] tells us how the upper and lower envelope of [x i ](t) at any time step t can be obtained. Proposition 1. For any t, the bounds of [x i ](t) can be computed as follows:

• x + i (t) is the solution of the system Ẋ = f (X, P E , P D ,t) with the configuration

H i ∈ H such that H i =        P = p + for all P ∈ P ,i P = p - for all P ∈ P ,i X(0) = x + (0) for all X ∈ X ,i X(0) = x -(0) for all X ∈ X ,i (5) • x - i (t)
is the solution of the system Ẋ = f (X, P E , P D ,t) with the configuration x - 1 (t) : This means that if the monotonic properties of f i are known, then one can solve the system with configuration H i to get x + i (t), and with H i to get x - i (t). This means that to retrieve the bounds of [x i ](t) for any variable X i , i = 1, . . . , m, the system Ẋ = f (X, P E , P D ,t) has to be solved twice with usual techniques not involving intervals. For such systems, propagating uncertainties modelled by intervals then becomes much easier than propagating probabilistic uncertainties (it comes down to solve at most 2•m different systems). Let us continue our example.

H i ∈ H such that H i =        P = p + for all P ∈ P ,i P = p - for all P ∈ P ,i X(0) = x + (0) for all X ∈ X ,i X(0) = x -(0) for all X ∈ X ,i (6)
p 1 = 0.5, x 1 (0) = 1 x + 1 (t) : p 1 = 1, x 1 (0) = 3
Example 2. Consider the system of Example 1. We have (assuming a positive P 1 ) that ∂ f 1 ∂ P 1 ≥ 0, hence f 1 is dynamically increasing w.r.t. P 1 and X 1 . Hence P ,i = P 1 , P ,i = / 0, X ,i = X 1 and X ,i = / 0. Hence the upper envelope x + 1 (t) is given by the configuration H 1 = {p + 1 , x + (0)}, and the lower one by x - 1 (t) by

H 1 = {p - 1 , x -(0)}.

Parameters optimisation

We now consider the problem of searching optimal values p D of design parameters P D when initial conditions and environmental parameters are intervalvalued, and when objectives (constraints) x(t) on the state variables are intervalvalued as well. We consider that, for each variable X i , the objectives x i : T → IR can be given over some subset T ⊆ T of the whole time domain (for instance, constraints on values may only be specified for the steady state only, that is after a time t * such that T = [t * ,t], . . . ). Note that objectives can be intervals, as their exact values can themselves be uncertainly known (or there may be many values that appear as optimal in the context). This way of formulating an optimisation problem is not usual, even in interval-analysis literature [START_REF] Raissi | Set membership state and parameter estimation for systems described by nonlinear differential equations[END_REF], where the objective function is usually precisely valued (e.g., corresponds to precise outputs). Again, such a problem may be in general difficult to solve. We propose a general way of formulating the problem, before proposing an easy-to-apply algorithm in the case of monotonic systems (i.e., systems satisfying Definition 1).

Consider some pre-defined objectives [ x i ](t) as well as some initial conditions [x i ](0), i = 1, . . . , m and some interval-valued uncertainty [p] for every environmental parameter P ∈ P E . Consider some (precise) value p D for design parameters P D . Design parameter values p D are said to form a guaranteed (resp. possible) solution if the solution of the system Ẋ = f (X, P E , P D ,t) with these values p D (and with given intervals on P E and X initial conditions) is such that

[x i ](t) ⊆ [ x i ](t) (resp. [x i ](t) ∩ [ x i ](t) = / 0) for every t ∈ T .
We denote the set of guaranteed solutions by G D and the set of possible solutions by S D . Both these sets may be empty, but we have the inclusion relationship G D ⊆ S D . A guaranteed solution is such that, despite of interval uncertainties, we are certain that with the given design parameter values, the true answer lies within the objective bounds, while a possible solution is such that, with the given design parameter values, the true answer may or may not lie within the objective bounds. Solutions that are totally outside the objective bounds are said non-admissible.

Example 3. Consider again the model of Example 2, except that this time parameter P 1 is considered as a design parameter that can be tuned through some control process (e.g., a speed of reaction controlled by some catalyser). Objective

x 1 (t) on the variable X 1 is specified for T = {10} and is such that [ x 1 ](10) = [START_REF] Dubois | The three semantics of fuzzy sets[END_REF][START_REF] Fortin | Gradual numbers and their application to fuzzy interval analysis[END_REF]. Figure 4 illustrates the notions of guaranteed, possible and non-admissible solutions for this case. For example, when P 1 = 0.75, we have x - 1 (t)

x i (0) ∈ [x i ](0), i ∈ {1, . . . , m}, p j ∈ [p j ] ∀P j ∈ P E ∀t ∈ T , i ∈ {1, . . . , m}, x + i (t) ≤ x i (t) + and x - i (t) ≥ x i (t) -, X 1 t x(t)
x + 1 (t) x - 1 (t)

P 1 = 0.75
x + 1 (t)

P 1 = 0.9 x - 1 (t)

x + 1 (t)

P 1 = 1.1
x i (0) ∈ [x i ](0), i ∈ {1, . . . , m}, p j ∈ [p j ] ∀P j ∈ P E ∀t ∈ T , i ∈ {1, . . . , m}, x - i (t) ≤ x i (t) + or x + i (t) ≥ x i (t) -.
Again, solving exactly such problems is in general difficult. However, when the model1 Ẋ = f (X, P E , p D ,t) is dynamically monotonic (either increasing or decreasing) in each variable and each environmental parameter, we propose a simple heuristic method to identify [G P ] and [S P ]. Algorithm 1 suggests some means to find them when each component of Ẋ = f (X, P E , p D ,t) is dynamically monotonic w.r.t. a given design parameter P. In this algorithm, H i P E ,X and H i P E ,X are the configurations of Proposition 1 reduced to environmental parameters and initial conditions, the values of design parameters being left unspecified.

The algorithm simply uses the known monotonic properties of the model to compute boundary values. For instance, consider Lines 2-3 and the case P ∈ Algorithm 1: Approximations of [G P ] and [S P ] of a design parameter P 10 G P = [max j=1,...,m g P, j -, min j=1,...,m g P, j + ] with G P = / 0 if max j=1,...,m g P, j

Input: model Ẋ = f (X, P E , p D ,
-≥ min j=1,...,m g P, j + ;

11 S P = [max j=1,...,m s P, j -, min j=1,...,m s P, j + ] with S P = / 0 if max j=1,...,m s P, j -≥ min j=1,...,m s P, j + ; P ,i , then considering the configuration of environmental parameters and initial conditions giving x + i (t) is the most constraining configuration we can have w.r.t. the objective upper bound

x i (t) + and such that there is still a chance that p is a guaranteed solution. We still have [S P ] ⊆ [G P ]. Note that Algorithm 1 works only when f satisfies Definition 1. In more complex cases, approximating sets G D and S D can be done by sampling different values of design parameters and then performing interval analysis (with propagation methods adapted to more general models [START_REF] Ramdani | A hybrid bounding method for computing an over-approximation for the reachable set of uncertain nonlinear systems[END_REF]) to check whether the sampled values are (guaranteed or possible) solutions.

Example 4. Consider again the model of Example 2, with the constraints of Example 3. Running Algorithm 1 to identify [G P 1 ] and [S P 1 ] gives the following solutions (note that here, initial conditions and objectives [ x 1 ](10) provides each time two points, hence totally determining P 1 in the equation x 1 (t) = P 1 •t +x 1 (0)):

• Line 2 of Alg. 1, i=1: Identification with

x + 1 (10) = 11 and H i P E ,X = {x 1 (0) = 3}, giving p = 0.8 as solution and g P 1 ,1 + = 0.8 (since P 1 ∈ P ,1 ).

• Line 4 of Alg. 1, i=1: Identification with x - 1 (10) = 8 and H i P E ,X = {x 1 (0) = 1}, giving p = 0.7 as solution and g P 1 ,1 -= 0.7 The obtained sets [S P ] and [G P ] for each design parameter P ∈ P D can then be transformed into fuzzy sets [START_REF] Zadeh | The concept of a linguistic variable and its application to approximate reasoning-i[END_REF]. Indeed, while guaranteed solutions all provide the same satisfaction to the designer or the decision maker (they all ensure that the true solution is within the objective boundaries), possible solutions can be seen as having gradual satisfaction degrees, as some of them will have a more significant overlap with the objective bounds than others.

•
First recall that a fuzzy set µ : V → [0, 1] is a mapping from a space V (here the real line) to the unit interval, where µ(v) is the membership degree of element 0 0.2 0.4 0.6 0.

µ ã(x) =    x-a 1 a 2 -a 1 if a 1 ≤ x ≤ a 2 1 if a 2 ≤ x ≤ a 3 x-a 4 a 3 -a 4 if a 3 ≤ x ≤ a 4 (7)
In our case, the fuzzy set degree expresses some satisfaction degree [START_REF] Dubois | The three semantics of fuzzy sets[END_REF] provided by a parameter value w.r.t. an interval-valued objective. We propose, for a design parameter P, to build the trapezoidal fuzzy number µ P such that P = {s P -, g P -, g P + , s P + } if G p i = / 0 and such that P = {s P -, (s P -+s P + ) /2, (s P -+s P + ) /2, s P + } otherwise. Figure 4 illustrates the fuzzy number obtained in Example 4.

Sensitivity analysis

When performing an uncertainty analysis, it is usual to perform it along with a sensitivity analysis [START_REF] Helton | Sensitivity analysis in conjunction with evidence theory representations of epistemic uncertainty[END_REF][START_REF] Saltelli | Global Sensitivity Analysis: The Primer[END_REF][START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. Sensitivity analysis consists in searching what input parameters or variables most contribute to the output uncertainty. Its results indicate the parameter or variable on which experimental efforts should concentrate in order to reduce the output uncertainty. Note that there are only very few works dealing with sensitivity analysis of interval approaches [START_REF] Moens | Interval sensitivity theory and its application to frequency response envelope analysis of uncertain structures[END_REF], contrarily to probabilistic approaches.

In this paper, we propose a very simple method to perform this sensitivity analysis for a given function f i . Let [x i ](t) be the interval-valued output resulting from initial interval uncertainty. Then, if we denote by L(x i (t j )) := (x

( i t j ) + -x - i (t j ))
the length of [x i ](t j ) for any time-step t j , we can define the overall imprecision

I(x i (t)) of x i (t) as I(x i (t)) := T ∑ j=1 L(x i (t j )), (8)
that is, the sum of interval lengths obtained at each time step. Now, to quantify the impact of each parameter and variable uncertainty on the output imprecision, we propose the following procedure (similar to some existing propositions in imprecise probability literature [START_REF] Ferson | Sensitivity analysis using probability bounding[END_REF]): reduce, for each parameter and variable, its uncertainty by a given fraction r ∈ [0, 1], such a reduction coming down to trans-

form an interval [a, b] into an interval [a , b ] such that a = M([a, b]) - L([a,b]) 2 (1 -r) , b = M([a, b]) + L([a,b]) 2 (1 -r) , (9
(x i (t)) ≤ I(x i (t)). If [x i ](t)
is the output obtained after reducing the imprecision of parameter P (or variable X) by r, the gain in the output precision generated by this reduction can then be defined as

G(P, r) := (I(x i (t)) -I(x i (t)) I(x i (t)) . ( 10
)
Let us illustrate this notion in our example.

Example 5. Let us consider again the model of Example 1, with initial uncertainty given by p 1 ∈ [0.5, 1] and x 1 (0) ∈ [1, 3], the resulting x 1 (t) being pictured in Figure 3. We have I(x i (t)) = 70 (in this case, the area between the two lines of Figure 3 can be computed analytically).

If we reduce the uncertainty of P 1 and X 1 (0) by 50 % (r = 0.5), we get p 1 ∈ [0.625, 0.875] and x 1 (0) ∈ [1.5, 2.5]. We then have G(P 1 , 0.5) = 70-45 /70 0.36, G(X 1 , 0.5) = 70-60 /70 0.14.

From these results, it appears that reducing the uncertainty of P 1 by 50 % has a more important impact on the output than reducing X 1 (0) uncertainty by the same amount. Hence, further experiments should focus on reducing the uncertainty around P 1 .

Finally, note that a reduction of 100 % (r = 1) of a parameter or variable comes down to consider the middle point of the interval it belongs to.

Application to food packaging

In this section, we apply the above framework to the design and optimisation of food packaging material, with the aim to maximise shelf life of packaged food (here, fresh fruits and vegetables). Note that this model will be integrated to an online decision support system (http://www.tailorpack.com/), hence quick computations and approximate (but reliable) answers should be privileged over more exact but computationally more demanding answers.

Problem presentation and model analysis

Preserving fresh food after harvest and on the shop shelves is an important issue in the food industry. Providing food with optimised packaging is a sure way to avoid premature decay and prolong food edibility. Such optimisation means that the packaging must be permeable enough to oxygen and carbon dioxide so as to allow food respiration, but not too permeable to these two gases, so that maturation and decaying process are slowed down.

In modified atmosphere packaging, oxygen and carbon dioxide partial pressures in packaging head-space are modified and settle to steady values after a transient phase. This modification in the internal gas partial pressures is achieved due to the mass balance between oxygen and carbon dioxide flux through the packaging material and O2 and CO2 consumption/production due to the product respiration. For a given fruit or vegetable, it is possible to experimentally determine (up to some uncertainties) oxygen and carbon dioxide partial pressures that will result in an optimal preservation. From these optimal partial pressure values, it is then possible to determine optimal packaging permeances, by using a mathematical model predicting the dynamic evolution of internal gases partial pressures. Several environmental parameters must be specified in this model such as respiration rate or respiratory quotient, packaging geometry, or environment variables such as temperature. The mass balance between gas transfer and respiration can be written as:

ṗ pkg O 2 = Pe O 2 • S e (p ext O 2 -p pkg O 2 ) -RR O 2 • m = f 1 (11) ṗpkg CO 2 = Pe CO 2 • S e (p ext CO 2 -p pkg CO 2 ) + RR O 2 • m • QR = f 2 (12)
Parameter Name Units with

Pe O 2 O 2 permeance mol.m -1 .s -1 .Pa -1 Pe CO 2 CO 2 permeance mol.m -1 .s -1 .
RR O 2 = RR O 2 max • p pkg O 2 (Km appO 2 + p pkg O 2 ) • 1 + p pkg CO 2
Ki CO 2 [START_REF] Helton | Sensitivity analysis in conjunction with evidence theory representations of epistemic uncertainty[END_REF] where the first part of the right-hand side describes gas flux per time unit through the packaging material, while the second part describes gas consumption (and emission) by the vegetable or fruit (modelled using a Michaëlis-Menten-type equation, see [START_REF] Helton | Sensitivity analysis in conjunction with evidence theory representations of epistemic uncertainty[END_REF]). There are two variables X = {p pkg O 2 , p pkg CO 2 }. Table 1 summarises the parameters with their names and units, while Table 2

Case study: chicory

Note that in our case study, the shape of the packaging (S and e) and the mass (m) of food inside it have already been fixed, so that the only parameters that can be adjusted to optimise the packaging are the material permeability properties (Pe O 2 and Pe CO 2 ). Note that, if a probabilistic approach was used, no less than 9 probability distributions would have to be specified, and independence between variables would have to be assumed, with some of the parameters sometimes measured only 3 times. Also, estimating boundaries by interval analysis for p pkg O 2 and

Parameter P D , P E , X p pkg CO 2 here requires only 4 simulations (far less than for a Monte-Carlo simulation [START_REF] Helton | Sensitivity analysis in conjunction with evidence theory representations of epistemic uncertainty[END_REF]).

f 1 f 2 Pe O 2 ∈ P D P ,
The considered food is chicory, which has been previously studied, but without any proper uncertainty analysis [START_REF] Charles | Modeling of active modified atmosphere packaging of endives exposed to several postharvest temperatures[END_REF].

Uncertainty propagation step

Knowing the permeability of the packaging material classically used to pack chicory, the mathematical model described by Eq (11)-( 13) can be used to simulate the evolution with time of the internal O2 and CO2 partial pressures. Parameters with their interval uncertainty are summarised in Table 3. Maximal shelf life is estimated to be 200 hours, therefore the simulation time domain is T = [0, 200]h and the domain is discretised so that each time step is one minute.

The results of interval analysis are displayed in Figure 6, and we can see that, at the steady states, p pkg O 2 ∈ [1.6, 8.1]% and p pkg CO 2 ∈ [1.9, 5.3]%. Using Proposition 1, only four simulations were needed (two for each gas) to estimate the bounds.

Also note that, in the current case, initial conditions p pkg O 2 (0) and p pkg CO 2 (0) for the two variables are perfectly known and only parameters are uncertain, as initial partial pressures are the same as in free atmosphere (21% for oxygen and 0% for carbon dioxide). Therefore, only uncertainty pertaining to environmental and design parameters P E and P D have to be taken into account here (contrary to Example 1, where x 1 (0) value is also uncertain).

Parameter

!"#$%&'()(*+$#&$%,",$-*,-*./*

Optimisation step

In the case where best oxygen and carbon dioxide permeabilities suiting a particular fruit or vegetable are not known a priori, it is possible to use Equations ( 11)-( 13) to perform a reverse engineering task. In this case, optimal oxygen and carbon dioxide concentrations in the packaging modified atmosphere have to be specified. Once this is done, Algorithm 1 can be run to find optimal permeabilities that allow reaching the specified goal.

In this study, Algorithm 1 was run with the same uncertainty on environmental parameters as the one of (values are expressed in mol.m -1 .s -1 .Pa -1 ) are shown in Figures 7 and8 (Levenberg-Marquardt Algorithm [START_REF] Nocedal | Numerical Optimization[END_REF] was used to identify parameters). Figure 9 shows the result of the interval analysis done with optimal parameters belonging to [G Pe O 2 ] × [G Pe CO 2 ]. We see that the resulting imprecise oxygen and carbon dioxide partial pressures well lie in the objective bounds at the steady state.

sensitivity analysis

Table 4 summarises the sensitivity analysis performed according to the method presented in Section 3.4. For each variable p pkg O 2 (t) and p pkg CO 2 (t), we have evaluated the precision gain after an uncertainty reduction of r = 0.5 of each (environmental and design) parameter. The last column (All) indicates the gain in the output precision when all parameters uncertainty is reduced by r = 0.5.

The results indicate, among other things, that while O 2 permeability uncertainty has an important impact on both O 2 and CO 2 internal partial pressures uncertainty (the reduction resulting in a gain of about 0.1 for each), CO 2 permeability uncertainty only impact the CO 2 internal pressures, and have almost none effects on O 2 internal pressure. As could be expected, the respiration rate QR uncertainty only impact on the CO 2 internal pressure, while Ki CO 2 uncertainty, due to the large value of Ki CO 2 , has (almost) no impact on the resulting uncertainty for both variables.

Finally, to reduce the most efficiently the uncertainty on O 2 and CO 2 internal partial pressures, one should focus on a better characterisation of O 2 permeability and thickness, and on O 2 /CO 2 permeability, respectively.

Conclusion

Interval analysis is a useful alternative to probabilistic analysis that requires far less information to be applicable, thus avoiding the need to introduce hypothesis unsupported by available information. In some situations where the models satisfy some monotonic properties, using intervals is also more computationally efficient, as sampling methods are not needed to achieve computations. In this paper, we have used interval analysis for two purposes: classical uncertainty propagation for monotonic dynamical models and robust design under uncertainties. We have also proposed an easy method to perform some first sensitivity analysis.

The proposed robust design optimisation method, although approximate, quickly produces optimal values for design parameters. It takes account of interval uncertainty and can cope with imprecisely specified goals, distinguishing between possible solutions (potentially satisfying the goal) and guaranteed solutions (certainly satisfying the goal). The difference between the two kinds of solutions is represented by the means of fuzzy sets describing optimal solutions. Methodological perspectives to this work include the development of more precise methods, eventually ending up with better and multi-dimensional approximations of the sets G D and S D . Also, it may be desirable to extend the current approach to hybrid uncertainty models mixing interval and probabilistic uncertainty [START_REF] Baudrit | Joint treatment of imprecision and variability in food engineering: Application to cheese mass loss during ripening[END_REF][START_REF] Baudrit | Joint propagation and exploitation of probabilistic and possibilistic information in risk assessment[END_REF][START_REF] Fuchs | Potential based clouds in robust design optimization[END_REF][START_REF] Nassreddine | State estimation using interval analysis and belief-function theory: Application to dynamic vehicle localization[END_REF], as information concerning different parameters may vary in quantity and quality. However, these two perspectives would involve more computationally demanding procedures, thus reducing the number of models one could work with.

More applied perspectives include the combination of the optimisation system with a decision support system where a user can search a database for optimal packaging [START_REF] Destercke | A flexible bipolar querying approach with imprecise data and guaranteed results[END_REF] (fuzzy sets describing sets of optimal solutions will then be used as user preferences) and the application of the system to other food products (such as mushrooms).

Figure 1 :

 1 Figure 1: Illustration of interval analysis with inclusion function (Equation (2)).

Figure 2 :

 2 Figure 2: Example 1 set of (extreme) configurations.

Figure 3 :

 3 Figure 3: Interval analysis with Example 2 model

  Figure 3 represents the two envelopes when [p 1 ] = [0.5, 1] and [x 1 ](0) = [1, 3] for t ∈ [0, 10].

  [x 1 ](10) = [x 1 -, x 1 + ](10) = [8.5, 10.5] and [x 1 ](10) ⊆ [ x 1 ](10). The exact bounds of sets G D and S D are, in practice, hard to find. However one may search, for each P ∈ P D , intervals [G P ] = [g P -, g P + ] and [S P ] = [s P -, s P + ] approximating G D and S D on each design dimension. Those values are the solutions of the following optimisation problems (provided such solutions exist): g P -= min P, g P + = max P, under the constraints

Fig 4

 4 Fig 4.A guaranteed solution

Fig 4

 4 Fig 4.B possible solution x

Fig 4 .

 4 Fig 4.C non-admissible solution

Figure 4 :

 4 Figure 4: Illustration of solutions with Example 2 model

1 • 5 Finally(

 15 Line 6 of Alg. 1, i=1: Identification with x + 1 (10) = 11 and H i P E ,X = {x 1 (0) = 1}, giving p = 1 as solution and s P 1 ,1 + = Line 8 of Alg. 1, i=1: Identification with x - 1 (10) = 8 and H i P E ,X = {x 1 (0) = 3}, giving p = 0.5 as solution and s p 1 ,1 -= 0.Lines 10 and 11 of Alg. 1), we obtain intervals [G P 1 ] = [0.7, 0.8] and [S P 1 ] = [0.5, 1] providing approximate values of P 1 that give guaranteed and possible solutions, respectively.

Figure 5 :

 5 Figure 5: Fuzzy set of Example 4

  ) with M([a, b]), L([a, b]) the middle and length of the interval, respectively. Such a reduction gives a new interval-valued solution [x i ](t) included in the previous one, and with an overall imprecision I

Figure 6 :Figure 8 :

 68 Figure 6: Uncertainty propagation.

  !$%&'%()*%+!,-./0/%!1,/ !2"#!$%&'%()*%+!3-###,%!1,/

  t), objectives x i (t) + and x i (t) -on T , Initial conditions [x i ](0) for i = 1, . . . , m and [p] for all P ∈ P E , Alg id : identification algorithm returning design paramter values Output: Intervals [G P ] and [S P ] 1 for i = 1, . . . , m do

	5 6	-= x i (t) + with configuration H i p, else set g P, j + = p; If P ∈ P ,i , set g P, j Run Alg id on P E ,X , get optimal value	p of
		P;
	7	If P ∈ P ,i , set s P, j	+ =	p, else set s P, j	-=	p; i P E ,X , get optimal value	p of
		P;
	9	If P ∈ P ,i , set s P, j	-=	p, else set s P, j	+ =	p;

2 Run Alg id on x i (t) + with configuration H i P E ,X , get optimal value p of P; 3 If P ∈ P ,i , set g P, j + = p, else set g P, j -= p; 4 Run Alg id on x i (t) -with configuration H i P E ,X , get optimal value p of P; 8 Run Alg id on x i (t) -with configuration H

Table 1 :

 1 Parameters names and units. ×= non-relevant.

	Pa -1

  summarises whether they are variable, environmental or design parameters, and their monotonic w.r.t. each function f 1 and f 2 . Note that (p ext O 2p pkg O 2 ) is always positive, while (p ext CO 2p pkg CO 2 ) is always negative, that p ext i are constant and that all parameters are positive (otherwise model would not be monotonic w.r.t. every parameter).

Table 2 :

 2 Parameters type and monotonic w.r.t. f 1 , f 2 . ×= non-relevant.

Table 3 :

 3 Parameters uncertainty for propagation.

0-1*2./*!*2./*2.345676785*79:7;767.9

	O2Max
	O2Min
	CO2Max
	CO2Min

Table 3

 3 (i.e., all parameters except packaging material permeabilities to O 2 and CO 2 , which are the design parameters). Given that T = [0, 200]h and the fact that there exists a transient phase, we have chosen oxygen partial pressure to be between 4 and 10 % and carbon dioxide partial pressure to be between 2 and 5 % at the steady state. Fuzzy sets obtained from Algorithm 1 and from the intervals [G Pe O 2 ] = [1.59, 1.92]10 -15 , [G Pe CO 2 ] = [4.94, 7.50]10 -15 and [S Pe O 2 ] = [4.34, 6.51]10 -15 , [S Pe CO 2 ] = [1.78, 19.96]10 -15

	T = [150, 200], is we want	p pkg O 2 (t) = [4, 10]% and	p pkg CO 2 (t) = [2, 5]% for every t ∈	T , that

Table 4 :

 4 Values G(P, 0.5) of sensitivity analysis on the parameters P after a reduction r = 0.5, for the two variables p pkg O 2 (t) and p pkg CO 2 (t)

	O2Max
	O2Min
	CO2Max
	CO2Min
	O2MaxObj
	O2MinObj
	CO2MaxObj
	CO2MinObj

Note that design parameters p D are no longer uncertain.