Claude Manté
email: claude.mante@univ-amu.fr

Guillaume Bernard

Patrick Bonhomme

D Nerini

Application of Ordinal Correspondence Analysis for Submerged Aquatic Vegetation Monitoring

Keywords: ordinal correspondence analysis, tests, Chi-squared partitioning, ecological monitoring, Macrophytes, Berre lagoon

The European Water Framework states that macrophyte communities (seaweeds and seagrass) are key indicators of the ecological health of lagoons. Furthermore the restoration of these communities, especially the Zostera meadows, is one of the main objectives of the Berre lagoon restoration plan. Consequently a monitoring program of the main macrophyte species still present in the lagoon was initiated in 1996. This monitoring resulted in a sequence of eleven spatially structured annual tables consisting of the observed density of these species. These tables are processed in this study.

First, we specify the principles of Beh's Ordinal Correspondence Analysis (OCA), designed for ordered row/column categories, and compare this method to classical Correspondence Analysis (CA). Then, we show that OCA is straightforwardly adaptable for processing a sequence of ordered contingency tables like ours. Both OCA and CA are afterwards used to reveal and test the main patterns of spatiotemporal changes of two macrophyte species in the Berre lagoon: Ulva and Zostera. The results we obtained are compared and discussed.

Introduction

Since the early to mid-20 th century coastal lagoons and estuaries have become among the most disturbed coastal ecosystems in the world. The initial causes are relatively well identified: eutrophication and organic pollution as a result of increasing agriculture, urbanisation in river catchments and port facilities, aquaculture, turbidity and oversedimentation [START_REF] Giesen | Eelgrass condition and turbidity in the Dutch Wadden Sea[END_REF][START_REF] Cardoso | Dynamic exchanges in seagrass assemblages under eutrophication and implications for recovery[END_REF]. Seagrass beds, which are common dwellers of these habitats, are declining throughout the world [START_REF] Short | Natural and human-induced disturbance of seagrasses[END_REF] at the expense of macroalgae. Indeed, eutrophication created suitable conditions for macroalgae blooms in several coastal lagoons, bays and estuaries over the last decades [START_REF] Valiela | Nitrogen loading from coastal watersheds to receiving estuaries: Review of methods and calculation of loading to Waquoit Bay[END_REF][START_REF] Bowen | The ecological effects of urbanization of coastal watersheds: Historicalincreases in nitrogen loads and eutrophication of WaquoitBay estuaries[END_REF][START_REF] Hauxwell | Light dependence of Zostera marina annual growth dynamics in estuaries subject to different degrees of eutrophication[END_REF]. In general, the process of eutrophication leads to a shift in the macrophytobenthic community from slow-growing seagrasses to phytoplankton and fast-growing macroalgae such as Ulva. The high surface-area to volume ratio of these macroalgae favours rapid nutrient uptake, high production and rapid growth rates, which enables them to outcompete the original vegetation [START_REF] Pedersen | Nitrogen dynamics of eelgrass Zostera marina during a latesummer period of high growth and low nutrient availability[END_REF][START_REF] Malta | Regulation of spatial and seasonal variation of macroalgal biomass in a brackish, eutrophic lake, Helgol[END_REF].

Monitoring of submerged aquatic vegetation distribution and abundance ranges from coarse assessments of presence/absence or area distribution in large areas -based on remotely sensed data and presented as macroscale maps -to fine-scale diver assessments of depth limits and cover (in % of the bottom surface area), biomass or shoot density along depth gradients [START_REF] Krause-Jensen | How are seagrass distribution and abundance monitored?[END_REF][START_REF] Cf | Préservation et conservation des herbiers à Posidonia oceanica[END_REF]. The cover is indeed the most frequently used measure for vegetation surveys, since it is not destructive and requires relatively little effort compared to other measurements of vegetation health.

The Berre lagoon (Provence, France) is one of the largest Mediterranean brackish lagoons (155 km²) and was occupied at the turn of the 20 th century by extensive seagrass meadows (perhaps over 6000 ha). Subsequently the lagoon was disturbed by urban and industrial pollution. The monitoring of macrophytes in this lagoon is based on visual censuses in scuba-diving along permanent transects, each compounded of 20 spatial units. A coding in six classes (intervals of cover) was adopted for studying the temporal variations (1996 to 2008) of the cover by eight species (or groups of species) of macrophytes. In this paper, we will focus on only two species of cardinal importance for evaluating the health of this ecosystem: Zostera and Ulva. Our goal is double:

1. from the spatio-temporal side, we want to investigate whether or not the distance to the shoreline (indirectly associated with coastal topography, depth, currents, etc.) is related to the abundance of each species and, when the answer is positive, we will try to identify homogeneous regions where its dynamics is important (regions of interest) 2. from the purely temporal side, it is capital to investigate the ways the new policy for freshwater discharge established since 2004 (the 8 th observation) impacted the abundance of each species in its regions of interest.

We chose to investigate these points through the Ordinal Correspondence Analysis (OCA) of Beh [START_REF] Beh | Simple Correspondence Analysis of ordinal cross-classifications using orthogonal polynomials[END_REF][START_REF] Beh | A comparative study of Scores for Correspondence Analysis with ordered categories[END_REF][START_REF] Beh | Simple Correspondence Analysis : a bibliographic review[END_REF] and classical Correspondence Analysis (CA). Because of the above decision-making concerns, we will focus on the different decompositions of the Peason chi-squared statistics associated with CA and OCA. For that purpose, we first compare decomposition and reconstitution formulas associated with both of these methods and propose original representations of the chi-squared statistics apportionment. We also propose original representations of the interaction terms in OCA, and a variant of this method suited for analyzing a sequence of doubly-ordered contingency tables.

Material and methods

Study site and sampling procedures

In 1966 the diversion of the Durance River towards a hydroelectric power plant and then into the Berre lagoon resulted in a heavy input of freshwater. Surface water salinity declined from 24-36 to 1-22 resulting in haline stratification, eutrophication and unstable ecological conditions [START_REF] Nérini | Analyse statistique de processus physiques et chimiques en océanologie côtière à l'aide d'une méthode de régression et de classification par arbre décisionnel[END_REF][START_REF] Nérini | Analysis of oxygen rate time series in a strongly polluted lagoon using a regression tree method[END_REF][START_REF] Stora | Time and spatial dynamics of Mediterranean lagoon macrobenthos during an exceptionally prolonged interruption of freshwater inputs[END_REF]. The seagrass beds underwent a dramatic decline. In fact they are functionally extinct since 1998; in 2004 they occupied a total surface area of 1.5 ha [START_REF] Bernard | Long term changes in Zostera meadows in the Berre lagoon (Provence, Mediterranean Sea), Estuar[END_REF]. During that time, large macroalgae blooms mainly constituted of Ulva were observed in the lagoon. Consequently enormous amounts of drift algae were washed up on the shore [START_REF] Pergent-Martini | Localisation et évolution des peuplements de phanérogames aquatiques de l'étang de Berre (Bouches du Rhône-France)[END_REF][START_REF] Soltan | Monitoring system of benthic macrophytes communities in Berre and Vaıne ponds: development of a new strategy[END_REF] (G. Bernard and P. Bonhomme pers. obs.).

A new policy for freshwater discharges aiming at restoring the lagoon started up in 2004. The policy induced strong changes in the lagoon: freshwater and silt inputs from theSt-Chamas hydroelectrical plant were heavily reduced. As a result, the global salinity of the lagoon increased while its stratification and the frequency of anoxia phenomena in its deepest part decreased.

[Figure 1 near here] Thirty one permanent transects [START_REF] Soltan | Monitoring system of benthic macrophytes communities in Berre and Vaıne ponds: development of a new strategy[END_REF] (G. Bernard & P. Bonhomme, pers. obs.), 100 m long by 1 m width, and regularly distributed along the shore (see Fig. 1), were annually explored by scuba-diving, in July, from 1996 to 1998 and from 2001 to 2008 (11 years of sampling). Each transect is divided in 20 successive individual units (5m² each). The cover of Zostera and Ulva (in % of the bottom surface area) has been evaluated by visual census for each individual unit and has been coded with an ordinal scale in 6 classes, ranging from 0 to 5 (see table 1). We used a "à la Braun-Blanquet'' [START_REF] Podany | Braun-Blanquet's legacy and data analysis in vegetation science[END_REF] cover score adapted for marine ecology, because the field of vision of divers is reduced.

[Table 1 near here] Thus, for each species, the data consists in 11 annual contingency tables, of common size 20 6 × ; the i th row of each table lists the overall frequency ( )

31 ≤ of the abundance code 1 - i
at each one of the 20 successive units composing transects. Thus, both the variables associated with rows and column of these tables are ordinal.

Statistical methods

Notice first that the data set associated with each macrophyte consist in a "cube" of size

K J I × ×
, with I=6, J=20, K=11. There are in the literature a lot of statistical methods suited for such three-ways data, which can be essentially divided into two large families. The first one is a whole family of genuine three-way methods:

CANDECOMP/PARAFAC, TUCKER, etc. [START_REF] Kroonenberg | Three-mode Principal Component Analysis[END_REF], which are extensions of PCA to three-ways quantitative data. These methods were not designed for contingency tables, but Carlier and Kroonenberg [START_REF] Carlier | Decompositions and biplots in three-way correspondence analysis[END_REF] proposed an extension of CA to three-way contingency tables. Nevertheless, to our knowledge, none of these methods has been adapted to ordinal data.

The second family was designed for processing a sequence of partial tables.

Roughly speaking, these methods consist in analyzing a "compromise" table built from this sequence. For instance, in the setting of Multiple Factor Analysis and related methods, the compromise is obtained by weighting the partial tables according to their dominant eigenvalues [START_REF] Zarraga | Simultaneous analysis and multiple factor analysis for contingency tables: two methods for the joint study of contingency tables[END_REF], while in the setting of STATIS [START_REF] Lavit | Analyse Conjointe de Tableaux Quantitatif[END_REF] the compromise is obtained from a preliminary PCA of operators associated with the partial tables. In the same line, Foucart [START_REF] Foucart | Sur les suites de tableaux de contingence indexés par le temps[END_REF] proposed a method consisting in performing CA of the average table, the rows and columns of the partial tables being represented as supplementary elements. This method is frequently used by ecologists [START_REF] Pavoine | A new technique for ordering asymmetrical three-dimensional data sets in[END_REF], and we will show that it is straightforwardly adaptable for ordinal data. The relationships between both these families were investigated by Kiers [START_REF] Kiers | Comparison of "Anglo-saxon" and "French" three-mode methods[END_REF][START_REF] Kiers | Hierarchical relations among three-way methods[END_REF].

Let us now introduce the notation used. Depending on the context, J I × will denote either the dimensions of a table of I rows and J columns, or the Cartesian product of the sets { }

I , , 1 L and { } J , , 1 L . We shall denote K 1 the constant vector t ) 1 , , 1 , 1 ( L of dimension K, ( ) ( ) J I j i j i V U V U V U × ∈ = = ⊗ , ' :
the tensor product of two vectors of respective length I and J, and o the element-wise (Hadamard) product of matrices: (

) j i j i j i B A B A , , , := o
, where j i A , denotes the entry of the matrix A at the i th row and j th column (for further details, see for instance [START_REF] Horn | Topics in Matrix Analysis[END_REF]). At last, the expression

" l ≅ V
" will mean that the random variable V obeys the distribution l , and p m δ will denote the Kronecker symbol ( 0 (Benzécri, 1967;[START_REF] Beh | Simple Correspondence Analysis : a bibliographic review[END_REF][START_REF] Greenacre | Theory and Applications of Correspondence Analysis[END_REF] It is time now to remind the reader of essential principles of CA. Let T be some

= p m δ if p m ≠ and 1 = p p δ ).

Correspondence Analysis in a nutshell

( ) M m B A m m m ≤ ≤ 0 : , ; λ , with ] 1 , 1 [ : - - = J I Min M (remember that 1 0 = λ , I A 1 0 = and J B 1 0 = ).
We can now write:

j m m M m i m j i B A , 0 , , λ θ ∑ = =
or, using tensor products [4,[START_REF] Horn | Topics in Matrix Analysis[END_REF] (see also the appendix):

) ( ) ( ) ) 2 ( 1 ∑ = ⊗ ⊗ + ⊗ = Θ ⊗ = M m m m J I
) 3 ( , ) , ( , 1 , , 1
,

p m j p j J j j m i p i I i i m B P B A P A M M p m δ = = × ∈ ∀ ∑ ∑ = = .
They also verify:

∑ × ∈ = J I j i j p j i i m p m m B P A ) , ( , , , ) 4 ( δ λ
.

Consider now the Pearson chi-squared statistics

( ) 2 ) , ( , 2 1 :
∑ × ∈ - × = Χ J I j i j i j i P P n θ .
Because of formula (1), we obtain an alternative decomposition of this statistics:

∑ = = Χ M m m n 1 2 2 ) 5 ( λ .

Correspondence Analysis with ordered categories

Since the principal axes issued from CA do not take into account the order of the categories associated with the rows and the column of P , Beh [START_REF] Beh | Simple Correspondence Analysis of ordinal cross-classifications using orthogonal polynomials[END_REF] proposed an alternative decomposition of 2 Χ , where the role of principal axes is played by orthonormal polynomials. To compare his Ordinal CA (OCA) with CA, we will use the same notations as in the previous section, denoting { }

1 1 : - ≤ ≤ I m A m
the system of row polynomials, and { }

1 1 : - ≤ ≤ J p B p
, the one of column polynomials. Each one of these systems only depends on a score [START_REF] Beh | A comparative study of Scores for Correspondence Analysis with ordered categories[END_REF], and on the associated marginal probability, either I P or J P . The chosen score reflects the ordered structure of the categories. The polynomials, resulting from Emerson's recurrence formula [START_REF] Emerson | Numerical construction of orthogonal polynomials from a general recurrence formula[END_REF][START_REF] Beh | A comparative study of Scores for Correspondence Analysis with ordered categories[END_REF], are centred and verify:

) 6 ( ), 1 ( ) 1 ( ) , ( , 1 , , 1 , p m j p j J j j m i p i I i i m B P B A P A J I p m δ = = - × - ∈ ∀ ∑ ∑ = =
There are substantial differences between principal axes of CA and these polynomials: as singular vectors, axes are paired, as shown in formula (1); moreover they are naturally sorted according to their common variance. Emerson's polynomials don't possess such properties, and the decomposition of the variance associated with OCA includes all the bivariate moments between row and column polynomials, giving

rise to a whole ) 1 ( ) 1 ( - × - J I matrix Columns Rows ,

Λ

, of general entry:

∑ × ∈ = J I j i j p j i i m p m B P A ) , ( , , , , ) 7 ( : λ .
We have also a reconstitution formula [START_REF] Beh | Simple Correspondence Analysis of ordinal cross-classifications using orthogonal polynomials[END_REF]:

( ) ( ) . ) 1 ( ) 1 ( ) , ( ,
∑ - × - ∈ ⊗ ⊗ + ⊗ = J I p m p m J I p m J I B A P P P P P o λ
While the interest of this saturated formula is only theoretical, it can be coupled with tests in order to get some "de-noised table". More precisely, suppose a set

J I S × ⊂
of relevant interactions has been determined. The associated de-noised table is given by the unsaturated reconstitution formula [START_REF] Beh | Simple Correspondence Analysis of ordinal cross-classifications using orthogonal polynomials[END_REF]: Notice the similarity between formulas ( 7) and ( 4), ( 3) and ( 6), and ( 8) and [START_REF] Beh | A comparative study of Scores for Correspondence Analysis with ordered categories[END_REF].

( ) ( ) ).

Representation of a sequence of ordered tables

For this purpose, we will straightforwardly adapt Foucart's Correspondence Analysis [START_REF] Foucart | Sur les suites de tableaux de contingence indexés par le temps[END_REF][START_REF] Pavoine | A new technique for ordering asymmetrical three-dimensional data sets in[END_REF]. Remember that the processed table is the sum of 11 annual partial tables, with identical rows and columns ("slices") of the data cube

∑ = = K k k T T 1 :
, and let us denote [START_REF] Foucart | Sur les suites de tableaux de contingence indexés par le temps[END_REF] analyzed T through CA, we will use OCA instead.

∑ × ∈ = J I j i k j i k T ) ,
∑ = = K k k k P P 1 α , with ∑ = = K k k k k 1 / π π α . While Foucart

Tests

The unique rigorous independence test in CA is based on 2 Χ which, under the hypothesis (H) of independence of the columns and rows of T, asymptotically obeys a chi-square distribution:

( ) [4] derived from this statistics and ( 5) a heuristic for determining the number of significant principal axes in CA of contingency tables. It consists in seeking for the smaller m such that, for some

) 1 ( ) 1 ( 2 2 - - ≅ Χ J I χ . Benzécri
fixed [ 5 . 0 , 0 ] ∈ τ : ( ) ( ) ) 9 ( , ) 1 ( ) 1 ( 1 2 2 ∑ + = - - < M m k k J I Q n τ χ λ
where the threshold ( ) τ , l Q denotes the quantile of order τ of the probability law l (Benzécri [4] merely proposed to use the mean

) 1 ( ) 1 ( - -J I
as a threshold). This heuristic is useful for data representation and reconstitution [4,[START_REF] Beh | Simple Correspondence Analysis of ordinal cross-classifications using orthogonal polynomials[END_REF][START_REF] Beh | Simple Correspondence Analysis : a bibliographic review[END_REF] from the set of significant axes.

Contrary to CA, OCA gives rise to a three-level battery of tests. Firstly, at the cell-level,

Rayner & Best [START_REF] Rayner | Smooth extensions of Pearson's product moment correlation and Spearman's Rho[END_REF] proved that under (H), we have asymptotically:

) 10 ( ) 1 , 0 ( ), 1 ( ) 1 ( ) , ( , N n J I p m p m ≅ - × - ∈ ∀ λ
and that all such bivariate moments are asymptotically independents. Thus, since every polynomial { }

1 1 : - ≤ ≤ I m A m can be correlated with every { } 1 1 : - ≤ ≤ J p B p ,
we will use [START_REF] Cambridge | Retrospective analysis of epiphyte assemblage in relation to seagrass loss in a eutrophic coastal embayment[END_REF] to test whether or not each interaction can be considered as significant.

At the intermediary level we can also test whether or not each row polynomial m A (respectively column polynomial p B ) significantly contributes to the Pearson chisquared statistics, since we have under (H): and these quantities are relevant criteria for selecting best projection planes for the columns (resp. the rows) profiles (because of the duality between rows and columns).

∑ - = - ≅ = 1 1 2 ) 1 ( 2 ,

Data analysis

Scores and weights

As we said in the previous section, the system { } (see Table 1). Notice that, due to the sampling conditions,

J J J P 1 1 =
does not depend on the data, while I P does.

Spatial and temporal dynamics of Ulva

We will now analyze the 20 6 × contingency table This will be done first through CA, and afterwards through OCA; issues from the second analysis will be explained in greater detail.

∑ = = 11

Results of CA of the Ulva data

The overall Pearson statistics was 346.05: this is highly significant since, under . In both cases, formula (9) indicated that only two axes are significant, associated respectively with 67.2% and 13.5% of the total inertia. The corresponding components are plotted on Figure 2; the individual units (in grey) are labelled by their distance to the shoreline along the sampling transects (from 0-5m to 95-100m), while the cover is labelled by the original code (from 0 to 5).

[Figure 2 near here]

We can infer from Figure 2 that the first 20 meters of the transects are equally associated with the codes {2,3,4,5} and concentrates most of the Ulva abundances trough the time series, while the last 50 meters are associated with the absence of Ulva (code 0).

Results of OCA

We will first use formulas [START_REF] Cardoso | Dynamic exchanges in seagrass assemblages under eutrophication and implications for recovery[END_REF] and [START_REF] Carlier | Decompositions and biplots in three-way correspondence analysis[END_REF] to determine the relevant polynomials.

On Figure 3 are plotted the contributions { }

1 1 : ˆ-≤ ≤ J j n j µ
on the left, and

{ } 1 1 : ˆ-≤ ≤ I i n i
η on the right; on these plots, grey and black lines respectively mark the quantiles of order 0.5 and 0.95 of the associated chi-squared distributions. It is noteworthy that while 1 μ corresponds to about 60.5% of total variance (close to the first principal axis: 67.2%), 1 η corresponds to only 29.4% of the same quantity. This stresses the complexity of column profiles of this table.

[Figure 3 near here]

On Figure 3, the only significant column polynomials are 1 B , 2 B and, at the very outside, 3

B , while all the { }

5 1 : ≤ ≤ m A m
seem statistically significant. Since in both cases both the first polynomials correspond to higher variances, it is reasonable to simultaneously display on Figure 4 the row profiles projected on the plane ( ) [Figure 4 near here]

Note that, even if Figure 4 is quite different from Figure 2, its interpretation is similar: the same clusters of rows and columns are evidenced in both plots.

Issues of the tests associated with formula [START_REF] Cambridge | Retrospective analysis of epiphyte assemblage in relation to seagrass loss in a eutrophic coastal embayment[END_REF] λ is associated with some cell whose grey level depends on its estimated value.

Precisely, each statistics

( ) ] 1 , 0 [ , 1 : ) , ( 2 , 2 1 ∈             > = p m n P Min p m Gl λ χ τ τ
has been computed and coded in gray levels for producing Figure 5 (and Figure 12 too).

[Figure 5 λ would be low; since it is not the case, we must be cautious! Furthermore, one should also be careful with polynomials of degree greater than 4, because the estimation variance of these moments can be excessive [START_REF] Rayner | Smooth extensions of Pearson's product moment correlation and Spearman's Rho[END_REF]. We can indeed see on Figure 5 , corresponding to about 20.5% of the total inertia. Thus, the salient feature of T is that the more we move away from the shoreline, the more the cover by Ulva decreases (see the left panel of Figure 6).

[Figure 6 near here]

The second dominant interaction is the quadratic by linear one 35 . 7 ˆ1 , 2 ≈ λ . In this case, when we move away from the shoreline (see the right panel of Figure 6), the cover, which is initially medium, gradually concentrates on extreme values: the dispersion increases. More insights about interactions are provided in the appendix.

The set S of significant interactions, displayed on Figure 5, enabled us to reconstitute the table by using formula (8); the result is displayed on Figure 7.

[Figure 7 near here]

Finally, we display on Figure 8 the spatio-temporal variations of cover by Ulva, projected onto the plane ( ) To save place, we discarded the units from 50 to 100m, where the dynamics is less important.

[Figure 8 near here]

Note that until now, the time was passive (projection of supplementary elements). But, since the units 0 to 20m consist in an homogeneous region, we can make the time play an active part by considering the "transposed cube" of type

J K I × × (I=6, J=20, K=11),
∑ = = 4 1 : k k U U
in the same way as T ; notice that while the row polynomials remain the same as in the previous analysis, the column polynomials

{ } 1 1 : - ≤ ≤ J p B p
are now associated with time, which is thus active. To save place, we will only mention the main issues of OCA and CA of U .

The overall chi-squared statistics associated with U is highly significant (under (H), X 2 =527.8 should obey ( )

50 2 χ
). All the row polynomials are significant, and the highest contributions of polynomials, 1 μ and 1 η correspond to only 44% of total variance, while 2 μ and 2 η correspond to respectively 15% and 9.8% of total variance.

Thus, CA seems more efficient than OCA for displaying the data: the first eigenvalues of CA corresponded to 83.6% and 10.12% of the total variance while, according to the heuristic ( 9), the next ones were not significant. We can see on Figure 9 [Figure 9 near here]

From the inference point of view, complementary insights can be derived from OCA, thanks to the table of bivariate moments (of type λ is also very important, and corresponds to a progressive regularization of the cover in the region of interest: extreme values (cover around 20% or 80%) got rarer, while medium values (about 50%) became more frequent.

Spatial and temporal dynamics of Zostera

We proceeded the same way for Zostera as for Ulva, processing first the table

∑ = = 11 1 : k k T T
successively through CA and OCA.

Results of CA of the Zostera data

Notice first that the overall Pearson statistics is now 228.4: this is significant for some ( )

95 2 χ
, but remember it was 346.05 for Ulva. Correlatively, we found from ( 9) that a single principal axis is significant, which is associated with 75.3% of the total inertia. The first principal plane is displayed on Figure 10.

[Figure 10 near here]

This plane essentially highlights that Zostera was mostly absent in the lagoon (the code 0 represent 97.3% of the observations) and that the highest modalities of cover were extremely rare for this species; it stresses that most Zostera abundances were observed in the first part of the sampling transects from 5 to 25 m, from the shore, in very shallow waters.

Results of OCA of the Zostera data

One can see on Figure 11 that four column polynomials { } [Figure 11 near here]

We separately plotted on Figure 12 the row and column profiles because, as one can see on Figure 13, simultaneous representation could be dubious: 5 [Figure 12 near here]

The reader can also see on Figure 13 that only a few of the polynomials of degree greater than 4 are important for Zostera. Nevertheless, the very dominant term in the decomposition is 4

. 7 ˆ1 , 1 - ≈ λ
, corresponding to about 23.8% of the total inertia. Thus, the salient feature of T is that the more we move away from the shoreline the more the cover by Zostera decreases. It is noteworthy that the same phenomenon appears for Ulva and Zostera with similar proportions of variance.

[Figure 13 near here]

The second important interaction in this Space Density,

Λ

table is the quintic by linear . We noted that in these two cases, interactions are dominated by low values of the cover. In these cases, the probability of not observing Zostera is highly oscillating along the transect, without any trend. This erratic spatial structure is is probably linked with the rarity of Zostera, and its patchiness.

On the other hand, notice that 1 μ corresponds to about 38.7% of total variance (much less than the first principal axis of CA: 75.3%), and 1 η corresponds to 53.6% of the same quantity.

Thanks to formulas ( 8) and ( 10), the table could be reconstituted from the set S of significant interactions, displayed on Figure 13. Since the "independence structure"

J I P P ⊗
is very strong, we preferred to plot the discrepancy of P from independence,

J I P P P ⊗ -
, and its approximation ( ) ( )

∑ ∈ ⊗ ⊗ = ⊗ - S m m m J I m J I S B A P P P P P o λ ~.
One can see on Figure 14 how the original discrepancy is filtered this way, and how much the position of null densities on the transect is important.

[Figure 14 near here]

Next, we display on Figure 15 the temporal variations of cover by Zostera, projected onto the optimal plane ( )

5 1 , A A .
We discarded the segment 50-100m, where there is practically no Zostera.

[Figure 15 near here]

We can see on Figure 15 that for Zostera too, through the whole time series, most part of the dynamics took place in the 20 first meters of the transect, very close to the shoreline, and that these four units were those where this species was most abundant (see Figures 10&12). Thus we will proceed the same way as for Ulva, keeping only the first four slices of the transposed cube for investigating further the temporal variations of Zostera. To save place, we will only mention the main results of OCA and CA of U .

The overall chi-squared statistics associated with this table is hardly significant (X 2 =70.16 while, under (H), it would obey ( )

50 2 χ
); remember that in the case of Ulva it was 527.8! The decomposition of U through OCA is also very different from that obtained for Ulva: only the two first row an column polynomials are significant. They correspond to moderate percentages of the total variance: 25.2 and 29.1% for the column polynomials, 26.5 and 27.6% for the row polynomials. In the case of CA, only the first component (62.7%) is significant; we will use OCA for displaying the data on Figure 16, where the observations posterior to 2004 are clearly opposed to the previous ones.

[Figure 16 near here]

The very dominant bivariate moment in the λ are significant too, but with much bigger Pvalues (0.037 and 0.045); each of them is associated with about 6% of total variance. Thus, through the whole time series, there is a significant decrease of Zostera abundance, even if it remains in low levels of abundances, close to 0. On Figure 16, the years 2004 to 2008 are associated with the score 0, indicating a slight but indubitable rarefaction of Zostera.

Ecological significant results for submerged aquatic vegetation monitoring

Among submerged aquatic vegetation, Ulva and Zostera are identified as biological indicators of ecological status in the Water Framework Directive for transition waters. As a consequence, their monitoring concerns most of the lagoon's managers worldwide. We have shown that the abundances of Ulva and Zostera in the Berre lagoon are statistically dependent of the distance to the shore, and that this spatial structuring is stronger but more complicated for Ulva than for Zostera. Most of the cover of these two species is located in the first 20 meters of the sampling transects, very close to the shore in very shallow water. The main part of the variability is also concentrated along the first 20 meters of the sampling transects, and remain high up to 50m from the shore. This is probably due to the Ulva way of growing, which spends most of its life-cycle as a free drift, transported by the currents and accumulated on the shore depending on the hydrodynamic conditions. As a consequence, for each year, bathymetry and coast-line topology could be related to Ulva density in a point.

Concerning Zostera, the present day extent of the seagrass meadows in the lagoon, restricted to very shallow waters, suggest light as the limiting factor [START_REF] Valiela | Nitrogen loading from coastal watersheds to receiving estuaries: Review of methods and calculation of loading to Waquoit Bay[END_REF][START_REF] Bernard | Long term changes in Zostera meadows in the Berre lagoon (Provence, Mediterranean Sea), Estuar[END_REF]. Light reduction could be due either to the silt input being still too high or to sediment resuspension. It could also be related to eutrophication of the lagoon which results in micro-and macro-algae blooms together with high levels of colonization of Zostera leaves by epiphytes and light reduction for the seagrass meadows.

Trough time, Ulva cover shows a clear increase from 1998 to 2008 despite a high interannual variability: years with high Ulva abundances (2002) can be distinguished from years with weak abundances (1996)(1997)(1998). From the first years of the time-series, the Zostera cover remained very weak and the abundances decrease with time. There is no clear sign of re-colonization since 2005, when the environmental pressure (inputs of freshwater, silts and nutrients) have been strongly reduced. The reasons could be, in addition to continuing nutrient inputs, the resuspension of silt, no longer trapped under the seagrass canopy, during wind episodes, which are frequent in the area, and/or the release of nutrients from the bare silt habitat, which would constitute an indication of a possible hysteresis of the system. Competition between the two species could also be one of the constraining factors for Zostera growth, the suffocation by an enteromorph-like Ulva of a Zostera meadow has been described in the literature [START_REF] Silberstein | The loss of seagrasses in cockburn sound,Western Australia. III. The effect of epiphytes on productivity of Posidonia australis Hook f[END_REF][START_REF] Brun | Resource allocation and proteolyticactivity in the seagrass Zostera noltii Hornem. subjected to light reduction by Ulva canopies[END_REF][START_REF] Brun | Effect of shading byUlva canopies on growth and carbon balance of the seagrass Zostera noltii[END_REF][START_REF] De Casabianca | Preliminary results on eelgrass regression and redseaweed dominance under increasing eutrophication (Thau Lagoon, France[END_REF][START_REF] Cambridge | Retrospective analysis of epiphyte assemblage in relation to seagrass loss in a eutrophic coastal embayment[END_REF].

Conclusion

We studied an 11 year long sequence of the abundance of Ulva and Zostera in a heavily disturbed Mediterranean lagoon, the Berre lagoon, through the Ordinal Correspondence Analysis of Beh. The results are compared with those issued from classical Correspondence Analysis. Both of these methods consist in distinct decompositions of the Pearson chi-squared statistics and of the corresponding contingency table, according to different systems.

In the case of CA, since the systems consist in singular vectors, there is no interaction between factors of different orders, which is optimal in a purely descriptive perspective.

The price to pay lies in the complexity of the factors.

The systems associated with OCA are very simple (orthogonal polynomials), but it takes into account numerous interactions between row and column systems, giving rise to representation difficulties (simultaneous projections, for instance). Consequently, we propose a number of graphical guidelines ( 2 χ plots, interactions plots, partial table plots) to help the user to select the relevant components and interpret the associated interactions.

Indeed the main advantages of OCA over CA (in addition to simplicity of the factors) is that it takes into account the ordinal nature of row and/or column variables, and that it is possible to test the significance of all the pairs of factors. We could thus obtain denoised tables given by the unsaturated reconstitution formula, and determine significant patterns of spatial (firstly) and temporal (secondly) changes of Ulva and Zostera abundance.

The spatial location of Ulva and Zostera seems to be strongly dependant on the distance to the shore, and indirectly to the depth. Main abundances and variability are observed very close to the shore since this is where free-floating Ulva drift can accumulate and limit light for the deeper growing Zostera. Through time Ulva shows a clear increase of its abundances while Zostera remained in very weak abundances; this suggests a high level of eutrophication despite recent reduction of nutrient inputs in the lagoon.

and it is straightforwardly extensible to three (or more) ways tables (see for instance [START_REF] Carlier | Decompositions and biplots in three-way correspondence analysis[END_REF] or most papers about PARAFAC) ; we will use such products in the next paragraph.

About bivariate moments

Remember the aim of OCA is to highlight the ways P possibly differs from the matrix

×

. Proceeding the same way with the sampled probability, we can finally describe each interaction with the bi-dimensional function Table 1. Classes and scores corresponding to the cover of the bottom surface area by Ulva and Zostera for each 5m² successive units of the 31 permanent transects along the Berre lagoon shore. See 3.1 for the definition of the scores.

  the reconstitution formula:

(

  profiles obtained from the OCA of P are situated at the barycentre of the partial probabilities projected as supplementary elements.

  the global level, we find again the classical test, associated with a third decomposition: η (resp. p µ ) is the part of the total information bore by m A (resp. p B ),

  the number of principal axes to keep, we fixed in formula

  e. the location and dispersion components in the terminology of Beh[START_REF] Beh | Simple Correspondence Analysis of ordinal cross-classifications using orthogonal polynomials[END_REF].

2 1 ,Figure 8

 18 Figure 8 that most of the dynamics of Ulva took place in the 20 first meters of the

  that through the whole time series, the Ulva abundance increased trough time, from 1996 to the years 2001 and followings (unfortunately, 1999 and 2000 are missing). The year 2005 is singular, while 2002 is associated with extreme abundances of Ulva.

  . Consequently, the best representation of the column profiles is obtained in the plane ( )

.

  Plotting the associated function, we saw that this interaction is dominated by very low values of the cover due to the rarity of Zostera. Indeed, the probability of not observing it decreases when we move away from the shoreline.

  domain. As an illustration, we have shown on Figure6the dominant interactions obtained in the analysis of the spatiotemporal distribution of Ulva (see section 3.2.2).

Figure 1 .

 1 Figure 1. location of the permanent transects along shore for monitoring the submerged aquatic vegetation in the Berre lagoon.

Figure 2

 2 Figure 2 Ulva: the first principal plane of CA. The individual sampling units clearly constitute three groups: 0 to 20m, 20 to 50m and 50 to 100m.

Figure 3 :

 3 Figure 3: Chi-square tests on row and column polynomials for Ulva.

Figure 4 :

 4 Figure 4: location and dispersion components for Ulva.

Figure 5 :

 5 Figure 5: significance of bivariate moments at the threshold 0.05 for Ulva.

Figure 6

 6 Figure 6 Ulva: the two main bivariate moments of the spatio-temporal OCA.

Figure 7 :

 7 Figure 7: approximation of P by S P ~, from formula (8), for Ulva.

Figure 8 :

 8 Figure 8: projection of the partial tables onto the plane (

Figure 9 :

 9 Figure 9: first plane of temporal CA for Ulva (83,6%, 10.1 %).

Figure 10 :

 10 Figure 10: Zostera: the first principal plane of CA.

Figure 11 :

 11 Figure 11: Chi-square tests on row and column polynomials for Zostera.

Figure 12 :

 12 Figure 12: row and column profiles projected onto the plane (1,5) of Zostera.

Figure 13 :

 13 Figure 13: significance of bivariate moments at the threshold 0.05, for Zostera.

Figure 14 :

 14 Figure 14: the discrepancy, and its approximation for Zostera.

Figure 15 :

 15 Figure 15: projection of the partial tables onto the plane (

Figure 16 :

 16 Figure 16: first plane of temporal OCA for Zostera.

  near here]

	(	Gl	(	m	,	p	)	≈	0	) correspond to black ones, and moderate interactions correspond to grey
	cells. It is noteworthy that issues of these tests have implications for simultaneous
	representation of the row and column. For instance, a glance on Figure 5 indicates that
	the simultaneous representation of rows and columns on Figure 4 would undoubtedly
	make sense if the bivariate moments 1 , 1 λ and 2 , 2 λ were high, while 2 , 1 λ and 1 , 2
						On this figure, non-significant bivariate moments (bilateral test at the threshold
	0.05) correspond to white cells (	Gl	(	m	,	p	)	≈	1	), while highly significant ones

  and keeping only its first four slices { }

		1 U L , U ,	4	. We will then analyse
	the	11 6 × contingency table

  indicates an increasing linear by linear interaction, the significance of this moment confirms that the abundance of Ulva increased with time in the lagoon. The moment

															Λ	Density,	Time	). The most important
	moments are:	ˆ1 , 5 ≈ λ	10	.	2	,	ˆ1 , 1 ≈ λ	8	.	5	and	ˆ10 , 1 λ	≈	. 8	8	. Thus, the most important and
	reliable interactions are associated with 1 B , i.e. the flow of time. Since	ˆ1 , 1 > λ	0
	5	1 ,

  2. suppose E is a vector space equipped with the basis { }

																					e 1 L ,	,	e	m	and F is another
	vector space, equipped with the basis {	f	, 1 L	,	f	p	}	; then	F E ⊗ is a third vector
	space, equipped with the basis { e	i	⊗	f	j		:	1	≤	i	≤	m	,	1	≤	j	≤	} p : this perfectly fit
	the bivariate moments tables structure (see figures 5 and 12)
	3. the tensor product of two vectors is straightforwardly extensible to functional
	spaces: if	f	(	y	)	∈	H	1	and	2 g ∈ , the tensor product ) ( H x	f	⊗	g	∈	H	1	⊗	H	2	of
	these functions is	f	⊗	g	( ) , x y	: =	f	(	y	)	g	(	x	)

  As we saw in Section 2.2.2, this difference can be decomposed (see formula[START_REF] Brun | Resource allocation and proteolyticactivity in the seagrass Zostera noltii Hornem. subjected to light reduction by Ulva canopies[END_REF]) in a sum of contributions of general form:

	recurrence formula [14], while the score	s J	( j	)	plays the part of the ordinate for
	obtaining the	A m	( y	)	. Consequently, it is natural to define the polynomial	A m	⊗	B	p	( ) x y ,
	on the rectangle [ s	J	), 1 (	s	J	(	J	)	] [ s	I	), 1 (	s	I	(	I	)	]
	J P P ⊗ . ( I p m , λ	I P	⊗	J P	) ( o	m A	⊗	B	p	)	.
	Notice that in that expression, the vectors m A and p B actually stem from the sampling
	of othonormal polynomials. The column polynomial p B is obtained from a recurrence
	formula [2], where the score	s I	(i	)	plays the part of the abscissa in Emerson's

Acknowledgements

The authors wish to thank the referees for their stimulating comments and suggestions, and Starrlight Augustine for help with the English text.

Appendix

Why tensor products?

, one should think that only regular matrix-vector algebra is necessary for our purpose. Nevertheless, we preferred to use the tensor product notation for three reasons:

1. as Benzécri noticed [4, p. 58], this is a natural notation for bilinear applications,


Mercedes Bonet

MD PhD Mélanie Durox

MSc Béatrice Blondel

PhD Pascal Boileau

Véronique Pierrat

MD PhD Jennifer Zeitlin

Management of mother's own milk for very preterm infants in tertiary-levelneonatal units in the Ile-de-France region in France

Keywords: breast milk, neonatal unit, France, Europe

Dear Editor:

Breastfeeding rates among very preterm infants at discharge from hospital are low in the French Ile-de-France region (24%) when compared to those in other European regionswhere rates varied from 19% to 70%. [START_REF] Bonet | Variations in breastfeeding rates for very preterm infants between regions and neonatal units in Europe: results from the MOSAIC cohort[END_REF] Variations in practices related to the management of breast milk in European neonatal units may explain some of these differences.We sought to describe procedures for management of mother's own milk for infants born very preterm in the Ile-de-France region and to assess how these relate to existing French regulations and guidelines.

Two regulations and one guideline related to the use of mother's own milk for preterm or hospitalized infants were in force during the study period (2009-2010) [START_REF]relative au don de lait personnalisé d'une mère à son enfant hospitalisé et rappel des dispositions en vigueur en matière d'allaitement maternel[END_REF][START_REF]de l'article L. 2323-1 du code de la santé publique : Règles de bonnes pratiques de collecte, de préparation, de qualification, de traitement, de conservation, de distribution et de délivrance sur prescription médicale du lait humain par les lactariums[END_REF][START_REF]Recommandations d'hygiène pour la préparation et la conservation des biberons[END_REF] .The content of these texts differswith respect to their target populations, time recommended for breast milk storage and recommendations for preventing infection transmission.One regulation specifically regulates breast milk administration in the neonatal unit, [START_REF]relative au don de lait personnalisé d'une mère à son enfant hospitalisé et rappel des dispositions en vigueur en matière d'allaitement maternel[END_REF] and allows the use of fresh breast milk without bacteriological test if milk is stored for less than 12 hours. The two otherstarget human milk banks [START_REF]de l'article L. 2323-1 du code de la santé publique : Règles de bonnes pratiques de collecte, de préparation, de qualification, de traitement, de conservation, de distribution et de délivrance sur prescription médicale du lait humain par les lactariums[END_REF] or handling and storage of milk bottles in hospitals, childcare centers and homes. [START_REF]Recommandations d'hygiène pour la préparation et la conservation des biberons[END_REF] They recommendedthe storage of fresh breast milk up to a maximum of 48 hours, [START_REF]de l'article L. 2323-1 du code de la santé publique : Règles de bonnes pratiques de collecte, de préparation, de qualification, de traitement, de conservation, de distribution et de délivrance sur prescription médicale du lait humain par les lactariums[END_REF][START_REF]Recommandations d'hygiène pour la préparation et la conservation des biberons[END_REF] performance of bacteriological tests the first time infants receive their mothers' own milk in the neonatal unit anduse of mother's own pasteurized milk or donor milk if maternal serology is positive forhuman cytomegalovirus(HCMV). [START_REF]Recommandations d'hygiène pour la préparation et la conservation des biberons[END_REF] Data on neonatal unit practices came from semi-structured face-to-face interviewswith the health care professionaldesignated by the head of unit as the person with the most knowledge about breastfeeding,in eight tertiary-level neonatal units in 2009-2010.Theseunits accounted for approximately half of the very preterm admissions in the 17 tertiary-level neonatal units in the region. Units were purposely chosen to ensure geographical variability (Paris and suburbs), diverse socio-cultural characteristics of the mothers (origin and social class) and variation of breastfeeding rates at discharge from the neonatal unit (high, mediumor low). [START_REF] Bonet | Variations in breastfeeding rates for very preterm infants between regions and neonatal units in Europe: results from the MOSAIC cohort[END_REF] No unit refused to participate in the study.Interviewees were two lactation consultants, two nurses in charge of the pasteurization unit, two head nurses, two nurses and one pediatrician. Written summaries were made of all interviews. We extracted information about the procedures for mother's own milk management: type of milk (fresh, frozen-thawed or pasteurized), storage time, bacteriological analysis, restrictions on the use of breast milk and availability of written protocols. This study was exempt of ethical approval according to French regulations.

Table 1 shows wide variability in procedures implemented by the units. Half used only pasteurized milk (mother's own or donor milk) before infants reached 32 weeksand no unit used frozen-thawed milk. All units required maternal HCMV serology before giving fresh milk for very preterm infants and used pasteurized milk if the mother had a positive HCMV serology, although limits for use of pasteurized milk were not the same. The majority of the units performed systematic and regular bacteriological analyses of fresh milk, but at different frequencies. Two units performed bacteriological analysis only if the mother or the infant presented signs of infection. Maximum storage time for fresh milk varied from 12 to 48 hours between units. Seven units authorized breast milk expression at home or at the neonatal unit, but one unit used fresh breast milk only if the mother expressed her milk on-site.Four units had an on-site facility, in the hospital or the neonatal unit,for pasteurizing breast milk.

Two units had a written protocol for the use of breast milk for infants born before 32 weeks at the moment of the study and the latest versions dated from August 2008 and April 2009.

This study contributes to the growing literature showing wide variations in the procedures for managing mother's own milk in European neonatal units. [START_REF] Buxmann | Feeding of very low birth weight infants born to HCMV-seropositive mothers in Germany, Austria and Switzerland[END_REF][START_REF] Cossey | The use of human milk in the neonatal intensive care unit: practices in Belgium and Luxembourg[END_REF][START_REF] Omarsdottir | Breastmilk handling routines for preterm infants in Sweden: a national cross-sectional study[END_REF] Variations in practices among neonatal units in Ile-de-Francemay be due to the fact that the guidelines did not always give the same recommendations, as for instance,for storage time for refrigerated fresh milk. These differences also relate to units' interpretation ofregulations and guidelines. Most units implemented more restrictive practicesrelated to fresh breast milk management, such as milk pasteurization, regular bacteriological analysis and non-use of frozen-thawed breast milk, than required by current guidelines.

The units' more conservative approach could have different explanations. Units might prefer to apply the principal of precaution because of concern about risks associated with infection transmission through breast milk. Indeed, evidence is not clear on bacteriological analysis for the use of breast milk from the mother to her own infant and the optimal time for milk storage.It could be an indication of the lack of updated knowledge of health care professionals on research about breast milk use and the latest recommendations published in France. It might also reflect the absence of active breastfeeding promotion policies in these units. Few units in our study had lactation consultants or written protocols for the use of breast milk for very preterm infants.

The impact of these practices on breastfeeding for preterm and hospitalized infants requires further study. Transfers between units are common for very preterm infants and significant changes in the procedures between units might constitute a barrier to successful lactation and breastfeeding. Indeed, inconsistencies in the information given by health care professionals to mothers might affect their capacity to express or breastfeed in the hospital.

More restrictive practices might also have a negative impact because of delays linked to the pasteurization of milk and of costs related to the need for special equipment and trained staff.

Including data about how these practices affect maternal motivation seems also important.

This study also raises more general questions about adherence to recommendations and highlights the need forstrategies for dissemination and implementation of regulations and guidelines and regular evaluation of clinical practices.

Table 1 . Use of mother's own milk for infants born very preterm in tertiary-level neonatal units in Ile-de-France (2009-2010)

 1

	HCMV: human cytomegalovirus
	a: One unit stored milk for less than 12 hours if diagnosis of any intestinal morbidity
	(necrotising enterocolitis, gastroschisis, etc)
	b: For infants between 30 and 32 weeks
	c: Performed only if infection suspected

Acknowledgments:

We thank the neonatal units, and all neonatal doctors and nurses who participated in the interviews.

Disclosure Statement:

Authors declared no conflict of interest.

Funding: This study was funded by INSERM.

Author's contribution:

MB and MD performed the literature review, developed the interview guide, conducted interviews and wrote the first draft of the paper. BB, PB, VP, JZ contributed to the development of the interview guide and the manuscript. All authors participated in revisions and approved the final manuscript.


Stefan Janaqi

Jorge Aguilera

Meriam Chèbre

Robust real-time optimization for the linear oil blending

Keywords: Blending, Robustness, Polyhedra, RTO, Linear Programing

Robust real-time optimization for the linear oil blending

Stefan Janaqi, Jorge Aguilera, Mèriam Chèbre

Introduction

The oil blending process consist in determining the optimal proportions to blend from a set of available feedstocks or components such that the final product obtained fulfills a set of specifications on their properties. A blending system is typically constituted by three functional subsystems: the scheduling subsystem, the online optimizer and the control subsystem.

The scheduling subsystem is the one in charge of the general refinery production planning. The scheduling subsystem uses a linear program to calculate for a given period (up to a month) the recipes , properties and volumes of reference for a sequence of blends , , , = 1, … , (typically ∈ 10,15 ). These calculations are based on the mean characteristics of components.

Several sources of uncertainty (see below) perturb the process. The online optimizer is then required to update the target recipe which may became sub-optimal, or even infeasible, due to these disturbances in the process. For the first blends of the sequence, the online process fits well with the forecasted , , . Nevertheless, after a number of blends, it happens sometimes that the blending environment differs significantly from the mean characteristics considered by the scheduling subsystem. In these cases (the example presented at the end of the paper is one of them), it is very difficult for the online optimizer to take the right decision. The blending environment is represented in this paper by the blending polytopes (see section 2) that change with time and are different from mean polytopes used by scheduling subsystem.

The feedback is based on measures of the blends' and components' properties gathered by online analyzers. Finally, the control subsystem is in charge to adjust the component's flow rates (the recipe ) in order to attain the current target recipe. Here we focus on the Real Time Optimization (RTO) system formed by the online optimizer and the control subsystem. An RTO loop consist of calculating the optimal recipe * for the new constraints and adjusting the actual recipe by * -. Typically, our control subsystem adjusts the recipe every 5 minutes by calling the online optimizer with different sets of constraints (up to 100 calls in 5 minutes). The total time for a blend is several hours (10-20 hours). A RTO is one that meets this time constraint. This motivates our approach for robust blend models that can be calculated efficiently. Our uncertainty model gives rise to linear models that are efficiently solved.

An important characteristic in the oil blending process is the strict requirements over some properties controlled by environmental, legal and technological standards. If these properties are not satisfied, one can correct the blend to a certain limit by dumping the appropriate additives in it. Otherwise, the blend must be reblended. The additives are too expensive whereas reblending reduces the refinery capacity, so one should consider these blends' bounds as hard constraints which must be satisfied.

On the other hand, oil blending is a complex process where several unknown and uncertain factors affect the blend's properties. In addition to the plantmodel error produced by the linearization of blending laws, there are other sources of uncertainty: measure errors on components and properties caused by instruments' precision; uncertain knowledge of the components' properties due to upstream process variations and uncontrolled blending conditions such as air temperature, humidity, etc. All these are typical uncertainties arising in any RTO problem (see [START_REF] Zhang | Results analysis for trust constrained real-time optimization[END_REF]).

Hard constraints on the blends' properties combined to these sources of uncertainty are the ideal characteristics to apply the Robust Optimization (RO) techniques (see [START_REF] Ben-Tal | Robust Optimization[END_REF]), which are an approach to optimize under uncertainty. Roughly speaking, robust optimization looks for optimal solutions that are strictly inside the optimization domain. The robust optimization models differ in the way they insure and measure the "how much inside the optimization domain?". A strength of the RO techniques is to consider a deterministic and set-based uncertainty model if one has no information on the random process(es). No probabilistic assumption is made over the uncertainty and the solution obtained is optimal for any realization of the uncertainty in a given set. RO techniques seem to be ad hoc to address an uncertain RTO problem where the feasibility is the primary concern. Some previous works have been devoted to find robust solutions to the blending problem. In [START_REF] Hendrix | Finding robust solutions for product design problems[END_REF], a geometric approach is presented for the product and mixture design problems but only considering the uncertainty due to measurement imprecision. To deal with components' properties uncertainty, [START_REF] Singh | Model-based real-time optimization of automotive gasoline blending operations[END_REF] proposed a nonlinear blend RTO system based on predictions of feedstocks properties whereas [START_REF] Wang | Online optimization model design of gasoline blending system under parametric uncertainty[END_REF] presented a chance constraint model and a hybrid neural networks-genetic algorithm solution. More recently, [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF] introduced a linear blending control algorithm which handles this type of uncertainties via an estimator of the components' properties. A more general method based on stochastic programming which covers various types of uncertainty is presented in [START_REF] Zhang | Real-time optimization under parametric uncertainty: a probability constrained approach[END_REF].

Our contribution begins by highlighting the polytopes underlying the blending process. This and the form of uncertainty on components (see below) permit to conceive robust models without the need of convex programming tools. Based on the data, we develop robust models that give rise to linear models. Consequently, a linear programming solver is all we need. With linear programming solvers we meet the real time constraint that is "less than 3s per solver call". Actually, we obtain 0.25s per solver call. This point is of greatest interest in the real time processes under investigation. Another result of interest is the calculation of a feasible volume's interval. As we'll see in the following, a blend may be infeasible because of its volume. It is then useless to try to modify the recipe while the blend volume is out of the feasible interval.

In this paper we study the blending process with linear blending laws and the components stocked in running tanks. The paper's structure follows. Section 2 presents the Blending Polytopes inherent to the Blending Problem and the Basic Equation which will be used to introduce the updated data in the RTO model. Section 3 analyzes the general case when the heel of a previous blend is used in the production of the new blend and a RTO approach is proposed to solve the Blending Problem. In Section 4 we construct uncertainty sets for various types of uncertainty. These sets determine the robust regions for the blending polytopes in which the new robust RTO (RRTO) method is based. Finally, the price of robustness and reblending is introduced and explained through a real numerical example of the RRTO method in Section 5.

Blending Polytopes

Let's denote by and the number of components and properties respectively. We will represent a blend's recipe (recipe in short) by a vector ∈ such that is the percentage of component present in the blend whereas the blend's properties (properties) will be denoted by a vector ∈ . Then the set of components is described by the × matrix where , is the i th physical and / or chemical property of the component j. It is customary when blending to use the heel of a previous blend to produce a new one. The process starts then with a volume of a previous blend with properties ∈ and continues by adding gradually a volume of the new blend to obtain a final product's volume ! "# = + . We suppose that a blend's properties are a linear combination (or more precisely, a convex combination) of its component's properties. This is not true in general, but a number of invertible transformations exist such that the transformed properties blend linearly. So, for any % ∈ 0, the blend's properties corresponding to the recipe are determined by the basic equation:

= & ' + ( ' (2.1)
where, & ' = ) * ) * +, , ( ' = , ) * +, .

Having a blend's volume with properties , equation (2.1) can be used at any time in the process to plan the production of a fixed volume % ∈ 0, . We just update and the matrix and take -% as the remaining volume to pour. For this planning horizon, the optimal recipe is obtained by solving the Blending Problem (BP):

min 1 2 3 ≤ ≤ (2.2) ≤ ≤ (2.
3)

1 3 = 1, 0 ≤ ≤ 1 (2.4)
Here the components' costs are represented by the vector 2 ∈ . A recipe is subject to components availability and hydraulic bounds which impose the component's constraints (2.2). Similarly, the properties constraints (2.3) are determined by properties lower and upper ( ) bounds. If = then we say that property is regulated, otherwise it is controlled. Additionally, properties' bounds are labeled as hard 5 6 , 6 7 or soft 5 8 , 8 7.

Hard bounds are related to legal, commercial and environmental specifications which must be satisfied whereas soft bounds can be violated but incurring into quality giveaway. This means that the blend is better than needed. Finally, (2.4) are the percentage constraints.

We note 9 : = ; ∈ | 1 3 = 1= and the simplex

> : = ; ∈ | 1 3 = 1, 0 ≤ ≤ 1=.
Now, there are 4 distinct polytopes participating in blending process. Geometrically, the intersection of the regions determined by the components constraints and > : defines the recipe's polytope:

? = @ ∈ > : | ≤ ≤ A (2.5)
which represents the set of all recipes that can be produced. On the other hand, the properties constraints define the polyhedron (as it may be unbounded):

B C = D ∈ | ≤ ≤ E (2.6)
Both are defined in the component's space.

In properties space, the image of ? under the basic equation is the possible blend's polytope:

F = ; ∈ | ∈ ?= (2.7)
Finally, the target polytope is defined by:

B = D ∈ | ≤ ≤ E (2.8)
The polytope F depends on the volume %. Let's define F ∞ = ; | ∈ ?= . Then, ∈ F if and only if -= ( ' -. This shows that F = F % is a homothetic transformation of F ∞ with homothetic center and homothetic ratio ( ' . When % grows continuously from 0 to , F % "moves" continuously from to F . Consequently, if B has a non empty intersection with the cone 2H I , F ∞ defined by the apex and the base F ∞ then, there is an interval % ∈ , "J such that F % ∩ B ≠ ∅ . In the following we suppose that 2H I , F ∞ ∩ B ≠ ∅. If this were not the case, one has to change the components (and F ∞ ) in order to assure feasibility.

Any point in F represents the properties of a blend issued from one or more recipes in ? whereas B depicts the target properties to be attained by the blend produced.

We notice that, among the four preceding polytopes, only F and B C depend on %. We omit the % when it is considered to be fixed.

At this point we have identified the polytopes underlying the BP and we observe that for a fixed volume %, a feasible solution of BP is a recipe

∈ ? ∩ B C such that ∈ F ∩ B.
Equivalently, F NO N PIQON RI ⟺ F ∩ B = ∅ ⟺ ? ∩ B C = ∅.

Blending from a previous blend: a RTO method

In any RTO loop, the online optimizer is required to produce a recipe that permits the controller to guide the blending process. In Section 2 we stated that for a fixed volume % to pour, the optimal recipe can be obtained by solving the BP. However, the BP can result to be infeasible because of the uncertainty perturbing the process or for instance, when a previous blend with properties out of specification is used in the production of the new blend. Therefore, a method to compute at any time the best recipe for the planning horizon is required. Later we present this method but first we discuss the BP feasibility when the heel of a previous blend is used.

If the blend is produced from scratch or the volume to pour is too big compared with the heel's volume, that is, if = 0 or % → ∞ then the basic equation reduces to = and the polytope F % does not depend on . However, the fact of using the actual blend's properties in the basic equation to update the model inside the RTO loop, makes blending from a previous blend the regular case. As explained in the preceding section, we suppose that there is an interval of volumes U , U"J such that F % ∩ B ≠ ∅, ∀ % ∈ U , U"J . The BP is feasible only for these volumes.

Having the values U and U"J helps to select the appropriate planning horizon at each RTO loop. Choosing a planning horizon by taking a volume % < U produces an infeasible BP and makes it necessary to generate an alternative recipe with a possible deterioration of the overall performance of the RTO method. On the other hand, U"J is the maximum volume to blend when looking for a blend within specifications. After U"J has been poured, F % moves away from B and the blend's properties deteriorate progressively. Furthermore, knowing U and U"J helps the control system to reduce unnecessary and inefficient interventions. The latter have as purpose to correct the blend to be within specifications but often they are based on a limited view of the problem. In order to compute the interval U , U"J , let's rewrite the basic equation (with = At any RTO loop, we proceed through the steps I. to VI. as follows:

I. Incorporate the newly available information , , , in the basic equation. If = 0, solve problems 3.3, 3.6 and 3.7 in this order until having a feasible solution, then STOP. Otherwise, go to step II. II. Compute U and U"J by solving the linear programming problem:

max min 1 * ,1 , ∈ ? ∩ B C , ≥ 0 (3.2)
If Problem 3.2 is infeasible or if ∉ U , U"J go to step IV. Otherwise, a blend without quality giveaway exists for any volume % ∈ U , U"J . Choose a blending horizon % ∈ U , U"J and go to step III.

III. Compute the optimal recipe:

min ^23 , ∈ ? ∩ B C (3.3)
Notice that solving Problem 3.3 with % free and the additional constraint:

U ≤ ≤ U"J (3.4)
we obtain at the same time the optimal blending horizon % * ∈ U , U"J and the optimal recipe for this volume. This is the best choice if it's not imperative to produce a particular volume. If this is the case STOP. Otherwise go to step IV.

IV. Forget the components costs and focus only on the hard constraints. Compute the interval U , U"J by solving the problem:

max min 1 * ,1 , ∈ ? ∩ B 6C , ≥ 0 (3.5)
If this problem is infeasible or if ∉ U , U"J go to step VI. Otherwise, go to step V.

V. Look for a recipe ∈ ? ∩ B 6C and a volume % * ∈ U , U"J producing the blend with minimal quality giveaway. To do this, solve the following problem with fixed % from a finite set of values, % ∈ ;% : , … , % _ = ⊂ U , U"J :

min 1, ‖ -b ‖ : , ∈ ? ∩ B 6C , b ∈ B C (3.6)
VI. Find ∈ ? to produce a volume of blend with properties as near as possible to the ideal blend 3 :

min 1 ‖ -3 ‖ : , ∈ ? (3.7)
Following these steps, the RTO method always produces a recipe that guides the control process. In step I, = 0 and solving problems 3.3, 3.6 and 3.7 we obtain the minimal cost recipe, the blend with minimal quality giveaway or the closest blend to 3 . If BP is feasible for some volume % ∈ U , U"J then a blend without quality giveaway exists and the method generates the one of minimal cost (step III). Otherwise, the hard constraints become the priority and the method searches for a blend satisfying them while minimizing the blend's quality giveaway (step V). Finally, if there is no blend satisfying the hard constraints, the recipe producing the blend with properties as close as possible to 3 is proposed (step VI). We finish this section by stressing the possible application of the RTO method to determine the appropriate heel's volume to use in the blend. In our analysis, is considered as fixed but we can solve the Problem 3.3 with the additional constraint 3.4 with and free in order to obtain the cheapest recipe * for some * ∈ U , U"J . Then we can find (by means of the relation

= X * '
) a suitable pair , . The choice of norm c : , in problems 3.6 and 3.7 permits to obtain blends violating a minimum number of properties. Moreover, it yields LP problems which can be solved efficiently.

Robust RTO

In the previous section we proposed a blending RTO method based on blending polytopes and its evolution with the blended volume. Looking to reduce the model deviations produced by some uncertainties, the method uses blend's and components' properties updates to feedback the model via the basic equation. However, model updating may fail to guarantee even a feasible solution. A main reason of this failure is the implicit assumption that data remains unchanged inside each RTO loop.

For instance, when online blending is used, values fluctuate with time because components are issued from different process presenting also perturbations. To address this problem, [START_REF] Singh | Model-based real-time optimization of automotive gasoline blending operations[END_REF] proposed a blending RTO method which updates the model with predictions of the components' properties. Although this method improves the model accuracy, it continues to be non robust as it depends on the quality of the predictions. Moreover, uncertainty in the blending process affects other factors than values as we will see below.

BP Uncertainty

In accordance with [START_REF] Zhang | Results analysis for trust constrained real-time optimization[END_REF], uncertainty in any RTO system may be of four types:

1. Process uncertainty: components properties, temperature, humidity, etc. 2. Measurement uncertainty. 3. Model uncertainty. 4. Market uncertainty: components availabilities and prices, blends demands, etc.

In this work we consider the uncertainties arising from components measurement and blend's properties measurement (type 2) and the uncertainty caused by imprecise knowledge of the components properties (type 1).

Measurement and components properties uncertainties manifest geometrically in different ways: for the first type, the real recipe (the real percentages of each component used in the blend) and its properties may differ from the nominal ones. So, the real recipe and its properties are located in neighborhoods of the nominal recipe and its properties respectively. For the second type, the real matrix d differs from the nominal matrix and hence the real polytope B e C is different from the nominal polytope B C . In both cases, when a nominal feasible recipe ∈ ? ∩ B C is computed, the real recipe may lie outside the polytope B C and the real blend's properties may be outside the polytope B. Then the real recipe results to be infeasible.

An intuitive idea to fix our RTO method against measure and uncertainties follows from the previous geometric information. The idea consists in computing for polytope B C its convex robust region B C such that any point in this region resists to uncertainty and to measurement errors. That is, any point in B C is guaranteed to remain inside the real polytope for all possible realizations of and any measurement error, within reasonable limits. So, the robust RTO method will consist in replacing polytope B C by B C in the RTO method from Section 3. To develop this idea, first we need to model and measure the uncertainty we would like to be protected against and then to compute the appropriate robust region of B C . This development follows the Robust Optimization (RO) theory developed by Ben-Tal et al (see [START_REF] Ben-Tal | Robust Optimization[END_REF]).

4.2

Robust Regions and Robust RTO Let be a nominal recipe and its properties. To model the components and properties measurement uncertainties, we suppose that the real recipe f lies in the ball > , g of radius g and center whereas the real blend's properties f lies in the ball:

> , g = ; ∈ | -g ≤ ≤ + g =
We are given the minimal and maximal values of . In order to model matrix uncertainty, we use interval sets. That is, we suppose that each real value d , is comprised in the interval h , = i , -j , k , , + j , Y l around its nominal value , for some positive values j , k , j , Y .

Here we could use different sets and any norm to model the uncertainty. The level of conservatism (how much we want to be protected against uncertainty) and the problem complexity depend on these selections. Taking interval sets and the max norm, the robust regions obtained are polytopes and the complexity in the model is preserved at the expense of being probably too conservative (we are protected from the worst deviations of all values and from the biggest measurement errors occurring all at the same time). Now we proceed to construct the robust region of polytope

B C . Any point ∈ B C is robust regarding uncertainty iff ≤ f ≤ . That is, iff for any d such that d , ∈ h , , ≤ & ' , + ( ' d ≤ , N = 1, … , (4.1)
holds. Notice that any row d can be expressed parametrically as

d m = + m k3 n k + m Y3 n Y , with, n k = oNQpqj ,: k , … , j , k r, n Y = oNQpqj ,: Y , … , j , Y r, m k = min 0, m , m Y =
max 0, m for some m ∈ such that ‖m‖ s ≤ 1. Therefore we can deduce that ∈ B C is robust regarding uncertainty iff

+ ( ' j k ≤ ≤ -( ' j Y , N = 1, … , (4.2)
In order to derive equation (4.2) let's define for a given , P m = m k3 n k + m Y3 n Y . We then resolve the following optimization problem:

max t P m , ‖m‖ s ≤ 1.
The optimal m saturates the constraints. Thus, for any = 1, … , one has:

either m = 1 = m Y or m = -1 = m k . The optimal value is P m = -∑ j , k ,v * wk: + ∑ j , Y ,v * w:
. This is clearly less than or equal to ‖n Y ‖ : = j Y (because of positivity of j Y and ) and this gives the right side inequality of (4.2). By resolving min t P m , ‖m‖ s ≤ 1 we deduce the left side inequality of (4.2) with a similar reasoning. This is a particular case of a duality result: when resolving max t Qm , ‖m‖ x ≤ 1 , the optimum is ‖Q‖ y , where ‖ * ‖ x , ‖ * ‖ y are dual norms with z { Y z | w:.

From now on we denote by } k = ( ' j k and } Y = ( ' j Y .

In addition to uncertainty, ∈ B C resists also to components' measure uncertainty if and only if any point in > , g satisfies (4.2). That is, iff

+ } k ≤ + b ≤ -} Y , N = 1, … , (4.3)
holds for all b ∈ such that ‖b‖ s ≤ g . Then, computing the minimum and maximum on b we get that (4.3) holds iff

+ } k + g ≤ + b ≤ -} Y -g , N = 1, … , (4.4)
with g = ( ' g ‖ ‖ : . The reasoning follows the same lines as the deduction of inequalities (4.2).

Equivalently, any recipe ∈ B C is robust regarding uncertainty and properties' measure uncertainty iff

+ } k ≤ + ~ ≤ -} Y , N = 1, … , (4.5)
holds for all ~∈ such that |~ | ≤ g . As previously, computing the minimum and maximum on ~, (4.5) holds if

+ } k + g ≤ ≤ -} Y -g , N = 1, … , (4.6)
Finally, letting ∆ = max g , g we obtain the robust region B C of polytope B C :

B C = D ∈ | + } k + ∆ ≤ ≤ -} Y -∆ , N = 1, … , E (4.7)
Any ∈ B C resists to uncertainty and to components and properties' measurement uncertainties.

To summarize, let's consider a nominal feasible recipe ∈ ? ∩ B C . If f ∈ > , g and f ∈ > , g for some g ∈ Y and g ∈ Y and if , ∈ h , , then the recipe f will be feasible in reality: f ∈ ? € ∩ B e C .

The RTO method proposed in Section 3 transforms then in a robust RTO method by a simple substitution of B C with B C . This reduces the impact in implementation as the structure of RTO remains the same. As and may be considered as free variables in the RTO method, we can describe the polytope B C by using explicitly Equations 4.4 and 4.6 in association with the identity ( ' = z z+• * instead of Equation 4.7 where a Max is involved. This robust RTO method depends completely on the robust region B C and to obtain it we only need to determine the values of j k , j Y , g , g . This study lies upon the fact that this information is available. It's worthwhile to note that while g and g are considered as fixed values independent of the RTO loop's length, j k and j Y depend on it. As fluctuations on may accumulate over time, the longest the loop's length is, the biggest these fluctuations can be.

Here we limit the analysis to measurement and components properties uncertainties. However, other types of uncertainty manifest geometrically in the same way and thus can be treated identically. For instance, when the uncertainty is due to uncontrolled factors like temperature or humidity, the real blend's properties are located in a ball around the nominal blend's properties. To model the uncertainty in the components prices, we can transform the optimization problem to one with certain objective function and such that uncertainty appears as a constraint (2 3 ≤ 2). Then we can construct the uncertainty sets and determine the way it affects the robust regions. We can proceed similarly for the uncertainty in the components availabilities which affects the robust region ? of polytope ?.

The price of robustness and reblending

In this section we present a case study based on real data to illustrate and compare some key aspects of the RTO method and its robust counterpart. The BP consists in producing a fixed volume ! "# = 5000 ' of blend from 8 components and = 2000 ' of the heel's volume from a previous blend. Each component and the previous blend have 7 properties to be controlled during the process and they are represented by the 7 × 8 matrix and vector respectively. Vectors and "J stand for the properties bounds while vector 2 denotes the components' cost. In order to produce a robust recipe, we assume that components and properties measurement errors are bounded by g = 0.12, 0.0003, 4.5, 0.0826, 0.028, 0.000049, 0.2 and g = 0.01 respectively. Regarding uncertainty, we dispose of k and Y the absolute lower and upper bounds of matrix . Let's define, ‹ k = -k , ‹ Y = Yand ‹ = Yk . We suppose that there are j k and j Y such that -j k ≤ d ≤ + j Y , with j k = min ‹ k , OE‹ and j Y = min ‹ Y , OE‹ for some 0 < OE < 1. As we stated in Section 4.2, the values of g and g are fixed during the process whereas OE is directly related to the RTO loops' length. We take OE = 0.01.

Solving the robust version of problem 3.2 with = 2000 ' , we get U • , U"J• = 0, 0.0362 and the corresponding robust feasible volumes interval: U • , U"J• = 55304, ∞ . The subscript " " will indicate the result of a robust version of an optimization problem. The corresponding interval for the nominal (non robust) case is U , U"J = 34587, ∞ . This means that we need to produce at least 34587 ' (55304 ' ) in order to get a (robust) blend within specifications which uses completely . Taking only the hard constraints, we obtain similar intervals.

If we decide to produce ! "# = 5000 ' using = 2000 ' then we solve the problem 3.7 and we obtain a recipe with cost 90.72 ŽI• ' but producing a blend out of specifications. Actually, this blend violates only one property's bound but by more than 80%. Instead of this, We compared this performance (in terms of CPU time) with the classical quadratic model of robustness [START_REF] Ben-Tal | Robust Optimization[END_REF]. This model gives rise to convex programming. We used the CVX toolbox with Matlab (see [START_REF] Grant | CVX User's Guide[END_REF]). The CVX solver took about 15s at each call and this is much above our time constraint (no more than 3s per call).

Then we generate the optimal nominal recipes when = 168.44 ( U"J for OE = 0.01) (with reblending) and = 0 (without reblending) are used in the production of 5000 ' of blend. The cost of these recipes are 88.29 and 87.04 respectively. This is normal in general as a "bad" heel is difficult to correct. We observe that for this case, the recipe's cost with reblending is ''.'"k'". • '". • =1.43% greater than the nominal recipe's cost whereas the relative price of robustness is of only 0.35% for a significant value of OE = 0.

Conclusions

In this paper a robust real time optimization method for the linear oil blending process has been introduced. The method is based on the RO techniques and it is intended to avoid reblending while minimizing the blend's cost and the quality giveaway. We constructed a set of models for different types of uncertainty arising in the blending process. The simplicity of these models may produce over conservative solutions (blends too expensive) but we showed via an example the convenience of the RO techniques even for these simple models.

A main characteristic of our RRTO method is the integration of the case when the heel of a previous blend has to be incorporated in the new blend. This feature provides meaningful information for the control system, for instance, to determine the appropriate heel's volume to use in the blend or to have an estimate of the volume to pour before getting a blend within specifications.

There are many factors associated with the cost of reblending a previous blend which failed to be within specifications: the tank use, the process time, the inventory costs. Another one is the recipe's cost increase for using the heel's volume of the previous blend. This cost depends (obviously) on the particular blend in which it is used. In this work we compared this cost of reblending with the price of robustness to stress the convenience of the RO techniques.

The results obtained comfort the idea that if reblending or additives cost is expensive enough, then using a more conservative technique like RO may improve the global performance of the blending process. The recipe cost rises by taking the robust recipe in place of the nominal one but reblending and additives expenses cancel. More experiments need to be realized in order to estimate the impact of RTO with a large panel of scenarios.

It appears as a good perspective (kindly suggested by referees) to deal with variable g and g especially when longer blending horizons are to be considered. Also, Monte Carlo simulation with variable could be a way to compare both methods.

  constraints remain linear. Now we present the RTO method:Blending RTO MethodLet's denote by B 6C = D ∈ | 6 ≤ ≤ 6 E the polyhedron defined by the hard constraints and let 3 ∈ B 6C an "ideal" target blend without quality giveaway provided by the scheduling subsystem.

Table 1 :

 1 Relative price of robustness for different levels of uncertainty.

	1.

Acknowledgements

We thank the referees for their suggestions and remarks.

we can compute the biggest heel's volume allowing us to produce 5000 ' of robust blend. This is the decision taken in practice. From relations U"J• = ) *

) and ! "# = + we obtain:

These 174.5 ' are far below the nominal = 2000 ' forecasted by the scheduling system. The scheduling system uses a linear program to calculate every month the recipes, properties and volumes of reference for a sequence of blends , , , = 1, … , (typically ∈ 10,15 ). For the first blends of the sequence, the online process fits well with the forecasted , , . Nevertheless, after a number of blends, it happens sometimes that the overall blending environment differs significantly from the mean characteristics considered by the scheduling system. In these cases (the present example is one of them), it is very difficult for the online optimizer to take the right decision. Our decision was based on the criterion "minimize reblending".

Next, fixing = 174.5 ' and = 4285.5 ' we solve the robust version of problem 3.3 to obtain the optimal robust recipe Therefore, producing a robust recipe induces a cost increase of 2.68%. However, we observe that if we take = 0 (no reblending), then the robust recipe cost is 87.70 and the nominal recipe cost is 87.04 producing a cost increase of only 0.76%.

From these results we are interested in comparing the price of robustness with the reblending cost (the cost difference between the recipes obtained when the heel's volume is used and when it is not). In order to provide a fair comparison, the price of robustness is obtained by taking = 0 (no reblending involved) and the reblending cost from the nominal recipes (no robustness involved).

To compute the price of robustness we conduct a blending simulation over 36 RTO loops of 2-hours length (10% of the total blend time). In Table 1 we show the average recipe's cost over the 36 periods and the percentage increase in cost (∆2) from the nominal recipe to the robust recipe for different OE values and taking = 0. We used the GLPK solver on an INTEL, 2 CPU, 32bits, 2GHz, 4Gb computer. The execution time for each call was about 0.25s.