A Dietz

M Schirlin

C Azzaro-Pantel

L Pibouleau

S Domenech

MULTICRITERIA OPTIMIZATION OF BATCH PLANT DESIGN AND RETROFIT: APPLICATION TO AN INDUSTRIAL FINE CHEMISTRY CASE STUDY

Keywords: multicriteria optimization, batch plant design, batch plant retrofitting

. The case study is based on an industrial multiproduct batch plant. Both the monocriterion and multicriteria design and retrofit were implemented. The results showed that the actual batch plant configuration works nearly at its production capacity limit. Besides, the monocriterion study reveals that a batch plant of lower cost operating;:ç in the same conditions can obtain the same production rate. Concerning the retrofitting case, the aim was to adapt the actual batch plant configuration to an increased demand while minimizing investment and adding as much as possible identical equipment items to the actual configuration. The previous developed methodology was applied and several compromise solutions were found.

INTRODUCTION

The design of multiproduct and multipurpose batch plants is a key problem in chemical engineering [START_REF] Biegler | Systematic Methods for Chemical Process Design[END_REF]) and much work has been presented in this area in the past few years, for instance [START_REF] De-Ming | Optimal Design of Batch Plants with Uncertain Demands Considering Switch Over of Operating Modes of Parallel Units[END_REF] [START_REF] Chunfeng | Ants Foraging Mechanism in the Design of Multiproduct Batch Chemical Process[END_REF], [START_REF] Harish | Optimal Reliable Retrofit Design of Multiproduct Batch Plants[END_REF], [START_REF] Heo | A New Algorithm for Cyclic Scheduling and Design of Multipurpose Batch Plants[END_REF], [START_REF] Montagna | The optimal retrofit of multiproduct batch plants[END_REF], [START_REF] Cavin | Multi-objective process design in multi-purpose batch plants using a Tabu Search optimization algorithm[END_REF], [START_REF] Tânia | Optimal design and retrofit of batch plants with a periodic mode of operation[END_REF]. The formulation of batch plant design generally involves mathematical programming methods, such as LP (Linear Programming), NLP (Non-Linear Programming), MILP (Mixed-Integer Linear Programming) or MINLP (Mixed-Integer Non-Linear Programming). The main drawback of this methodology is the difficulty, even impossibility, to describe with a high degree of sophistication, the real constraints (various storage policies or human operators shift, for instance …). In other cases, the number of equations to take as constraints often renders the problem impossible to solve. An alternative was proposed by Bérard et al. (2000), consisting in coupling a Discrete Event Simulator (DES) in order to evaluate the feasibility of the production at medium term scheduling, with a master optimization procedure based on a Genetic Algorithm (GA). The optimization variables take discrete values and the problem presents a marked combinatorial feature (the equipment sizes are considered as discrete values). These ideas were generalized in [START_REF] Dietz | Integrating Environmental Impact Minimization Into Batch Plant Design: Application To Protein Production ESCAPE-14[END_REF] within a framework for multicriteria design and retrofitting taking into account both the batch plant structure and the operating conditions as optimization variables. The choice of a hybrid method GA/DES was then all the more justified as several criteria were simultaneously taken into account. This paper is divided into three sections. The first is dedicated to the presentation of the multiproduct batch plant and to the results obtained with the discrete event simulator. In the second section, the multicriteria optimization procedure is presented. Then, the results of the optimal design and retrofitting are analyzed. Finally, the conclusions and perspectives of this work are given.

SIMULATION OF A MULTIPRODUCT BATCH PLANT

The case study is an industrial multiproduct batch plant producing four fine chemical products. For confidentiality reasons, only two of them are considered here. A simplified formulation was obtained from the detailed information about the production recipes in order to be implemented in the discrete event simulator. Figure 1 presents the manufacturing sequence for the two considered products, respectively A1 and B2. Each recipe involves five treatments steps. Only two of them (R3 and R6) are shared by the two products considered. The so-called time and size factors of each stage are presented in Table 1. Let us recall here that a size factor is defined as the real batch size at the i stage by kg of product obtained at the end of batch plant treatment. The time factor T i represents the treatment time at the treatment stage I and can be considered as constant in the case of pure discontinuous stages (i.e. reactor) or can be a function of the batch size and/or the equipment item size that performs the treatment task, as it is the case for semi-continuous devices (i.e. filter). It must be emphasized that the size factors have the same order of magnitude for all the treatment stages and are very similar between the two products to manufacture. Concerning time factors, the limiting value corresponds to the synthesis stage in both cases. Since the products to manufacture have similar production recipes both from time and size factors viewpoint, it is interesting to produce them in the same batch plant sharing the available resources. According to this approach, the batch plant design problem has been formulated as a multiproduct batch plant sharing two treatment stages, R3 and R6, manufacturing two products, A1 and B2 (as seen in Figure 2). The hypothesis of an equipment item at each treatment stage has been imposed, which is usually considered at the batch plant design phase. These recipes and data were thus embedded into the Discrete Event Simulator previously developed (see [START_REF] Dietz | Integrating Environmental Impact Minimization Into Batch Plant Design: Application To Protein Production ESCAPE-14[END_REF]. The classical expression C[€] = a + b V[m 3 ] c was used to compute the equipment item costs, where C is the equipment item cost, V the equipment item volume and a, b et c are constants. The c coefficient was set at 0,6 as it is usually assumed in investment cost computing. The a constant was fixed as 10000 €, which represents the fixed cost of the equipment item. A value of 20000 € was adopted for the b constant , taking as reference the cost of a 16 m 3 equipment item, that is around 122000 €. Table 1 presents the required production for each product, A1 and B2 respectively. Two production policies were validated with plant data, i.e. mono and multiproduct cases. Figure 3 presents the Gantt chart corresponding to a monoproduct campaign policy for four batches of both products A1 and B2. The equipment items show a regular use and a relatively high use rate. This behaviour comes from the design stage of the batch plant where only this production policy was considered. The Gantt diagram corresponding to a multiproduct policy is presented in Figure 4. As the time factors corresponding to the reaction exhibit the highest values for both products, this stage becomes the limiting one. The other stages have a relatively low use rate, even though regular, since they treat only one of the two products, which suggests that oversized equipment is involved as far as the multiproduct policy is considered.

Figure 3 -Gantt chart for a monoproduct campaign production policy for product A1 and B2 respectively.

Figure 4 -Gantt chart for a multiproduct campaign production policy.

A TWO-STAGE METHODOLOGY FOR PLANT DESIGN AND RETROFIT

-Basic concepts

The approach used consists in coupling a stochastic algorithm, indeed a Genetic Algorithm (GA) with the Discrete Event Simulator (DES). The management of production constraints and the computation of the objective functions are carried out by the DES. The master MultiObjective Genetic Algorithm (MOGA) procedure aims at improving the performance criteria and involves a Pareto sort procedure [START_REF] Dietz | Integrating Environmental Impact Minimization Into Batch Plant Design: Application To Protein Production ESCAPE-14[END_REF].

Figure 5 -Batch plant design framework.

References to multiobjective optimization could be found in [START_REF] Bhaskar | Applications of Multiobjective Optimization in Chemical Engineering[END_REF], [START_REF] Coello | An Updated Survey of GA-based Multiobjective Optimization Techniques[END_REF], [START_REF] Ehrgott | Lecture Notes in Economics and Mathematical Systems -Multicriteria Optimization[END_REF]. Lately, there has been a large development of different types of multiobjective genetic algorithms, which are reflected in the literature. For an overview on genetic algorithms in multiobjective optimization, see [START_REF] Fonseca | An overview of evolutionary algorithms in multiobjective optimization[END_REF]. Literature surveys and comparative studies on multiobjective genetic algorithms are also given in [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF]) [START_REF] Bhaskar | Applications of Multiobjective Optimization in Chemical Engineering[END_REF], [START_REF] Coello | An Updated Survey of GA-based Multiobjective Optimization Techniques[END_REF]. Fonseca and Fleming have also divided multiobjective genetic algorithms in non-Pareto (Schaffer, 1985) and Pareto based approaches [START_REF] Goldberg | Algorithmes Génétiques[END_REF].

The main advantage of genetic algorithms over other methods, particularly over other stochastic procedures such as Simulated Annealing, is that a GA manipulates a population of individuals. It is therefore tempting to develop a strategy in which the population captures the whole Pareto front in one single optimization run. This approach was adopted here.

-Multiobjective Genetic Algorithm

The flowchart presented in Figure 6 illustrates the cycle {Evaluation, Selection, Crossover and Mutation} which is repeated until a stop criterion is reached. After this cycle, the Pareto sort is applied. The multiobjective aspects are taken into account at the selection stage and the compromise solution search at the cross-over stage. The selection procedure is carried out by the biased Goldberg's roulette, with a roulette for each criterion to optimize. An equal number of individuals for each criterion was selected to complete the total number of individuals passing by the survival procedure to the next population. The cross-over procedure that proposes compromise solutions was not modified. It must be pointed that the population was composed of "good" individuals for each selected criterion thanks to the Goldberg's roulette and that the individuals were chosen through a random procedure. This allows the crossover of "good" solutions for a criterion with "good" solutions for another with a strong probability to generate a compromise between both criteria. In the case that both "good" chosen solutions correspond to the same criterion, the cross-over procedure will carry out the traditional function of generating a better solution than the two previous ones. The mutation procedure is not modified and its aim is, as usual, to diversify the search and to avoid local optimum solutions.

The aim here is to propose a generic multiobjective genetic algorithm able to evolve naturally towards the whole set of optimal Pareto solution. This evolution must be done from an initial population generally randomly generated composed of individuals not adapted to the considered criteria.

The genetic algorithm includes the following steps (see Figure 6):

Step 1: an initial population is generated randomly. This procedure guarantees a diversified initial population covering the complete space search.

The following steps are performed in order to pass from the actual population (k) to the next one (k+1). First an intermediate population is generated in steps 2 to 4.

Step 2: an elitism procedure selects the best individual for each criterion and is placed in intermediate population.

Step 3: an equal number of individuals is chosen for each criterion using the biased Goldberg's roulette.

Step 4: the intermediate population is completed with individuals generated by the cross-over procedure. Let us note that the individuals to whom this procedure is applied are chosen randomly from the population k. The simplest cross-over procedure is used with only one cutting point.

Step 5: The mutation procedure is applied to a fixed number of randomly chosen individuals. Only one point of the chromosome is modified, changing its value from 0 to 1 or the opposite.

Step 6: the new population becomes the current one and Steps 2 to 5 are repeated until the maximal number of generations is reached.

Step 7: A Pareto sort procedure is carried out over all the individuals evaluated over generations. At the end of the algorithm, the set of Pareto non dominated solutions is obtained.

BATCH PLANT DESIGN

Monocriterion Batch Plant Design

First, the batch plant design was carried out in order to validate the approach. The results obtained showed that a batch plant with an investment cost of around 75% of the actual one can reach the same production rate (Figure 7). The investment reduction is achieved over all stages except the reaction stage (R3), which has the longest operating time. Both solutions propose the same capacity for this stage, using yet two equipment items half sized instead of one. This gives more flexibility to the batch plant. Since the first stage size is reduced, the batch size is also reduced, and, consequently, the downstream equipment item sizes, R21, R6, R5, F, S.

The configuration obtained with a multiproduct policy is 5% lower than the monoproduct one. The difference between both configurations lies in the spin-dryer stage sizing. When a monoproduct policy is proposed, a bigger spin-dryer is necessary for treating the batch in a reduced time, thus releasing the filter (F) and R4 equipment items available for production purposes and reducing the batch storage time in these equipment items.

It is interesting to visualize the weight of each treatment stage on the total investment cost criterion. Figure 7 presents the contribution of each treatment stage to the total investment cost for the current batch plant configuration and the two configurations obtained by the batch plant design framework. In all three cases, the contribution to the global cost for the different stages is more or less equivalent. The most expensive equipment item is R5, i.e., a distillation column, which presents the highest values for cost coefficients. The main difference between the actual configuration and the two proposed configurations lies in the reaction stage, R3, which contribution percentage is doubled for both new configurations. The decrease in investment cost for batch plant design framework was observed in most of the treatment stages, except the reaction stage, which has the highest time factor. This behaviour was already highlighted in the simulation phase that was carried out to validate the model.

Figure 7 -Contribution of all stages to the total investment cost.

The batch plant configuration for each solution is illustrated in Figure 8. The first remark concerns the synthesis stage, R3: both proposed configurations have the same capacity as the actual configuration, which means that it is the limiting stage; but they propose to use two smaller size reactors instead of one middle-sized item, giving some flexibility to the batch plant. The decrease in the size factor of the first treatment stage involves in turn a decrease in the batch size and in the equipment item size for the downstream process, i.e., R21, R6, R5, F, S. In the actual configuration, let us note that an equipment item, R6, was oversized because it was initially designed to manufacture a product that was finally discarded in the industrial production step.

Figure 8 -Monocriterion batch plant design results.

Bicriteria Batch Plant Design

The bicriteria batch plant design problem based on an {investment cost -production capacity} approach was then investigated while keeping the previous design formulation with a fixed time horizon and production. The production capacity criterion is taken into account through the necessary time to manufacture the imposed production. If the production is manufactured in half time horizon, the batch plant production capacity in the complete horizon is therefore twice the fixed production. The investment cost is a criterion to minimize while the production capacity must be maximized. Parameter setting of the MOGA results from a preliminary study and is presented in Table 3. The results obtained for both production policies are presented in Figure 9. Three optimisation runs were carried out for each case where it can be seen that the multiproduct policy is slightly better, which assesses once more the interest of a multiobjective batch plant framework. At each implementation, the framework was able to find the best solution for each criterion considered as well as a complete set of compromise solutions.

BATCH PLANT RETROFIT

-Monocriterion Batch Plant Retrofit

In the retrofit application, the current configuration is now supposed to adapt to a new market demand, i.e., 2.5 times the actual demand for both products, A1 and B2. The design problem was then extended to treat the retrofit case and can be formulated easily as the number of equipment items and their corresponding size to find, while optimizing the various criteria. Figure 10 shows schematically the equipment items that have to be added to the actual configuration for the two production policies considered. The results obtained have the same investment cost for the shared stages, R3 and R6. The main difference concerns the stages that are used by only one of the products. The advantages of using a multiproduct policy are highlighted, because no additional equipment items are needed at the non-shared stages. In order to explain briefly these differences, the shared stages must treat a higher demand in the same time, which leads to a need in resources. With a multiproduct policy, the non-shared stages, are used over all the time horizon; with a monoproduct policy, they are used only by the product being manufactured, i.e. over only part of the time horizon and, consequently, more resources are necessary to treat the same production in a shorter time. Figures 11 also presents the weight contribution of the additional investment cost for each treatment stage for both production policies, mono and multiproduct respectively. The monoproduct production policy needs a higher additional investment of about 50% than the one corresponding to the multiproduct policy. An investment of the same order of magnitude is necessary for the non-shared treatment stages for a monoproduct policy. Concerning the shared stages, R3 and R6, a higher level of investment is needed for R3 treatment stage for both policies, which can be attributed to the difference in product time factor.

-Bicriteria Batch Plant Retrofit

In the previous configurations, the added equipment items are not in all cases of the same size as the existing ones. However, it is required, for flexibility and maintenance reasons, that all the equipment items of a given stage be identical, which drives us to a compromise in investment cost. In order to explore the compromise configurations which respect these two criteria, investment cost and different equipment items, the bicriteria retrofitting is carried out. Given the stochastic aspects of the genetic algorithm, three optimization runs were carried out. For the monoproduct policy, the results obtained are illustrated in Figure 13, and the corresponding configurations are presented in Table 4. It is interesting to visualize the dominated Pareto solutions obtained in the second optimization run (S6, S7, S8 of Table 2), which correspond to a local optimum where the main investment is carried out over the A1 production chain with R5, ES and R4 instead of the B2 production chain (ES, F, S) which corresponds to the optimal solutions. Given the few solutions obtained, it is advised to keep them for the final decision procedure, because they could be interesting for other criteria, such as maximal production capacity. Concerning the multiproduct policy, the two expected solutions were found (see Table 4).

CONCLUSIONS AND PERSPECTIVES

In this work, an integrated multicriteria optimization framework is implemented for batch plant design and retrofit. The different stages of the methodology were presented. At first, the DES was implemented and validated with plant data. The AG-DES coupling was validated with the monocriterion batch plant design and the results obtained agree with the actual configuration. Then, the batch plant retrofit was tackled and both the mono and multicriteria cases were treated. Concerning the multicriteria case, the method was able to find the set of non-dominated Pareto optimal solutions as well as some other interesting configurations. Generally, the results highlighted the interest of a multiproduct production policy. Some perspectives could now be to include aspects such as maintenance and uncertainty on the market demand in the batch plant model.

Figure 1 -

 1 Figure 1 -Production recipes.

Figure 2 -

 2 Figure 2 -Multiproduct batch plant formulation.

Figure 6 -

 6 Figure 6 -MuliObjective genetic algorithm (MOGA).

Figure 9 -

 9 Figure 9 -Bicriteria investment cost -production capacity batch plant design: global results.5. BATCH PLANT RETROFIT

Figure 10 -

 10 Figure 10 -Monocriterion retrofitting results with a monoproduct production policy.

Figure 11 -

 11 Figure 11 -Monocriterion retrofitting results with a multiproduct production policy.

Figure 13 -

 13 Figure 13 -Monoproduct multicriteria retrofitting results

Table 1 :

 1 Time and size factors.

	Product A1	S i .10 3 [m 3 /kg]	T i [h.]	Product B2	S i .10 3 [m 3 /kg]	T i [h.]
	R3	5.91	16.5	R3	5.68	33
	R6	5.91	10	R21-2	6.70	8.5
	R5	5.45	11.5	R21-1	5.01	8.5
	ES	4.26	1.5	R6	5.01	20
	R4	2.48	6	F	2.66	--
				S	1.25	23

Table 2

 2

		: Production data.
	Product	Production (kg/year)
	A1	300000
	B2	200000

Table 3 :

 3 GA parameters.

	Type of parameter
	Population size	400
	Generation number	1000
	Survival rate	0.5
	Mutation rate	0.4
	Elitism	1

Table 4 :

 4 Multicriteria batch plant retrofit results.

	MonoP Criteria	R3	R6	R5	ES	R4	R21	F	S
	S1	768240/4	2-16	1-16	0	1-2	0	0	1-8/4	1-3
	S2	804157/3	2-16	1-16	0	1-2	0	0	1-16/4	1-3
	S3	856303/2	2-16	1-16	0	1-2	0	0	1-16/8	1-3
	S4	910743/1	2-16	1-32	0	1-2	0	0	1-16/8	1-3
	S5	970561/0	2-16	1-32	0	1-2	0	0	1-16/8	1-6
	S6	1014890/2	2-16	1-16	1-16	2-2	1-8	0	0	0
	S7	1050810/1	2-16	1-16	1-16	2-2	1-16	0	0	0
	S8	1105250/0	2-16	1-32	1-16	2-2	1-16	0	0	0
	MultiP Criteria	R3	R6	R5	ES	R4	R21	F	S
	S1 346690 €/1 2-16	1-16	0	0	0	0	0	0
	S2	401130/0	2-16	1-32	0	0	0	0	0	0


Sofiane Hocine

Luc Pibouleau

Catherine Azzaro-Pantel

Serge Domenech

MODELLING A FLUIDIZED BED FOR WATER TREATMENT VIA GAMS PACKAGE

) by using a combined stochastic/NLP approach. The problem is solved for two values of fluidization velocity (50L/h and 90L/h), which represent two experimental behaviours, and with two GAMS tools (DICOPT and SBB). For the two flowrates considered separately, the solutions of Montastruc and SBB are quite equivalent. However, when both flowrates are considered together, SBB easily reaches an optimal solution with a minimum number of elementary units. Due to crisp constraints, the Montastruc's method and DICOPT need "good" initializations for giving only worse solutions than SBB. From this study, it appears that the SBB procedure is more efficient than the DICOPT one and the hybrid method of Montastruc for modelling a fluidized bed through a superstructure.

INTRODUCTION

Phosphorus recovery from wastewater accords with the increasing demand for sustainable development of the phosphate industry and the associated stringent environmental quality standards. During the last decade, numerous engineering solutions addressing phosphorus recovery from wastewater by precipitation of calcium phosphate in a recyclable form have been proposed (Morse et al., 1998). An advanced alternative is to implement the so-called pellet reactor approach first presented by Seckler et al., (1996), and recently used by [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF]. In this work, a two-step procedure is proposed with Ca-PO 4 -H 2 O as a support system. The first modelling level was previously presented by [START_REF] Montastruc | Modélisation et optimisation d'un réacteur en lit fluidise de déphosphatation d'effluents aqueux[END_REF], and leads to the development of a thermodynamic model for predicting phosphate conversion. The second step is related to the computation of the pellet reactor's efficiency. It must be emphasized that a phenomenological way to treat the problem would imply the precise knowledge of some key parameters of the induced agglomeration process (such as coating), which are in fact rather difficult to obtain practically, such as calcium phosphate density and thickness. For these reasons, in the work of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF], another alternative based on a network-oriented reactor model is proposed: the pellet reactor model is identified as a combination of elementary systems, representing ideal flow patterns, as perfect mixed flows, plug flows, continuous stirred tank reactors, etc. A superstructure involving the set of all possible solutions corresponding to the physical reactor is defined. Then, potential solutions are extracted from this superstructure and evaluated according to a given performance function. In the work of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF], the MINLP problem (Mixed NonLinear Programming Problem) arising from the superstructure formulation is solved by a combination of a simulated annealing (SA) and a quadratic programming (QP) method from the IMSL library. Each potential structure proposed by the simulated annealing procedure is optimised with the QP. The major difficulty of the approach of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF], lies in the initialisation of the QP procedure at each iteration of the SA, probably leading in some cases to quite poor solutions. The present work is a continuation of the one of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF], where the MINLP problem of the model identification is formulated using the high level input language of the GAMS software. According to two values of fluidization velocity (50L/h and 90L/h), which represent two observed experimental behaviours, and several objective functions, the problem is solved by means of two GAMS tools (DICOPT and SBB).

The DICOPT and SBB software packages can be used for solving mixed-integer nonlinear programming (MINLP) problems that involve linear binary or integer variables as well as linear and nonlinear continuous variables. Both solvers are based on different solution algorithms, leading to two different ways of solving the problem, in terms of efficiency. This article is organized as follows. First, the fluidized bed is described for the two experimental fluidization velocities. Then, a brief survey of MINLP methods is recalled, and followed by the presentation of the numerical tools chosen within the GAMS package. Then, the fluidised bed is modelled by means of a superstructure and the problem constraints are listed. The problem data and the various objective functions are described in the following part. The final section concerns the results presentation and their comparison with those obtained by [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF].

PROCESS DESCRIPTION

The process is based on calcium phosphate precipitation obtained by mixing a phosphate solution with calcium ions and a base. More precisely, it involves a fluidized bed of sand continuously fed with aqueous solutions (see Figure 1). Calcium phosphate precipitates upon the surface of sand grains; at the same time, small particles, also called "fines," leave the bed with the remaining phosphate not recovered in the reactor. Two different regimes were experimentally observed (see Figure 2). On the one hand, for high values of the fluidization velocity (>90L/h), only one zone is observed at the bed top, in which fines leave the bed with the liquid effluent. Process efficiency is thus mainly attributed to fine coating on sand grains. On the other hand, for low values of the fluidization velocity (<50 L/h), an additional layer is observed at the upper zone of the fluidized bed where fines stagnate and agglomerate. Because of the stratified form and size of the agglomerates, the fines remain at the surface of the bed. Inside the fluidized bed, fines totally cover the sand grains and constitute large particles of complex structure. In this paper, only the second step of modelling approach is considered, that is to say that the pellet reactor efficiency is computed from the identification of the pellet reactor model as a reactor network, represented by a superstructure involving a combination of ideal flow patterns.

NUMERICAL SOLUTION PROCEDURES

Survey of MINLP methods

Numerical procedures for solving MINLP problems have appeared only quite recently, compared to other mathematical programming approaches, like classical linear and nonlinear ones. Indeed, MINLP procedures are based on efficient NLP and MIP or IP (Mixed Integer Programming or Integer Programming) packages. Two approaches can be distinguished for solving MINLP problems: the stochastic method and the deterministic one, as presented in figure 3. A third solution technique, combining the two previous methods, is the hybrid approach. The stochastic approach mimics natural phenomena, like annealing of solids (Pibouleau et al,2005), ant colony behaviour, or natural evolution of species (Dedieu et al,2003). It gives birth to evolutionary algorithms, which can reach a "good" solution, not necessary an optimal one in terms of Kuhn-Tucker conditions, starting from an initial and feasible solution according to the constraints, through a set of local improvements. Insofar as they do not require the knowledge of mathematical properties of the problem, they are generally easy to implement. In deterministic methods, the optimal solution (at least a local one) is guaranteed, but they are more complex to implement, due to strict mathematical properties required for the problem under consideration. A recent alternative [START_REF] Hocine | A hybrid optimization method (GA/SQP) for methyl acetate production process by reactive batch distillation[END_REF][START_REF] Montastruc | Modélisation et optimisation d'un réacteur en lit fluidise de déphosphatation d'effluents aqueux[END_REF] consists in using hybrid methods, combining a stochastic procedure for solving the master (MIP or IP) problem, and a deterministic one by implementing a NLP package for solving the slave (continuous) problem.

Because on the one hand, in the stochastic approach an optimal solution is not guaranteed, and on the other hand in the hybrid method, the initialisation of the NLP procedure at each iteration on the master problem is a rather difficult task, a deterministic approach using the GAMS environment has been adopted in this study.

GAMS optimisation tools

GAMS is a modelling system for optimization that provides an interface with a variety of different algorithms. Models are supplied by the user to GAMS in an input file in the form of algebraic equations using a high level language. GAMS then compiles the model and interfaces automatically with a "solver" (i.e. optimisation algorithm, SBB and DICOPT in this study). The compiled model as well as the solution found by the solver are then reported back to the user through an output file [START_REF] Gams | A user's guide[END_REF]. Two procedures, presenting different characteristics in terms of numerical methods used and variable initialisation phase, are used in this study.

DICOPT (Discrete Continuous Optimization Package)

The solver is based on the OA/ER optimization algorithm (Kocis and Grossmann, 1987). The MINLP algorithm solves a series of NLP and MIP sub-problems by using the solver CONOPT for the NLP part, and the solver CPLEX for the MIP sub-problem.

SBB (Simple Branch & Bound)

This more recent solution procedure is based on Branch and Bound method, where some integer variables are considered as continuous ones at some nodes of the search space tree, providing lower bounds, by solving NLP problems (with the CONOPT solver), for guiding the search [START_REF] Allen | How to use mixed integer programming[END_REF][START_REF] Gupta | Branch and Bound experiments in nonlinear integer programming[END_REF]. The main difference between the two solvers lies in the variable initialisation phase, which is not required for the SBB solver. Let us recall that the initialisation phase was a crucial feature [START_REF] Montastruc | Modélisation et optimisation d'un réacteur en lit fluidise de déphosphatation d'effluents aqueux[END_REF].

FLUIDIZED BED MODELING

In order to get a sufficient degree of accuracy, a variant of the superstructure of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF] was proposed, where each plug flow reactor was replaced by a series of four well-stirred tank reactors and additional flows were added, as shown in Figure 4. The superstructure, includes the one of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF], and involves a dead zone (volume V1), a by-pass (number 15), a recycle stream (number 2), and 12 continuous stirred tank reactors (CSTR), (volumes V2 to V13). The binary variables of the MINLP problem represent the existence or absence of elementary units, and the continuous ones the corresponding volumes. Since it involves 15 binary variables, the problem combinatorics is 2 15 =32,768 potential solutions.

For the two flowrates, taken separately or simultaneously, the problem consists in extracting from the superstructure, the model allowing to reach the specified output concentration, with a minimum number of streams and elementary units. The problem formulation comprises the following equations (where Vi, Fi, and Ci represent respectively the reactor volumes, the total flowrates, and the molar concentrations in fines), whereas the term yi is related to the presence or absence of element i.

Volume Constraint

total i i V V = ∑ = 13 1 (1) Input Node Constraints 14 11 1 13 12 10 0 F F F F F F F + + = + + + (2a) 1 14 11 1 13 12 17 10 9 0 0 ) ( C F F F F C F C F C F C f + + = + + + (2b) Output Node Constraints 0 3 12 6 14 11 5 F F F F F F F + + + = + + (3a) f C F F F F F C F C F C
) ( 0 13 12 6 14 1 11 13 5 5

+ + + = + + (3b)
Balance on reactors -these equations have the following form (example for reactor 7):

2 6 10 9 2 2 9 8 1 1 C V k F C F C F C F C + + = + (4a) 10 7 2 9 7 1 F y F F y F + = + (4b) 10 7 9 7 ) 1 ( ) 1 ( F y F y - = - (4c) Concentration Equations 8 7 2 7 9 ) 1 ( C y C y C - + = (5a) 7 8 3 8 8 ) 1 ( C y C y C - + = (5b) 6 9 4 9 7 ) 1 ( C y C y C - + = (5c) f C y C y C ) 1 ( 10 5 10 6 - + = (5d)
The optimisation problem involves the following items. Linear Constraints: 1 equation for the volumes, 8 total balance equations on the reactors, 2 total balance equations on nodes, 4 concentration equations Bilinear Constraints: 2 partial balance equations on nodes, 12 partial balance equations on the reactors Bounded Variables: 13 volumes, 14 flowrates, 13 concentrations

PROBLEM DATA AND SOLUTION

The following sets of data were used for solving the model identification problem.

Case 1 (low flowrate)

Case 2 (high flowrate) F 0 = 50 l/h C 0 = 1 mol/l (inlet concentration) C out = 0.258 mol/l (inlet concentration) K = 91 h -1 V tot = 1.3 l (total volume) F 0 = 90 l/h C 0 = 1 mol/l (inlet concentration) C out = 0.477 mol/l (outlet concentration) K = 68.5 h -1 V to t = 1.9 l (total volume) For high (90 l/h) and low (50 l/h) input flowrates considered separately, the problem is solved according to three objective functions to be minimized: ) represents the output concentration given by the model (respectively the given experimental output).

3)

∑ = + - 15 1 2 exp , mod , ) ( i i out out y C C ρ
, in order to obtain the best model with regard to the concentration, but also the simplest one with respect to the number of elementary units (ρ is a positive penalty coefficient adjusted so that both terms of the objective function have the same order of magnitude). Let us note that the search for the model involving the smallest possible number of elementary units, constitutes a key point for future control purposes of the process. When both flowrates (90 and 50 l/h) are considered together, the goal is to extract from the superstructure a model that fits well with the two experimental outputs. In that case, the first objective function cannot be used because the problem would involve two conflicting equality constraints with the same left hand side and two different right hand sides (the two output experimental concentrations); so the problem is solved according to the two objective functions to be minimized:

2 2 exp , mod , 1 2 exp , mod , ) ( ) ( case out out case out out C C C C - + -
, where case 1 and case 2 are related to the two studied flowrates.

∑ = + - + - 15 1 2 2 exp , mod , 1 2 exp , mod , ) ( ) ( i i case out out case out out y C C C C ρ
, in order, as in the previous case, to obtain the best, but also the simplest model.

NUMERICAL RESULTS

The results obtained by the two methods SBB and Dicopt, for the three flowrate value cases and for the three objective functions are presented. The results obtained previously by [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF] The study presented in this paper compared the performances of two different approaches for tackling the problem of model identification, involving an MINLP formulation : a purely deterministic framework via the GAMS package and a hybrid SA/QP approaches within the GAMS environment, two solvers were used, i.e., SBB and DICOPT for MINLP problem. The solutions obtained from DICOPT strongly depend on initial values, and are only local optima. This is not the case with SBB where initial values for optimization values are not necessary. For the two flowrates considered separately, the solutions of SA/QP and SBB are quite similar. However, when two flowrates are considered together, SBB easily reaches an optimal solution with a minimum number of elementary units, while due to crisp constraints, SA/QP and DICOPT need "good" initializations for giving only worse solutions than SBB. As above mentioned, the high level input language of GAMS and the efficiency of its incorporated solvers (particularly SBB), make easy the solution of large MINLP problems, arising from model identification of complex processes.

For the modeling of the considered fluidized bed, the approach based on the GAMS environment (particularly SBB) has proved to be more efficient than the hybrid approach of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF] using a combination of the Simulated Annealing procedure with a nonlinear programming package. Current works are related on the model identification of more complex processes, such as the residence time distribution in a ventilated room, by introducing some dynamic (time dependent) variables.

FigureFigure 2 .

 2 Figure 1. Schematic representation of the pellet reactor

Figure 3 .

 3 Figure 3. The two approaches for solving MINLP problems.

Figure 4 .

 4 Figure 4 . Superstructure for the fluidized bed reactor

  that case the given output concentration corresponding to each flowrate is introduced as an

  are also recalled for comparison purposes.

		( C	out	,	mod	-	out C	,	exp	)	2													0.6
																																			NONOPTIMAL STRUCTURE
																																			0.82	0.18
																																			1.6
																																			( C	out	,	mod	-	C out	,	exp	)	2	=	6	84 .	e	-	14
	( C	out	,	mod	-	C	out	,	exp	)	2	+	ρ	∑ 15	i y										0.10	0.90	1.01	0.10	0.90	1.06
																					i	= 1											( C	out	,	mod	-out C	,	exp	)	2	=	0	0.01	( C	out	,	mod	-out C	,	exp	)	2	=	0 0.061
	Literature result																					Montastruc (SA/QP)
	( C	out	,	mod	-	C	out	,	exp	)	2	+	ρ	∑ = 15 1 i	i y									1.33	0.50	2.74	0.50	1.33	( C out	,	mod	-	C out	,	exp	)	2	=	1	e	-	20
																																	0.33	1.74	0.33
	Case 2 (high flowrate)
	Objective function																		SBB solver	Dicopt solver
								15																								0.82
								∑	i y
								i	= 1																								0.24	0.76
																																			0.18
																																( C out	,	mod	-	C out	,	exp	)	2	=	6	6 .	e	-	21	( C out	,	mod	-	C out	,	exp	)	2	=	4	1 .	e	-	22
																																		0.13	0.45
		( C	out	,	mod	-	out C	,	exp	)	2
																																			0.55	0.24	0.76
																																	0.87
																																( C	out	,	mod	-	C out	,	exp	)	2	=	5	36 .	e	-	23	( C out	,	mod	-	C out	,	exp	)	2	=	4	1 .	e	-	22
																																			0.82	0.21
	( C	out	,	mod	-	C	out	,	exp	)	2	+	ρ	∑ 15	i y										1.79
																					i	= 1											0.93	0.07
																																			0.18	0.79
																																( C	out	,	mod	-	C	out	,	exp	)	2	=	. 2	5	e	-	21	( C out	,	mod	-	C	out	,	exp	)	2	=	1	e	-	19
	Literature result																					Montastruc (SA/QP)
	Case 1 (low flowrate) + -exp , mod , ) ( out out C C 2	ρ	= 1 ∑ i 15	i y										0.6 0.4	0.36	( C out	,	mod	-	C	out	,	exp	)	2	=	1	e	-	20
	Objective function																		SBB solver 1.36	Dicopt solver
								15																								1.01	1.06
	∑ = 1 i Case 3 (two flowrates considered together) i y 0.10	0.90	0.10	0.90
																																			0.01	0.061
	Objective function =	( C	out	,	mod	-	C	out	,	exp	)	2	( C 1 out case +	mod ( , out -out C mod , C	exp C -,	) out 2	0 exp = ,	)	2	case	2	( C Quadratic mod , -out out C deviation	,	exp	)	2	=	0
																																			( C	out	,	mod	-	out C	,	exp	)	2

In case 1 (low flowrate), the results obtained by the solver Dicopt for objective functions 1 and 3 are close to those obtained by the solver SBB. On the other hand, the structure obtained with the second objective function by using the solver Dicopt is nonoptimal (i.e. this solution comprises a number of items greater or equal than 5). This difference is explained by the initialization of variables, which is a crucial point for this kind of problem implying strong interactions between the variables. In case 2 (strong flowrate), in spite of the difference between the obtained structures, the results are equivalent. In both cases, the number of elementary structures does not exceed two. One can think that the problem comprises several local optima because several models including the same minimal number of elementary operations (here two) can give the same exit concentration. This explains the fact that the three procedures lead to relatively different results, while remaining rather close. The same remark as in case 1 can be made for the use of Dicopt where the obtained results may vary according to the initialization of the variables (in particular for the discrete variables). Finally, in case 3 (two flowrates together), one notices that SBB finds the best solution for the two objective functions, whereas Dicopt does not find an optimal solution. The need for a good initialization for Dicopt is clearly confirmed in this last case. It can be noted that the structures obtained are rather comparable except in certain cases for the solutions provided by Dicopt, where the need for a good initialization appears to be a basic point for avoiding convergence towards a non optimal solution. This kind of situation is frequently encountered in the treatment of problems of large size with the Dicopt solver. Even if the results obtained using the GAMS environment are close to the results presented in (Monstatruc et al, 2004), some issues such as the computation time and initialization (generally not required for the SBB solver) show that the use of the deterministic methods implemented in this study appear more efficient than the hybrid method (SA/QP) of [START_REF] Montastruc | A Systemic Approach for Pellet Reactor Modeling: Application to Water Treatement[END_REF].


F De Belleville

L Dairaine

M Gineste

C Fraboul

Reducing Satellite Communication Cost Using Terrestrial Peer-to-Peer for Lost Recovery

A practical solution to implement IP multicast service may consist in using a geostationary satellite. The broadcast nature and the large coverage zone of such systems make it possible for a source to reach a potentially very large number of receivers with only one hop. In the context of reliable multicast communications, a hybrid satellite/terrestrial approach based on communication costs is described. Due to particular data dissemination resulting from the satellite communication phase, ad-hoc lost recovery peer-to-peer mechanisms are specially designed and evaluated through simulations.

I. Introduction

The support for point to multipoint communications is a useful and performing service that is not sufficiently available in nowadays Internet networks. Numerous applications like multimedia streaming or software deployment and update would indeed take a great advantage of such a service. IP Multicast [START_REF] Deering | Host Extensions for IP Multicasting[END_REF] offers an appropriate multipoint service in the context of the ubiquitous Internet Protocol (IP). Nevertheless, its availability for large communication groups over the whole internet is still very limited for the time being. This lack of deployment is mainly due to technical concerns and economical issues [START_REF] Diot | Deployment issues for the IP multicast service and architecture[END_REF].

A practical solution to implement IP multicast service may consist in using a geostationary satellite. The broadcast nature and the large coverage zone of such systems make possible for a source to reach a huge number of receivers with only one hop. Satellite broadcasting may seems expensive at first sight, but per-receiver cost becomes less than using terrestrial network when the number of receivers increases. For this reason this study is focused on large scale reliable multipoint communications. Moreover we consider applications with no time constraints, because the long transmission delay of satellite links is not really compatible with such applications. Software updates or cache feeding in content delivery networks are examples of possible target applications.

The main objective of the proposed approach, named Hybrid Satellite Terrestrial Reliable Multicast (HSTRM), is to provide a low cost fully reliable multipoint service. Thus the satellite link may be used only to transmit packets useful for a large group of receivers. A way to partially achieve this goal is to use Hybrid ARQ type II [START_REF] Nonnenmacher | Parity-based loss recovery for reliable multicast transmission[END_REF] for retransmissions: with this technique, any packet transmitted can be exploited by any receiver waiting for information. The remaining issue is to ensure that the satellite link is not used for a small group containing few receivers only. During a satellite transmission the receivers which experienced losses have indeed to wait for more information to recover those losses. In consequence more time is required to transmit data to them. Thus the number of receivers which continue to really use the satellite broadcast will decrease all along the transmission. To solve this problem, our approach consists in estimating regularly the broadcast audience during the transmission [START_REF] De Belleville | Group size estimation for hybrid satellite/terrestrial reliable multicast[END_REF]: data are then broadcasted via satellite only when the number of receivers is sufficiently high (i.e. when per-receiver cost is less than with terrestrial networks) and via terrestrial transmissions otherwise.

The present paper focuses on issues related to the terrestrial phase of the communication. In our context, peer-topeer appeared to be the most suitable technique for seeking and exchanging data over the terrestrial network. Nevertheless, because of the atypical context (hybrid communication, satellite broadcasting, use of Forward Error Correction) and the specific objective of the approach (minimize the overall communication cost), a set of ad-hoc mechanisms has to be proposed. This led us to the proposition of a specific mechanism for terrestrial error recovery which takes advantage of both the satellite and the terrestrial networks. Some of its parameters are set according to theoretical computation results. The remaining ones are set thanks to simulation results. Performances of the proposed mechanism are eventually studied through simulations.

The paper is structured as follow. In the second section, the context of the study is described. The specific communication environment and the multicast transport service are then presented. Then, the principles of the Hybrid Satellite Terrestrial Reliable Multicast (HSTRM) protocol are given. A statistical study comparing costs for terrestrial, satellite and hybrid terrestrial/satellite multipoint is finally proposed. In the third section, the paper focuses on the terrestrial phase of the HSTRM protocol with the complete definition of an efficient low-cost peer-topeer mechanism. In the last section, parameter settings and performance study of the mechanism are provided.

II. A Reliable Transport Protocol Designed for Hybrid Satellite-Terrestrial Network

In this section, the context of the study is specified. Then the target transport service is described, as well as a proposition of a hybrid satellite/terrestrial approach to achieve a reliable multicast transport service. The last part presents a statistical study of this approach.

A. Communication System

The considered system is a hybrid satellite/terrestrial network, i.e. end-users are connected to both satellite and terrestrial networks. Assumptions have been made in the context of the French DIPCAST project [5]: end-users are connected to terrestrial network via a high speed access network (e.g. xDSL or LAN), and to satellite system either via a high speed access network, or directly with a Very Small Aperture Terminal (VSAT).

The satellite system uses a geostationary satellite, and proposes a best effort multipoint communication service based on IP Multicast. This supposes that a protocol which manages joining and leaving procedures is integrated to the satellite system, as well as tree establishment algorithm. This problem of integration has addressed out during the DIPCAST project but is out of the scope of this paper.

In order to be representative of today Internet, the terrestrial network is not supposed to support multipoint transmissions. Thus any terrestrial communication in this hybrid system is a point-to-point transmission. Finally the following assumptions are made on the application:

-The application transmits data from one source to a large group of receivers (one-to-many communications) -This application does not have strong delay constraints. -A service transmitting session characteristics (file properties, start time, associated group address, etc.) is available. Those transmissions can be done via out of band means (e.g. e-mail), or via session management tools.

-No receiver can join the session after the beginning (late joining is not supported).

B. Multicast Transport Service

Since satellites are really advantageous for large scale data transmissions, applications which transfer files to a large number of receivers (several hundreds or more) are particularly considered in the present paper. Furthermore we consider applications requiring a full reliability, i.e. which must be assured that the whole group has received transmitted information. An example of such applications is the transport of multimedia files (e.g., video, music, games, etc.) towards a large set of users. Another considered issue concerns the overall communication cost. The utilisation of satellite links is indeed quite expensive. Nevertheless when the number of receivers increases, the peruser cost decreases. Thus any application being charged according to its bandwidth utilisation may prefer protocols which carefully watch communication cost.

Considering file transport using the best effort protocol IP Multicast over satellite, a fully reliable transport protocol must be used. In the following paragraphs some features of this protocol required in the above-mentioned context are exposed.

The transport protocol must guarantee that all receivers in the group receive transmitted information. Several multicast transport protocols propose a statistically reliable service [START_REF] Koyabe | Reliable Multiast by Satellite: A Comparison Survey and Taxonomy[END_REF]. Although it allows designing transport protocols with no return channel, full reliability is not ensured because there is no adaptation to losses that really happened. This technique is then not convenient for the aimed purpose.

Since satellite bandwidth is expensive, the transport protocol should ensure that any useless satellite transmission is avoided. This confirms that statistical reliability is not recommended in the presented context because systematic coding of information implies a potential waste of bandwidth [START_REF] Jung | Uni-directional versus Bi-directional Communication[END_REF]. Finally, as considered applications are designed for transmissions towards very large groups, underlying protocols (and specifically transport protocol) must scale very well. In particular for the transport layer, mechanisms of feedback suppression like [START_REF] Nonnenmacher | Optimal Multicast Feedback[END_REF] must be studied and configured for a satellite link.

C. A Hybrid Satellite/Terrestrial Approach

Numerous multicast transport protocols have been designed in the last few years to achieve efficient and scalable multicast transmissions [START_REF] Koyabe | Reliable Multiast by Satellite: A Comparison Survey and Taxonomy[END_REF]. From all researches on reliable multicast transport, a technique referenced as Hybrid ARQ type II [START_REF] Nonnenmacher | Parity-based loss recovery for reliable multicast transmission[END_REF] has emerged. It is an efficient way to reduce the used bandwidth and to improve the scalability. It consists in using FEC combined with Automatic Repeat reQuest (ARQ) in the following way: after a transmission, the source asks for the maximum number M of missing packets. Then it generates and transmits M new encoded packets. As those M packets have not already been transmitted, they are useful for any receiver which experienced losses. Hybrid ARQ type II allows to radically reduce the amount of retransmitted information in multipoint transmissions, and is then particularly interesting for any large scale full reliable transport protocol [START_REF] Nonnenmacher | Parity-based loss recovery for reliable multicast transmission[END_REF].

Most of the protocols presented in [START_REF] Koyabe | Reliable Multiast by Satellite: A Comparison Survey and Taxonomy[END_REF] are usable with satellite links because they support asymmetric transmissions. Nevertheless no previous work considers the overall communication cost. With the used reliability mechanisms, either packets are systematically retransmitted to the whole group, or each missing packet is retrieved with a point-to-point connection. In our context when point-to-point retransmissions are used, if a packet is requested by numerous receivers, satellite broadcasting may be profitable. On the contrary, if every missing packet is transmitted to the group via satellite, when only few receivers are concerned, few terrestrial point-to-point connections can be cheaper. According to this simple statement, the trade-off between satellite and terrestrial bandwidth use may be studied, and it is our belief that optimisation of this trade-off is an interesting way to reduce the overall communication cost.

An interesting approach may be to define a threshold R min representing the minimum number of receivers for satellite broadcasting to be advantageous (taking economical costs into account) [START_REF] De Belleville | Une Approche Hybride Satellite/Terrestre pour le transport fiable multipoint à grande échelle[END_REF]. A session would then behave as follows: at a predefined time, the satellite starts broadcasting over the IP multicast group. During this transmission, the source periodically estimates the session size. A receiver is considered to belong to the session as long as it has not received the whole transmitted information. Several papers have addressed the question of estimating multicast session size [START_REF] Alouf | Optimal on-line estimation of the size of a dynamic multicast group[END_REF]. An appropriate mechanism has been proposed in the specific context of HSTRM [START_REF] De Belleville | Group size estimation for hybrid satellite/terrestrial reliable multicast[END_REF] and is assume to be available in this paper. Thus once the entire initial information has been transmitted, session size is likely to decrease (all receivers which experienced no losses quit the session). The source then goes on estimating session size while it transmits encoded redundancy packets to repair losses. When estimated session size goes below R min , the satellite transmission stops. Receivers which do not have enough information to decode received data (i.e. receivers experiencing high loss rates) then contact other receivers to recover missing data. This terrestrial recovery can be done using peer-to-peer services. When all receivers have fully received the information, the session stops.

D. Statistical Study

The present paragraph gives an illustration of the gain generated by the hybrid satellite/terrestrial approach. For that purpose a statistical study is presented, with the following assumptions:

Errors and Fades

For terrestrial transmissions, packet losses are only supposed to be produced by network congestions. These congestions cause at the transport layer a Packet Loss Rate (PLR) of 5% for point-to-point communications [START_REF] Paxson | End-to-end Internet packet dynamics[END_REF] and 10% for multipoint communications [START_REF] Yajnik | Packet loss correlation in MBone multicast network[END_REF].

For satellite transmissions, three categories of receivers are considered. The first one corresponds to all receivers under a clear sky. Assumption is made that these ones do not experience any losses (possible issues due to scintillation are not considered) and that this category encompass 90% of the receivers. The second category is supposed to include 9.9% of the receivers, and corresponds to end-users under a rainy sky. Fades due to light rain are supposed to cause a PLR of 20% at transport level. Finally, the last category encompasses receivers under a stormy weather. This category is supposed to include 0.1% of the receivers which experience a PLR of 60%. According to [START_REF] Pinder | Four years of Experimental Results from the New Mexico ACTS Propagation Terminal at 20.185 and 27.505 GHz[END_REF] those values are realistic for a satellite communication using the Ka Band.

Note that all losses are supposed to be independent and uniformly distributed among packets. This assumption is not realistic, but since transport layer is supposed to implement ARQ type II technique, only the amount of lost packets is important for our computations. Assumption was also made that no packet are lost on the return channel.

Cost Function

In order to compare costs generated by hybrid satellite/terrestrial multipoint communications with pure terrestrial and satellite multipoint communications, it is necessary to first define a cost function. We choose to adopt a perpacket cost approach, and then define a cost function as follows:

( ) ( ) K R C K R F X X X , , × = α , { } SM TM TU X , , ∈ , ( 1
)
where K is the number of packets to transmit, R the number of receivers and α x the per-packet transmission cost. C X (R,K) represents the average number of packets passing through the network in order to transmit K packets to R receivers (taking losses into account). The indices TU, TM and SM correspond to Terrestrial Unicast, Terrestrial Multicast and Satellite Multicast communications.

For terrestrial point-to-point communications R connections experiencing independent and uniformly distributed losses (represented by the PLR) are used to transmit K packets to R receivers. Then:

( ) ( ) PLR R K PLR PLR i R K K R C i i TU - × = × - × × × = ∑ ∞ = - 1 1 , 1 1 .
(

) 2
For multipoint communications, the probability P(N,PLR,R) that exactly N packets are needed so that all endusers receive K packets can be expressed as: the probability that each end-user receive K packets among N, minus the probability that R receivers receive K packets among (N-1). Then P(N,PLR,R) can be calculated as:

( ) = R PLR N P , , . (3)
For terrestrial multicast communications, according to [START_REF] Chuang | Sirbu: Pricing Multicast Communication: A Cost-Based Approach[END_REF] transmitting data to R receivers is equivalent to 8 . 0 R point-to-point connections. The average number of packets passing through the network is then:

( ) ( ) 8 . 0 , , . , R R PLR N P N K R C K N TM × ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ∑ ∞ = . ( 4
)
For satellite multipoint communications, since we defined three categories of receivers and since no multicast tree is established, we have:

( ) ( ) [ ] ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∈ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ∑ ∞ = 2 , 0 , , , . max , i R PLR N P N K N C K N i i SM β (5)
where PLR 0 =0, PLR 1 =20%, PLR 2 =60%, β 0 =90%, β 1 =9.9% and β 2 =0.1%. For hybrid satellite/terrestrial transmissions, the cost is defined as the sum of the costs generated by terrestrial network utilization, and satellite system utilization.

Results

Using the cost function defined above, R min has been computed for group sizes ranging from 100 to 600,000. Then the hybrid satellite/terrestrial approach has been compared to terrestrial point-to-point and multipoint communications as well as to pure multipoint satellite communications. Figure 1 shows the results for α U = α M = 1, α S = 100 and K = 100. The different levels perceptible on the curve representing the hybrid communication cost are due to the model definition: when group size increases, the number of receivers under a rainy sky or a stormy sky increases as well. Thus it becomes necessary to repair more and more losses using satellite link. In this example, the hybrid satellite/terrestrial approach induces a gain ranging from 10% to 50% compared with the most advantageous of the three classical approaches.

(

) The hybrid satellite/terrestrial approach relies on hybrid ARQ type II. As mentioned previously this technique implies to encode transmitted information. The remaining of the paper is focused on the design of an efficient low cost mechanism for the terrestrial phase of the hybrid approach.

[ ] K N when PLR R K = - , 1 ( ) ( ) K N when PLR PLR i N PLR PLR i N R i N i K N i R i N i K N i > ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - - = - - = ∑ ∑ , 1 1 1 ) 1 (

III. An Efficient Low-Cost Peer-to-Peer Mechanism for Terrestrial Communications

The terrestrial phase of transmission starts as soon as the number of receivers became lower than a cost-related threshold R min formally defined into [START_REF] De Belleville | Une Approche Hybride Satellite/Terrestre pour le transport fiable multipoint à grande échelle[END_REF]. When the number of receivers is lower than this threshold, which will be noted MinRcv, satellite transmission is supposed to be not advantageous any more. During this phase, receivers seek missing data in order to supplement those which they received during the satellite phase of transmission. One major need for this approach is a mechanism for retrieving packet losses using the terrestrial network. In our context, peerto-peer appeared to be the most suitable technique for exchanging data over the terrestrial network. Nevertheless, the analysis of existing peer-to-peer mechanisms reveals that none of these proposals take the cost criteria into account. Moreover, most the classical mechanisms are designed either for conditions where the localisation of data stored into only few peers or for the dissemination of a large amount of data.

A. Objectives and hypothesis

When the satellite transmission phase stops, the data repartition over the network is slightly different from the classical hypothesis taken in the context of peer to peer networking. HSTRM proposes a very large scale group communication service. As shown previously, if the satellite transmission stops, this means the number of receivers having not received the information has became relatively small in regards to the total number of receivers. Then, only a small fraction of the total number of receivers does not have all the information. Nevertheless, the protocol does not know how many data the lossy receivers got during the satellite transmission, with possible receivers having received no information at all. Moreover, considering the proximity in the network is not related to the geographical proximity, the system does not have information about the repartition of the information into the network. Recall that all receivers are connected to the terrestrial network by a high speed access network such as xDSL line.

HSTRM is designed to be used in the context of large group; then the terrestrial mechanism must be also suitable for large scale network. This means that the P2P mechanisms should not have an important impact on the underlying network or the receivers. This is particularly true for the P2P network establishment because all group members are involved into this phase. To allow an easy deployment of the protocol, the terrestrial mechanism should not impose specific requirements on the underlying network, such as IP multicast over wired network.

Finally, a very important constraint is related to the cost of use of the terrestrial communication. For the global proposal to be of interest, the proposed terrestrial recovery mechanism should not be more expensive than a simple satellite broadcast or terrestrial multicast.

B. Related work

In such a way to ensure HSTRM aims, different types of approaches have been analysed, and particularly the techniques of cooperative downloads, application-level multicast and peer-to-peer downloads.

The first type of mechanism is known under the appellation of cooperative downloads. In this case, the main objective is to relax the load imposed on a centralized server. This situation classically occurs when a server has a very popular content such as a software update. This problem has been addressed e.g. by Slurpie [START_REF] Sherwood | Slurpie: A cooperative bulk data transfer protocol[END_REF] and Bit Torrent [START_REF] Mikel | Dissecting BitTorrent : five months in a torrent's lifetime[END_REF]. The aims of these solutions are to replicate as quickest as possible a large file into a set of distributed hosts. To achieve this, new receivers that want to get the file are redirect towards other users which already have partially or totally the file. The load is then distributed on the various participants to this download. Nevertheless, the main objective of these contributions differs from HSTRM. During the terrestrial phase, the point is not to duplicate content on the whole group where nobody has the information. At this phase of communication, only few receivers are directly concerned and usually look after a small amount of data (i.e., the blocks of data that have not been received during the broadcasting phase).

Another classical solution to allow a group of receivers to communicate with other receivers consists in using an application-level multicast service that consists in managing the data replication to various receivers in end-system high layer. This approach has been addressed in many studies [START_REF] Costello | Search party: Using randomcast for reliable multicast with local recovery[END_REF][START_REF] Nonnenmacher | How bad is reliable multicast without local recovery[END_REF]. Nevertheless, as the previous type of problem, the objectives differs from ours in transmitting a file to a whole set of receivers and the context is a transmission from a source to the group.

Another way to disseminate data among a set of receivers make use of peer to peer techniques. CAN (Content Addressable Network [START_REF] Ratnasamy | A scalable content-addressable network[END_REF]) and CHORD [START_REF] Stoica | Chord: A scalable peer-to-peer lookup service for internet applications[END_REF] are two classical techniques proposed to localise efficiently information in a P2P network. They make use of hash functions: when a node looks for specific information, it calls the hash function that computes an identifier. This identifier can be considered as the information's address into the P2P network. This identifier is then searched into the network. These approaches allow localizing information among a set of receivers which are offering large information diversity. Recall to our context, the receiver looks after data that are already present in a large set of receivers. The localization then becomes far less complex than in a classical P2P context. Almost all receivers have a neighbour that already owns the requested information.

Local recovery mechanisms allow requesting retransmission to hosts located "near" the receiver. The problematic addressed into these works are then very related to our context. A proposal consists in dividing the receivers into several sub-groups, each sub-group having a host in charge [START_REF] Cao | Reliable multicast via satellites[END_REF]. When losses occur in a receiver, it seeks missing information within the sub-group. If nobody has this information, the host in charge contacts the source for that it retransmits the data by satellite. This approach targets a service very near HSTRM needs. Nevertheless, the choice to use a FEC coding at transport level can modify the properties of considered system. Moreover the study does not tackle at all the problems related to the P2P network setup, the choice of the host in charge, and the search for information. Finally, loss can lead to satellite retransmissions, while in our case all data must be transferred by means of the ground network.

The objective for a large number of works is to minimize the duration of the information downloading. To satisfy this need, several mechanisms are used to make the P2P network to converge towards a more powerful network. However, in our particular case study, among the whole number of peers, a large number of receivers will not make any requests (i.e. the number of receivers that entirely got the data is large from the main principle of the proposed approach). It is thus useless for the P2P network to arrange the peers in such a way to bring them closer between each other: they just must remain reachable.

Among all referred work, none takes into account the constraint of cost related to the use of the network. However, taking into account such constraint makes impossible for example the use of mechanisms allowing making the network of peers to converge (because these mechanisms imply a considerable additional use of the network). That also explains the little number of studies studying the setting of the P2P network in the error recovery context: the number of generated messages is not a criterion taken into account.

There is one considerable advantage related to the use of the hybrid approach: after the satellite diffusion phase, the FEC encoded information is naturally present in most of the peers. Consequently, there is no need of any replication and dissemination policies (see for example [START_REF] Rhea | Maintenance-free global data storage[END_REF][START_REF] Lancérica | Evaluation of content-acces QoS for various dissemination strategies in peer to peer networks[END_REF][START_REF] Cuenca-Acuna | PlanetP : Using gossiping and random replication to support reliable peerto-peer content search and retrieval[END_REF]). On the contrary, this natural dissemination of encoded data should be used to enhance the recovery phase. Nevertheless, the hybrid approach is associated to a certain number of problems specific which are not considered in the analyzed approaches. For example, none of the studied systems consider a large proportion of users already having the data (the opposite case is generally considered). This assumption allows simplifying the searching mechanisms. Finally the assumption of transmission in a (very) large groups, associated to the fact that this network is set up at particular instant (when satellite transmission stops), induces scalability problems. The problems presented above leads to propose a set of ad-hoc mechanisms tailored to adapt the studied context. These mechanisms are detailed in the following section.

C. Mechanisms for the terrestrial recovery phase

Principal objectives of the terrestrial mechanism for recovery errors come directly from the preceding remarks. They must be designed:

to allow information finding among the whole receiver group; to avoid an expensive use of the terrestrial network;

to be used in a large scale context, without implying unacceptable waiting times from the receivers point of view. To achieve these aims, a P2P network is established so that the receivers can communicate between them. This network consists in a regular logical v-ary tree (see Figure 2). the network creation is simple; this limits the number of known peers from a given receiver, which avoids the risks of message implosion and distributes the storage need to memorize the state over the tree; the tree traversal is relatively fast.

In the following mechanisms presentation, the number of logical trees is notedτ , the number of peers in each tree is i η , and the height of trees is i h (with

τ ≤ ≤ i 1
). To constitute this network, it is necessary to have roots (which are the base of the network). Then, the source collects the addresses of the set of receivers which answers to the Inquiry messages during the initial phase of HSTRM [START_REF] De Belleville | Transport multipoint fiable à très grande échelle : intégration de critères de coût en environnement Internet hybride satellite / terrestre[END_REF] (those answers are also used for the estimation of the size of the group [START_REF] De Belleville | Group size estimation for hybrid satellite/terrestrial reliable multicast[END_REF]). It then transmits these addresses to the entire group by means of messages P2PRootList. Those messages are sent every P2PRLPeriod units of time, allowing lossy receivers to get them. The source sends with this list a duration TMaxP2P, corresponding to the time the receivers have to be connected.

When the reception of the first message P2PRootList occurs, the receivers compute randomly a latency ranging between 0 and TMaxP2P. When the associated timer expires, the receivers choose randomly a root among the received list, and sends a message ConnectP2P to the selected root, specifying their IP address. If the selected root does not answer, another root is choosen.

The first reachable root will answer by a message PeerConnected, if it has less than v neighbours. In this case the considered receiver becomes a direct root neighbour. In the other case, the root returns a message Redirected which specifies a dependent peer, to which the newcomer receiver must contact. The receiver sends to the specified peer a message ConnectP2P. The peer answers by a message PeerConnected. The peer that the newcomer tries to reach can however be not reachable. In this case, the new peer will then contact its root which will retransmit the address of another peer.

These operations are illustrated in with an example where v=3 ( i @ represents the address of receiver i). Moreover, in such a way to avoid a partition of the network because of a node breakdown, the peers keep in memory a certain number of potential root addresses. Thus, if a receiver becomes isolated (i.e., its neighbours are not reachable), it can always contact the other roots and to have access to the remaining of the P2P network. Any peer thus has access to the whole network.

IV. Parameter settings and Performance study

The following sections presents the choices on v,τ and TMaxP2P which was determined either by the way of a theoretical study or by simulations. The cost of the P2P network setup, the memory required to represent the network state, and the temporal spacing of connections are studied through theoretical calculations. Then, the cost ) and the efficiency of information retrieval are studied by the way of simulations. Let us note that the results obtained by means of theoretical calculations were also confirmed by simulation.

A. Cost to setup the P2P network

It is possible to compute the number of messages exchanged during the establishment of the P2P network, when the process described in the paragraph III.C is used:

-For the τ roots, no message is used; -For the υ τ × peers directly connected too the τ roots, 2 messages are exchanged ; -For the ( )

1 + × - υ τ
( ) ( ) ( ) ( ) 1 4 2 , , 2 + × - × + × × × = υ τ υ τ α τ υ R R F U P P (6)
This value constraints the minimum size of the data that can be transmitted so that the hybrid communication remains advantageous. According to the paragraph II.D, the cost related to the transmission of K packets is ( )

K R F X , .
And when S bites are to be sent, ( )

⎡ ⎤ TDU K S × /
blocs of K packets are needed so that the whole file is transmitted. In consequence, the cost file C related to the transmission of the file is:

( ) ⎡ ⎤ ( ) K R F TDU K S C X file , / × × = (7)
It is then possible to compute the minimum size S so that the cost to setup the P2P network does not exceed a certain percentage of the file transmission cost. For example with R=600,000, 10 = υ and 166 = τ (those values are explained below), the condition ( )

file P P C R F 1 . 0 , , 2 ≤ τ υ
implies that more than 60 MB have to be sent.

B. Memory requirements

According to the process described in the paragraph III.C, all the root either accept the newcomers, or redirect them toward a peer located at the bottom of the logical tree. In consequence the roots have to keep the addresses of the peers located at the bottom of the tree (see Figure 4). This paragraph focuses on the evaluation of the memory required for that purpose.

With the notations introduced into the Figure 2, the following inequalities are respected:

∑ ∑ = - = ≤ < i h i i h i υ η υ 0 1 0
. In order to approximate the amount of memory needed, we only consider the case when the equality is verified (this is a conservative approach):

( ) ( ) 1 ln / 1 1 ln - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - ∝ υ υ τ υ R L
. The objective is then to find convenient values for τ and υ. When υ increases, the number of neighbours increases for each peer, and, as a consequence, the robustness of the overall network improves. Nevertheless, the number of addresses stored increases for each root. The parameter τ has equally a significant impact: when it increases L decreases. In exchange the network becomes more divided. The Figure 5 shows the amount of memory used for the addresses storage for different values of τ and υ.

As it can be shown in the Figure 5, the parameter τ allows reducing drastically the memory requirements. Once this parameter is fixed, v has only small effects on memory space. As a result, the choice of 30 = τ when R=600000 requests only less than 80Ko per root. In addition, each node of the tree has to keep also IP addresses of a number of peers to ensure a robust connectivity. As a result, the nodes needs to keep the addresses of known roots, the addresses of all intermediary nodes between them and the root and the addresses of their child. The number of addresses stored is v h + + τ . For realist values ofτ , either this size if negligible in regards to L, or R is sufficiently low for the number of stored addresses to be not a real problem. Only constraints on L have to be kept in mind.

C. Connections spacing

The choice of τ is not only dependant of memory considerations. In fact, the value of this parameter has an important impact over the number of requests sent to the roots. The principle of the connection spacing mechanism is fairly simple: consists to distribute over time the connection requests from the receivers, in such a way to avoid root congestions. The very classical Nack Suppression mechanism is used as follow: when a P2PRootList is received, each receiver that are not part of the list wait for a time between 0 and TMAXP2P before replying. The present section studies the issues related to the choice of τ and TMAXP2P. The average number of connections requests C A that each root receive is connected to these two parameters:

( )

P TMAXP R A C 2 / × = τ .

Several constraints must then be taken into account:

-C A must not be higher than a predefined threshold

Max C A .
-The time TMAXP2P which represents the P2P network setup time must not be higher than a maximum, denoted as TMAXP2P Max , so that the time to establish the P2P network remains acceptable.

-τ should not be more than the average number of answers collected for a single Inquiry message so that a simple exchange of message is sufficient to get all the addresses needed. Denote as moy τ this average number.

τ should be chosen in such a way that each logical tree contain at least η min peers (i.e. the network is not too divided in regards to the total number of peers).

Considering the previous constraints, the choice of the parameter τ and TMAXP2P can be done with the following algorithm: After the definition of this algorithm, our objective was to study the performances of the network established. For that purpose, several simulations were conducted.

D. Simulation results

A simulation study has been achieved to verify the effectiveness of the techniques described in the previous paragraphs. Simulations have been carried out also to determine a suitable value for the number of neighbors. The simulation software has been developed with Java language. All the simulations have been setup with the following parameters:

-30 = moy τ : This value is coherent with the group size estimation study provided in [START_REF] De Belleville | Group size estimation for hybrid satellite/terrestrial reliable multicast[END_REF]. Then, when the number of root is lesser than moy A large number of simulations have been achieved in this context, in a first time when all the receivers of the initial group are reachable during the terrestrial phase, and in a second time when a subset of the original receivers are unreachable. In this context, the studies cover various performance aspects of the terrestrial phase, particularly the memory requirements, the load over the network, the cost linked to the information localization and retrieval. The second set of simulations was also addressing robustness issues in evaluating the impact on the overall mechanism of a set of unreachable receivers. In this paper, we only consider a subset of the studies achieved when all receivers are reachable. Only the results concerning cost and efficiency of information search is provided. The complete results of study aiming at determining the impact of a set of non reachable receivers on the provided service is available in [START_REF] De Belleville | Transport multipoint fiable à très grande échelle : intégration de critères de coût en environnement Internet hybride satellite / terrestre[END_REF].

Information searching cost

The objective of this set of simulations is to ensure that the cost related to the search of data is not too high; so that the proposed hybrid approach remains advantageous. In the achieved simulations, MinRCV (the value of MinRCV is chosen according to [START_REF] De Belleville | Une Approche Hybride Satellite/Terrestre pour le transport fiable multipoint à grande échelle[END_REF]) are randomly chosen among all the receivers: they represent the receivers which do not receive the entire information. In consequence, they contact other peer in the overlay network to retrieve lacking information. The number of search messages is then computed. It appears in the experiences that when the group size increases, the number of messages decreases. This behaviour is due to the fact that the ratio of lossy receivers decreases when the group size increases (MinRCV is independent of the group size). In consequence the probability that a peer is connected to another peer who received the whole information increases. Figure 7 shows the number of generated messages for different values of υ with 600 receivers and 300 lossy receivers (each curve represents the average curve obtained from 300 Monte Carlo runs). According to the Figure 7 several cases may be distinguished:

-The receivers who are a leaf of the logical tree and whose father received all the data generate only one search message. -The receivers located at the middle of the tree and whose father and son received all the data generate ( )

1 + υ search messages.
-The remaining lossy receivers generate more messages because the request is propagated in the tree. As one can see on this figure the influence of the parameter υ decreases when the group size increases. Moreover, when the group is more than 6,000 receivers, the cost is almost independent ofυ . As a consequence, the value of υ should be fixed as 10 because it increases the robustness of the network and the efficiency of the mechanism.

Searching performance

We evaluate the performance as i) the number of peers who successfully retrieve the information and ii) the minimum number of hops realised in the overlay network before the information is found. For the first criterion, during all the simulations realized, no receiver was unable to retrieve the information. In consequence the mechanism can be considered as very efficient.

The second criterion characterize the minimum search time before the information is successfully retrieved. For the reasons evoked previously, results are less satisfying with small groups. Figure 9: represents the values obtained with a small group (R=600) and for different values of υ. According to the results when υ increases i) the number of peers who do not find the data after two hops decreases, and ii) the maximum number of hops required decreases. This is merely due to the fact that the number of peers contacted after two hops increases when υ increases.

Finally to characterize the overall efficiency of the mechanism, we considered the maximum of the minimum number of hops required to retrieve the information (i.e. the maximum of the curves represented in the Figure 9: ). The Figure 10 shows the evolution of the maxima obtained as a function of the group size. As expected for small groups, the maximum decreases as υ increases. On the other hand, the trend is inverted for large groups. Nevertheless, the influence of υ is really important for small group, and quite limited for large groups. As a consequence the value of υ should be chosen as high as possible. According to the results presented in this paper, υ=10 is a satisfying value because the mechanism is robust and efficient at low cost.

V. Conclusion

According to the proposed hybrid approach, only a small part of the initial receivers group has to retrieve data using the terrestrial network. Moreover, most of the receivers' just need a few part of the FEC encoded information. In this context we propose a mechanism which takes advantage of the satellite transmission to establish a logical network composed of T v-ary trees. In particular, this mechanism uses timers to avoid the saturation of the terrestrial network when large groups are involved in the establishment of the network. Then, peers can search for data contacting their neighbors in the logical structure. The mechanisms are completely defined, configured using theoretical results, and evaluated through simulations. According to the simulation results, the mechanisms allow to efficiently find data in the network at very low cost.

Figure 1 :

 1 Figure 1: Cost generated by terrestrial, satellite and hybrid communications. α TU = α TM = 1, α SM = 100 and K = 100

Figure 2 .

 2 Figure 2. Logical representation of the P2P network. Several advantages are related to such receiver's organization:the network creation is simple; this limits the number of known peers from a given receiver, which avoids the risks of message implosion

Figure 3 .

 3 Figure 3. Connection of a newcomer into the P2P network.

Figure 5 :

 5 Figure 5: Memory requirements (R = 600 000).

Figure 4

 4 Figure 4: addresses kept in memory ( 3 , 15 = = υ η

  setup the P2P network. According to the previous remarks, the expression of α TU is the terrestrial per-packet transmission cost):

  context, R receivers are distributed among τ logical trees. Thus each tree contains on average τ / R peers. In consequence each root have to keep L addresses in memory, where

Figure 6 :

 6 Figure 6: Algorithm for the choice of τ and TMAXP2P

τ,

  only one Inquiry message is necessary to obtain all the necessary addresses. -20 = cMAX A message per seconds corresponding approx. to a low throughput of 3kb/s

Figure 7 :Figure 8 :

 78 Figure 7: number of search messages generated for each lossy receivers. R=600, MinRCV=300 Figure 8: total number of search messages generated

Figure 8

 8 presents the evolution of the total number of search messages generated as a function of the group size for different values ofυ .

Figure 9 :

 9 Figure 9: Minimum number of hops for each receivers. R=600, MinRCV=300.

Figure 10 :

 10 Figure 10: Minimum number of hops for each receivers. R=600, MinRCV=300.