Sigal Portnoy

Nicolas Vuillerme

Yohan Payan

Amit Gefen
email: gefen@eng.tau.ac.il

Prof Amit Gefen

Clinically-Oriented Real-Time Monitoring of the Individual's Risk for Deep Tissue Injury

Keywords: Pressure ulcer, Rehabilitation, Sitting, Wheelchair, Biomechanical model

Spinal cord injury patients are under daily risk for developing deep tissue injurya severe pressure ulcer that initiates in soft tissues at the bones' proximity. We aimed to formulate a patient-specific biomechanical model that can continuously monitor internal tissue stresses in real-time. We adopted a formulation solving an axisymmetric contact problem of a finite-thickness, elastic layer (soft tissue) and a rigid spherical indentor (ischial tuberosity). We utilized finite element analyses to expand the formulation for large deformations. Sensitivity analyses showed that the soft tissue mechanical properties are the most influential factors in this modeling. We then used synthetic surface pressure data and actual surface pressures recorded under the buttocks of 5 paraplegic wheelchair-users to demonstrate clinical feasibility.

Output parameters were designed to be simple so that they can be easily interpreted by the user. Specifically, we calculated peak and average internal von Mises stress and stress dose, under each buttock, and also a time-dependent stress asymmetry index, to account for frequency of posture adjustments. Inter-subject variability was higher than the intra-subject variability. The heaviest subject had the highest maximal and average peak internal soft tissue stress. We believe that this method holds a high potential for clinical applications.

INTRODUCTION

Spinal cord injury (SCI) patients are under daily risk for developing deep tissue injury (DTI), which is a severe pressure ulcer that typically initiates in soft tissues at the proximity of the ischial tuberosity (IT) bones [START_REF] Black | Deep tissue injury: an evolving science[END_REF][START_REF] Ceelen | Compression-induced damage and internal tissue strains are related[END_REF][START_REF] Gefen | Deep tissue injury from a bioengineering point of view[END_REF][START_REF] Olesen | Missing links in pressure ulcer research -an interdisciplinary overview[END_REF][START_REF] Oomens | The importance of internal strain as opposed to interface pressure in the prevention of pressure related deep tissue injury[END_REF][START_REF] Stekelenburg | Deep tissue injury: how deep is our understanding?[END_REF] . DTI has recently become a consensus pathway of injury, which is reflected in the International Pressure Ulcer Guidelines published jointly by the European and American Pressure Ulcer Advisory Panels (see www.epuap.org and www.npuap.org, respectively). While surface pressure measurements are believed to be effective in alerting wheelchairusers against focal pressures that may cause skin injury [START_REF] Pipkin | Effect of model design, cushion construction, and interface pressure mats on interface pressure and immersion[END_REF], these measurements cannot predict dangerous internal tissue loading [START_REF] Agam | Toward real-time detection of deep tissue injury risk in wheelchair users using Hertz contact theory[END_REF][START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. For example, a similar pressure map may register under the buttocks of a heavy person with sharp ITs and a thin person with blunt ITs; however, their susceptibility to DTI depends on the IT curvature as well as the thickness of the soft tissues, as demonstrated by MRI-based finite element (FE) analyses conducted by Elsner and Gefen [START_REF] Elsner | Is obesity a risk factor for deep tissue injury in patients spinal cord injury?[END_REF] and Sopher et al. [START_REF] Sopher | Exposure to internal muscle tissue loads under the ischial tuberosities during sitting is elevated at abnormally high or low body mass indices[END_REF].

Patient-specific FE analyses demand relatively high computer power, at least for large tissue deformations, which is the case for the buttocks area, and the model creation process cannot be easily automated [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: application to 60 clinical cases[END_REF]. Furthermore, monitoring sitting pressures through pressure mats is mostly reserved to rehabilitation institutes due to high costs.

Therefore, internal loads are presently unmonitored -neither in the clinical setting nor in the patient's daily surroundings.

Previous studies [START_REF] Agam | Toward real-time detection of deep tissue injury risk in wheelchair users using Hertz contact theory[END_REF][START_REF] Linder-Ganz | How do normals move during prolonged wheelchair-sitting?[END_REF] showed that real-time patient-specific monitoring of internal soft tissue strains and stresses under the ITs is, in principle, feasible.

However, the Agam and Gefen [START_REF] Agam | Toward real-time detection of deep tissue injury risk in wheelchair users using Hertz contact theory[END_REF] model did not directly account for the individual's soft tissue thickness under the ITs. The Linder-Ganz et al. [START_REF] Linder-Ganz | How do normals move during prolonged wheelchair-sitting?[END_REF] model was more anatomically accurate but required the computer power of a laptop to run. Here we sought to develop a biomechanical model that is simple enough for running on low computer power, i.e. can be programmed on low-cost chips, in real-time, but that still holds essential individual anatomical characteristics. We further aimed to establish "user-friendly" output parameters that can be easily interpreted by users. The biomechanical model is presented herein and its clinical feasibility is demonstrated.

The present study applies to SCI patients, who are considered more prone to DTI with respect to other patient populations [START_REF] Black | Deep tissue injury: an evolving science[END_REF][START_REF] Ceelen | Compression-induced damage and internal tissue strains are related[END_REF][START_REF] Gefen | Deep tissue injury from a bioengineering point of view[END_REF][START_REF] Olesen | Missing links in pressure ulcer research -an interdisciplinary overview[END_REF][START_REF] Oomens | The importance of internal strain as opposed to interface pressure in the prevention of pressure related deep tissue injury[END_REF][START_REF] Stekelenburg | Deep tissue injury: how deep is our understanding?[END_REF]. Note that we do not aim at demonstrating prevention capabilities of the method at this stage, but rather, we present the algorithm and provide a proof-of-concept. Demonstrating prevention capacities will require further clinical studies that are underway.

METHODS

Biomechanical Model

In order to evaluate internal stresses in sub-dermal soft tissues at the proximity of the IT in real-time, using low computing power and low hardware costs, we adopted a formulation solving an axisymmetric contact problem of a finite-thickness, transversely isotropic elastic layer and a rigid spherical indentor [START_REF] Ning | Asymptotic solutions for axisymmetric contact of a thin, transversely isotropic elastic layer[END_REF]. The IT bone edge is modeled as a hemisphere with radius R (Fig. 1a). Due to the bodyweight, the IT hemisphere is indenting, to a depth of 0 z , against a soft tissue layer that effectively represents the mechanical contributions of the overlying muscle, fat and skin. A contact radius a is formed between the IT and soft tissue layer. The downward displacement of the IT is given by

2 2 0 a R R w z - - + =
(Fig. 1a). Assuming small deformations of the soft tissue at this stage of the modeling, R a w 2 2 ≈ [START_REF] Ning | Asymptotic solutions for axisymmetric contact of a thin, transversely isotropic elastic layer[END_REF]. The thickness of the soft tissue is h. The load applied to the IT by the bodyweight at time t is P(t), which induces an equal and measurable reaction force at the surface of the body (Fig. 1a). The contact between the IT and the soft tissue can be considered as either free-slip, i.e. the bone can slide freely on the overlying soft tissue or as no-slip, meaning that there is no movement between the bone and the soft tissue (see further below).

The formulation of the stress tensor in the soft tissue for the IT-soft tissue contact problem depicted in Fig. 1a is as follows. For a no-slip IT-soft tissue interface [START_REF] Ning | Asymptotic solutions for axisymmetric contact of a thin, transversely isotropic elastic layer[END_REF]:

( ) [ ] a r z h h R r c h R r a c h R r a c rz zz rr , 0 2 2 44 2 2 33 2 2 13 ∈ ∀          - - = - - = - - = σ σ σ where ( ) 4 1 33 4         = c t hP R a π (1)
where C ij are the constitutive parameters of the soft tissue (given later in eq. 3).

For free-slip IT-soft tissue interface, the stress tensor is [START_REF] Ning | Asymptotic solutions for axisymmetric contact of a thin, transversely isotropic elastic layer[END_REF]: Assuming now an isotropic soft tissue layer, the constitutive parameters of the soft tissue are:

( ) ( ) ( ) [ ] a r r a h R c c c c r a h R
( ) ( )( ) ( )( ) ( ) ν α ν ν να ν ν α ν + = - + = = - + - = = 1 2 2 1 1 2 1 1 1 44 13 12 33 11 E c E c c E c c (3)
where E and ν are the effective elastic modulus and Poisson's ratio of the soft tissue, respectively. A geometry function α(R/h) was introduced herein to expand the formulation for determining soft tissue stresses under large deformations as well. The α(R/h) functions were characterized separately for the no-slip and free-slip IT-soft tissue contact conditions, using large-deformation FE models (described later on) of the geometry depicted in Fig. 1a.

The von Mises stresses in the soft tissue around the IT at each time point t are calculated by:

( ) [ ] 2 6 , , 2 2 2 2 rz rr zz zz rr VM z r t σ σ σ σ σ σ + + + - = (4)
In the real-world situation, slipping of the IT upon the soft tissue is hindered by interlacing fibers between the IT and the surrounding soft tissue, and possibly by some degree of tissue adhesion [START_REF] Shacham | Measurements of the static friction coefficient between bone and muscle tissues[END_REF]. Since the real-world contact between the IT and soft tissue is neither frictionless nor completely bonded, we aimed at calculating the von Mises stress in the soft tissue assuming frictional IT-soft tissue contact. For this purpose, we used a weight-function with β proportional to the IT-soft tissue friction coefficient, as follows:

( ) ( ) ( ) ( ) ( ) ( ) z r t z r t z r t slip free VM slip no VM VM , , 1 , , , , - - - + = σ β σ β σ (5)
Substituting β=1 provides the no-slip solution and substituting β=0 provides the freeslip solution to the soft tissue stresses. An alternative approach was also explored, in which the weighted average of the two stress tensors (eq. 1 and eq. 2) was first calculated, followed by the calculation of the von Mises stress, according to eq.4.

Differences between this approach and the algorithm expressed by eq. 5 did not exceed 7.5% in internal stress values, hence we proceeded with the eq. ( 5) option, consistent with the separate free-slip and no-slip FE analyses conducted to determine α(R/h).

Focusing on spatial peaks of soft tissue stresses, which appear at the IT-soft tissue contact [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF], we limit eqs. ( 1) and (2) to the symmetry axis ( 0 = r ; Fig. 1a) where soft tissue deformations maximize, and do not depend on the z coordinate. The geometry functions ( )

slip no h R / α and ( ) slip free h R /
α mentioned above were determined by comparing von Mises stresses computed using eq. ( 5) for different h R / values, with peak von Mises stresses computed with large-deformation FE simulations of corresponding geometries (Fig. 1). Specifically, each axisymmetric FE model (ABAQUS v. 6.9, SIMULIA, RI, USA) consisted of a rigid IT with radius R, and a soft tissue layer with thickness h, effective elastic modulus E=25.3kPa (as measured for transversely indented flaccid muscle tissue by Palevski and colleagues) [START_REF] Palevski | Stress relaxation of porcine gluteus muscle subjected to sudden transverse deformation as related to pressure sore modeling[END_REF] and Poisson's ratio ν=0.49 [10]. The soft tissue layer was modeled with 8-node biquadratic axisymmetric quadrilateral elements ("CGAX8R" in ABAQUS) with higher mesh density at the contact region (Fig. 1). We selected the mesh densities after gradually condensing the meshes until the resulted stresses did not change by more than 2%. We analyzed internal soft tissue stresses where R/h ranged between 0.5 to 2.3 (R was set to 20mm and h values ranged between 8.7mm to 40mm), which represents a broad range of anatomical and pathoanatomical conditions characteristic to paraplegia, as previously measured by means of MRI [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. We collected data from a total of 6 FE model geometries (see R/h values in Fig. 1d,e), applied with downward displacement so that the ratios of IT indentation depth over the undeformed soft tissue thickness were 20% to 50% at 6%-steps. Each model was simulated with either no-slip or free-slip IT-soft tissue contact (Fig. 1). From the aforementioned FE simulations we obtained that α(R/h) behaves as a power law function (Fig. 1d,e):

B h R A       = α (6)
where the constants A are 0.9718 and 0.0686 for the free-slip and no-slip contact conditions, respectively, and the constants B are -0.277 and -0.285 for the free-slip and no-slip contact conditions, respectively (Fig. 1). This function was chosen since it provided better fit than polynomial functions and had just two parameters. Note that for the anatomically-relevant R/h range, and in a large deformation regime, the α(R/h) functions mostly act to lower ( )

t VM σ
. The criterion for determining the A and B constants was a minimal root mean square error (RMSE) between the von Mises stresses predicted from eq. ( 5) and the peak von Mises stresses computed by FE simulations with corresponding h R / geometry:

∑         - = 2 , , 1 FE VM FE VM VM N RMSE σ σ σ (7)
After introducing the α(R/h) functions in eq. ( 3), the mean RMSE were 5.4% for the no-slip and 3% for the free-slip IT-soft tissue contact conditions, respectively, for the entire range of R/h. The maximal RMSE were 7.2% (at R/h=0.8) for the no-slip and 5.5% (at R/h=2.3) for the free-slip IT-soft tissue contact conditions, respectively.

The correction factor has been calculated for R/h in the undeformed configuration, so that it decreases with a decrease in the applied load (i.e. when loads are small the correction factor tends to 1). Nevertheless, real-world conditions involve IT indentation depths over undeformed soft tissue thicknesses in the range of 20-50% (as applied herein) [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF], and hence, in practice, all calculations involved large loads/deformations, and R/h values that were higher than zero (at least 0.5).

In order to evaluate the history of soft tissue exposure to stresses over time, an internal stress dose (ISD) parameter was computed -as the numerical time integral of the peak internal von Mises stress in the soft tissue. At the occurrence of a load-relief, the ISD is reduced proportionally to the duration of the relief. Specifically, the ISD is multiplied by a factor of 1/t (time in seconds), so that if a person lifts his buttock for 1s, there is no reduction in accumulated ISD, if he or she lifts it for 2s, the stress dose reduces to half and so on.

Sensitivity Analyses

Subsequently, our analytical model was used to perform sensitivity analyses of intrinsic patient factors, i.e. anatomy and tissue mechanical properties. The elastic modulus of the soft tissue varied between 13.7kPa and 40kPa according to experimental data from loading flaccid muscle tissue perpendicularly to its fibers [START_REF] Palevski | Stress relaxation of porcine gluteus muscle subjected to sudden transverse deformation as related to pressure sore modeling[END_REF].

Poisson's ratios were varied between 0.49 and 0.499, as soft tissues of the buttocks contain ~80% water [10]. The geometry of the IT radii (12mm -25mm) and soft tissue thicknesses (8mm -40mm) were adopted from the MRI studies of Linder-Ganz and colleagues [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. The parameter β, which is proportional to the IT-soft tissue friction coefficient, was varied through its entire range (0 to 1; eq. 5). The results of the sensitivity analyses, presented in Table S1 (included with the supplementary material), revealed that the mechanical properties of the soft tissue (E, ν) are the most influential factors in this modeling, followed by the internal anatomy (R, h).

Accordingly, the impact of variation in soft tissue mechanical properties on the peak internal von Mises stress and ISD was tested via simulation of static and dynamic sitting conditions (see Section 2.3 below).

Static and Dynamic Synthetic Data

We simulated two possible wheelchair sitting conditions: static sitting and a dynamic loading during wheelchair propulsion. To simulate static sitting we created synthetic data consisting of 2 minutes of bilateral weight bearing, followed by load shifting to the monitored side for two minutes, then 5 seconds lift-off (zero load on the monitored buttock) and back to 2 minutes of bilateral weight bearing (Fig. 2).

Simulations of dynamic sitting were created according to force patterns measured under the buttocks of wheelchair users during wheelchair propulsion [START_REF] Dabnichki | Pressure variation under the ischial tuberosity during a push cycle[END_REF] (Fig. 3). In each simulation we varied either the elastic modulus E the Poisson's ratio ν for the ranges defined in Section 2.2. In each simulation we calculated the internal peak von Mises stress and the ISD as function of time (Figs. 2,3).

Determining the Optimal Pressure Sensing Area under the Ischial Tuberosity

At this stage of the formulation, all parameters of our model are well-defined except the area on the pressure mat to be used for calculating body loads transferred through the IT on each side, P(t). We therefore compared previously-published subject-specific MRI-based FE buttocks stress results [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF] to our present stress calculations, made for different surface area sizes and for the same subjects, in order to test which surface pressure area produces stress values that best fit the MRI-FE calculations by Linder-Ganz and colleagues [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. Specifically, we compared values of peak soft tissue von Mises stresses calculated by our present algorithm against patient-specific FE internal soft tissue von Mises stresses published by Linder-Ganz et al. [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF] for the same healthy (n = 6) and paraplegic (n = 5) subjects during static sitting. For our calculations we assumed that the surface load under each IT, was dispersed across a square of the buttock-sit contact area around the projection of the IT in the pressure map. The P(t) force was hence calculated from pressure mat measurements acquired from each subject (''Tactilus'', SensorProduct Co., NJ, USA), by multiplying the measured pressure [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF] and different sizes of the aforementioned square. We assumed that E muscle =25.3kPa [START_REF] Palevski | Stress relaxation of porcine gluteus muscle subjected to sudden transverse deformation as related to pressure sore modeling[END_REF], E fat =0.852kPa [START_REF] Gefen | Viscoelastic properties of ovine adipose tissue covering the gluteus muscles[END_REF], ν=0.49 [10] and β=0.2 [START_REF] Shacham | Measurements of the static friction coefficient between bone and muscle tissues[END_REF] for all patients. The effective elastic modulus assigned to the soft tissue bulk of each individual subject was weighed using their individual MRI-measured [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF] muscle-to-fat ratios (averaged for the left and right body sides; Table S2). A 5x5cm 2 area size provided the best outcome in terms of an X-Y plot (Fig. 4a) and an average-difference (Bland-Altman) plot (Fig. 4b), and it is also physiologically reasonable -as a projection of the IT with a sufficiently large soft tissue region under it. Stresses calculated using the present model were compared with effective values of von Mises stress in muscle and fat tissues in Linder-Ganz et al. [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF], which were correspondingly weighed together according to the subject-specific muscle-to-fat ratios. The average-difference plot comparing the present soft tissue stress data with the previous subject-specific MRI-based FE stress analyses of the same 11 subjects published by Linder-Ganz et al. [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF], as shown in Fig. 4b, overall demonstrates good agreement between the two methods, as all data points in the plot are within the ±2 standard deviation (SD) boundaries of the difference (Fig. 4b), and the absolute mean difference between the methods is ~4kPa.

Patient Studies

The present algorithm was tested with sitting pressure data recorded from paraplegic patients in order to demonstrate its clinical feasibility. Specifically, pressure measurements between the wheelchair seat and the buttocks of five paraplegic subjects were previously recorded at a frequency of 1Hz with a 256-sensors pressure mat (''Tactilus'') by Linder-Ganz et al. [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. An ethical committee approval was obtained for these studies and all patients gave their informed consent.

We have selected a convenience sample, which reflects the midrange of patient body weights; the body characteristics and relevant personal information of the subjects are provided in Table S2. For the purpose of assessing feasibility of clinical application, we assumed, at this stage, that E muscle =25.3kPa [START_REF] Palevski | Stress relaxation of porcine gluteus muscle subjected to sudden transverse deformation as related to pressure sore modeling[END_REF], E fat =0.852kPa [START_REF] Gefen | Viscoelastic properties of ovine adipose tissue covering the gluteus muscles[END_REF], ν=0.49 [10] and β=0.2 [START_REF] Shacham | Measurements of the static friction coefficient between bone and muscle tissues[END_REF] for all patients. The elastic modulus assigned to the soft tissue bulk of each subject was weighted using their individual MRI-measured [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF] muscle-to-fat ratios (Table S2). Our algorithm for calculation of internal peak von Mises stresses was employed separately for the left and right buttock of each subject over pressure data recorded for a duration of 2 minutes of static sitting. A code was written (LabView 8.0, National Instruments Co., TX, USA) which finds the peak sitting pressure from one sensor (sensing area of 2.5x2.5cm 2 ''Tactilus'') under the left or right buttock in order to identify the location of the projection of the IT at each side.

The code then chooses 4 sensors at this location that together form a square of 5x5cm which includes the peak pressure, and which also have the highest sum of pressure measurements. This code was verified by visually inspecting its selections of IT projection locations against the measured pressure maps. We computed the maximal internal peak von Mises stress in the soft tissue at each body side for each patient, as well as the averages of von Mises stresses over time. We then calculated the coefficient of variation (COV) for both intra-subject and inter-subject stress data.

A time-dependent soft tissue stress asymmetry index (SAI) was further defined, based on an analogue sitting pressure asymmetry index introduced in Linder-Ganz et al. [START_REF] Linder-Ganz | How do normals move during prolonged wheelchair-sitting?[END_REF]:

( )         = IT right VM IT left VM t SAI σ σ log (8)
This SAI is a measure of the extent of asymmetry in internal soft tissue load sharing between the two body sides. The SAI produces zero for ideally balanced sitting, is positive when the peak internal soft tissue stress is higher under the left IT and negative when the peak stress is higher under the right IT. We calculated the peak and average SAI per each patient and the COV for both intra-subject and inter-subject data. We further calculated the time integral of the SAI and the median SAI frequency as a measure of how often subjects change posture.

RESULTS

The synthetic data simulations (Figs. 2,3) demonstrated that consistent with Table S1, variations in ν overall had mild influence on the evaluated soft tissue stress and ISD, provided that ν was kept within the range of 0.49-0.496, which well the physiologically-relevant range of Poisson's ratios for soft tissues [10]. Hence E is the more important patient-specific mechanical property that needs to be considered in clinical application of the present method.

The maximal and average internal von Mises stresses calculated separately under the left and right IT of 5 paraplegic subjects (whose characteristics are provided in Table S2) are presented in Table 3. An example of the loads, peak internal von Mises stress and ISD under the left and right IT of subject D3 (Table S2) is plotted in Fig. 5. In a two-minutes period of static sitting, the soft tissues under the left IT of this particular subject consistently bore higher stresses and ISD than the soft tissues under his right IT, despite that at some periods, e.g. during the timeframe between 40s and 65s, the at the left and right sides were approximately equal (Fig. 5). This indicates the importance of weighing the effect of the internal anatomy as well, which indeed demonstrates a slightly sharper IT and a thinner soft tissue layer at the left buttock side of this subject (Table S2), which caused the higher internal stresses at his left buttock side.

The averages of maximal peak internal von Mises stresses over time, as well as corresponding COV are calculated for both intra-subject and inter-subject data (Table 1). The heaviest subject, W5 (Table S2), had the highest internal von Mises stresses under his left IT. High values of COV of soft tissue stress were calculated for subject K4, who had the thickest soft tissue bulk (Table S2). The inter-subject variability, as determined by 53.4% COV of the average of peak von Mises stresses across subjects, was higher than the intra-subject variability, for which the COV maximal value was 37.9% and which averaged below 12% (Table 1).

The SAI parameter (eq. 8) is plotted in Fig. 6 for all the paraplegic subjects (Table S2). Peak and average of peak internal SAI, COV of the SAI and median SAI frequency are presented in Table 2. Average SAI and corresponding COV are calculated for both intra-subject and inter-subject data (Table 2). The peak, average and time integral of SAI values calculated for subject K4 were at least 2-times higher than SAI values of the subjects (Table 2, Fig. 6). Consistent with that, the median SAI frequency calculated for subject K4 was the maximal among the five subjects, indicating that subject K4 often shifted her weight from side to side.

DISCUSSION

In this study, we aimed at designing a patient-specific modeling method which is clinically-oriented and that will monitor in real-time and alert wheelchair users against high deep tissue stresses that may cause DTI. We employed clinically-oriented intuitive measures of the risk for DTI in an individual patient, being the current peak internal tissue stress and the stress history (ISD), as well as the balance between stresses in the right and left buttock sides (SAI) and its frequency of change and history. This entire set of parameters is based on the von Mises stress, which was selected because it is a scalar that weighs the contributions of the different components of the stress tensor (compression, tension and shear) together. In view of the target user communities -clinicians and perhaps the patients themselves -this is practical and efficient in terms of data reduction. The SAI in particular is a dimensionless number that indicates how often patients shift weight between sides.

Weight shifting may temporarily relieve high loads under one IT, yet cause extreme loading on the contralateral IT. Asymmetry is a very important parameter used by to evaluate postures and load-bearing activities. Physicians, physical therapists and care-takers are accustomed to terminology concerning with asymmetry and hence it is a simple measure that provides data which clinicians can easily interpret. Quantification of internal stresses, as well as the weight shifting frequency, provide clinicians with simple quantitative data concerning the patient's sitting habits (posture, load-asymmetry and load shifting frequency). Following continuous monitoring with the proposed system, patients can learn to regulate their internal soft tissue loading and adapt to safer sitting habits in term of DTI prevention. In other words, the system is useful for clinical evaluation as well as biofeedback applications.

It should be noted that the concept of the SAI is considerably different than simply tracking the regulation of sitting pressures between the two body sides, because an inherent anatomical left-right asymmetry (demonstrated in most of the R,h pairs of parameters in Table S2) imposes a need to apply uneven sitting pressures between the two buttock sides, in order to compensate for the anatomical asymmetry.

Previous models of our group, designed to predict patient-specific internal soft tissue stresses under the IT in real-time, were limited by several assumptions, i.e. modeling the soft tissue bulk as an infinite half-space (assumed in Herz's theory in Agam and Gefen [START_REF] Agam | Toward real-time detection of deep tissue injury risk in wheelchair users using Hertz contact theory[END_REF]) or utilizing linear FE with no-slip between all tissues and requiring laptop computing power [START_REF] Linder-Ganz | Real-time finite element monitoring of sub-dermal tissue stresses in individuals with spinal cord injury: towards prevention of pressure ulcers[END_REF]. The modeling designed herein better approximates the actual anatomy of the buttocks by considering finite soft tissue thickness and frictional IT-soft tissue interaction, but still we were able to keep a closed-form analytical solution to the stress state, which is essential for real-time, low-computer power applications. In fact, the relative simplicity of the present modeling allows it to be programmed in a chip that can be embedded in a wheelchair.

There remain, of course, several limitations, mainly concerning the geometry of the IT and soft tissues, which are not an ideal sphere and flat layer. The soft tissues are modeled as a homogenous isotropic material with no differentiation between muscle, fat, skin and connective tissues. Additionally, we assigned mechanical properties taken from animal studies of transversal loading of ovine [START_REF] Gefen | Viscoelastic properties of ovine adipose tissue covering the gluteus muscles[END_REF] and porcine [START_REF] Palevski | Stress relaxation of porcine gluteus muscle subjected to sudden transverse deformation as related to pressure sore modeling[END_REF] fat and muscle tissues, respectively. We considered the muscle-to-fat ratio in the buttocks of subject and weighed the elastic modulus according to the patient's soft tissue anatomy to account, in some degree, for the inter-subject variability in soft tissue stiffness. We do foresee, however, that in clinical use of the present method, it will be necessary to feed patient models with individual, directly-measured tissue mechanical properties in addition to the individual anatomy. One very promising way for obtaining the individual E noninvasively and without ionizing radiation or expensive MRI is to use ultrasound elastography [START_REF] Deprez | 3D ultrasound elastography for early detection of lesions. Evaluation on a pressure ulcer mimicking phantom[END_REF] which examines internal anatomy and tissue stiffness in a single scan. The model could be automated by developing a dedicated image processing code for segmentation of bone contours, muscle and fat layer contours, and calculating a mean radius of curvature of bone and muscle/fat ratio from the segmentation. The spatial resolution of a commercial ultrasound is in the 0.1-1mm which should be sufficient for the present application.

The present algorithm was tested with actual sitting pressure data recorded at the buttocks-wheelchair interface of paraplegic wheelchair users to demonstrate clinical feasibility. It should be noted that a small sample size test group was tested here.

Once the system is put into large scale clinical trials it would be possible to collect population data on exposure of deep tissues to internal stresses, including time history of stresses and stress asymmetry. Through follow-up studies, this will ultimately allow to define threshold values indicating danger of DTI. In an example depicted in Fig. 5, loads measured with the pressure mat for patient D3 (Table S2) were similar the left and right IT during the 40-65s timeframe. Contrarily, internal von Mises soft tissue stresses at that time were evaluated to be higher under the left IT.

Undoubtedly, the higher tissue stresses within the left buttock in this particular subject were produced by a sharper left IT indenting a thinner soft tissue layer at that side (Table S2). In common clinical practice, the equally distributed surface pressures would have been regarded as optimal, although internal tissue stresses are not evenly distributed. This example clearly demonstrates the importance of using patientspecific models that calculate internal soft tissue stresses for DTI prevention.

The heaviest subject, W5 (Table S2), had the highest maximal and average peak internal soft tissue stresses (Table 1). This is consistent with findings of Elsner and Gefen [START_REF] Elsner | Is obesity a risk factor for deep tissue injury in patients spinal cord injury?[END_REF] and Sopher et al. [START_REF] Sopher | Exposure to internal muscle tissue loads under the ischial tuberosities during sitting is elevated at abnormally high or low body mass indices[END_REF] who showed, using FE, that abnormal body mass indices, either overweight or underweight, contribute to the state of elevated internal tissue loads, which consequently increases the risk for sitting-acquired-DTI.

In conclusion, we propose herein a practical algorithm for real-time patientspecific monitoring of internal tissue stresses under the IT of wheelchair users, in clinical or home settings. We believe that this method holds a high potential for clinical applications, bridging between sophisticated modeling efforts at the research and everyday clinical needs. As a next step, we hope to translate it to a medical device ready for large-scale human studies, with the aims of screening patients at risk for DTI, based on high ISD or SAI values, and also, for training patients by means of biofeedback to avoid DTI. S2). Fig. 6 The stress asymmetry index (SAI) defined in eq. 9, calculated for 5 paraplegic subjects (whose characteristics are detailed in Table S2) sitting in a wheelchair for 2 minutes.

Figure captions

Fig. 1 (

 1 Fig.1(a) The geometry of the analytical model. Examples of finite element simulations of a specific geometry (R=20mm, h=15.38mm) with interface conditions of (b) free-slip and (c) no-slip. The α(R/h) functions for (d) free-slip and (e) no-slip simulations and the power law functions fitted to each condition.

Fig. 2

 2 Fig. 2Synthetic simulation of static sitting where (a) the load is applied for 2 minutes of bilateral weight bearing, followed by load shifting to the monitored side for two minutes, then 5 seconds lift-off (zero load on the monitored buttock) and back to 2 minutes of bilateral weight bearing. The computed von Mises stresses for variations in (b) the elastic modulus and (c) the Poisson's ratio are presented as well as the stress doses for variations in (d) the elastic modulus and (e) Poisson's ratio.

Fig. 3

 3 Fig. 3Synthetic simulation of dynamic sitting during wheelchair propulsion where (a) a repetitive load is applied for 3 seconds. The computed von Mises stresses for variations in (b) the elastic modulus and (c) the Poisson's ratio are presented as well as the stress doses for variations (d) the elastic modulus and (e) Poisson's ratio.

Fig. 4

 4 Fig. 4The pressure sensing area taken into account for body load calculations on each side was determined by the best achieved (a) X-Y plot and (b) average-difference (Bland-Altman) plot comparing between von Mises stresses calculated by our algorithm against patient-specific finite element von Mises stresses published by Linder-Ganz et al (2008) for the same healthy (n = 6) and paraplegic (n = 5) subjects during static sitting. The average-difference plot shows good agreement between the two calculation methods as all data points are within ±2 standard deviation boundaries of the difference, and the absolute mean difference is ~4kPa.

Fig. 5

 5 Fig. 5An example of the (a) loads, (b) peak internal von Mises stress and (c) internal stress dose under the left and right ischial tuberosities of D3 (TableS2). Fig.6The stress asymmetry index (SAI) defined in eq. 9, calculated for 5 paraplegic subjects (whose characteristics are detailed in TableS2) sitting in a wheelchair for 2 minutes.

Acknowledgement

Funding for this joint Israel-France scientific collaboration was provided by the Universite Joseph Fourier (Grenoble, France) and the IDS-Arcan company (Montceau-les-Mines, France). The contribution of NV was supported by a grant from the Centre National de La Recherche Scientifique (France).


Abdessamad ; Grirrane

Eleuterio Álvarez

Antonio Pastor

Agustin Galindo

Pastor,* [a] ] Antonio

Agustín Galindo

Thiodipropionate Zn(II) complexes: Synthesis, DFT studies and X-ray structure of [{Zn(phen)(H 2 O)} 2 (µ-tdp) 2 ]•3H 2 O

Keywords: Zinc, Density functional theory, Coordination modes, Thiodipropionate, X-ray Zinc, Density functional theory, Coordination modes, Thiodipropionate, X-ray

published or not. The documents may come     Thiodipropionate Zn(II) complexes: Synthesis, DFT studies and X-ray structure of [Zn(phen)(H2O)2(µ-tdp)2]*3H2O

Introduction

Thiodicarboxylic acids are interesting ligands due to their capacity to coordinate to metal centres through both carboxylate groups and the sulphur atom. Continuing our line of investigation into applications of the thiodiacetate ligand [1], we planned to explore the coordinating properties of the higher homologue of thiodiacetic acid, thiodipropionic acid S(CH 2 CH 2 COOH) 2 , as investigations of structures containing the thiodipropionate ligand (tdp) are comparatively less extensive than those of the thiodiacetate [2]. Additionally, the tiodipropionate anion is more flexible than the thiodiacetate, possibly allowing the formation of novel topologies. Only a few examples of compounds containing tdp [3] or the hydrogenthiodipropionate anion has been X-ray characterised [4], although several structures are known for the parent thiodipropionic acid [START_REF]See for example: (a) V. R. Pedireddi[END_REF].

Results and Discussion

Synthesis and characterisation

Through reaction of aqueous zinc dichloride with a 1:1 mixture of Na 2 CO 3 and thiodipropionic acid colourless crystals of [Zn(tdp)] (1) were obtained. Analytical and IR data were consistent with the proposed formula. Absence of water in the complex is confirmed by the absence of the O-H absorption in the IR spectrum and by analytical data. A very strong and broad absorption band at 1538 cm -1 , corresponding to the antisymmetric vibrations of the carboxylate groups, confirmed the presence of tdp coordinated to the Zn 2+ through the O atoms [6]. The asymmetric unit contains one half of the dimeric complex, the other half being related by an inversion centre. The Zn atoms are each in a distorted octahedral coordination sphere defined by two N atoms of the phen ligand, three O atoms from the two carboxylate anions and the O atom of a water molecule trans to the O atom of the bidentate carboxylate group. Each thiodipropionate anion bridges the two Zn 2+ ions in a different coordination mode. One terminal carboxylate group binds in a bidentate mode through two O atoms, whilst the other is bonded to the metal centre through just one O atom. This bonding mode is similar to that observed in polymeric related complex [Zn(tdp)(bipy)•(H 2 O)] n [3]. The Zn-O(carboxylate) bond lengths fall in the region 2.048-2.272 Å, and the Zn-N bond distances are 2.110(3) and 2.162(3) Å respectively. These distances are comparable to those found in [Zn(tdp)(bipy)•(H 2 O)] n [3], where the Zn-O (tdp) lengths are in the region 2.042-2.214 Å and the Zn-N distances are 2.112(2) and 2.198(2) Å. In the Zn-tdp-bipy compound and in complex 2 the thiodipropionate acts as bridge between the Zn atoms with the same coordination mode µ 2 -κ 2 (O,O')-κ 1 (O''). However, the formation of a binuclear compound in 2 or a polymeric complex in [Zn(tdp)(bipy)•(H 2 O)] n [3] is almost certainly due to the presence of water molecules of solvation in the former and the existence of π-stacking in the latter (see Supporting Information). Other selected structural parameters are collected in Table 1. In the crystal structure, the aqua ligand forms an intermolecular hydrogen bond with the uncoordinated carboxylate oxygen atom (O4 … O5 length 2.618(4) Å). Complex 2 is isostructural to the related adipate Zn compound [7].

Theoretical calculations

Due to the paucity of structural data for thiodipropionate compounds of transition metals [START_REF]Cambridge Structural Database System[END_REF], we decided to explore theoretically the bonding preferences of this ligand using the DFT approach. The first approach was to select the simplest neutral fragment [Zn(tdp)] and to optimize, without symmetry restrictions, the two models I and II, schematically depicted in Scheme 1. Model I displays the thiodipropionate ligand coordinated in κ 3 (O,S,O')-fashion, which is the usual coordination mode for related thiodiacetate mononuclear compounds [1], with both carboxylate groups κ 1 -bonded to the metal. In model II the thiodipropionate group is again tricoordinated, but now the sulphur atom is not bonded to the metal, while the two carboxylate moieties are κ 2 and κ 1 coordinated to zinc. Scheme 1. Models of [Zn(tdp)] theoretically studied. Fig. 2 shows the resulting optimised structures. Coordination of the sulphur atom to the Zn metal, moving from model II to I, is an exergonic process of -3.7 kcal•mol -1 . This low value suggests that such coordination is probably not essential in this type of compounds. If the sulphur atom is not bonded, the coordination of both carboxylates to different metal centres should be favoured over coordination to a single one. Confirmation of this assumption is provided by pimalate compounds (analogous to thiodipropionate with the S atom substituted for CH 2 ). In all structurally characterised pimalate complexes of transition metals each carboxylate group is bonded to a different metal centre [START_REF]Cambridge Structural Database System[END_REF]. To confirm the preferentiality of the coordination mode we extended the DFT calculations to a third model, III (see Scheme 1), in which the thiodipropionate bridges two zinc ions. It is actually a binuclear species [Zn(tdp)] 2 where one carboxylate group is κ 2 -bonded to a Zn atom, with the other carboxylate acting as bridge in a syn,anti-µ 2 -κ 1 (O),κ 1 (O') coordination. Fig. 3 shows the optimised structure. From an energetic point of view, the dimerisation process from II to III is clearly favourable because the conversion is very exergonic (-72.4 kcal•mol -1 ). Consequently, the coordination of tdp as a bridge is evidently favoured over tricoordination to a single metal. This indicates that thiodipropionate most likely behaves as a bridge ligand in compound 1, the compound probably being polymeric in nature [Zn(tdp)] n .

Experimental Section

General

All preparations and other operations were carried out under aerobic conditions. Thiodipropionic acid and other chemicals were obtained from commercial sources and were used without further purification. Infrared spectra were recorded on Perkin-Elmer Model 883 spectrophotometer (KBr pellets). 1 H spectra were run on Bruker AMX-300 spectrometer. Microanalyses (C, H, N) were performed on an Elemental LECO CHNS 93 analyzer by the "Servicios Generales de Investigación" (CITIUS) of the University of Sevilla.

Synthesis of [Zn(tdp)] (1)

To a solution of S(CH 2 CH 2 COOH) 2 (0.53 g, 3 mmol) in water (15 ml) Na 2 CO 3 (0.32 g, 3 mmol) was added and the mixture stirred until evolution of CO 2 ceased. This was then mixed with ZnCl 2 (0.41 g, 3 mmol) dissolved in water (15 ml), the resulting solution stirred for 10 minutes. Acetone was added and the solution was kept at room temperature overnight. Colourless crystals of 1 formed, which were filtered off, washed with acetone and diethyl ether and air dried (0.55 g, 75 % yield). IR (cm

Synthesis of [{Zn(phen)(H 2 O)} 2 (µ µ µ µ-tdp) 2 ]•3H 2 O (2)

A solution of o-phen•H 2 O (0.099 g, 0.5 mmol) in 10 ml of water was added to a suspension of 1 (0.121 g, 0.5 mmol) in 20 ml of water. The mixture was heated at 80 ºC for 1 h, then concentrated to about 20 ml. Slow evaporation of the solvent resulted in the formation of colourless crystals of [{Zn(phen)(H 2 O)} 2 (µ-tdp) 2 ]•3H 2 O (2), which were collected by filtration, washed with water, acetone and diethyl ether and dried in air (0.160 g, 68 %). IR (cm

X-ray structural study of 2

A summary of the crystal data of compound 2 is given in Table 2. The X-ray data were collected at 100 K on a Bruker-Nonius X8APEX-II CCD diffractometer equipped with a graphite monochromator and Mo-K α radiation. The data were reduced [START_REF]SAINT 6.02[END_REF] and corrected for Lorentz polarisation effects and absorption by multiscan method applied by SADABS [START_REF] Sadabs | [END_REF]. The structure was solved by direct methods (SIR-2002 [START_REF] Sir2002 | [END_REF]) and refined against all F 2 data by full-matrix least-squares techniques (SHELXTL-6.12 [12]). All non-hydrogen atoms were refined with anisotropic displacement parameters. The hydrogen atoms were included from calculated positions and refined riding on their respective carbon atoms with isotropic displacement parameters.

Computational details

The electronic structure and geometries of the model complexes I-V were computed within the density functional theory at the B3LYP [13,14] level. The basis set used for all the atoms was 6-311G**. All the optimised geometries were characterised as local energy minima by diagonalisation of the analytically computed Hessian (vibrational frequencies calculations). The DFT calculations were performed using the Gaussian 98 suite of programs [15]. Cartesian coordinates for the optimised molecules are available from the authors upon request.

Supporting Information Crystallographic data for the structure reported in this paper have been deposited with the Cambridge Crystallographic Data Centre as CCDC-661939. Copies of the data can be obtained free of charge on application to CCDC, 12 Union Road, Cambridge CB2 1EZ, UK [fax.: (internat.) + 44 1223/336-033; e-mail: deposit@ccdc.cam.ac.uk].

____________*

  Professor Agustín Galindo, Dr. Antonio Pastor Fax: +34 954 557153 E-mail: galindo@us.es, anpastor@us.es [a] Departamento de Química Inorgánica, Facultad de Química, Universidad de Sevilla, Aptdo 1203, 41071 Sevilla, Spain [b] Instituto de Investigaciones Químicas, CSIC-Universidad de Sevilla, Avda. Américo Vespucio 49, 41092 Sevilla, Spain Supporting information for this article is available on the WWW under http://dx.doi.org/10.1002/zaac.201000xxx or from the author. The reaction of 1 with o-phenanthroline gives the new complex [{Zn(phen)(H 2 O)} 2 (µ-tdp) 2 ]•3H 2 O (2). The crystal structure of 2 was determined by X-ray diffraction (Fig. 1).

Fig. 1 .

 1 Fig. 1. ORTEP view of the dinuclear Zn(II) complex 2. The displacement ellipsoids are drawn at a 50 % probability level.

Table 1 .

 1 Selected bond distances (Å) and angles (º) for compound 2 Zn(1)-O(3) 2.048(3) N(1)-Zn(1)-N(2) 77.86(12) Zn(1)-O(5) 2.063(3) O(3)-Zn(1)-O(2) 97.75(11) Zn(1)-N(1) 2.110(3) O(5)-Zn(1)-O(2) 94.31(11) Zn(1)-N(2) 2.162(3) N(1)-Zn(1)-O(2) 150.05(11) Zn(1)-O(2) 2.187(3) N(2)-Zn(1)-O(2) 89.69(11) Zn(1)-O(1) 2.272(3) O(3)-Zn(1)-O(1) 83.28(11) O(1)-C(1) 1.253(5) O(5)-Zn(1)-O(1) 151.34(11) O(2)-C(1) 1.273(5) N(1)-Zn(1)-O(1) 95.28(11) O(3)-C(6) 1.280(4) N(2)-Zn(1)-O(1) 96.79(11) O(4)-C(6) 1.241(5) O(2)-Zn(1)-O(1) 58.88(10) Zn(1)-C(1) 2.566(4) O(3)-Zn(1)-C(1) 89.20(12) O(3)-Zn(1)-O(5) 91.35(11) O(5)-Zn(1)-C(1) 123.08(12) O(3)-Zn(1)-N(1) 93.43(12) N(1)-Zn(1)-C(1) 123.60(12) O(5)-Zn(1)-N(1) 113.17(12) N(2)-Zn(1)-C(1) 95.11(11) O(3)-Zn(1)-N(2) 171.27(12) O(2)-Zn(1)-C(1) 29.72(11)

Fig. 2 .

 2 Fig. 2. Optimised structures of models I-II.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Optimised structure of model III.

  -1 ): 1538 s, 1451 s, 1413 s, 1327 m, 1316 m, 1211 m, 1167 w, 1063 w, 955 m, 785 s, 742 m, 698 s, 576 w, 512 w. 1 H NMR (300 MHz, D 2 O): δ 2.59 (t, 3 J HH = 7.2 Hz, CH 2 ), 2.30 (t, 3 J HH = 7.2 Hz, CH 2 ). Elemental analysis for C 6 H 8 O 4 SZn (Mw = 241.4): calcd. C, 29.82; H, 3.31; found C, 29.87; H, 3.20 %.

Acknowledgement

This work was supported by Spanish Ministerio de Ciencia e Innovación and Junta de Andalucía.


Yazid ; Mati
email: matie@qu.edu.sa.

Chams Lahlou
email: chams.lahlou@emn.fr

Stephane Dauzere-Peres
email: dauzere-peres@emse.fr

Al-Qassim

Stéphane Dauzère-Pérès

Modelling and Solving a

Keywords: flexible, blocking Job-shop scheduling, genetic algorithm, flexible, blocking, design

HAL is a

1. Introduction Designing a manufacturing system involves complex decisions with conflicting objectives such as cost reduction, high equipment utilization and customer satisfaction. One of the important issues is the selection of the proper equipment type and number. In modern manufacturing systems, machines are more and more flexible, and choosing the adequate type is crucial to handle the variety of possible products, demand fluctuations and machine failures. Further, determining the right equipment number helps reducing the procurement and maintenance costs while increasing the equipment utilization.

This paper discusses the problem of selecting machines in a workshop where key components of aircrafts are produced. With the increase of demand and in order to satisfy customers while minimizing costs, the company under consideration was setting a new workshop that could make different part types on the same set of machines. The main challenge the company had to face was the determination of an adequate configuration of the workshop, in order to minimize the cycle times of products. After a first analysis carried out at the company, focusing on demand forecast, a first configuration of the workshop was proposed. In this configuration, products have to go through workstations, and there is one machine in each workstation. The company was interested in evaluating various choices before building the workshop such as: the number of entrant workstations, the functions of some workstations (i.e. the machine type in these workstations), etc.

In the classical literature on equipment selection (see [START_REF] Chen | A heuristic for solving manufacturing process and equipment selection problems[END_REF], for instance), aggregate models are proposed and solved in which the details of the actual processing of jobs is not considered. In our approach, and because the scheduling constraints are important, we needed to solve a more detailed model that integrates detailed scheduling features to be sure that the design choices were correct. Our problem is thus formulated as a scheduling problem with several possible scenarios for the design of the workshop. The resulting scheduling problem is a complex shop problem with a material handling system for transporting products, resource flexibility, lack of buffer space between machines, etc.

Flexible shop scheduling problems have attracted many researchers who mainly addressed the job-shop case with multi-purpose machines (also called flexible jobshop), in which an operation may be processed on several machines in a given set, and one of them must be selected. Complexity results are given by [START_REF] Brucker | Job-shop scheduling with multi-purpose machines[END_REF], [START_REF] Mati | The complexity of the two-job shop scheduling problems with multi-purpose unrelated machines[END_REF] for the two-job case. For the flexible job-shop with several jobs, various heuristics are proposed. Most of them use a disjunctive graph for representing feasible solutions and a tabu search heuristic to guide the search. Among these heuristics we may cite the hierarchical approach of [START_REF] Paulli | A hierarchical approach for the FMS scheduling problem[END_REF] and the integrated approach of Dauzère-Pérès and [START_REF] Dauzère-Pérès | An integrated approach for modelling and solving the multiprocessor job-shop scheduling problem using tabu search[END_REF]. A heuristic based on genetic algorithms is developed in [START_REF] Chan | Flexible job-shop scheduling problem under resource constraints[END_REF]. This approach is used in an iterative approach to evaluate different flexibility levels of assignment of operations to machines. A generalization to the case where an operation may be processed on several resources simultaneously (multi-resource) is done by [START_REF] Dauzère-Pérès | Multi-resource shop scheduling with resource flexibility[END_REF], Mati and Xie (2008a). In the first paper, a tabu search heuristic is proposed using the disjunctive graph representation. The second paper presents a greedy heuristic that schedules jobs one after another according to a job sequence. A genetic algorithm is used to optimize the job sequence.

The blocking constraint arises in production environment where buffer capacity between workstations is limited. Scheduling problems that take into account this constraint are not abundant and the majority of existing studies concentrates on the flow-shop and hybrid flow-shop scheduling problems. A survey on blocking con- straints is proposed by [START_REF] Hall | A survey of machine scheduling problems with blocking and no-wait in process[END_REF]. For the flow-shop problem, [START_REF] Grabowski | The permutation flow shop problem with blocking: A tabu search approach[END_REF] use the disjunctive graph model for modelling and propose tabu search heuristics for solving the problem. A genetic algorithm is proposed by [START_REF] Caraffa | Minimizing makespan in a blocking flowshop using genetic algorithms[END_REF] for minimizing the makespan. A branch and bound method is proposed by Ronconi (2008) and some complexity results are given by [START_REF] Martinez | Complexity of flow-shop with a new blocking constraint[END_REF] for a new type of blocking constraint. For the hybrid flowshop problem, some heuristic methods are proposed by Wang and Tang (2008), [START_REF] Wardono | A tabu search algorithm for the multi-stage parallel machine problem with limited buffer capacities[END_REF] and exact methods are investigated by [START_REF] Sawik | Scheduling flexible flow lines with no in-process buffers[END_REF][START_REF] Sawik | Mixed integer programming for scheduling flexible flow lines with limited intermediate buffers[END_REF].

Considering blocking constraints in job-shop environments leads to a control problem called deadlock situation. A system deadlock is a situation that arises when the flow of products is permanently inhibited and/or operations on products cannot be performed. For example, in a robotic cell, a deadlock arises if the robot carries a part to a machine which is working on another part. Under such a situation, the robot cannot unload the part since the machine is not free and needs the robot to download the part it is working on.

Various solutions have been proposed to handle deadlock situations in automated manufacturing systems [START_REF] Chu | Deadlock analysis of Petri nets using siphons and mathematical programming[END_REF]: deadlock prevention, detection and avoidance methods. The majority of previous studies that deal with scheduling and deadlock issues uses Petri nets and reachability graphs to obtain optimal solutions [START_REF] Lee | Scheduling flexible manufacturing systems using Petri nets and heuristic search[END_REF]DiCesare 1994, ElMaraghy and[START_REF] Elmaraghy | Deadlock-free Scheduling in Flexible Manufacturing Systems Using Petri Nets[END_REF]. Mati et al. (2001a,b) propose two heuristics for solving the job-shop scheduling problem with multi-resource operations. The first heuristic uses a greedy heuristic and genetic algorithms. The second heuristic makes use of the disjunctive graph and tabu search. Integer linear programming models are proposed in [START_REF] Gomes | Optimal scheduling REFERENCES for flexible job shop operation[END_REF] to schedule flexible job shops under limited intermediate buffers. These models are validated on realistic examples using commercial mathematical programming software. Various characteristics make our problem more complex. A deadlock-free reactive scheduling approach is proposed in [START_REF] Fahmy | A generic deadlock-free reactive scheduling approach[END_REF], which aims at answering in real time to the occurrence of a variety of possible disruptions in flexible job shops.

This paper formulates a practical equipment selection problem as a complex scheduling problem and presents a heuristic approach for solving it. The rest of the paper is organized as follows. Section 2 gives a detailed description of the problem encountered by the company, in particular the features of the production environment. Section 3 presents how the problem was modelled, with the assumptions made together with the company. Our resolution method is then described in Section 4. Numerical experiments are shown in Section 5, together with the analysis that was performed. Finally, Section 6 ends the paper with some conclusions and perspectives.

Description of the problem

The problem of selecting equipment for designing a production workshop is formulated as a scheduling problem that is solved for each design scenario of the workshop. This approach is possible because the number of alternatives is limited. This section describes the various features of the production environment and the connection between the resulting scheduling problem and those studied in the scheduling literature.

The company has a set of N products (8 in the specific setting we are considering) J = {J 1 , J 2 , . . . , J N } that must be performed in the planning horizon H. Each product J i must be performed in a number of units that differs from one product to another. There are 8 machines (i.e. workstations) in the workshop M = {M 1 , M 2 , . . . , M 8 }, where machines M 1 and M 3 (resp. M 2 and M 4 ) are of the same type (i.e. can process the same set of operations). Each machine can perform no more than one product at a time and is assumed to be continuously available during the planning horizon. Five different tasks T ∈ {T 1 , . . . , T 5 } are performed in the production workshop. Task T 1 can be carried out by machines M 2 or M 4 . Task T 2 can be performed by any machine in the set {M 1 , M 3 , M 5 }, where M 5 is an optional machine that the company could acquire or not. Tasks T 3 , T 4 and T 5 are performed on machines M 6 , M 7 and M 8 , respectively. The layout of the workshop is described in Figure 1.

The products enter the workshop via a loading workstation (called entrant workstation in this paper) where they are put on pallets. Task T 1 is the first operation to be performed for each product. The remaining sequences of operations slightly vary from one product type to another, but the estimated processing times can be very different. Some tasks are not needed for some products and others are repeated several times. The sequence of tasks is (T

1 , T 2 , T 3 , T 2 , T 1 , T 2 , T 3 , T 4 , T 5 ) for products J 1 , J 2 and J 3 , (T 1 , T 2 , T 1 , T 2 , T 3 , T 4 , T 5 ) for product J 4 , and (T 1 , T 2 , T 1 , T 2 , T 4 , T 5 ) for products J 5 through J 8 .
Due to the size and weight of the products, there is no buffer space for storing products between machines. Thus, if a product finishes its processing on a machine and the machine on which the product has to be processed next is occupied, then the product has to wait on the current machine, thus blocking it for performing other products. This feature is referred in the scheduling literature to as the blocking constraint [START_REF] Hall | A survey of machine scheduling problems with blocking and no-wait in process[END_REF]. It is also called the hold-while-wait constraint by authors who address the scheduling and deadlock avoidance problems (Mati et al. 2001b).

With the above features, and according to the scheduling literature, the problem can be seen as a flexible job-shop scheduling problem with blocking constraints, for which some operations have several possibilities to be performed and the repetition of operations is allowed, i.e. the same operation may appear several times in the sequence of a product. An important feature that makes this scheduling problem different from other problems with blocking investigated in the literature is that poor sequencing of machines may lead to undesirable deadlocks (Mati et al. 2001a,b). More precisely, let us consider two products P 1 and P 2 with two operations each. The first operation of P 1 can be performed by M 1 or M 2 and the second operation by M 3 . The first operation of P 2 is performed by M 3 and the second by M 2 . Thus, if one decides to assign M 2 to the first operation of P 1 and to allow the first operations of the two products to be performed simultaneously, then a deadlock situation occurs since P 1 cannot move to its second operation (M 3 is occupied by P 2 ) and P 2 cannot progress to its second operation (M 2 is occupied by P 1 ).

Moreover, the manufacturing environment is characterized by the following features that lead to an even more complex scheduling problem.

• Limited number of pallets. Products to be produced are put on pallets at the loading (or entrant) workstation. After having visited all necessary machines in the manufacturing process, a product is removed from the pallet on which it is attached, and the pallet is moved to the entrant workstation for processing a new product. Due to the cost of a pallet, only a limited number is allowed limiting the number of products in the workshop at the same time. on the machine on which it is processed next, this machine must be available. This is a no-wait constraint that prevents the AMHS to hold products.

• Storage on machines. If the machine that processes the next operation of a product is occupied then the product can be stored on any empty machine. This alternative minimizes the blocking times and improves the utilization of critical machines. • Negative time lags. Some operations can be started even though the preceding operation on the routing is not yet completed. This time corresponds to the preparation of the machine by an operator, and represents about 5% of the processing time.

The company was interested in evaluating various choices that had to be made before building the workshop. The first choice concerns the number of entrant workstations from which products are brought into the workshop. Since the processing of products begins with task T 1 , it could be possible, at an additional investment cost, to also use M 4 as an entrant workstation, instead of only M 2 . The second choice concerns the need of purchasing a new workstation (i.e. a new machine M 5 ) and thus changing the functions of some workstations by modifying manufacturing processes. The last choice concerns the determination of the right number of pallets to be used in the workshop. A smaller number of pallets means less products in the workshop, which may decrease the utilization of machines and therefore increases the total time to process all products. On the other hand, more pallets may lead to more products in the workshop and thus increased congestion due to blocking constraints.

We shall show in the next sections how we modelled and solved the scheduling problem to help the company answer these three questions.

Modelling the problem

Taking into account all the features of the production environment results in a complex scheduling model. In order to solve the problem, some constraints were relaxed, validated by the company. This section first discusses how the constraints are modelled and then specifies those that are relaxed. The manufacturing processes of jobs are modified. In what follows, we use the following notation to describe an operation O i,j of a product J i : (S 1 , . . . , S t , p i,j ), where S k ⊂ M is the k th set of resources that can be used to perform O i,j , and p i,j is the processing time which does not depend on the chosen set of resources. As it will be shown below, the resources in our problem are the machines, the Automated Material Handling System (AMHS) and the pallets.

Let us illustrate our modelling choices in the following subsections on two consecutive operations O i,j and O i,j+1 of a product J i described by the routing (1):

J i ≡ . . . ; ({M 1 }, p i,j ); ({M 2 }, {M 3 }, p i,j+1 ); . . . (1)
This means that, at operation O i,j , the product J i must be processed on machine M 1 with a processing time of p i,j and, at operation O i,j+1 , can be processed either on machine M 2 or on machine M 3 with a processing time of p i,j+1 .

Storage on machines

To integrate the possibility of storing J i between O i,j and O i,j+1 on machines M 4 or M 5 , a waiting operation O i,j + = ({M 1 }, {M 4 }, {M 5 }, 0) is added between O i,j and O i,j+1 . Machines M 1 is integrated in the set of alternative machines to express the fact that J i can also wait on this machine. If one of the machines that can perform O i,j+1 is empty, then J i can move to O i,j+1 by advancing instantaneously to O i,j + on machine M 1 and then to O i,j+1 . The new manufacturing process of J i is given by the routing (2):

J i ≡ . . . ; ({M 1 }, p i,j ); ({M 1 }, {M 4 }, {M 5 }, 0); ({M 2 }, {M 3 }, p i,j+1 ); . . . (2)

Material handling

The AMHS can be explicitly integrated by adding a transportation operation between any two consecutive operations. The added operation (noted j + in the example below) needs the AMHS and the processing time corresponds to the time needed for moving products. Adding this information to the manufacturing process of job J i leads to the one given by the routing (3) where A denotes the AMHS:

J i ≡ . . . ; ({M 1 }, p i,j ); ({A}, t j,j + ); ({M 1 }, {M 4 }, {M 5 }, 0); ({A}, t j + ,j+1 ); . . . (3)
where t j,k is the time needed for the AMHS to bring product J i from the workstation selected for operation O i,j to the workstation selected for operation O i,k . This time depends on the distance between these workstations in the workshop.

Pallets

Up to this point, every operation needs only one resource to be processed. However, it is necessary to take into account the fact that a pallet is needed for each product from the beginning until the end of its process. Hence, the pallets are modelled as a resource R which is available in a given number of units. The manufacturing process is then modelled by the routing (4):

J i ≡ . . . ; ({M 1 , R}, p i,j ); ({A, R}, t j,j + ); ({M 1 , R}, {M 4 , R}, {M 5 , R}, 0); . . . (4)
The fact that the pallet moves to the entrant workstation after completing the product is modelled by adding an operation at the end of the routing of each product. This operation needs the pallet and the AMHS with a processing time equal to the time needed to move to the entrant workstation.

Note that if the number of pallets is equal to one, the scheduling problem becomes trivial and the solution consists of bringing products one after another in the workshop. If the number of pallets is equal to the number of machines, the problem is equivalent to a flexible job-shop scheduling problem with blocking and thus pallets can be omitted.

The relaxed constraints

The no-wait constraint that prevents the AMHS to hold the products is always fulfilled since the schedules generated by the solving method are feasible. Otherwise there may be a case where the AMHS moves a product from a given machine M i and finds the next machine M j already occupied. This is a deadlock situation since M j can only be released by the AMHS.

The first constraint that is relaxed in our modelling and solving method is the transportation time that depends on the distance between machines. We use an average time of the transportation times between all workstations of the workshop. After solving the scheduling problem, we recalculate the starting and finishing times of operations by integrating the actual values of the transportation times. Of course, if a waiting operation uses the same previous machine in the manufacturing process (for instance M 1 in the routing (4)) then the moving time is equal to zero.

The second relaxed constraint is the time lags between two consecutive operations on the routing. We solve the problem without these times and, after finding a solution, a shifting procedure is called to find the actual starting and finishing times of operations.

The third relaxed constraint is related to the limited number of pallets n p . Integrating this feature in our scheduling model is difficult for the resolution method of (Mati and Xie 2008b), which will be used for evaluating the fitness of chromosomes in the proposed algorithm. One way would be to model each pallet as a separate resource and to integrate flexibility each time an operation needs a pallet. This solution is not convenient as it leads to a huge number of alternative sets (Dauzère-Pérès and Paulli 1997). In our case, we assume that the number of pallets is infinite and, during the search process, we save the best schedule with a given maximum number of pallets (n p ∈ {3, 4, 5, 6, 7} in our numerical experiments). The number of pallets required for a given schedule is determined by finding, over the entire schedule, the maximum number of jobs progressing in parallel in the workshop.

With the above modelling and simplification in mind, the resulting scheduling problem that is solved is a multi-resource job-shop scheduling problem with resource flexibility and blocking, in which an operation may need more than one resource. The AMHS is modelled as a real resource that must be taken into account when building a schedule.

Solving the problem

One way for tackling the problem is to model it as a mixed integer linear programming model and to solve it using a commercial solver. Computational experiments carried out by Mati and Xie (2008b) show that commercial solvers are not able to provide reasonably good solutions for the flexible job-shop scheduling problem with blocking constraints, even for instances with less than 10 jobs and 5 machines.

The method used in this paper to answer the questions of the company is an efficient genetic algorithm that uses the greedy procedure proposed by Mati and Xie (2008b) for evaluating the fitness of chromosomes. In this section, genetic algorithms are first briefly described. Then the basic parameters of our genetic algorithm are discussed, including the crossover and mutation operators and the strategy for population management. The method used for evaluating the fitness of the chromosomes is also recalled.

A brief description of genetic algorithms

Genetic algorithms are population-based procedures that start from an initial population containing a set of feasible solutions represented by chromosomes. The first question to answer when using genetic algorithms is the chromosome encoding (or encoding scheme). There are two types of encoding, the direct and indirect encoding, and the selection of the suitable type depends on the characteristics of the problem. In the direct encoding, a chromosome describes a solution whereas, in the indirect encoding, the chromosome contains information that is used to build a solution. Genetic algorithms iteratively apply genetic operators to replace the current population by a new and hopefully improved population. The first operator is the crossover operator that combines the characteristics of two selected chromo- somes to generate one or several offsprings. Chromosomes may also be submitted to a mutation operator that consists of a small perturbation of the genes. The purpose of mutation is to prevent the search process from being trapped in local optima if chromosomes in the population become too similar. Genetic algorithms also require a suitable mechanism for population management that decides which chromosomes are selected at each generation for applying the crossover operator. It also requires a fitness function that represents a measure of the quality of each chromosome.

To describe our genetic algorithm, we present below the different parameters of the algorithm, which are fairly standard except the fitness evaluation function that is based on a greedy procedure proposed by Mati and Xie (2008b). somes, mapping the portion of one parent located between these positions onto the other parent, and then exchanging the remaining information (see [START_REF] Goldberg | Alleles, loci, and the traveling salesman problem[END_REF]). After using the crossover operator, a mutation operator is applied that deletes a randomly selected job from its current position in the chromosome and inserts it at another randomly selected position. One of the advantages of this mutation operator is its good ability to modify the chromosomes and thus to diversify the population. The crossover operator is applied with a probability equal to 0.9, and the mutation probability is equal to 0.05. These values are standard for scheduling problems and, as for the population size, have been validated on academic test instances.

Stopping criterion. The genetic algorithm terminates after a fixed number of generations. This number has been set to 1,000 through experiments on the industrial data sets. Moreover, if no improvement of the fitness function is done after 200 generations, the algorithm is also stopped.

Fitness evaluation

The fitness function evaluates the quality of the chromosomes and plays a major role in the selection process. The solution method proposed by Mati and Xie (2008b) is used to measure the fitness of the chromosomes. This method is a greedy heuristic that schedules jobs consecutively according to the job sequence provided by the chromosome. It starts by optimally scheduling the two first jobs using a polynomial algorithm, called geometric approach, that reduces the scheduling problem to a shortest path problem in a two-dimensional plane with obstacles. After obtaining the optimal schedule for the two jobs, a combined job is defined by combining into a single operation the operations that are performed in parallel. Resource assignment and operation sequencing decisions of scheduled jobs will not be modified. The geometric approach then solves the scheduling problem of the combined job and the next job in the order of the chromosome. The resulting schedule is again used to define a combined job, and the process continues until all jobs are scheduled. The makespan of the final schedule corresponds to the fitness evaluation of the chromosome. For a detailed description of the greedy heuristic for solving shop scheduling problems, the reader is referred to Mati et al. (2001a), Mati and Xie (2008a,b).

Numerical experiments and analysis

We only present and discuss the results obtained on scenarios with 5 or 6 pallets since the remaining scenarios (with 3, 4 and 7 pallets) give solutions with larger total completion times and waiting times. The computational experiments were conducted on a one-month horizon, during which the workshop must produce in total 30 units of products (see the units per product in Table 1), and 323 hours of work per month (the workshop is open 17 hours per day and 19 days per month). The genetic algorithm was implemented in C++, and the experiments were ran on a Pentium IV (2 GHz). Each scenario took about 10 minutes of CPU time.

The combination of the three parameters (i.e. without or with one new workstation, one or two entrant workstations, five or six pallets) led us to define 8 scenarios. The results for the scenarios are detailed in Tables 2 through 5 and discussed below. For each scenario, the tables show the mean cycle time per product, the total waiting time per machine and the total completion time of all products. The cycle time is the time the unit of a product stayed in the workshop (i.e. completion time minus arrival time), and the mean cycle time is provided for the number of units of each product.

Scenarios with no new workstation

With one entrant workstation, the total completion time is 374.51 hours (respectively 375.09 hours) if 5 pallets (resp. 6 pallets) are used. The solutions are detailed in Table 2. It is interesting to note that, in this configuration, allowing one more pallet only slightly improves the total waiting time on the machines.

With two entrant workstations, the total completion time is reduced to 367.70 hours (resp. 361.28 hours) if 5 pallets (resp. 6 pallets) are used. As shown in the solutions of Table 3, adding one more pallet improves more the total waiting time than in the previous case. The gain is equal to 3.87 (i.e. 378 -374.13) in the previous case and equal to 35.79 (i.e. 460.96 -425.17) in this case.

It must be noted that using two entrant workstations instead of one leads to an improvement of the total completion time (ranging from 7 to 14 hours, i.e. less than one working day), but also to much larger waiting times, in particular when only 5 pallets are used. This can be explained by the fact that, when entrant workstations are not at all bottlenecks, arriving products are quickly pushed forward in the workshop.

Scenarios with a new workstation

With one entrant workstation, the total completion time is 369.75 hours (resp. 354.91 hours) if 5 pallets (resp. 6 pallets) are used. Hence, using 6 pallets is clearly much better and, as shown in Table 4, also leads to drastic improvements of the total waiting time (from 517.35 hours to 430.63 hours).

With two entrant workstations, the total completion time is reduced to 345.38 hours (resp. 332.16 hours) if 5 pallets (resp. 6 pallets) are used. Again, the solutions in Table 5 show that using 6 pallets also improves the total waiting time.

Using two entrant workstations reduces the total waiting time and the total completion time : the gain on the total completion time ranges from 22 hours to 24 hours (369.75 -345.38 = 24.37 hours in case of 5 pallets, and 354.91 -332.16 = 22.75 hours in case of 6 pallets), i.e. more than 1 working day. It is interesting to note that, contrary to the case without a new workstation, the total waiting time does not increase with two entrant workstations. This is because the capacity of the workshop has increased, and thus it can handle the products that are sent by the entrant workstations.

Answers to the questions of the company

According to the results presented in the previous subsections, our answers to the questions of the company were the following:

• It is always better to have two entrant workstations, although the improvement is really interesting, given the investment, only if one new workstation is added. If one considers the possible configurations, the results of the experiments show that the worst configuration is to use one entrant workstation without a new workstation. An improved configuration is to use either two entrant workstations but no new workstation, or one new workstation but only one entrant workstation. However, in both cases, the investment might be hard to justify given the improvement. Finally, the best configuration is to use two entrant workstations combined with a new workstation.

Finally, the experiments show that 323 working hours per month are not sufficient to handle the production, and hence the company should increase the number of working hours.

Conclusion

We have shown how a complex practical scheduling problem could be solved by making some realistic assumptions validated by the company. The various scenarios that have been run helped the company to take important decisions on key features of the new workshop. The next step will probably be, when the workshop is operational, to use the proposed resolution method to design a Decision Support System for its weekly planning. Other characteristics will most likely have to be taken into account, such as the availability of workers.

F

•

  Material handling. The workstations are linked by an Automated Material Handling System (AMHS) that moves the products between the machines, and be-

	February 11, 2010 16:37 International Journal of Production Research Mati˙Lahlou˙Dauzere-
	Peres-IJPR-revision2
	International Journal of Production Research	5
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

tween the machines and the loading/unloading (or entrant) workstations. The AMHS can handle only one product at a time. When the AMHS moves a product

Table 1 .

 1 Product J 1 J 2 J 3 J 4 J 5 J 6 J 7 J 8 Number of units per product.

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

13

Page 13 of 19 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk

Table 2 .

 2 Results with no new workstation and one entrant workstation.

	Page 14 of 19

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk

Table 3 .

 3 Results with no new workstation and two entrant workstations.

	Page 15 of 19

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	February 11, 2010 16:37 International Journal of Production Research Mati˙Lahlou˙Dauzere-
				Peres-IJPR-revision2
	Product Mean cycle time (hours) Machine	Waiting time (hours)
		5 pallets	6 pallets		5 pallets 6 pallets
	J 1	66.93	66.16	M 1	45	45.21
	J 2	65.64	80.78	M 2	110.64	79.7
	J 3	65.98	66.9	M 3	17.06	10.78
	J 4	47.32	41.81	M 4	95.85	60.35
	J 5	37.98	30.51	M 5 (new)	81.33	87.11
	J 6	35.06	44.9	M 6	55.1	62.24
	J 7	67.17	75.31	M 7	41.82	41.82
	J 8	38.29	32.04	M 8	70.55	43.42
				Total	517.35	430.63
		Total completion time (hours)	369.75	354.91
	F o r
		P
		e
			e r
			R
			e
				v i e
				w
					O n l
						y

Table 4 .

 4 Results with a new workstation and one entrant workstation.

16

Page 16 of 19 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	February 11, 2010 16:37 International Journal of Production Research Mati˙Lahlou˙Dauzere-
				Peres-IJPR-revision2
	Product Mean cycle time (hours) Machine	Waiting time (hours)
		5 pallets	6 pallets		5 pallets 6 pallets
	J 1	60.06	52.7	M 1	14.8	15.35
	J 2	61.89	70.37	M 2	46.72	33.39
	J 3	61.06	64.81	M 3	16.57	24.38
	J 4	55.01	45.51	M 4	82.72	90.88
	J 5	30.91	40.61	M 5 (new)	50	68.12
	J 6	28.85	34.07	M 6	48.1	23.46
	J 7	69.7	64.46	M 7	38.66	19
	J 8	38.33	36.06	M 8	70.88	49.41
				Total	368.45	323.99
		Total completion time (hours)	345.38	332.16
	F o r
		P
		e
			e r
			R
			e
				v i e
				w
					O n l
						y

Table 5 .

 5 Results with a new workstation and two entrant workstations.

	Page 17 of 19

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	February 11, 2010 16:37 International Journal of Production Research Mati˙Lahlou˙Dauzere-
	Peres-IJPR-revision2
	List of figure captions
	1. Layout of the production workshop.
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Parameters of our genetic algorithm

Chromosome encoding. In order to exploit standard crossover operators and to keep the feasibility of the solutions in the reproduction process, we use an indirect encoding in which a chromosome contains an n-string corresponding to a permutation of all integers from 1 to n. Each gene in the chromosome represents a job (product). The chromosome is then represented by a vector that gives the preference list for the treatment of all jobs. This list is used, by the fitness function described bellow, to build the actual schedule associated to the chromosome.

Initial population. In many implementations of genetic algorithms, a good starting population is important for the success of the algorithm. In our case, this population is constructed in a very simple way by randomly generating the chromosomes. However, we have found that this generally leads to a population with a very poor quality. Since it is not easy to generate a good initial population we have focused on improving it, mainly by increasing the diversity of chromosomes during the population management phase, as explained below. The size of the population is set to 20 chromosomes. This population size was fixed through numerical experiments conducted on academic test instances.

Chromosome evaluation. A greedy heuristic proposed by Mati and Xie (2008b) is employed to evaluate the chromosomes. This heuristic is recalled below.

Population management. There are two common strategies to manage the population: population replacement and incremental replacement. In the first case, the new population contains only chromosomes created by crossover and mutations operators, whereas, in the second case, only one chromosome of the parents is replaced with a generated chromosome. The population replacement has been preferred in our application since it allows a high population diversity, allowing the greedy heuristic to start from diverse chromosomes. A large part of the solution space can then be explored and hence the quality of the initial population can be improved. For selecting pairs of parents to be crossed, the simplest and most used roulette wheel selection is applied, which uses a distribution on the current population that depends on the fitness of the solutions. The probability of each chromosome to be selected is directly proportional to its fitness value.

Crossover and mutation operators. A general and suitable operator for crossover, in case of representing chromosomes by permutations, is the partially mapped crossover. This crossover has been applied successfully on many scheduling problems, and more specifically on shop problems. The crossover consists of randomly selecting two cut positions in parents chromo-