Nicolas Chevaugeon
email: nicolas.chevaugeon@ec-nantes.fr

Erwan Verron

Bernard Peseux

FINITE ELEMENT ANALYSIS OF NONLINEAR TRANSVERSELY ISOTROPIC HYPERELASTIC MEMBRANES FOR THERMOFORMING APPLICATIONS

Keywords: Hyperelasticity, Thermoforming, Fibers-reinforced material, Transversely isotropic constitutive equation

come

Finite element analysis of nonlinear transversely isotropic hyperelastic membranes for thermoforming applications

INTRODUCTION

The finite element simulation of thermoforming process for isotropic materials has been extensively examined in the past. Both axisymmetrical and three-dimensional models have been developed. In 1989, Zamani et al. [START_REF] Zamani | Status of the finite element method in the thermoforming process[END_REF] reviewed publications dealing with this subject. In the present introduction, we only focus on the more recent works on the subject. deLorenzi and Nied [START_REF] Delorenzi | Blow molding and thermoforming of plastics: finite element modeling[END_REF][START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF] examine both blow-moulding and thermoforming processes simulation. Using experimental observations, they assume that the thermoplastic parison can be modeled by a rubberlike hyperelastic membrane. In the same period, Charrier and coworkers [START_REF] Charrier | Free and constrained inflation of elastic membranes in relation to thermoforming -non-axisymmetric problems[END_REF][START_REF] Shrivastava | Large deformation finite element analysis of non-linear viscoelastic membranes with reference to thermoforming[END_REF] used the same approach for viscoelastic membranes. All these papers report on a quasi-static treatment of the problem. More recently, Bourgin et al. [START_REF] Bourgin | A first step towards the modelling of the thermoforming of plastic sheets[END_REF] and Verron et al. [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF][START_REF] Marckmann | Finite element analysis of blow-moulding and thermoforming process using a dynamic explicit procedure[END_REF] prefer to consider inertial effects and use a dynamic explicit numerical procedure to solve the problem.

All these works are focused on isotropic materials. There are just few studies concerning the thermoforming of fibers-reinforced thermoplastics. Some authors use the approach developed in injection molding process of fiber loaded material [START_REF] Bhattacharyya | Composite sheet forming[END_REF]. The material is described at a micromechanical scale (unit cell model) and is considered as a viscous fluid loaded by rigid fibers. From this microscopic point of view and using homogeneization technics, authors build macroscopic constitutive equations, formulated in terms of strain rate. The corresponding behaviour is anisotropic viscous fluid like [START_REF] Okine | Numerical modeling of sheet forming processes for thermoplastic composites[END_REF][START_REF] De Luca | An explicit finite element solution for the forming prediction of continuous fibre-reinforced thermoplastic sheets[END_REF][START_REF] Pickett | Numerical simulation of the forming process in long fibre reinforced thermoplastics[END_REF].

The main goal of our work consists in using the classical developments relative to thermoforming of isotropic materials and in adapting them to the general case of anisotropic materials. Following the work of Spencer [START_REF] Spencer | Continiuum Theory of the mechanics of fibre reinforced Composite[END_REF], we choose to describe the anisotropy induced by aligned fibers in the rubberlike thermoplastics at the macromechanical scale using material symmetry properties. The present paper reports on the use of transversely hyperelastic constitutive models in the thermoforming simulation. Such theory was successfully used by Kyriacou et al. [START_REF] Kyriacou | Finite element analysis of nonlinear orthotropic hyperelastic membranes[END_REF] in membrane inflation problems.

In the next section the theoretical background used to describe stress and strain for finite deformation in the Lagrangian point of view is recalled, and the general form of the strain energy function for hyperelastic transversely isotropic materials is established. Section 3 is devoted to the presentation of the finite element formulation suitable to describe the inflation of membrane under an inflating pressure. The numerical implementation of the constitutive equation is highlighted. Some numerical results obtained by our method for free and confined inflation cases are presented in section 4. Last, concluding remarks and perspectives are proposed.

HYPERELASTIC TRANSVERSELY ISOTROPIC CONSTITUTIVE EQUA-TION

In this paper, fibers-reinforced thermoplastic materials are considered. Fibers are assumed to be aligned and materials are studied far above their glass transition temperature. Therefore materials obey hyperelastic transversely isotropic constitutive equations.

Stress and strain definitions

Consider a material point P of a body B. P occupies the position X relative to a frame of reference R 0 : {O, e I } in the undeformed initial configuration C 0 . In the deformed configuration C, x is the position of P relative to an observer frame of reference R : {o, e i }. The transformation between C 0 and C is defined by:

x = ξ(X, t) ( 1
)
where t is the present time.

In the neighborhood of P, the deformation gradient F, the right Cauchy-Green deformation tensor C and the Green-Lagrange strain tensor E are respectively defined by:

F(X, t) = ∂x(X, t) ∂X (2) C(X, t) = F T F (3) E(X, t) = 1 2 [C(X, t) -I] ( 4
)
where I is the second-order identity tensor and • T denotes the transposition.

Let σ and S be the Cauchy "true" stress tensor defined in the deformed configuration and the second Piola-Kirchhoff stress tensor relative to the undeformed configuration. They are respectively defined by: df = σ . n ds (5) dF = S . N dS [START_REF] Bourgin | A first step towards the modelling of the thermoforming of plastic sheets[END_REF] in which df is the force in the deformed configuration acting on a deformed oriented area n ds and dF is the force in the undeformed configuration acting on an undeformed oriented area N dS. df and n ds are the transformed vectors of dF and N dS by the deformation gradient. With some standard manipulations the two stress tensors can be related by:

S = J F -1 σ F -T ( 7
)
where J is the Jacobian of the transformation (for a volume-preserving deformation: J = 1).

Construction of a Strain energy function

Hyperelastic materials are defined by the existence of a scalar function W (F) called stored energy or strain energy function, from which stresses can be derived at each point. In order to satisfy the objectivity principle, the strain energy function must be invariant under changes of observer frame of reference. It is well-known that the Cauchy-Green deformation tensor is invariant under changes of observer frame of reference [START_REF] Beatty | Topics in finite elasticity: hyperelasticity of rubber, elastomers and biological tissues -with examples[END_REF]. Thus, if the strain energy function W can be written as a function of C, it automatically satisfies the objectivity principle.

We now have to take into account the material symmetries that restrict the way W depends on C. Every orthogonal transformation member of the material symmetry group must keep the strain energy unchanged. For isotropic material, the material symmetry group consists in the entire group of orthogonal transformations. In the case of unidirectional fibers-reinforced composite, the symmetry group has to be determined. The fibers orientation is defined by a local unit vector denoted a 0 (X). Therefore, the material symmetry group is the group of orthogonal transformations which keep a 0 (X) unchanged. Consequently the material is isotropic in all planes orthogonal to the fibers direction and is called transversely isotropic. The general form of the corresponding strain energy function must satisfy:

W (C) = W (Q.C.Q T ) ∀Q orthogonal with Q.a 0 (X) = a 0 (X) ( 8 )
Another way of writing W consists in introducing explicitly the dependence on a 0 [START_REF] Spencer | Continiuum Theory of the mechanics of fibre reinforced Composite[END_REF][START_REF] Govindjee | Finite implementation of incompressible, transversaly isotropic hyperelasticity[END_REF].

The invariance with respect to the symmetry group is now expressed by:

W (C, a 0 ) = W (Q.C.Q T , Q.a 0 ) ∀Q orthogonal (9)
As the sign of a 0 has no physical meaning, we can introduce the orientation tensor A = a 0 ⊗ a 0 with no loss of information (⊗ stands for the tensor outer product). Then the previous equation ( 9) becomes:

W (C, A) = W (Q.C.Q T , Q.A.Q T ) ∀Q orthogonal ( 10
)
Thus W is an isotropic function of two symmetric tensors C and A. Using representation theorems [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], W can be written as a function of the invariants of C, A and their products. Following equations give a possible set of ten independent scalar invariants on which W can depend [START_REF] Sidoroff | Cours sur les grandes transformations[END_REF]:

I 1 = trC I 2 = 1 2 (trC) 2 -trC 2 I 3 = detC (11) A 1 = trA A 2 = 1 2 (trA) 2 -trA 2 A 3 = detA ( 12
)
I 4 = tr(C : A) I 5 = tr(C 2 : A) I 6 = tr(C : A 2 ) ( 13
)
I 7 = tr(C 2 : A 2 ) ( 14
)
where (I i ) i=1,3 are the invariants of C, (A i ) i=1,3 those of A and (I i ) i=4,7 are product invariants.

In the present case of transversely isotropic behaviour, A = A 2 and the invariants of A become:

A 1 = 1 A 2 = 0 A 3 = 0 (15)
Moreover, some relations can be written for the invariants of the products:

I 5 = I 7 I 4 = I 6 (16)
Therefore, the number of relevant independent invariants on which W depends is reduce to five [START_REF] Burton | A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations[END_REF]:

W = W (I 1 , I 2 , I 3 , I 4 , I 5 ) (17)
Note that in the isotropic case, the number of invariants is reduce to three [START_REF] Treloar | The mechanics of rubber elasticity[END_REF]:

I 1 , I 2 and I 3 .
Taking into account the incompressibility assumption, i.e.: I 3 = 1, the number of independent invariants is finally reduced to four.

The 2nd Piola-Kirchhoff stress tensor can be calculated by differentiation of the strain energy function with respect to the Cauchy-Green dilatation tensor C [START_REF] Beatty | Topics in finite elasticity: hyperelasticity of rubber, elastomers and biological tissues -with examples[END_REF]:

S = -pC -1 + 2 ∂W ∂C ( 18
)
where p is an arbitrary hydrostatic pressure preventing for any volume change. This pressure has no physical meaning and is determined from equilibrium equations. Using the differentiation chain rule, the stress tensor can be expressed as:

S = -pC -1 + 2 i=1,2,4,5 ∂W ∂I i ∂I i ∂C ( 19
)
with:

∂I 1 ∂C = I ( 20
)
∂I 2 ∂C = I 1 • I -C ( 21
)
∂I 3 ∂C = I 3 • C -1 (22) ∂I 4 ∂C = a 0 ⊗ a 0 ( 23
)
∂I 5 ∂C = a 0 ⊗ C • a 0 + a 0 • C ⊗ a 0 (24)

FINITE ELEMENT FORMULATION

In this section, we briefly present the finite element formulation adopted to simulate the thermoforming process. For more details, reader can refer to [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF][START_REF] Marckmann | Finite element analysis of blow-moulding and thermoforming process using a dynamic explicit procedure[END_REF].

Governing equations

The polymeric sheet is supposed to be a membrane and is represented by a twodimensional continuum. Its thickness becomes a function of the position on the midsurface. Assuming that there are no body forces, the Principle of Virtual Work, is expressed on the undeformed configuration for the inertial effects and internal work, and on the deformed configuration for the pressure follower force work:

B 0 δu • ρ 0 üdV 0 = - B 0 δE : SdV 0 + δB δu.TdS ∀δu ( 25
)
where B 0 and δB stand respectively for the volume of the undeformed membrane and the surface of the deformed membrane, ρ 0 is the mass density, ü is the acceleration vector, T is the external surface force due to pressure, and δu is a virtual displacement vector compatible with displacement boundary conditions. The contact with the mold is assumed sticky: the heated parison cools down and becomes stiffer after contact [START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF]. Thus the blowing pressure is not sufficient to deform the membrane. Every material points which enter in contact is fixed until the end of the simulation. This assumption simplifies the numerical procedure [START_REF] Khayat | Multiple contact and axisymmetric inflation of a hyperelastic cylindrical membrane[END_REF].

Spatial discretization

The previous equation ( 25) is discretized by the finite element method. In the membrane context, 3-nodes isoparametrical triangular elements are chosen. These elements only deform in their plane, remain flat and triangular The use of these simple elements provides analytical integration which is time saving. But it has a major flow: lot of elements are needed to obtain a good stiffness in high curvature regions of the parison. To avoid this problem, an automatic remeshing technique described in [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF] is implemented. After discretization and integration on each element surface, the problem reduces to a system of ordinary second order differential equations:

[M]{ Ü (t)} = F ext (t) -F int (t) (26)
where [M] stands for the mass matrix, U(t) is the nodal displacement vector, F ext (t),and F int (t) are respectively the external and internal nodal force vectors. By the use of the special lumping technique [START_REF] Taylor | The finite element method. Volume I: Basic formulation and linear problems[END_REF], [M] can be reduced to a diagonal matrix, splitting the system in n independent differential equations where n is the number of degrees of freedom of the model.

Time-integration scheme

In order to solve the previous system (26), we use the classical explicit second-order central difference method. The velocity of the degree of freedom i, Ui (t), and its acceleration, Üi (t), are given by:

Ui (t) = U i (t + ∆t) -U i (t -∆t) 2∆t (27) Üi (t) = U i (t + ∆t) -2U i (t) + U i (t -∆t) ∆t 2 ( 28
)
where ∆t is the time step. Consequently, the problem ( 26) is reduced to the following system:

U i (t + ∆t) = ∆t 2 M d ii F ext i (t) -F int i (t) + 2U i (t) -U i (t -∆t) (29)
for each degree of freedom i.

It is well-known that the stability of this scheme depends on the time step ∆t. ∆t must be less than the critical time step ∆t cr . ∆t cr is the smallest time step required for a wave to cross an element. As the present system is highly non-linear (large strains and material behaviour), stiffness properties change during calculations. To ensure the convergence of the scheme, ∆t cr is periodically estimated by determining the stiffest element in the whole mesh.

Constitutive equation implementation

In order to calculate the internal force term in Eq. ( 29), the relation between deformations and stresses must be computed in each element. Using the plane stress state assumption, deformation and stress matrices in a single triangular element have the following forms:

[C] =   C 11 C 12 0 C 12 C 22 0 0 0 C 33   (30)
and:

[S] =   S 11 S 12 0 S 12 S 22 0 0 0 0   (31)
where 1, 2, 3 represent the local axis in the element plane. The direction 3 is the normal to the surface.

Using nodal initial positions and displacements in the global referential axis, C 11 , C 12 and C 22 are easily computed in each element. From the incompressibility constraint detC = 1, we have:

C 33 = 1 C 11 C 22 -C 2 12 (32)
By definition, C 33 is the square of the principal stretch ratio in the membrane thickness direction:

C 33 = h h 0 ( 33
)
where h and h 0 are respectively the deformed and undeformed thicknesses of the membrane. By now, all terms of the strain tensor can be computed in each element. The components of A in element local axis are computed using the fibers direction a 0 previously defined in each undeformed element:

[A] =   a 2 01 a 01 a 02 0 a 01 a 02 a 2 02 0 0 0 0   (34)
All terms corresponding to the third direction are null because the fibers always stay in the element plane. For a given expression of the strain energy function, S can now be calculated in the element local frame. The hydrostatic pressure p is determined by using:

S 33 = 0 (35)
In order to illustrate this general method, a very simple strain energy function is chosen:

W = Φ(I 1 -3) + Ψ(I 2 -3) + ΘI 4 (36)
where Φ, Ψ and Θ are the material parameters. This constitutive equation is a tranversely isotropic generalization of the widely used Mooney-Rivlin model [START_REF] Treloar | The mechanics of rubber elasticity[END_REF]. Using Eqs ( 18), (32), ( 35) and (36), the hydrostatic pressure p is given by:

p = 2 Φ + Ψ(C 11 + C 22 ) C 11 C 22 -C 2 12 ( 37
)
And the stress components becomes:

S 11 = 2(-C 22 C 2 33 [Φ + Ψ(C 11 + C 22 )] + Φ + Ψ(C 22 + C 33 ) + ΘA 11 ) (38) S 22 = 2(-C 11 C 2 33 [Φ + Ψ(C 11 + C 22 )] + Φ + Ψ(C 11 + C 33 ) + ΘA 22 ) (39) S 12 = 2(C 12 C 2 33 [Φ + Ψ(C 11 + C 22 )] + Φ + Ψ(-C 12 ) + ΘA 12 ) (40)
where C 33 is given by Eq. (32).

NUMERICAL EXAMPLES

Validation: free inflation

In order to validate the implementation of the previous constitutive equation, our results are compared with those obtained by Kyriacou et al. [START_REF] Kyriacou | Finite element analysis of nonlinear orthotropic hyperelastic membranes[END_REF]. They use a quasi-static approach with four nodes elements to solve membrane inflation problems.

The free inflation of an initially flat square membrane submitted to an internal pressure is examined. The length of the undeformed membrane is denoted l 0 and its thickness h 0 . The inflating pressure is noted P . The material obey the modified Mooney-Rivlin model previously defined. To simplify the discussion the dimensionless inflating pressure P * is defined by:

P * = P l 0 Φh 0 (41)
Φ is the first material parameter. The square membrane is first biaxially pre-stretched in its undeformed plane until the deformed square length l reaches 1.1 l 0 . Kyriacou et al. use this pre-stretching to avoid any singularity in the tangent stiffness matrix at the first loading step. It is to note that no pre-streching is needed to start the computation in our dynamical approach. Numerical values are set to:

Φ = 1. Ψ = 0. Θ = 1. H = 10 -3 L = 1. ( 42
)
In this example, the square membrane is meshed with 400 triangular elements. Kyriacou et al. use 100 four nodes finite elements. Moreover, when the fibers direction coincides with one of the global axis directions, i.e.: a 0 = X or a 0 = Y, only a quarter of the membrane is meshed using symmetry properties. First, the inflating pressure P * is increased from 1.0 to 9.5. Figure 1 shows the height of the center node of the square sheet versus the inflating pressure. At the maximum pressure, the height is equal to 0.52. Our results are identical to the previously published data. Figure 2 shows the height of the center node at P * = 6. for different fibers directions. In the figure, a relative angle between one side of the square and the fibers direction equal to zero corresponds to a 0 = X. The curve has a sinus-like shape and is symmetrical with respect to the 45 o angle because of the geometrical symmetry of the problem. The membrane is less stiff when fibers are aligned with the diagonal direction of the square. The relative difference in height between our curve and Kyriacou's one does not exceed 0.03%.

It is well-known that the pressure-displacement curve corresponding to the free inflation of a rubberlike membrane exhibits limit points [START_REF] Khayat | Inflation of an elastic cylindrical membrane: non-linear deformation and instability[END_REF]. In order to be able to reach higher deformation ratios, the load is imposed using a constant gaz flow rate. The pressure at time t is computed step by step by using the Perfect Gaz Law and the volume inside the membrane at t [START_REF] Verron | Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes[END_REF]. Using this continuation technique, the limit point can be overcomed as shown in Figure 3. Thus the membrane reaches higher deformations.

The deformed geometry of a transversely isotropic inflated square membrane for two different fibers directions is presented in Fig. 4. The lines represent the fibers direction. The classical geometrical symmetry of the isotropic case disappears: membranes tend to exhibit a bulge perpendicularly to the fibers direction.

Finally, this example demonstrates that our results show good agreement with those of Kyriacou et al. Moreover, by the use of a gaz flow rate loading instead of a pressure loading, our method is able to go beyond limit points.

Thermoforming application

In the present example, some results concerning a simple thermoforming application are examined. The inflation of an initial circular flat sheet in a cylindrical mold is considered. Here automatic remeshing is used in order to be able to fill the corners of the mold.

As mentioned elsewhere [START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF], one of the main goal of the numerical simulation of such processes is to predict the final thickness distribution of molded parts. In order to evaluate the influence of fibers-reinforcement on the final geometry, comparisons are made between the membrane thickness distributions in the isotropic and the transversely isotropic cases. Material parameters are set to:

Φ = 1. Ψ = 0. Θ = 3. ( 43
)
Figures 5 and6 present respectively four inflation stage of the moulding of a cylinder in the isotropic and the anisotropic cases. Due to geometrical and material symmetries, only a quarter of the membrane is considered. In the isotropic case, the thickness distribution is axially symmetric (see Fig. 5). Figure 6 highlights the loss of symmetry due to the presence of fibers. The fibers direction is represented by the lines on the parison. It is to note that in this case axisymmetrical calculations would be unable to predict the final thickness distribution. In the transversely isotropic case, fibers are oriented as in Figure 6 and distributions differ for the two lines. The results corresponding to lines XZ and Y Z are respectively plotted with a pointed line and a dashed line. In both cases, isotropic and transversely isotropic (for the two lines), the highest thickness is observed near the clamped edge and the smaller thickness corresponds to the corner of the mold. An important remark found in the bibliography [START_REF] Delorenzi | Finite element simulation of thermoforming and blow molding[END_REF] relative to the thermoforming of isotropic material concerns the membrane thickness distribution at the end of the process. It is said that this distribution does not depend on the material model. But as pointed out in the present example, the anisotropic material behaviour due to the presence of fibers highly influence thickness distribution in thermoformed parts.

CONCLUDING REMARKS

In this paper a general form of the strain energy function for transversely isotropic hyperelastic materials has been developed. This function depends on a limited set of tensor invariants and ensure the objectivity principle. Moreover as shown by Kyriacou et al., the formulation of the constitutive equation is easy implemented in a quasi-static finite element program. Here it is shown that the implementation in a dynamic explicit finite element code is easier. Efficiency of the method is demonstrated on both free and confined inflation problems.

As a perspective, we are working on a viscoelastic counterpart of the present hyperelastic model and on models for misaligned thermoplastic composites using a more general orientation tensor. In order to determinate the material parameters, a biaxial rheometer for high strain rates is in development.

Figure 1 :Figure 2 :

 12 Figure 1: Height of the center node versus inflating pressure for a square membrane

Figure 3 :Figure 4 :

 34 Figure 3: Height of the center node versus inflating pressure for a square membrane loaded with a constant gaz flow rate

Figure 5 :

 5 Figure 5: Thermoforming stages of a circular plane sheet in a cylindrical box: isotropic case. Gray scale: relative membrane thickness

Figure 6 :

 6 Figure 6: Thermoforming stages of a circular plane sheet in a cylindrical box: transversely isotropic case. Gray scale: relative membrane thickness

Figure 7 Figure 7 :

 77 Figure 7 presents the membrane thickness evolution along two edges shown in the small figure in the corner. The results correspond to the previous simulations. In the


Keywords: Robust Linear Programming, duality, robust PERT scheduling

quite different model, bears some resemblance to the well-known Bertsimas-Sim approach to robustness). We show that in both cases, the resulting robust optimization problem can be efficiently solved in polynomial time.

Duality, Robustness, and 2-stage robust LP decision models. Application to Robust PERT Scheduling

Michel Minoux *

Résumé

Les modèles de programmation linéaire robuste étudiés dans la littérature relèvent principalement, soit de l'incertitude par colonnes, soit de l'incertitude par lignes. Ils supposent essentiellement que les colonnes (resp : les lignes) de la matrice des contraintes ont des coefficients pouvant fluctuer, tout en restant inclus dans des ensembles d'incertitude bien définis. On commence par étudier la question de savoir si la théorie de la dualité en programmation linéaire peut être utilisée pour transformer un problème de PL robuste avec incertitude disons par lignes, en un problème de PL robuste avec incertitude par colonnes. On exhibe des exemples simples montrant clairement que le dual d'un programme linéaire robuste donné n'est en général PAS équivalent à la version robuste du dual (en considérant bien sûr que les ensembles d'incertitude pour les lignes -resp : les colonnes -du primal sont définis exactement de la même façon que les ensembles d'incertitude pour les colonnesresp : les lignes -du dual). On poursuit l'étude de cette question de la dualité dans le contexte de la robustesse en considérant la sous-classe des programmes linéaires robustes avec incertitude sur le second membre seulement (cette sous-classe ne semble pas, jusqu'ici, avoir fait l'objet d'études significatives). Dans ce contexte nous introduisons le concept de 'modèle de programmation linéaire robuste sur deux étapes', par contraste avec le cas standard (qui pourrait être désigné sous le nom de 'modèle de robustesse sur une étape'), et nous formulons successivement (a) le dual du programme linéaire robuste ; (b) la version robuste du dual. Les expressions de la fonction objectif a optimiser dans les deux cas, bien que présentant des similitudes formelles, apparaissent clairement distinctes. De plus, du point de vue de la difficulté de résolution, un des problèmes est susceptible de résolution efficace (il s'agit d'un problème d'optimisation convexe), tandis que l'autre, en tant que problème non convexe, peut être difficile à résoudre dans le cas général. Comme application du modèle de programmation linéaire robuste à deux étapes proposé ici, nous discutons finalement le problème d'ordonnancement PERT robuste en considérant deux définitions naturelles possibles de l'ensemble d'incertitude pour les durées des tâches : le cas où l'ensemble d'incertitude est, à une mise à l'échelle près, une boule relativement à la norme Linfini ; et le cas où l'ensemble d'incertitude est, à une mise à l'échelle près, une boule de rayon maximum fixé relativement à la distance de Hamming. On montre que, dans les deux cas, le problème d'optimisation robuste peut être résolu efficacement en temps polynomial. On explique enfin pourquoi le modèle et l'approche antérieurement proposés par Bertsimas et Sim (pour le problème du plus court / plus long chemin robuste) ne pourraient s'appliquer au problème d'ordonnancement PERT robuste étudié ici, du fait de structures sous-jacentes essentiellement différentes.

Mots-clefs : Programmation Linéaire robuste, dualité, ordonnancement PERT robuste

Abstract

The various robust linear programming models investigated so far in the literature essentially appear to be based either on what is referred to as 'rowwise' uncertainty models or on 'columnwise' uncertainty models (these basically assume that the rows -resp: the columns -of the constraint matrix are subject to changes within a well specified uncertainty set). We first address the question of whether LP duality can be used to convert a robust LP problem with, say, rowwise uncertainty, into a robust LP with columnwise uncertainty. We provide simple examples supporting the general statement that the dual of a given robust linear programming model is NOT equivalent to the robust version of the dual (assuming of course that uncertainty for the rows -resp: the columns -of the primal, is specified in exactly the same way as uncertainty for the corresponding columns -resp: rows -of the dual). Next, we further investigate this issue of duality in the context of robustness by considering the subclass of robust LP models with uncertainty limited to the right handside only. (this subclass does not appear to have been significantly investigated so far). In this context we introduce the concept of ' two-stage robust LP model' as opposed to the standard case (which might be referred to as 'single-stage robust LP model') and we show how to derive both statements of (a) the dual to the robust model and (b) the robust version of the dual. The resulting expressions of the objective function to be optimized in both cases, though bearing some formal resemblance, appear to be clearly distinct. Moreover, from a complexity point-of-view, one appears to be efficiently solvable ( it reduces to a convex optimization problem) whereas the other, as a nonconvex optimization problem, is expected to be computationally difficult in the general case. As an application of the 2-stage robust LP model introduced here, we next investigate the robust PERT scheduling problem, considering two possible natural ways of specifying the uncertainty set for the task durations : the case where the uncertainty set is a scaled ball with respect to the L ∞ norm; the case where the uncertainty set is a scaled Hamming ball of bounded radius (which, though leading to a

Introduction

Various models for handling robustness objectives with respect to uncertainties on some specified coefficients in linear programming models have been proposed in the literature. We can mention [START_REF] Soyster | Convex Programming with Set-Inclusive Constraints and Applications to Inexact Linear Programming[END_REF], Ben-Tal and Nemirovski (1998,2000), Bertsimas andSim (2003, 2004).

The various approaches proposed can roughly be divided into two distinct categories, depending on whether the underlying uncertainty model refers to possible fluctuations on the row vectors of the constraint matrix (we call this 'rowwise uncertainty'), or on column vectors (we call this 'columnwise uncertainty').

Columnwise uncertainty was first considered by [START_REF] Soyster | Convex Programming with Set-Inclusive Constraints and Applications to Inexact Linear Programming[END_REF]. In this model each column A j of the m×n constraint matrix is either supposed to be exactly known, or is only known to belong to a given subset K j ⊂ R m ('uncertainty set'). The cost vector and the right handside are supposed to be certain. A robust solution is a solution which is feasible for all possible choices of the uncertain column vectors in their respective uncertainty sets. With this definition, assuming nonnegativity constraints on all variables of the LP, it can easily be shown that the problem of finding an optimal robust solution reduces to solving an ordinary LP with constraint matrix A = (a i,j ) where, ∀i, j, the coefficient a i,j is defined as:

• a i,j = max v∈K j {v i } in case of a ith constraint of the form ≤ • a i,j = min v∈K j {v i } in case of a ith constraint of the form ≥ Note that the above maximization (or minimization) can easily be carried out if we assume the uncertainty sets either of finite cardinality (and not too big !), or closed convex.

As observed by many authors, a drawback of Soyster's model is that it usually leads to rather conservative solutions, in other words the price to pay for robustness in the above sense is often too high. Contrasting with the above, rowwise uncertainty has attracted more interest and has been studied, among others, by Ben-Tal and Nemirovski (1998,2000), and more recently by Bertsimas andSim (2003, 2004).

Ben-Tal and Nemirovski start with the assumption that each row A i of the constraint matrix belongs to a known uncertainty set consisting of an ellipsoid E i ⊂ R n , and a solution x ∈ R n , x ≥ 0 is said to be robust in this context iff it satisfies:

for all i : A i x ≤ b i , ∀A i ∈ E i .
Ben-Tal and Nemirovski then show that finding an optimal robust solution reduces to solving a conic quadratic problem, which can be done in polynomial time. Obviously a drawback of this approach is that, due to nonlinearity of the resulting robust model, it is not practically applicable to very large problems.

A way to obviate nonlinearity, while retaining the idea of rowwise uncertainty, was proposed by Bertsimas andSim (2003, 2004), considering a slightly different model of uncertainty. More precisely, they assume that each uncertain coefficient a i,j can take values in a given interval [a i,j -α i,j , a i,j + α i,j ] and, for each row i, a positive parameter Γ i > 0 (not larger than the total number of uncertain coefficients in row i) is considered. A solution x is then qualified as Γ-robust (in the sense of Bertsimas and Sim) iff for all i = 1, ..., m this solution satisfies the ith constraint for all possible choices of the coefficients in row i such that at most Γ i of the uncertain coefficients in the row are allowed to deviate from the nominal values a i,j . (note that the above statement implicitly assumes the Γ i parameters to be integers, but Bertsimas and Sim show that a slightly more general definition, allowing for nonintegral values of the Γ i 's can be handled in the same way). With this model of uncertainty, Bertsimas and Sim show that finding an optimal Γ-robust solution can be reduced to solving an ordinary linear program only moderately increased in size, thus opening the way to large scale applications. Moreover the approach readily extends to optimization problems including integrality constraints on all or part of the variables, in that case the robust version of the problem is a MIP, but, again, the resulting robust model is only moderately increased in size as compared with the original model.

Duality and columnwise/rowwise robustness

We address in this section the question of whether duality can prove in any way useful to convert a columnwise uncertain linear program into a rowwise uncertain linear program, assuming of course the same uncertainty model for the columns of the given linear program and for the corresponding rows in the dual. Intuitively, no nice (i.e. strong) duality result is to be expected when taking into account robustness constraints since, in both the primal and the dual, there is a price to pay for uncertainty, therefore: if we maximize in the primal, the robust primal optimal solution value will be (in general strictly) less than the optimal solution value of the 'nominal' primal LP; and minimizing in the dual will lead to a robust dual optimal solution value (in general) strictly larger than the same value.

Let us illustrate the phenomenon on a small typical example. Consider the following LP (a continuous knapsack problem actually) with two uncertain coefficients a 1 and a 2 in the constraint matrix:

Maximize 4x 1 + 3x 2 (P) s.t: a 1 x 1 + a 2 x 2 ≤ 4 x 1 ≥ 0, x 2 ≥ 0 the standard LP dual of which reads: Minimize 4u (D) s.t: a 1 u ≥ 4 a 2 u ≥ 3 u ≥ 0
Let us assume that the uncertainty set for a 1 is the real interval [START_REF] Ben-Tal A | Robust Solution of Linear Programming Problems Contaminated with Uncertain Data[END_REF][START_REF] Bertsimas | Robust Discrete Optimization and Network Flows[END_REF] , the uncertainty set for a 2 is the real interval [START_REF] Ben-Tal A | Robust Convex Optimization[END_REF][START_REF] Ben-Tal A | Robust Solution of Linear Programming Problems Contaminated with Uncertain Data[END_REF], and let us take as definition of a robust solution in both (P) and (D) a solution which is feasible for any possible values of a 1 and a 2 in their respective uncertainty sets. Then it is easily seen that the optimal robust primal solution is x 0 = [0, 2] with corresponding primal objective function value 6; and the optimal robust dual solution is u 0 = 3 with corresponding dual objective function value 12. This example thus clearly shows that no natural extension of the usual properties related to LP duality is to be expected in the context of robust Linear Programming. The results discussed in the following section will provide additional evidence of this.

Duality and robustness for LP's with uncertainty in the RHS

A special case of columnwise uncertainty in Linear Programming is when uncertainty only concerns the coefficients of the right handside (RHS). Such problems frequently arise in practical applications. As a typical example, we mention the robust PERT scheduling problem with uncertainty on the durations of (some of) the tasks, assuming that a robust earliest termination date has to be determined. More precisely, we want to determine the minimum total duration of the project under any possible assignment of task durations, taken in a given uncertainty set. The 2-stage robust model discussed in § 3.2 below will appear to be relevant to such applications.

Consider the following LP:

max c T x Ax ≤ b x ≥ 0
and assume that the right handside b is not known exactly, but only known to belong to some uncertainty set B ⊂ R m . The set B may be finite or infinite (we will introduce additional assumptions on this set when necessary).

Two distinct robustness models for LP's with uncertain RHS will be successively discussed in this section, namely single-stage robust decision models ( § 3.1) and two-stage robust decision models ( § 3.2). For both cases it will be shown that, even in the restricted situation addressed here (uncertainty in the RHS only) one cannot use standard duality theory to convert a columnwise uncertain linear program into a rowwise uncertain linear program while preserving equivalence. This will thus provide additional evidence of the fact that strong duality does not hold in general for uncertain LP's.

Single-stage robust (LP) decision model

We first consider the simplest case where the values of all the decision variables x have to be fixed (taking into account uncertainty) before we get any kind of information on the actual realization of the uncertain parameters. In such a simple model (indeed a special case of Soyster's model) feasibility has to be ensured for any b ∈ B , and the problem to be solved simply reduces to:

max c T x Ax ≤ b x ≥ 0 where, ∀i, b i = min b∈B {b i }
The (standard) LP dual to the above problem reads:

(D1) min u T b u T A ≥ c u ≥ 0
On the other hand, if we consider the dual to:

max c T x Ax ≤ b x ≥ 0
we get:

min u T b u T A ≥ c u ≥ 0
Now consider the robust version of this dual problem where the cost vector b is uncertain and can take any value in B. A simple and natural objective in this context is to find u achieving a minimum value of max u T b over all possible b ∈ B, thus leading to:

(D2) min u max b∈B {u T b} s.t: u T A ≥ c u ≥ 0
It is clearly realized that (D1) and (D2) are completely different optimization problems on the same solution sets since, ∀u ≥ 0, u T b ≤ max b∈B {u T b} with strict inequality holding in the general case.

If the set B is either finite, or closed convex, it is observed that (D2) can be efficiently solved via standard convex optimization techniques since the function:

u → max b∈B {u T b}
is convex. Also, in this case, b can be efficiently computed since min b∈B {b i } too is a convex optimization problem, therefore problem (D1), too, can be efficiently solved.

Two-stage robust (LP) decision model

It frequently arises in applications that the process of decision-making under uncertainty can be decomposed in two successive steps (two-stage decision making) or more (multi-stage decision making). For simplicity of presentation, we restrict here to the case of two-stage decision making. In this case, the set of decision variables x is decomposed (partitioned) into two distinct sets of variables which we denote y and z. The y variables concern the decisions to be taken in the first stage (before knowing anything about which realization of uncertainty will arise) and the z variables concern the decisions to be taken in the second stage (after realization of uncertainty).

Limiting ourselves, to make the discussion easier to follow, to the case where the objective function only depends on the decision variables of the first stage, our decision problem can thus be rewritten:

(I) max γ T y s.t: F y + Gz ≥ b y ≥ 0, z ≥ 0
where γ, b are vectors and F and G are matrices of appropriate dimensions. (Observe that the reason for restricting to the case of an objective not depending on the z variables is only for the sake of simplicity in the presentation, our two-stage robust decision model would readily handle the general case of an objective depending on both the y and z variables). Now, since the RHS b in (I) is uncertain, we have to make our robustness objectives precise. In the sequel, we consider robustness for a solution y by requiring that feasibility can be ensured for any possible RHS b ∈ B by using the second stage decision variables z (by analogy with the terminology used in stochastic programming, the z variables might be referred to as 'recourse' variables). So, if we define Y = {y/y ≥ 0 and ∀b ∈ B, ∃z ≥ 0 : Gz ≤ b -F y}, we want to solve:

max y∈Y {γ T y}.
Note that in the above, for any given robust solution y, the value taken by the z variables depends on which b ∈ B is actually realized. This is an important feature of our model which explains why it can produce less conservative solutions as compared with Soyster's model (see example given in Remark 1 below).

According to Farkas' Lemma, we know that, for fixed b ∈ B, a necessary and sufficient condition for the existence of z ≥ 0 verifying Gz ≤ b -F y is that:

u T (b -F y) ≥ 0 for all u in the polyhedral cone: C = {u/u T G ≥ 0 and u ≥ 0}.
Denoting u 1 , u 2 , ..., u p , the extreme rays of the above cone, the set Y can equivalently be represented as the system of linear inequalities:

(u j ) T F y ≤ (u j ) T b, ∀b ∈ B, ∀j = 1, ...p which is equivalent to: (u j ) T F y ≤ w j = min b∈B (u j ) T b.
The robust 2-stage decision problem is then reformulated as:

(I)' max γ T y s.t: (u j ) T F y ≤ w j , ∀j = 1, ..., p y ≥ 0

Since we are interested in investigating duality in the context of robustness, let us state the (standard) dual to (I)'. So, introducing dual variables λ j , (j = 1, ...p), the dual to (I)' can be written as:

min j λ j w j = j λ j min b∈B (u j ) T b s.t: j λ j F T u j ≥ γ λ j ≥ 0, ∀j = 1, ..., p
and since: {u/u = λ j u j , λ j ≥ 0} = {u/G T u ≥ 0, u ≥ 0}, this can be rewritten as:

(DI)' min u min b∈B u T b s.t:

F T u ≥ γ G T u ≥ 0 u ≥ 0
or , equivalently if we denote W denote the set of solutions to (DI)'

min u∈W min b∈B {u T b} (1)
An a priori different way of using duality in our context would be to take the (standard) LP dual to (I) for fixed b, and then to carry out robustness analysis with respect to the coefficients b of the objective in the dual, allowing b to take all possible values in B . The LP dual to (I) reads:

(DI) min u T b s.t: F T u ≥ γ G T u ≥ 0 u ≥ 0
A natural robust version of (DI) consists in finding the dual solution u minimizing the value of u T b produced by the worst possible b, reads:

min u∈W max b∈B {u T b} (2)
which is to be contrasted with (1): indeed, it is seen that the robust version of the dual significantly differs from the dual of the robust version of the initial (primal) problem (I) because the function of u to be minimized is min b∈B {u T b} in one case, and max b∈B {u T b} in the other case.

It is worth observing that this structural difference between the two functions also implies a difference with respect to the practical solvability of the corresponding problems. The objective function in ( 2) is convex in u, making the robust version of (DI) efficiently solvable, whereas the objective in ( 1) is concave in u, making (DI)' and thus (by standard LP duality) (I)' too, difficult problems in the general case.

Remark 1: As already suggested above, the two-stage robust decision model proposed here is capable of producing less conservative solutions as compared with Soyster's model. The reason for this is that if, for a given uncertainty set B, we consider Soyster's model for problem (I), the problem to be solved is

max y∈Ys {γ T y}
where the set Y s is defined as {y/y ≥ 0 and ∃z ≥ 0 : Gz ≤ b -F y} with b defined as:

∀i, b i = min b∈B {b i }
It is easily seen that Y s ⊆ Y = {y/y ≥ 0 and ∀b ∈ B, ∃z ≥ 0 : Gz ≤ b -F y} and cases where strict inclusion holds (leading to an improved robust optimal solution value over the optimal value of Soyster's model) can easily be found, as illustrated by the following example.

In this example we consider 3 variables y ≥ 0, z 1 ≥ 0 and z 2 ≥ 0, and 3 constraints:

y -z 1 ≤ b 1 y -z 2 ≤ b 2 z 1 + z 2 ≤ b 3
and the uncertainty set B is taken as the set containing the two vectors (1, 0, 1) T and (0, 1, 1) T . The objective function is to maximize y. It is easily checked that the set Y S corresponding to Soyster's model is in this case the real interval [0, 1/2] leading to an optimal robust solution value 0.5. On the other hand, the set Y corresponding to our two-stage model is the real interval [0, 1] leading to the (less conservative) optimal robust solution value 1. Indeed, the value y = 1 is feasible in our model because, in case b = (1, 0, 1) T occurs, we can take z 1 = 0 and z 2 = 1; and, in case b = (0, 1, 0) T occurs, we can take z 1 = 1 and z 2 = 0. (Observe, as already pointed out above, that the value taken by the z variables indeed depends on which b ∈ B is actually realized). Of course, this example does not rule out the possibility of having Y = Y S for some special instances. As will be seen in the next section, this possibility will arise in connection with the robust PERT scheduling problem, in the special case (referred to there as Case 1) where the uncertainty set on the task durations is the cartesian product of a family of real intervals.

An application of the 2-stage robust LP model with RHS uncertainty: robust PERT scheduling

In this section we specialize the general two-stage robust LP decision model investigated above to robust PERT scheduling, with an uncertainty set D on the durations of the tasks, supposed to be given as a (finite or infinite) list of 'scenarios'. More precisely, we want to determine an earliest termination date which can be achieved for any realization of the task durations d in a given uncertainty set D.

Formulation as a 2-stage robust LP model

Consider a PERT network represented as a directed circuitless graph N in which the nodes correspond to tasks (the tasks are numbered i = 1, 2, ..n, the set of tasks is denoted I), and there is an arc (i, j) with length (duration) d i whenever there is a precedence constraint stating that processing of task j should not start before completion of task i. The set of arcs is denoted U . We assume that node 1 has no immediate predecessor (it thus represents the initial task) and node n has no direct successor (it thus represents the terminal task of the project). Denoting y j (j = 1, ..n) the starting date for each task j, and assuming first that the task durations d i are exactly known, we want to minimize the total duration of the project while satisfying all precedence constraints, in other words: Maximize -y n s.t:

y 1 = 0 y i -y j ≤ -d i , ∀(i, j) ∈ U
Indeed, it is easy to check that in the above, y 1 = 0 can be replaced by y 1 ≥ 0, or equivalently -y 1 ≤ 0, thus the problem can be rewritten:

Maximize -y n s.t: -y 1 ≤ 0 y i -y j ≤ -d i , ∀(i, j) ∈ U
This model is recognized as a special case of (I), the constraint matrix [F, G] being formed by the transpose of the node-arc incidence matrix of N with an additional row involving variable y 1 only (with associated coefficient -1). F is reduced in this case to a single column (the column corresponding to node n in the transpose of the incidence matrix of N ). The right handside vector b is the vector with coefficients equal to the opposite of the task durations (more specifically, the right handside coefficient for the constraint corresponding to arc (i, j) ∈ U is equal to -d i ). Note that we do not state explicitly the nonnegativity conditions on y, since they are implied by the precedence constraints and nonnegativity of the d i coefficients.

Thus the problem is cast in a form very similar to (I) the only difference being that the nonnegativity conditions on y and z are dropped. The consequence of this on the analysis of § 3.2 is just that we have to consider the polyhedral cone C ′ = {u/u T G = 0 and u ≥ 0} instead of the polyhedral cone: C = {u/u T G ≥ 0 and u ≥ 0.} Due to the special structure of the G matrix arising in the PERT scheduling problem, we have the following result: Proposition 1: The extreme rays of the polyhedral cone C' are in 1-1 correspondence with the characteristic vectors of the various paths between node 1 and node n in N .

Proof: By observing that G T is the node-arc incidence matrix of the graph N without the row associated with node n but with an extra column with coefficient -1 in the first row and all other coefficients 0, it is realized that u satisfying G T u = 0 and u ≥ 0 corresponds to a nonnegative flow between node 1 and node n in N with value equal to u 1,1 , the component of u corresponding to the extra column with coefficient -1 in the first row and all other coefficients 0. Therefore the extreme rays of the cone C ′ correspond to the incidence vectors of the various paths connecting node 1 to node n in N .

Let us denote P = {π 1 , π 2 , ..., π K } the set of all paths between 1 and n in N , u 1 , u 2 , ..., u K the corresponding characteristic vectors, the condition (u k ) T (b -F y) ≥ 0 specializes to:

-i∈π k d i + y n ≥ 0 (this is because, in that case: (u k ) T b = -i∈π k d i and (u k ) T F y = -y n ). So the condition for feasibility is that for each path π k : y n ≥ i∈π k d i .
In view of this, the robust PERT scheduling problem can be reformulated as:

max -y n s.t: y n ≥ i∈π k d i , ∀π k ∈ P, ∀d ∈ D
where we recall that D denotes the uncertainty set for the task durations.

This problem therefore reduces to determining the path π k maximizing, over the set P of all possible paths in N , the objective function:

max d∈D { i∈π k d i }
in other words we want to solve:

max π∈P max d∈D { i∈π d i } (RPS) ('Robust Pert Scheduling' problem)
Now, if we want to go further into the analysis of (RPS), we have to specify how the uncertainty set D is defined. Of course there are many possible ways for this; we content ourselves below to examine two among the most natural possible definitions, and show that, for each of them, the above robust optimization problem (RPS) can be efficiently solved .

Case 1: D is a scaled ball w.r.t. the L ∞ norm.

The first easy special case is when, for each task i, the duration d i can take any value in a given real interval

[d - i , d + i ] with 0 ≤ d - i ≤ d + i . In this case D is the cartesian product: [d - 1 , d + 1 ] × [d - 2 , d + 2 ] × ... × [d - n , d + n ],
which may be viewed as a scaled ball w.r.t. the L ∞ norm (using component-wise scaling to have all intervals of equal width). It is easily seen that an optimal robust solution for problem (RPS) can be obtained in this case by looking at a longest path (critical path) in N when each of the tasks i is assigned the longest possible duration d + i . As an illustration of the above, consider the following example with n = 7 tasks, where the graph of precedence constraints has the following arcs: (1, 2), (1, 3), (2, 3), (2, 5), [START_REF] Ben-Tal A | Robust Solution of Linear Programming Problems Contaminated with Uncertain Data[END_REF][START_REF] Mulvey | Robust Optimization of Large Scale Systems[END_REF], [START_REF] Bertsimas | Robust Discrete Optimization and Network Flows[END_REF][START_REF] Bertsimas | The Price of Robustness[END_REF], [START_REF] Bertsimas | Robust Discrete Optimization and Network Flows[END_REF][START_REF] Soyster | Convex Programming with Set-Inclusive Constraints and Applications to Inexact Linear Programming[END_REF], (4, 5), (5, 7) and (6, 7). Thus task 2 cannot be started before completion of task 1, etc. Also note that the tasks are numbered according to a topological ordering of the graph, since there is no arc (i, j) with i > j. The associated intervals Task 7 is not shown in the above table because it is a dummy task (of duration 0, without uncertainty) representing the end of the schedule. It is easy to see that in this example, the optimal solution to the (RPS) problem has duration 34 and corresponds to the critical path (1, 2, 3, 4, 5, 7). Indeed 34 is the earliest achievable termination date if we require that the schedule remains feasible for any possible choice of the task durations in the cartesian product

[d - 1 , d + 1 ] × [d - 2 , d + 2 ] × ... × [d - 6 , d + 6 ]
. This corresponds to the situation where the duration of each task i is d + i . Case 2: D is a scaled Hamming ball of bounded radius Γ.

Here again we assume that, for each task i, the duration d i can take any value in a given real interval [d i , d i + ∆ i ] with d i ≥ 0. d i is called the nominal value of the duration for task i, d i +∆ i being the possible extreme (or worst-case) value for the task duration. As is actually the case in many practical applications, it is unlikely that all tasks simultaneously take on worst-case values. To take this observation into account, we will impose an upper bound Γ on the number of task durations which are allowed to take on a worst-case value, given that all task durations which do not take on a worst-case value are assumed to be equal to their nominal value. More formally, associating with each task i a 0-1 integer variable u i , the uncertainty set D corresponding to this definition is:

D = {θ = (θ) i=1,...,n /θ i = d i + ∆ i u i (i = 1, ..., n) such that: n i=1 u i ≤ Γ, u i ∈ {0, 1}, ∀i}
As can be seen from the above definition, in the special case where all ∆ i are equal to 1, D is recognized as the Hamming ball of radius Γ centered at d = (d i ), i = 1, ..n (in other words, θ -d can be any 0-1 vector with at most Γ components equal to 1). When the ∆ i 's take on arbitrary positive values, the Hamming ball structure is still present after applying scaling to each component i with respect to the corresponding ∆ i value.

We note here that, in spite of the fact that the definition above is close in spirit to the concept of uncertainty suggested by Bertsimas andSim (2003, 2004), our model is fairly different from the one studied by these authors, since they restrict themselves to rowwise uncertainty, whereas in our robust PERT scheduling problem, we have uncertainty on the RHS only (a special case of columnwise uncertainty). For more detailed discussion of this issue, see §4.2 below.

We now show that, with the above definition of the uncertainty set D, problem (RPS) can be efficiently solved via a dynamic programming recursion. To that aim, we will consider the problem with parameter Γ as only one representative of the class of problems (RPS[i, k]) for i running from 1 to n (the number of tasks) and k running from 0 to Γ. More precisely, assuming that the tasks are numbered according to a topological ordering of the circuitless graph N , (RPS[i, k]) consists of the robust PERT scheduling subproblem corresponding to the subset of tasks 1, 2,...,i, the durations of at most k of which are allowed to take on their worst-case values. The case k=0 (no deviation allowed) corresponds to the usual PERT scheduling problem in terms of the nominal values for the task durations. We denote v * [i, k] the optimal objective function value for problem (RPS[i, k]), and for any task i, we denote P red[i] the set of tasks j such that (j, i) is an arc of N (the set of direct predecessors of node i). Bellman's optimality principle then leads to the following dynamic programming recursion:

∀i ∈ [1, n], ∀k = 0, 1, ...Γ : v * [i, k] = max j∈P red[i] max{v * [j, k] + d j ; v * [j, k -1] + d j + ∆ j } (3)
The optimal value of the robust PERT scheduling problem we are interested in is then: v * [n, Γ]. The rationale behind Eq. ( 3) can easily be explained as follows. Consider the set of all paths from 1 to j ∈ P red [i]. The duration of arc (j, i) has nominal value d j and worst-case value d j + ∆ j . The maximum duration of a path from 1 to i through j with at most k tasks allowed to deviate from their nominal values can be obtained: either by allowing for at most k deviations on the subset of tasks {1, 2, ..j} and taking the nominal duration for arc (j, i); or by allowing at most k -1 deviations on the subset of tasks {1, 2, ..j} and taking the worst-case duration for arc (j, i). Thus the optimal value for node i via node j is the maximum value among these two alternatives, and the optimal value for node i is the maximum taken on the set of all direct predecessors of i. Obviously, solving the recursion (3) is achieved in polynomial time O(m×n), where m is the number of arcs and n the number of nodes of the PERT network; more precisely the complexity is O(m × Γ) where Γ, the parameter defining the uncertainty set, is at most n, the number of tasks (but often significantly smaller than n in practical applications).

Let us illustrate the above on the same 7 task example as the one considered to illustrate case 1. We thus consider the same intervals for the task durations, the lower bound of each interval representing the nominal task duration, and the upper bound representing the worst-case duration. For Γ = 3, application of the recursion (3) leads to the v * [i, k] values shown in the following table. The maximum above is obtained for j = 6 and the corresponding optimal path is (1, 2, 6, 7).

i = 2 i = 3 i = 4 i = 5 i = 6 i = 7 k = 0 2
Similarly we have:

v * [7, 3] = max{v * [3, 3] + 3; v * [3, 2] + 6; v * [5, 3] + 4; v * [5, 2] + 8;
v * [6, 3] + 8; v * [6, 2] + 16} = max{15, 18, 28, 29, 20, 28} = 29.

The maximum above is obtained for j = 5 and the corresponding optimal path is (1, 2, 3, 4, 5, 7). It is thus seen that, depending on the choice of the control parameter Γ, various optimal paths are obtained which, of course, may differ from the optimal solution to the non-robust PERT scheduling problem (considering only the nominal values for the task durations). Also, observe that the value v * [7, 3] = 29 corresponds to a less conservative robust situation as compared with the one obtained in case 1 above.

Differences with Bertsimas and Sim's approach

We now turn to show that, in spite of the similarity in the definition of the uncertainty sets, the robust version of the PERT scheduling problem investigated here is essentially different from the model proposed by Bertsimas and Sim [START_REF] Bertsimas | Robust Discrete Optimization and Network Flows[END_REF] for the robust version of the shortest path problem. From an abstract point-of-view, the difference basically stems from the fact that, in our case, we are faced with a LP problem with uncertainty on the right handside, whereas Bertsimas and Sim address a LP problem with uncertainty on the cost coefficients. However, to further understand the source of this difference, we show below which difficulties would arise if we wanted to apply the Bertsimas-Sim approach to the robust longest (critical) path problem on a directed circuitless graph G.

Following these authors, the robust shortest s-t path problem in G with uncertainty parameter Γ (assuming Γ ∈ N) is formulated as:

(RSP) min x∈X { (i,j)∈U c i,j x i,j + max S⊆U,|S|≤Γ (i,j)∈S ∆ i,j x i,j }
where X denotes the set of incidence vectors of all s-t paths in G; c i,j denotes the nominal cost of arc (i, j) and c i,j + ∆ i,j is the worst-case cost of arc (i, j). After transformation of (RSP) using the duality theorem to convert the second term in the brackets into a minimization, the problem is reformulated as a standard LP , the solution of which reduces to m + 1 aplications of a standard shortest path algorithm. Observe that one of the reasons for all the above to work so nicely is that the second term in the brackets, as a function of x, is convex in x , since it is the pointwise maximum of a finite number of linear functions.

The above approach is still valid if, instead of looking for an optimum robust minimum cost path, we were looking for an optimum robust maximum benefit path: c i,j > 0 being interpreted as a reward associated with the use of arc (i, j), the effect of uncertainty being to reduce the nominal reward c i,j by the amount i,j . The problem would then take the form: max x∈X { (i,j)∈U c i,j x i,j -max S⊆U,|S|≤Γ (i,j)∈S ∆ i,j x i,j } which is essentially analogous to the above robust minimum cost path, up to a change in the signs of the coefficients in the objective (still assuming, of course, that the graph G under consideration is circuitless). In particular, we note that the function to be maximized is concave, so we still have a convex optimization problem.

By contrast, the robust PERT scheduling problem addressed in the present paper is formulated as: max x∈X { (i,j)∈U c i,j x i,j + max S⊆U,|S|≤Γ (i,j)∈S ∆ i,j x i,j } (with c i,j > 0 and ∆ i,j > 0).

It is then readily observed that this problem consists in maximizing a convex function of x on {0, 1} m , and it is well-known that this cannot be simply reduced to ordinary linear programming as is the case for Bertsimas and Sim's approach. Thus robust PERT scheduling may be viewed a typical illustration of the big differences between models featuring rowwise uncertainty and models featuring columnwise uncertainty in robust Linear Programming.

Conclusions

In this paper, various Robust Linear Programming problems have been investigated, and the question of whether LP duality can still be used to help in solving such problems has been addressed and answered negatively. Among the problems considered, Robust Linear Programming problems with uncertainty in the right handsides only, have been recognized as an interesting sub-class of problems, for which the solution techniques should not confine themselves to the classical approach proposed by [START_REF] Soyster | Inexact Linear Programming with Generalized Resource Sets[END_REF]. In this respect we have been lead to propose a new class of robust LP models referred to here as 'Two-Stage Robust Decision Models' which can be expected to lead to less 'conservative' optimal robust solutions than those usually obtained from Soyster's model. In order to show the practical usefulness of this 2-Stage model, a specialization robust PERT scheduling has been discussed, leading , under two natural ways of defining the uncertainty set w.r.t. the task durations, to efficient solution methods. Also some fundamental difference between our approach to robust PERT scheduling and the one proposed by Bertsimas and Sim [START_REF] Bertsimas | Robust Discrete Optimization and Network Flows[END_REF] in the context of the robust shortest path problem has been pointed out. We think that many other possible applications of this 2-Stage robust modeling approach would deserve further investigations, for instance in dynamic inventory management, optimal resource allocation problems, telecommunication problems, etc. This will be the subject of future research.

  [d - i , d + i ]for the durations of the tasks are the following:

Acknowledgments

We thank the referees for their careful reading and constructive comments.