Ghailen Laouini

Olivier Meste

member IEEE Marianna Meo

Analysis of Heart Rate Variability Using Time-Varying Filtering of Heart Transplanted Patients

In this paper, we analyze the heart rate variability (HRV), obtained by using the time-varying integral pulse frequency modulation (TVIPFM) which is well adapted to the exercise stress testing. We consider that the mean heart period is varying function of time, during exercise. This technique allows the estimation of the autonomic nervous system modulation (ANS) from the beat occurrences. The estimated respiratory sinus arrhythmia is then filtered in the time-frequency domain around the respiration using a time-varying filter. It is proven that the Spectrogram is a convenient time-frequency representation that allows the implementation of such filter. The recorded data comes from exercise test performed by ten heart transplant patients. The magnitude of the filtered modulation of the heart rate due to respiration is compared to the date of transplantation taking into account the volume of respiration. It reveals that the normalized magnitude of the filtered variability, is significantly increased as the age of transplantation is higher with a high correlation coefficient (R=0.74, p=0.01). This correlation raised to 0.82 when considering dynamic behavior of the parameters.

Applied to our dataset, standard parameter fails to exhibit such correlation.

I. INTRODUCTION

Several works [START_REF] Bayly | Spectral analysis of pulse frequency modulation in the nervous system[END_REF] have used the Integral pulse frequency modulation (IPFM) model to explain the regulation of the heart rate by the ANS, in different physiological situations, even during stress testing exercise [START_REF] Bailon | Analysis of heart rate variability during exercise stress testing using respiratory information[END_REF]. However the IPFM model assumes a constant threshold [START_REF] Mateo | Improved hart rate variability time-domain signal construction from the beat occurence time according to IPFM model[END_REF][START_REF] Mateo | Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal[END_REF] considering a constant mean heart period. During dynamic protocol the use of the time-varying threshold IPFM (TVIPFM) model to correct the heart rate variability is relevant for the estimation of the modulating signal [START_REF] Bailon | The integral pulse frequency modulation model with time-varying threshold: application to heart rate variability analysis during exercise stress testing[END_REF][START_REF] Laouini | Evidence of the influence of respiration on the heart rate variability after human heart transplantation: Role of observation model[END_REF]. Many authors have emphasized the close relationship between respiration and HRV [START_REF] Meste | Some considerations on the IPFM Model for the heart rate variability analysis[END_REF][START_REF] Meste | Time-varying analysis methods and models for the respiratory and cardiac system coupling in graded exercise[END_REF][START_REF] Blain | Assessment of ventilatory thresholds during graded and maximal exercise test using time-varying analysis of respiratory sinus arrythmia[END_REF]. This component appears to be the principal non stationary modulation of high frequency. Since it is embedded in a mixture of modulating signals, its analysis requires a filtering that is time-varying in that case. The performance of stress test exercises allows the use of the respiratory sinus arrhythmia (RSA) information for the ANS tone assessment. This modulation is effective with respects to the innervation of the heart. After orthotopic heart transplantation the loss of the innervation should suppress this modulation. During long-term follow-up it has been shown [START_REF] Folino | Heart rate variability in patients with orthotopic heart transplantation: Long-term follow-up[END_REF] that the HRV is higher as long is the age of the transplantation, with a strong correlation. In the study [START_REF] Folino | Heart rate variability in patients with orthotopic heart transplantation: Long-term follow-up[END_REF] the variability measurement is addressed by using the standard deviation of the data during resting conditions. G. Laouini is with Laboratory I3S, University of Nice and CNRS, Sophia Antipolis 06903, France, (e-mail:laouini@i3s.unice.fr) O. Meste is with Laboratory I3S, University of Nice and CNRS, Sophia Antipolis 06903, France, (e-mail:meste@i3s.unice.fr) Thus a correlation between this parameter and the age of the transplantation is exhibited under controlled protocol, limiting the impact of such analysis. It would be of greater interest to show this correlation with relaxed and dynamic conditions corresponding to increasing physical exercise.

The paper is organized as follows. The methods and materials are presented in Section II, where we use a mathematical development to demonstrate the feasibility of filtering in the time-frequency domain, using the Spectrogram. The TVIPFM model is described as well in this section to obtain a corrected HRV signal. Section III presents the results and Section IV the discussion.

II. METHODS

A. Mathematical development

In this study the signal to be filtered is non stationary. It means that the characteristics of the bandpass filter must vary with time. The time-frequency transformation that will permit this particular filtering is the Spectrogram defined by:

S(t, f ) = | x(s)h * (s -t)e -i2πf s ds| 2
where h(t) stands for the smoothing window. In the following it will be demonstrated that by using this representation in addition with a matched template G(t, f ), an efficient time-varying filter can be applied. Note that continuous time is considered in the entire development but results are still valid with sampled time. The spectrogram being the squared magnitude of the short Fourier Transform, it is linked to the Wigner-Ville (WV) distribution by the convolution products in time and frequency:

S(t, f ) = (W h * * W x )(t, f )
Assuming that the WV transform of h(t) is separable such that W h (t, f ) = g(t)l(f ) (verified when h(t) is a gaussian function), when x(t) is a linear frequency modulation, i.e. x(t) = exp(iπαt 2 ), the Spectrogram is then:

S(t, f ) = g(t -s)l(f -ν)δ(s -αν)dsdν S(t, f ) = g(t -αν)l(f -ν)dν
Note that the signal x(t) is locally approximated by this model whose modulation slope is α. This assumption is valid with respect to the respiration component where even during exercise the frequency varies smoothly. This quantity 34th Annual International Conference of the IEEE EMBS San Diego, California USA, 28 August -1 September, 2012 3436 978-1-4577-1787-1/12/$26.00 ©2012 IEEE is nonzero when the two windows g and l are simultaneously nonzero.

g(t) = 0 f or -T < t < T l(f ) = 0 f or -F < f < F (1)
where g(t) and l(f ) are assumed hermitian symmetric.Then:

g(t -αν) = 0 f or -T +t α < ν < +T +t α l(f -ν) = 0 f or -F + ν < f < F + ν (2)
It can be concluded that for a given time t 0

-F + -T + t 0 α < ν < F + +T + t 0 α So [-F + -T +t0 α ; F + +T +t0 α
] is the support of S(t 0 , f ) along the frequency axis. The mean value of this interval is t0 α that corresponds to the instantaneous frequency of x(t) for t = t 0 .As shown in the sequel, because this result means that the support is centered on the instantaneous frequency of x(t), the integration around this interval will provide the magnitude without bias. This allows us to impose to the template G(t, f ), that will be multiplied with S(t, f ), to be nonzero around f = t0 α for t = t 0 in order to get rid off the other components. It is well known that for the Spectrogram, the marginal is:

S(t, f )df = (p x * p h )(t)
with p x (t) being the squared envelop |x(t)| 2 of x(t) . If this envelop is locally a constant A, then the computation of this integration will produce A 2 if |h(t)| 2 = 1. This result imposes to the smoothing window h(t) to be energy normalized. The slow variation of the respiration component in the heart period perfectly match the local assumption.

This demonstration proves the role of the time-varying filtering around the respiration to retrieve the envelope of the respiration component observed within the HRV signal. This modulation will be corrected or not by the TVIPFM model. The template G(t, f ) that conveys the knowledge of the respiration frequency is defined as:

G(t, f ) = 1 if f ∈ B R ∀t G(t, f ) = 0 otherwise (3)
With B R is the band frequency of respiration. Details on the computation of the time-varying filter using G(t, f ) is fully provided in [START_REF] Meste | Time-varying analysis methods and models for the respiratory and cardiac system coupling in graded exercise[END_REF].

B. The time-varying threshold Integral pulse frequency modulation model

The TVIPFM model is based on the hypothesis that the ANS influence on the sinus atrial SA node can be represented by the modulating signal m(t). In the TVIPFM model the integral of 1+m(t) is compared to a time-varying threshold T (t), representing the time-varying mean heart period. Assuming that the first beat occurs at time t 0 = 0, the beat occurrence time series can be generated as solution of [START_REF] Rompleman | The measurement of heart rate variability spectra with the help of a personal computer[END_REF]:

k = t k 0 1 + m(t) T (t) dt
with its extension to continuous case:

k(t) = t 0 1 + m(s) T (s) ds
the instantaneous heart rate is defined as below

d HR (t) = 1 + m(t) T (t) = dk(t) dt
Supposing the variations of the term 1 T (t) are slower than those of the term m(t) T , and their spectral components don't overlap, we define the time-varying heart rate as:

d HRM (t) = 1 T (t)
and the heart rate variability signal is defined as

d HRV (t) = d HR (t) -d HRM (t) = m(t) T (t)
So, to obtain the modulating signal, we have to correct the HRV signal d HRV (t) by the time-varying mean heart rate d HRM (t).

In the TVIPFM, the variations of the time-varying threshold T AC (t), are small compared to its mean value T DC . Rewriting the instantaneous heart rate

d HR (t) = 1 + m(t) T (t) = 1 + m(t) T DC + T AC (t) = 1 + m(t) T DC (1 + T AC (t) T DC ) ≈ 1 + m(t) T DC (1 + T AC (t) T DC )
with the assumption m(t) < 1, we can neglect the second order term, then we obtain:

d HR (t) = 1 + m(t) -T AC (t) T DC T DC we assume d HRM (t) = 1 -T AC (t) T DC T DC = 1 T DC
To obtain the modulating signal m(t), first we estimate d HR (t), where k(t) is estimated from pairs (t k , k), then we estimate the d HRM (t) term by low-pass filtering d HR (t) and we calculate the d HRV (t) term. Finally, the modulating signal is obtained as:

m(t) = d HRV (t) d HRM (t)

III. MATERIALS AND RESULTS

Recorded signals have been obtained from surface electrodes that produce interferences such as noise and baseline wander. So, to exploit all signals we must cancel the disturbances caused by breathing. First we filter the signal using 500-th order high-pass finite impulse response filter in order to remove the baseline generated by the respiration. After using a threshold technique we demodulate the filtered signal to detect all the time occurrences of R-waves, namely the t k . The differences of consecutive t k provide the heart period. Finally, the HRV is computed by removing the trend of the heart period. In addition to the ECG, the respiration signal is recorded to extract the frequency and the volume. Ten heart transplant subjects with different ages of transplantation months are analyzed by using different parameters. The parameter introduced in [START_REF] Folino | Heart rate variability in patients with orthotopic heart transplantation: Long-term follow-up[END_REF] is indeed a standard deviation (stand) evaluated on the HRV at rest. It will be compared to the quantity mag(t) corresponding to the magnitude of the time-varying filtered respiration component extracted from the HRV. When corrected by the TVIPFM model the subscript will be tv. Because the respiration volume has not been controlled during the protocol, the time-varying quantities will be normalized by the respiration volume computed from the respiration signal. This normalization denoted by the subscript resp will be exclusively applied to mag because it contains only the respiration influence. All these quantities are computed at rest and at the maximum of the exercise. An interval of 100 seconds will be used to compute parameter stand and to average the quantities mag(t). In summary, the analysis will provide the parameters stand, mag, mag tv and mag tv,resp . Figure [START_REF] Bayly | Spectral analysis of pulse frequency modulation in the nervous system[END_REF], reveals that the RSA modulation is high at the beginning of exercise, then the magnitude of this modulation decreases abruptly after the beginning of exercise due to the exercise pressor reflex causing ventilatory responses to exercise. Then period from second 1400 to 1800 corresponds to the increasing exercise. The intensity of the RSA is at lower values during this interval. During the recovery, there is an abrupt increase which gradually decreases as the heart rate decreases reaching values close to those observed at the beginning of the test. The corrected variability with the TVIPFM model (top), the magnitude of the envelop filtered around the respiration frequency (middle) and the breathing volume (bottom), the red vertical lines correspond to the begining and the maximum of the exercise, respectively. Note that the magnitude unit is meaningless because of the correction procedure, the volume has been normalized by the greater value within all the subjects Figure [START_REF] Mateo | Improved hart rate variability time-domain signal construction from the beat occurence time according to IPFM model[END_REF] shows first the HRV corrected with the TVIPFM model (upper). This variability is then filtered in the timefrequency domain, producing an estimation of the RSA magnitude (middle). If needed, the volume of the respiration is computed from the respiration signal (bottom). From this function of time (middle) an average over 100 seconds is computed in the resting period and at the maximum of the exercise. Table I, displays all the parameters extracted from the ten patients. Its clear that mag tv,resp outperforms the other parameters whereas the stand parameter fails in exhibiting a linear correlation. The difference (rest)-(max) is provided in order to investigate a dynamic behavior. Once again, mag tv,resp shows a clear correlation. Figure [START_REF] Mateo | Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal[END_REF], represents the magnitude of the heart rate variability mag tv,resp (rest), filtered around the frequency band of respiration and normalized with the respiratory volume, function of the age of transplantation. It shows the high correlation (R=0.74, P=0.01), between the age of the transplantation and the RSA with TVIPFM correction and respiration normalization.

One of the important results in this paper is that from the parameters obtained from the HRV, only the variability corrected by the TVIPFM and normalized is highly correlated with the age of the transplantation. The result was expected from the literature but with a parameter that fails in that case (stand). However, since a higher variability is related to a reinnervation, our result is in line with the physiology knowledge.

IV. DISCUSSION

In this paper, a method of time-varying filtering is proposed, centering the filter around the respiration frequency band. The relevance of using the Spectrogram for this aim is proven. In combination with this filter we use the TVIPFM model which has a time-varying mean heart period adapted to non stationary conditions. The time-varying mean heart rate is estimated by low-pass filtering d HR (t). Thus, the modulating signal is filtered in the time-frequency domain and the instantaneous power is then calculated by integrating over frequencies to obtain the mag tv .

Probably because no respiration profiles have been imposed in the protocol, the normalization with the respiration volume has increased the correlation. In contrast to standard parameter, namely the standard deviation stand, the proposed one corroborates the physiologic expectation on our data set. The reinnervation should be higher as the age of the transplantation is longer, inducing a higher variability. Note that the global variability addressed by the parameter stand didn't convey any positive result. The high correlation between the RSA and the age of transplantation is due to the TVIPFM correction and the time-varying filtering. This significant correlation is revealed at rest and not during exercise. The dynamic behavior corresponding to (rest)-(max) is even more correlated. This result is of great interest because it highlights the gain of the HRV dynamics due to the innervation process.

Fig. 1 .

 1 Fig. 1. The spectrogram of the heart rate variability with the template G(t, f ) boundaries (in black) defined by the respiratory frequency band

  Fig. 2.The corrected variability with the TVIPFM model (top), the magnitude of the envelop filtered around the respiration frequency (middle) and the breathing volume (bottom), the red vertical lines correspond to the begining and the maximum of the exercise, respectively. Note that the magnitude unit is meaningless because of the correction procedure, the volume has been normalized by the greater value within all the subjects

TABLE I THE

 I CORRELATION OF PARAMETERS WITH THE AGE OF TRANSPLANTATION At rest. The magnitude of the heart rate variability corrected with TVIPFM model filtered around the frequency band of respiration and normalized with the respiratory volume, function of the age of transplantations

											(rest)	(max)	(rest)-(max)
								stand		R=0.25, p=0.48	R=0.27 , p=0.45	R=0.11 , p=0.76
								mag		R=0.43, p=0.22	R=-0.21, p=0.56	R=0.29 , p=0.40
								magtv		R=0.61, p=0.06	R=-0.14, p=0.69	R=0.67 , p=0.03
							magtv,resp	R=0.74, p=0.01	R=0.09 , p=0.80 R=0.82, p=0.003
		6	x 10 -4
		5
	Arbitrary Unit	3 4
		2
		1
		0	0	20	40	60	80	100	120	140	160
							months
	Fig. 3.


Jean-Marc Roussel

Jean-Jacques Lesage

Algebraic synthesis of logical controllers despite inconsistencies in specifications

Keywords: Controller synthesis, Algebraic approaches, Boolean algebra, Dependable system

This paper deals with the problem of consistency of the requirements which are the starting point of controller synthesis methods for Discrete Event Systems (DES). In previous papers, an algebraic synthesis method for logical systems has been proposed. This method includes a theorem allowing the proof of consistency of the set of specifications. In this paper, we show that if inconsistencies are detected, the conditions of these inconsistencies can be given to the designer. It is therefore possible for the designer to propose priority rules between the involved requirements allowing the correction of incoherent specifications. Both the mathematical and the methodological aspects of this work are presented. In an illustrative purpose, the case study of a hydraulic press is developed.

INTRODUCTION

Programmable Logic Controllers (PLCs) are frequently used for the control of safety-critical systems. Since failure of such systems may have disastrous effects, the use of formal methods for the design of control algorithms is often mandatory (IEC 61508 (1998)). Automatic synthesis methods are among the best solutions to handle this problem and are one of the most active subjects of research in the field of Discrete Event Systems (DES) since the end of 80's (cf. numerous sessions organized in this field in the past WODES editions 1 ). Most part of recent works in this area are still based onto the Supervisory Control Theory (SCT) [START_REF] Ramadge | The control of discrete event systems[END_REF]) and are aiming to the synthesis of a Supervisor, and not directly to the Controller of an automated system. Furthermore, the use of state models (Finite Automata, Petri Nets...) and their composition for the construction of the models of the plant and of the specifications generates a complexity which remains problematic for the synthesis of a supervisor for complex systems [START_REF] Gohari | On the complexity of supervisory control design in the RW framework[END_REF]). It is therefore interesting to explore other ways for performing synthesis, such as algebraic approaches. In previous works, we proposed a method specifically developed to get the control laws that can be directly implemented into the controller [START_REF] Hietter | Synthèse algébrique de lois de commande pour les systèmes à évènements discrets logiques[END_REF]). We have chosen to synthesize these control laws under the form of recurrent Boolean equations because of the wide possibilities they offer for the formalization of safety requirements and for implementation.

Nevertheless, whatever is the synthesis method used, one of the weak links of the automatic generation of the control laws is the step of formal transcription (within state models or algebraic expressions) by the designer of the informal requirements and safety properties the controller has to satisfy. In the case of SCT, some authors 1 available at http://www.diee.unica.it/giua/WODES/ have proposed more or less generic approaches for the construction of the models of the plant [START_REF] Hanisch | A modular plant modeling technique and related controller synthesis problems[END_REF]) or of the specifications [START_REF] Roussel | Designing dependable logic controllers using the supervisory control theory[END_REF]). But in any case, the hypothesis that requirements can be inconsistent has never been taken into account, except in [START_REF] Zowghi | A logical framework for modeling and reasoning about the evolution of requirements[END_REF]) where authors propose a theorization of requirements evolution. Unfortunately in the framework of industrial collaborations we have been able to verify that it is often the case. In this paper we show how, in consideration of specific hypotheses, it is possible to install a correction loop for helping the designer to formalize these requirements and so, improving the synthesis method robustness to the lack of precision of the specifications. This paper is organized as follows. Some basics of algebraic synthesis are given in Section 2, and Section 3 recalls the main steps of our method. Section 4 presents the mathematical framework of our approach and new results that allow us to accept inconsistencies in specifications. The strategy we developed for making the synthesis more robust to the lack of consistency of the specifications is described in section 5 thanks to a case study.

PROBLEM STATEMENT

Fig. 1 proposes a generic representation of a DES whose the controller has p Boolean inputs (u i ), q Boolean outputs (y j ) and r Boolean state variables (x l ). Plant and Controller are connected through a closed loop exchanging inputs and outputs signals. The state variables, needed for expressing sequential behaviors of the controller, are represented by internal variables.

The algebraic modeling of the control laws of the controller necessitates the definition of (q + r) switching functions of (p + r) variables. Even if this representation is very compact (the r Boolean state variables allow the representation of 2 r different states), the construction by hands of

u i y j x l Controller Plant / q / p / r y j [k] = F j (u 1 [k], • • • , u p [k], x 1 [k -1], • • • , x r [k -1]) x l [k] = F q+l (u 1 [k], • • • , u p [k], x 1 [k -1], • • • , x r [k -1])
Fig. 1. A sequential DES these switching functions is a very tedious and error-prone task [START_REF] Huffman | The synthesis of sequential switching circuits[END_REF]): the controller of Fig. 1, admits 2 p inputs combinations, can send 2 q outputs combinations and can express (22 (p+r) ) (q+r) sequential behaviors. That is the reason why algebraic modeling approaches have been replaced by methods based on state models since the middle of 50's [START_REF] Mealy | A method for synthesizing sequential circuits[END_REF], [START_REF] Moore | Gedanken Experiments on Sequential Machines[END_REF]). Nevertheless, thanks to recent mathematical results obtained onto Boolean algebras [START_REF] Rudeanu | Lattice Functions and Equations[END_REF], [START_REF] Brown | Boolean Reasoning: The Logic of Boolean Equations[END_REF]), the automatic algebraic synthesis of switching functions is now possible.

In [START_REF] Princeton | On the synthesis of a reactive module[END_REF] an interesting approach for the systematic construction of a reactive program from its formal specification is proposed. In this work, the program synthesis is considered as a theorem proving activity. A program with input x and output y, specified by the formula ϕ(x, y), is constructed as a by-product of proving the theorem (∀x)(∃y)ϕ(x, y). The specification ϕ(x, y) characterizes the expected relation between the input x and the output y computed by the program. This approach is based on the observation that the formula (∀x)(∃y)ϕ(x, y) is equivalent to the second order formula (∃f )(∀x)ϕ(x, f (x)), stating the existence of a function f , such that ϕ(x, f (x)) holds for every x.

This approach provides a conceptual framework for the rigorous derivation of a program from its formal specification.

It has also been used to synthesize specifications under the form of finite automata from their Linear Temporal Logic (LTL) description [START_REF] Filiot | Antichains and compositional algorithms for LTL synthesis[END_REF]).

The core of our approach is based on this strategy: we aim at deducing the (q + r) switching functions of (p + r) variables which define the behavior of the controller from a formula ϕ(u

i [k], x l [k-1], y j [k], x l [k]) that holds for every k, every u i [k] and every x l [k -1].
To cope with combinatorial explosion, switching functions will be handled through a symbolic representation (and not their truth-tables which contains 2 (p+r) Boolean values). Each input u i (res. output y j ) of the controller will be represented by a switching function U i (res. Y j ). To take into account the recursive aspect of state variables, each state variable x l will be represented by two switching functions: X l (for time [k]) and p X l (for time [k -1]).

According to this representation, the synthesis of control laws of a logical system from its specification can now be transformed into the search of the solution of the mathematical problem:

(∀U i )(∀ p X l )(∃Y j )(∃X l )ϕ(U i , p X l , Y j , X l )
where (U i , p X l , Y j , X l ) are (p + q + 2r) switching functions of (p + r) variables.

OVERVIEW OF OUR METHOD

The input data of the proposed method (Fig 2 ) are unformal functional and safety requirements given by the designer. In practice, these requirements are most often given in a textual form and/or by using technical Taylor-made languages (Gantt Diagrams, Function Blocks Diagrams, Grafcet...). All the steps of our synthesis method are implemented into a prototype software tool developed in Python 2 . The first step is the formalization of requirements within an algebraic description. Requirements expressed with a state model can directly be translated into recurrent Boolean equations, thanks to the algorithm proposed by [START_REF] Machado | Logic controllers dependability verification using a plant model[END_REF]. Generic formalization of some requirements (change of operation modes, emergency management...) are also proposed for helping the designer to formalize "standard" requirements given in textual form (examples are given in Section 5.2). In case where the know-how of the designer enables him to build a priori the global form of the solution (or of a part of the whole solution) it is also possible to give fragments of solution as requirements (Hietter et al. (2008c)).

The second step consists in checking the consistency of the set of requirements by symbolic calculation. In a previous paper (Hietter et al. (2008a)), a theorem allowing the proof of consistency of the set of requirements has been proposed. In this paper, we show that if inconsistencies are detected, the conditions of these inconsistencies can be given to the designer. It is therefore possible for the designer to propose priority between the involved requirements allowing the correction of inconsistent specifications. We will develop more especially the mathematical results on which this step of the method is based in Section 4.3; it will be applied to a case study in Section 5.3.

The core of the method is the third step, which consists in the synthesis of the control laws. This step is performed by solving the system of equations which represents the set of consistent requirements. The mathematical results we have obtained (Theorem 2 recalled in Section 4.2), allow to find a symbolic expression of the set of solutions under the form of a parametric expression [START_REF] Hietter | Algebraic synthesis of dependable logic controllers[END_REF]).

In the fourth step of the method, a particular solution has to be chosen among the set of solutions. For that, a specific value of each parameter of the general solution has to be fixed. In some cases, optimal solutions, according to given criteria, can be automatically found [START_REF] Leroux | Algebraic synthesis of logical controllers with optimization criteria[END_REF]).

For ergonomic reasons, the synthesized control laws can finally be displayed under the form of a finite automaton [START_REF] Guignard | Symbolic generation of the automaton representing an algebraic description of a logic system[END_REF]).

After the mathematical background of the method has been recalled, we are going to show how, in consideration of specific hypotheses, the second step of the method can be improved by a correction loop helping the designer to formalize the requirements and so improving the robustness of our synthesis method to the lack of precision of the specifications.

MATHEMATICAL BASICS

Section 4.1 recalls classical notations and results onto the Boolean algebra of n-variable switching functions. Section 4.2 recalls the principal results on which our synthesis method is based. Section 4.3 presents the new results we use for helping the designer to formalize requirements.

The Boolean algebra of n-variable switching functions

To avoid confusion between Boolean variables and Boolean functions of Boolean variables, each Boolean variable b i is denoted as b b i . The set of the two Boolean values b 0 and b 1 is denoted as:

B = { b 0, b 1}. The classical two-element Boolean Algebra is ({ b 0, b 1}, ∨, ∧, ¬, b 0, b 1).
An n-variable switching function is a mapping 3 f :

B n → B where B = { b 0, b 1} ( b b 1 , . . . , b b n ) → f ( b b 1 , . . . , b b n )
Let F n (B) be the set of the 2 2 n n-variable switching functions. F n (B) contains (n + 2) specific n-variable switching functions: the 2 constant functions (0, 1), and the n projection-functions (f i Pj ). These functions are defined by: 0 : [START_REF] Brown | Boolean Reasoning: The Logic of Boolean Equations[END_REF].

B n → B ( b b 1 , . . . , b b n ) → b 0 1 : B n → B ( b b 1 , . . . , b b n ) → b 1 f i Pj : B n → B ( b b 1 , . . . , b b n ) → b b i 3 From Section 3.11 of
F n (B) can be equipped with three closed operations in order to obtain the Boolean Algebra (F n (B), +, •, , 0, 1):

Op. + : F n (B) 2 → F n (B) (f, g) → f + g Op. • : F n (B) 2 → F n (B) (f, g) → f • g Op. : F n (B) → F n (B) f → f where ∀( b b 1 , . . . , b b n ) ∈ B n , (f + g)( b b 1 , . . . , b b n ) = f ( b b 1 , . . . , b b n ) ∨ g( b b 1 , . . . , b b n ) (f • g)( b b 1 , . . . , b b n ) = f ( b b 1 , . . . , b b n ) ∧ g( b b 1 , . . . , b b n ) f ( b b 1 , . . . , b b n ) = ¬f ( b b 1 , . . . , b b n ) A Boolean formula onto F n (B) is a composition F(α 1 , . . . , α n ) of (α 1 , .
. . , α n ) by operations +, • and . It represents one and only one element of F n (B) and can be expressed as:

F(α 1 , . . . , α n ) = F 0 (α 2 , . . . , α n ) • α 1 + F 1 (α 2 , . . . , α n ) • α 1
where F 0 (α 2 , . . . , α n ) and F 1 (α 2 , . . . , α n ) are Boolean formulae of α 2 , . . . , α n . These two formulae can be directly obtained from F(α 1 , . . . , α n ) as follows:

F 0 (α 2 , . . . , α n ) = F(α 1 , . . . , α n )| α1←0 = F(0, α 2 , . . . , α n ) F 1 (α 2 , . . . , α n ) = F(α 1 , . . . , α n )| α1←1 = F(1, α 2 , . . . , α n )
It is also possible to define a partial order relation, called

Inclusion-Relation, between elements of F n (B). Definition 1. (Inclusion-Relation). 4 ∀(x, y) ∈ F n (B) 2 , define x ≤ y if and only if x • y = x.
Relation Inclusion is a partial order: it is reflexive (x ≤ x), antisymmetric (if x ≤ y and y ≤ x, then x = y) and transitive (if x ≤ y and y ≤ z, then x ≤ z). Since x • y = x ⇔ x • y = 0, we also have x ≤ y ⇔ x • y = 0.

Solutions of Boolean equations over F n (B)

Consider the Boolean algebra of n-variable switching functions (F n (B), +, •, , 0, 1).

• Let (f 1 Pj , • • • , f n Pj ) be the projection-functions of F n (B). • Let (x 1 , • • • , x k ) be k elements of F n (B) considered
as unknowns.

For notational convenience, we note 'X k ' the vector

(x 1 , • • • , x k ) of the k unknowns and 'Pj' the vector (f 1 Pj , • • • , f n Pj ) of the n projection-functions of F n (B)
. Any set of simultaneously-asserted relations of n-variable switching functions can be reduced to a single equivalent relation as: F(X k , Pj) = 0. For that, it is sufficient:

• to rewrite each equality according to:

F 1 (X k , Pj) = F 2 (X k , Pj) ⇔ F 1 (X k , Pj) • F 2 (X k , Pj) + F 1 (X k , Pj) • F 2 (X k , Pj) = 0 • to rewrite each inclusion according to: F 1 (X k , Pj) ≤ F 2 (X k , Pj) ⇔ F 1 (X k , Pj) • F 2 (X k , Pj) = 0
• to group rewritten equalities as follows:

F 1 (X k , Pj) = 0 F 2 (X k , Pj) = 0 ⇔ F 1 (X k , Pj) + F 2 (X k , Pj) = 0
In order to express Eq(X k , Pj) = 0 in a canonic form, the following notation has to be introduced: for x ∈ F n (B) and a ∈ {0, 1}, x a is defined by

x 0 = x x 1 =
x This notation can be used for vectors as follows: for

X k = (x 1 , • • • , x k ) ∈ F n (B) k and A k = (a 1 , • • • , a k ) ∈ {0, 1} k , X A k k is defined by X A k k = i=k i=1 x ai i = x a1 1 • . . . • x a k k
Theorem 1. Canonic form of a Boolean equation Any Boolean equation Eq(X k , Pj) = 0 can be expressed in the canonic form

A k ∈{0,1} k Eq(A k , Pj) • X A k k = 0
where Eq(A k ,Pj) (with A k ∈ {0, 1} k ) are the 2 k discriminants5 of Eq(X k , Pj) according to X k .

This canonic form is obtained by expressing Eq(X k , Pj) according to the k unknowns. For example, we have: Eq(x 1 , x 2 , Pj) = Eq(0, 0, Pj)

• x 1 • x 2 + Eq(0, 1, Pj) • x 1 • x 2 + Eq(1, 0, Pj) • x 1 • x 2 + Eq(1, 1, Pj) • x 1 • x 2 Theorem 2. Solution of a k-unknown equation The Boolean equation over F n (B)
Eq 0 (X k , Pj) = 0 (1) is consistent (i.e. has at least one solution) if and only if the following condition is satisfied:

A k ∈{0,1} k Eq 0 (A k , Pj) = 0 (2)
When ( 2) is satisfied, Equation (1) admits one or more k-tuple solutions (S(x 1 ), • • • , S(x k )) where each S(x i ) is defined by

S(x i ) = A k-i ∈{0,1} k-i Eq i-1 (0, A k-i , Pj) + p i • A k-i ∈{0,1} k-i Eq i-1 (1, A k-i , Pj) (3)
where

• Eq i (x i+1 , • • • , x k , Pj) = Eq i-1 (x i , • • • , x k , Pj)| xi←S(xi)
• p i is an arbitrary parameter, i.e, a freely-chosen element of F n (B).

The general form of this theorem can be found in [START_REF] Rudeanu | Lattice Functions and Equations[END_REF] or [START_REF] Brown | Boolean Reasoning: The Logic of Boolean Equations[END_REF]. However, in these works, the solution is not expressed in a parametric form, but only with intervals. The proof of this parametric formulation, which is mandatory in our approach, can be found in [START_REF] Hietter | Synthèse algébrique de lois de commande pour les systèmes à évènements discrets logiques[END_REF].

Two new theorems for coping with inconsistencies of specifications

In practice, it is very difficult for a designer to specify the whole requirements of a complex system without inconsistencies. It is the reason why requirements given by the designer are often declared as inconsistent according to Theorem 2. Since the inconsistency condition is a Boolean

• Several requirements cannot be simultaneously respected. In this case, a hierarchy between requirements can be proposed in order to find a solution. The requirements which have the lower priority have to be corrected for becoming consistent with the requirements which have the higher priority. This strategy is based on Theorem 4. • The detected inconsistency refers to combinations of the projection-functions for which the designer knows that there is no solution. To avoid blocking the synthesis process, it is necessary to introduce new assumptions and to use Theorem 3. Theorem 3. Solution of a Boolean equation according to an assumption among the projection-functions.

The following problem    

Equation to solve: Eq 0 (X k , Pj) = 0 Assumptions:

A(Pj) = 0 (4)
admits the same solutions as the following equation:

Eq 0 (X k , Pj) ≤ A(Pj) (5)
Proof. According to A(Pj) = 0, Eq 0 (X k , Pj) = 0 can be rewritten as:

Eq 0 (X k , Pj) = 0 A(Pj) = 0 ⇔ A(Pj) + Eq 0 (X k , Pj) = 0 ⇔ A(Pj) + A(Pj) • Eq 0 (X k , Pj) = 0 ⇔ A(Pj) • Eq 0 (X k , Pj) = 0 A(Pj) = 0 ⇔ Eq 0 (X k , Pj) ≤ A(Pj) A(Pj) = 0
Equation A(Pj) • Eq 0 (X k , Pj) = 0 is consistent if and only if the following condition is true (Theorem 2):

A(Pj) • A k ∈{0,1} k Eq 0 (A k , Pj) = 0 (6)
By construction, this new condition is the subset of the initial condition ( A k ∈{0,1} k Eq 0 (A k , Pj) = 0) for which the proposed assumption is satisfied. All the others terms have been removed.

If ( 6) is satisfied, (5) admits one or more k-tuple solutions where each component S(x i ) is defined by

S(x i ) =A(Pj) •   A k-i ∈{0,1} k-i Eq i-1 (0, A k-i , Pj) +p i • A k-i ∈{0,1} k-i Eq i-1 (1, A k-i , Pj)   + A(Pj) • p i As A(Pj) = 0, S(x i
) can also be expressed as:

S(x i ) = A k-i ∈{0,1} k-i Eq i-1 (0, A k-i , Pj) + p i • A k-i ∈{0,1} k-i Eq i-1 (1, A k-i , Pj)
When A(Pj) = 0 is satisfied, the solutions of (5) are also solution of Equation Eq 0 (X k , Pj) = 0. Theorem 4. Solution of a Boolean equation system according to a priority rule between requirements. The following problem 

       
Equations system to solve: where:

   HR F H (X k , Pj) = 0 LR F L (X k , Pj) = 0 OR F O (X k , Pj) =
• F H (X k , Pj) = 0 is the formal expression of the requirements which have the higher priority (HR). • F L (X k , Pj) = 0 is the formal expression of the requirements which have the lower priority (LR). • F O (X k , Pj) = 0 is the formal expression of the others requirements (OR). • HR LR is the priority rule between inconsistent requirements. admits the same solutions as the system of equations:

   F H (X k , Pj) = 0 F L (X k , Pj) ≤ I(Pj) F O (X k , Pj) = 0 (8)
where I(Pj) is the inconsistency condition between requirements 'HR' and 'LR':

I(Pj) = A k ∈{0,1} k (F H (A k , Pj) + F L (A k , Pj))
Proof. Thanks to Theorem 2, the inconsistency condition I(Pj) between requirements 'HR' and 'LR' can be found by solving equation: F H (X k , Pj) + F L (X k , Pj) = 0. We have:

I(Pj) = A k ∈{0,1} k (F H (A k , Pj) + F L (A k , Pj))
To remove the inconsistency between Requirements 'HR' and 'LR' according to the priority rule 'HR LR', it is necessary to restrict the range of requirement 'LR' to the part for which there is no inconsistency, i.e. I(Pj) = 0. That is the case, when F L (X k , Pj) = 0 replaced by F L (X k , Pj) ≤ I(Pj).

Thanks to Theorem 2, (9) admits always one or more ktuple solutions and it is impossible to find a less restrictive condition over Requirement 'LR'.

F H (X k , Pj) = 0 F L (X k , Pj) ≤ I(Pj) (9) 2

SYNTHESIS OF CONTROLLERS DESPITE INCOHERENT REQUIREMENTS

Control system specifications

Let us consider a hydraulic press with a vertical ram (Fig 3). A safety-light curtain is used to safeguard operators during the movements of the ram. A control panel allows to select the mode of operation: Manual or Automatic mode.

• In Manual mode, all the operations are carried out by pressing the corresponding push-buttons. As soon as a push-button is released, the ram movement stops. • In Automatic mode, the cycle starts by pressing the 'Start' push-button: the ram is going down and comes back to the up position after the press operation has been done.

Fig. 3. A hydraulic press and its human-machine interfaces

Inputs and outputs of the controller (Fig. 4a) Each movement of the ram is controlled thanks to a Boolean output ('goUp' and 'goDown'). The controller is informed of the position of the ram thanks to inputs 'up' and 'down'. The safety light curtain is connected to input 'lc' (lc = b 1 when the operator is not in the detection zone of the light curtain). The control panel of the press is composed of an emergency stop button (input: 'esb'), a three position center-off switch for the operating mode selection ('sbA', off: no mode selected, 'sbM') and four push-buttons (inputs: 'vpb', 'spb', 'uppb' and 'dopb').

Control laws to synthetize It is not possible to identify automatically how many and which state variables must be used. The designer has to fix the state variables by expertise.

For this case study, we propose to use 5 state variables: one for each output; one for each mode of operation (Automatic, Manual)6 and one for characterizing a state where the press is in a failure mode (Fail). According to this choice, 5 15-variable switching functions (Auto, Manual, Fail, GoUp and GoDown) have to be synthesized (Fig. 4b). The 15 projection-functions of F 15 (B) are therefore:

• The 10 switching functions (Up, Down, Lc, Esb, SbA, SbM, Vpb, Spb, Uppb, and Dopb) which characterize the behavior of the inputs of the controller. They are defined as follows: Up :

B 15 → B (up[k], . . . , goDown[k -1]) → up[k]
• The 5 switching functions ( p Auto, p Manual, p Fail, p GoUp and p GoDown) which characterize the previous behavior of the state variables of the controller. They are defined as follows:

p Auto :

B 15 → B (up[k], . . . , goDown[k -1]) → auto[k -1]

Algebraic formalization of requirements

The complete formalization of the behavior of the hydraulic press is given Fig. 4c. In order to illustrate the power of expression of relations Equality and Inclusion,

             auto[k] = Auto(up[k], . . . , dopb[k], auto[k -1], manual[k -1], fail[k -1], goUp[k -1], goDown[k -1]) manual[k] = Manual(up[k], . . . , dopb[k], auto[k -1], manual[k -1], fail[k -1], goUp[k -1], goDown[k -1]) fail[k] = Fail(up[k], . . . , dopb[k], auto[k -1], manual[k -1], fail[k -1], goUp[k -1], goDown[k -1]) goUp[k] = GoUp(up[k], . . . , dopb[k], auto[k -1], manual[k -1], fail[k -1], goUp[k -1], goDown[k -1]) goDown[k] = GoDown(up[k], . . . , dopb[k], auto[k -1], manual[k -1], fail[k -1], goUp[k -1], goDown[k -1]) c) Complete formal specification                                                                   Requirements:                                            R1 Auto • Manual + Auto • Fail + Manual • Fail = 0 R2 Esb ≤ Fail R3 Up • Down ≤ Fail R4 Fail • p Fail ≤ Lc R5 Auto = SbA R6 Auto • p Auto ≤ Up • Lc R7 Auto • p Auto ≤ Vpb R8 Auto • Lc ≤ Up R9 Auto • p Auto • (Lc + Up) ≤ Fail R10 Manual = SbM R11 Manual • p Manual ≤ Lc R12 Manual • p Manual ≤ Vpb                                                  R20 GoUp • GoDown = 0 R21 Up ≤ GoUp R22 Down ≤ GoDown R23 (GoUp + GoDown) ≤ Lc R24 (GoUp + GoDown) ≤ (Auto + Manual) R25 Manual ≤ ((GoUp • Uppb) + (GoUp • Uppb)) R26 Manual ≤ ((GoDown • Dopb) + (GoDown • Dopb)) R27 Manual • (Uppb • Dopb) ≤ (GoUp + GoDown) R28 Auto • (GoUp • p GoUp) ≤ Up R29 Auto • (GoDown • p GoDown) ≤ Down R30 Auto ≤ (Up + GoUp + GoDown) R31 Auto ≤ ((GoDown • p GoDown) • (Up • Spb) +(GoDown • p GoDown) • ( p GoDown • Up • Spb)) Priority rules:                                    {R1,R2,R3,R4} R5
                Fail = Esb + Up • Down + Lc • p Fail + Up • Lc • p Auto Auto = Esb • SbA • (Up • Down • Lc • Vpb + p Auto • (Up • Down + Up • Lc)) Manual = Esb • SbM • (Up • Down) • (Lc • Vpb + p Manual) GoUp = Esb • Up • Lc • (SbA • p Auto • (Down + p GoDown) + SbM • Uppb • Dopb • (Vpb + p Manual)) GoDown = Esb • Down • Lc • (SbA • (Spb • Up • (Vpb + p Auto) + p GoDown • ( p Auto + (Vpb • Up))) +SbM • Dopb • Uppb • (Vpb + p Manual))
e) Control laws of the hydraulic press

                 fail[k] = esb[k] ∨ up[k] ∧ down[k] ∨ ¬lc[k] ∧ fail[k -1] ∨ ¬up[k] ∧ ¬lc[k] ∧ auto[k -1] auto[k] = ¬esb[k] ∧ sbA[k] ∧ (up[k] ∧ ¬down[k] ∧ lc[k] ∧ vpb[k] ∨ auto[k -1] ∧ (up[k] ∧ ¬down[k] ∨ ¬up[k] ∧ lc[k])) manual[k] = ¬esb[k] ∧ sbM[k] ∧ ¬(up[k] ∧ down[k]) ∧ (lc[k] ∧ vpb[k] ∨ manual[k -1]) goUp[k] = ¬esb[k] ∧ ¬up[k] ∧ lc[k] ∧ (sbA[k] ∧ auto[k -1] ∧ (down[k] ∨ ¬goDown[k -1]) ∨ sbM[k] ∧ uppb[k] ∧ ¬dopb[k] ∧ (vpb[k] ∨ manual[k -1])) goDown[k] = ¬esb[k] ∧ ¬down[k] ∧ lc[k] ∧ (sbA[k] ∧ (spb[k] ∧ up[k] ∧ (vpb[k] ∨ auto[k -1]) ∨ goDown[k -1] ∧ (auto[k -1] ∨ (vpb[k] ∧ up[k]))) ∨sbM[k] ∧ dopb[k] ∧ ¬uppb[k] ∧ (vpb[k] ∨ manual[k -1]))
Fig. 4. Details of the case study all the requirements used for the synthesis of the operation modes are expressed in textual form hereafter:

• R1 The three modes (Automatic, Manual and Fail) are exclusive. • R2 While the Emergency stop button 'esb' is pressed, the press is in Failure mode. • R3 If the observed position of the ram is both 'up' and 'down', the press is in Failure mode. • R4 For leaving the Failure mode, the operator must not be in the detection zone of the light curtain. • R5 The press is in Automatic mode if and only if the three position center-off switch is turned on 'sbA' position. • R6 For reaching the Automatic mode, the press ram must be in 'up' position and the operator must not be in the detection zone of the light curtain. • R7 For reaching the Automatic mode, the 'vpb' pushbutton must be pressed. • R11 For reaching the Manual mode, the press ram must be in 'up' position and the operator must not be in the detection zone of the light curtain. • R12 For reaching the Manual mode, the 'vpb' pushbutton must be pressed.

The relation "Inclusion" is the more used in our approach since it allows the expression of both sufficient (as R2) and necessary conditions (as R4).

Synthesis process

In traditional design methods, the design procedure of a logic controller is not a linear process, but an iterative one converging to an acceptable solution. At the beginning of the design, requirements are neither complete nor without errors. Most often, new requirements are added during the search of solutions, others are corrected. This complementary information is given by the designer after analysis of the partial solutions he found or when inconsistencies have been detected. If we do not make the hypothesis that the specifications are complete and consistent, designing a controller with a synthesis technique will also be an iterative process in which the designer plays an important role.

Synthesis of the operation modes For this case study, we have started with the synthesis of operation modes. The first four studied requirements have been: R1, R2, R5 and R10. For this subset of requirements, the result given by our software tool was the following inconsistency condition: I 0 = SbA • SbM + Esb • SbA + Esb • SbM Since requirements are declared as inconsistent, we have to give complementary information to precise our specifica-tion. By analyzing each term of this formula, it is possible to detect the origin of the inconsistency:

• SbA•SbM: What happens if the Automatic mode and the Manual mode are simultaneously selected? We consider that it is not possible (due to the technology of the three position center-off switch) and we have added Assumption A1. With this complementary information, Problem (10) admits as parametric solution (11):

                     Requirements:              R1 Auto • Manual + Auto • Fail +Manual • Fail = 0 R2 Esb ≤ Fail R5 Auto = SbA R10 Manual = SbM Priority rules: {R1,R2} R5 {R1,R2} R10 Assumptions: A1 SbA • SbM = 0 (10)    Fail = Esb + p 1 • SbA • SbA p 1 ∈ F 15 (B) Auto = Esb • SbA Manual = Esb • SbM (11)
By gradually adding all the requirements and selecting the solution which minimizes the Fail mode (the parameter p 1 is fixed to 0) the solution we obtain is proposed on Fig. 4d (the first three equations). In an illustrative purpose, a state model representation of the synthesized operation modes management, automatically built thanks to [START_REF] Guignard | Symbolic generation of the automaton representing an algebraic description of a logic system[END_REF], is given in (Fig. 5). The transition conditions are non-trivial Boolean expressions of inputs that take into account the whole set of specifications given in Section 5.2.

Synthesis of the reactive control laws

To obtain the complete control law of the hydraulic press, the solution previously obtained for the operation modes has to be completed thanks to requirements R20 to R31. All the priorities rules and assumptions used for this case study are given on Fig. 4c. The solution we obtain is proposed on Fig. 4d.

The control laws presented Fig. 4e have been obtained by translating the expression of the five unknowns according to the projection-functions into relations between recurrent Boolean equations. These control laws can be automatically translated in the syntax of the Ladder Diagram language (IEC 61131-3 (2003)) before being implemented 1 : ∅ 2 : Manual 3 : Auto 4 : Fail

                                                   E 1-2 = ¬esb ∧

CONCLUSION

Many research works in the field of DES aim at formalizing steps of the systems life cycle. Since 20 years, significant progresses have been obtained for the synthesis, verification, performance evaluation, diagnosis ... of DESs. Nevertheless, one of the common difficulties of these works is that there comes a time when it is necessary to pass from the informal expression of the knowledge of a system to its formalization. Few works have paid attention to this important task which is very error prone. In this paper, we proposed an iterative process that allows coping with inconsistencies of the requirements during the synthesis of the controller. The framework in which we proposed this approach is an algebraic synthesis method. Since the problem is located to the frontier between formal and informal, intervention of the designer is necessary. Nevertheless, we have shown that this intervention can be guided by the results the formal method provides.

Fig. 2 .

 2 Fig. 2. The algebraic synthesis method step by step

(

  Press ram in Up position) up (Press ram in Down position) down (Safety light curtain not broken) lc (Emergency stop button) esb (Three position switch in Automatic Position) sbA (Three position switch in Manual Position) sbM ('Validation' Push-button) vpb ('Start' Push-button) spb ('Up' Push-button) uppb ('Down' Push-button) dopb goUp (To move up the press ram) goDown (To move up the press ram) b) General form of the expected control laws

  (* Fail mode has priority on Automatic mode. *) {R6,R7,R8} R5 (* Starting rules have priority on the switch position rule. *) {R1,R2,R3,R4,R8} R10 (* Fail mode has priority on Manual mode. *) {R11,R12} R10 (* Starting rules have priority on the switch position rule. *) {R21,R23,R27} R25 (* Safety requirements have priority on functional requirements. *) {R21,R23,R27} R26 (* Safety requirements have priority on functional requirements. *) R23 R28 (* Safety requirements have priority on functional requirements. *) R23 R29 (* Safety requirements have priority on functional requirements. *) R23 R30 (* Safety requirements have priority on functional requirements. *) {R22,R23} R31 (* Safety requirements have priority on functional requirements. *) Assumptions:    A1 SbA • SbM = 0 (* The 2 positions of the switch button are exclusive. *) A2 p Auto • p Manual + p Auto • p Fail + p Manual • p Fail = 0 (* Consequence of R1. *) A3 p GoUp • p GoDown = 0 (* Consequence of R20. *) Optimization criteria: Minimisation of: Fail d) Solution obtained by symbolic calculation 

  lc ∧ sbM ∧ vpb ∧ (¬up ∨ ¬down) E 1-3 = ¬esb ∧ lc ∧ sbA ∧ vpb ∧ up ∧ ¬down E 1-4 = esb ∨ up ∧ down E 2-1 = ¬esb ∧ ¬sbM ∧ (¬up ∨ ¬down ∧ ¬(lc ∧ sbA ∧ vpb)) E 2-3 = ¬esb ∧ lc ∧ sbA ∧ vpb ∧ up ∧ ¬down E 2-4 = esb ∨ up ∧ down E 3-1 = ¬esb ∧ ¬sbA ∧ (¬lc ∧ up ∧ ¬down∨ (¬sbM ∨ ¬vpb) ∧ (up ∧ ¬down ∨ lc ∧ ¬up)) E 3-2 = ¬esb ∧ lc ∧ sbM ∧ vpb ∧ (¬up ∨ ¬down) E 3-4 = esb ∨ up ∧ down ∨ ¬up ∧ ¬lc E 4-1 = ¬esb ∧ lc ∧ (¬up ∨ ¬down)∧ (¬sbA ∧ ¬sbM ∨ ¬vpb ∨ sbA ∧ ¬up) E 4-2 = ¬esb ∧ lc ∧ sbM ∧ vpb ∧ (¬up ∨ ¬down) E 4-3 = ¬esb ∧ lc ∧ sbA ∧ vpb ∧ up ∧ ¬down

Fig. 5 .

 5 Fig. 5. State model of modes of operation into a PLC. The complete code is composed of only 10 rungs.

  • R8 During the Automatic mode, the operator can be in the detection zone of the light curtain without to be in danger only if the press ram is in 'up' position. • R9 During the Automatic mode, if the operator is detected by the light curtain while the press ram is not in 'up' position, one has to switch into the Failure mode. • R10 The press is in Manual mode if and only if the three position center-off switch is turned on 'sbM' position.

Case studies are available online: http://www.lurpa.enscachan.fr/isa/asc/case studies.html

The term of 'discriminant' comes from[START_REF] Brown | Boolean Reasoning: The Logic of Boolean Equations[END_REF].formula, it is possible to use it for the detection of the origin of inconsistencies. Two cases have to be considered:

These mode variables must not be confused with inputs 'sbA' or 'sbM' which are the demand of operator for reaching these modes.


Marin Bougeret

Pierre-Franc ¸ois Dutot

Klaus Jansen

Christina Robenek

Denis Trystram

Tight approximation for scheduling parallel jobs on identical clusters

We consider the Multiple Cluster Scheduling Problem (M CSP ), where the objective is to schedule n parallel rigid jobs on N identical clusters, minimizing the maximum completion time (makespan). M CSP is 2-inapproximable (unless P = N P), and several approximation algorithms have already been proposed. However, ratio 2 has only been reached by algorithms that use extremely costly and complex subroutines as "black boxes" which are polynomial and yet impractical due to prohibitive constants.

Our objective within this work is to determine a reasonable restriction of M CSP where the inapproximability lower bound could be tightened in almost linear time. Thus, we consider a restriction of M CSP where jobs do not require strictly more than half of the processors of a cluster, and we provide a 2-approximation running in O(log(nh max )n(N + log(n))), where h max is the processing time of the longest job. This approximation is the best possible, as this restriction (and even simpler ones) remains 2-inapproximable.

Introduction

Problem statement

In new parallel computing platforms, several clusters share their computing resources in order to distribute the workload. Each cluster is composed of a set of identical processors connected by a local interconnection network. Jobs are submitted in successive packets called batches. The objective is to minimize the time when all the jobs of a batch are completed, to start the following batch as soon as possible. Many such computational systems are available all over the world, and the efficient management of the resources is a crucial problem. Let us start by defining formally the Multiple Cluster Scheduling Problem (M CSP ). Definition 1 (M CSP ) We are given n parallel rigid jobs J j , 1 ≤ j ≤ n, and N clusters. A job J j requires q j processors during p j units of time, and each cluster owns m identical processors. The objective is to schedule all the jobs in the clusters, minimizing the maximum completion time (makespan) under the following constraints:

1. the q j processors allocated to job J j must belong to the same cluster (as proposed in [START_REF] Dutot | Scheduling on hierarchical clusters using malleable tasks[END_REF])

2. at any time, the total number of processors used in any cluster must be lower or equal to m M CSP is closely related to the Multiple Strip Packing Problem (M SP P ) where the objective is to pack n rectangles in N strips, minimizing the maximum reached height. The only difference between these two problems is the "contiguous" constraint. Indeed, in M SP P rectangles must be allocated contiguously, which in terms of job scheduling amounts to force jobs to use consecutive indexes of processors (see Figure 1).

Of course, the results for M CSP generally do not apply to M SP P , as algorithms may schedule jobs in a non con-tiguous way. The converse is also not clear, as ratio of approximation algorithms for M SP P may not be preserved when considering M CSP (optimal value of an M CSP instance may be strictly better than the corresponding one for M SP P ). However, as we can notice in Figure 2, many results for M SP P directly apply to M CSP , as the proposed algorithms build contiguous schedules that are compared to non-contiguous optimal solutions.

We refer to the packing context in this paper, so that the solutions can be described using the classical vocabulary of packing problems. Thus, the problem treated in this paper (M CSP ) consists in packing n rectangles in N strips of width 1, minimizing the maximum reached height, and without contiguous constraints 1 .

Related Work

As shown in [START_REF] Zhuk | Approximate algorithms to pack rectangles into several strips[END_REF] using a gap reduction from 2 partition, M CSP (and M SP P ) are 2-inapproximable in polynomial time unless P = N P, even for N = 2. The main positive results for M CSP are summarized below, in Figure 2

Fast algorithm this paper

Figure 2. Main results

We distinguish the 3-approximation [START_REF] Schwiegelshohn | Online scheduling in grids[END_REF] and the 5 2approximation [START_REF] Bougeret | Approximating the non-contiguous multiple organization packing problem[END_REF] that have a low computational complexity (these algorithms are usable on real size instances) from the 2-approximation [START_REF] Bougeret | Approximation algorithm for multiple strip packing[END_REF] and the 2 + ǫ-approximation [START_REF] Ye | On-Line Multiple-Strip Packing[END_REF].

At a first sight, the best results seem to be the 2approximation in [START_REF] Bougeret | Approximation algorithm for multiple strip packing[END_REF] and the 2 + ǫ-approximation in [START_REF] Ye | On-Line Multiple-Strip Packing[END_REF]. However, the 2-approximation requires using a high running time algorithm when the number of clusters is lower than a huge constant N 0 . Thus, any exponential dependency in N 0 is hidden, and the value of this constant (N 0 ≈ 10 4 ) makes this algorithm practically impossible to use. Moreover, the 2 + ǫ-approximation requires solving the famous P ||C max problem (which is makespan minimization when scheduling sequential jobs on identical machines) with a ratio 1 + ǫ 2 . Thus, to give a rough idea, applying this technique with ǫ = 1 3 would lead to a Ω(n 36 ) algorithm, using the PTAS of [START_REF] Hochbaum | Using dual approximation algorithms for scheduling problems theoretical and practical results[END_REF]. 2

Motivations and contributions

As explained above, the two previous results of [START_REF] Bougeret | Approximation algorithm for multiple strip packing[END_REF][START_REF] Ye | On-Line Multiple-Strip Packing[END_REF] require both extremely high running time, and do not provide any insight on M CSP because of these black boxes subroutines. Thus, we followed in [START_REF] Bougeret | Approximating the non-contiguous multiple organization packing problem[END_REF] another approach by looking for fast and direct algorithms for M CSP . Our previous 5 2 -approximation in [START_REF] Bougeret | Approximating the non-contiguous multiple organization packing problem[END_REF] is based on the discarding technique presented in Section 2.2. What we call discarding technique is a classical framework in scheduling problems. The idea is to define properly a set of negligeable items (items are rectangles here), and to prove that it is possible to add these items only at the end of the algorithm without degrading the approximation ratio. Thus, the effort can be focused on the set I ′ of remaining large items, that are generally more structured.

The 5 2 -approximation was obtained through a basic application of this technique (i.e. with a set I ′ containing only really huge rectangles, in this case rectangles whose width is larger than 1 2 ). As we believe that the discarding technique of Section 2.2 is well suited for M CSP , we apply it again using a more challenging set I ′ . A natural direction would be to improve the 5 2 ratio for M CSP by targeting a fixed ratio ρ < 5 2 . Typically, one could target ρ = 7 3 by defining the small jobs as those whose length is lower than 1 3 (instead of 1 2 ). However, as the relative performance improvement is getting smaller, and the difficulty of these "ratio tailored" proofs is likely to increase rapidly, we consider here a different approach.

Our objective is to find a reasonable restriction of M CSP where the inapproximability lower bound could be tightened in almost linear time. In this spirit, we study a restriction of M CSP where all rectangles have a width lower than 1 2 , meaning that jobs submitted to the clusters do not require strictly more than half of the processors. We provide for this problem a very fast 2-approximation running in O(log(nh max )n(N + log(n))), where h max is the maximum height of any rectangle. It turns out that this result is the best possible approximation, as this restriction of M CSP (and even simpler ones, where the width of rectangles is lower than 1 c , c ∈ N, c ≥ 2) remains 2inapproximable unless P = N P.

General principles

In this section, we generalize the framework used in the 5 2 -approximation of [START_REF] Bougeret | Approximating the non-contiguous multiple organization packing problem[END_REF]. This framework will be applied in Section 3 to get the 2-approximation.

Preliminaries

Recall that our objective is to (non contiguously) pack a set I of n rectangles r j into N strips of width 1. Rectangle r j has a height h j and a width w j . We denote by s(r j ) = w j h j the surface of r j . These notations are extended to W (X), H(X) and S(X) (where X is a set of rectangles), which denote the sum of the widths (resp. heights, surfaces) of rectangles in X.

A layer is a set of rectangles packed one on top of the other in the same strip (as depicted Figure 3). The height of layer Lay is H(Lay), the sum of the height of all the rectangles in Lay. A shelf is a set of rectangles that are packed in the same strip, such as the bottom level of all the rectangles is the same. Even if it is not relevant for the noncontiguous case, we consider for the sake of simplicity that in a shelf, the right side of any rectangle (except the right most one) is adjacent to the left side of the next rectangle in the shelf. Given a shelf sh (sh denotes the set of rectangles in the shelf), the value W (sh) is called the width of sh.

Packing a shelf at level l means that all the rectangles of the shelf have their bottom at level l. A bin is a rectangular area that can be seen as a reserved space in a particular strip for packing rectangles. As the width of a bin is always 1, we define a bin by giving its height h b , its bottom level l b and the index i b of the strip it belongs to. Packing a shelf sh in a bin b means that sh is packed in strip S i b at level l b . Moreover we always guarantee that the height of any rectangle of sh is lower than h b .

The utilization u π i (l) of a packing π in strip S i at level l (sometimes simply denoted by u(l) or u i (l)) is the sum of the width of all the rectangles packed in S i that cut the horizontal line-level l (see Figure 3). Of course we have 0 ≤ u π i (l) ≤ 1 for any l and i. Let us now describe three useful basic procedures. The CreateLayer(X, h) procedure creates a layer Lay (using rectangles of X) of height at most h, using a Best Fit (according to the height) policy (BFH). Thus, CreateLayer(X, h) adds at each step the highest rectangle that fits. Of course, the layer produced by the procedure is such that H(Lay) ≤ h. Moreover, notice that we will always pack the layers in the strips with the narrowest rectangles on the top. CreateShelf (X, w) creates a shelf sh (using rectangles of X) of width at most w, using the Best Fit (according to the width) policy (BFW). Thus, CreateShelf (X, w) adds at each step the widest rectangle that fits. Of course, the shelf produced by the procedure is such that W (sh) ≤ w. Throughout the paper, we consider that the sets of jobs used as parameters in the algorithms are modified after the calls.

Let us now state a standard lemma about the efficiency of the "best fit" policies.

Lemma 2 Let Sh denote the shelf created by

CreateShelf (X, w). If the k widest rectangles of X are added to Sh, then W (Sh) > k k+1 w.

Proof Let x be the cardinality of X. Let us assume that

w i ≥ w i+1 for 1 ≤ i < x. Let i 0 ≥ k + 1 be the first index such that r i0 is not in Sh. Let a = Σ i0-1 i=1 w i . We have W (Sh) ≥ a ≥ (i 0 -1)w i0 > (i 0 -1)(w -a) leading to a > i0-1 i0 w ≥ k k+1 w.

Discarding technique

How to pack all rectangles in three steps

Discarding techniques are common for solving packing and scheduling problem. As mentioned before, the idea is to define properly a set of small items (rectangles here), and to prove that adding these small items only at the end of the algorithm will not degrade the approximation ratio. Thus, the effort can be focused on the remaining large items. In this section we present an adaptation of this general technique to the context of non-contiguous multiple strip packing. Thus, we will define a set of big rectangles I ′ ⊂ I, and the larger the set I ′ , the better the approximation ratio (as the remaining small rectangles become really negligible). In order to partition rectangles according to their height, we use the well-known dual approximation technique [START_REF] Hochbaum | A polynomial approximation scheme for scheduling on uniform processors: Using the dual approximation approach[END_REF].

We denote by v the guess of the optimal value. Given an instance I, let L W D = {w j > α} be the set of wide rectangles, L H = {h j > βv} be the set of high rectangles, and I ′ = L W D ∪ L H be the set of the big rectangles, with 0 < α < 1 and 0 < β < 1. Let r(α, β) = 1 1-α + β be the approximation ratio we target (the origin of this formula will be explained in Section 2.2.2). We also need the following definition. Definition 3 A packing is x-compact if and only if for every strip S i there exists a level l i such that for all l ≤ l i , u i (l) > x and u i restricted to l > l i is non-increasing. Figure 4 provides an example showing a (1-α) compact packing, and why step c) is simple. Indeed, adding as soon as possible a small rectangle r j (h j ≤ βv and w j ≤ α) to a (1α) compact packing cannot exceed v( 11-α + β). Let us now describe the three main steps of our approach. Notice that what we call a preallocation is a normal packing (i.e. that defines the bottom level of each rectangle, which is sufficient) that is based on simple structures like shelves and layers. We will prove that to get a r(α, β) = ( 11-α + β) ratio, it is sufficient to: a) construct a preallocation π 0 of I ′ that fits in r(α, β)v, and such that rectangles of L W D ⊂ I ′ are already packed in a (1α)-compact way b) turn π 0 into a (1-α)-compact packing π 1 by repacking rectangles of I ′ \ L W D using the list algorithm LS π0 of Lemma 4 c) add the small remaining rectangles (I \ I ′ ) using algorithm LS (see Lemma 5)

l 3 = l 4 l 1 l 2 α ≤ α ≤ βv r j S 1 S 2 S 3 S 4
Step a) is the most difficult one. Thus, Section 3 is entirely devoted to the construction of π 0 (targeting α = 1 3 and β = 1 2 ). Of course, building the preallocation becomes harder when α and β are small, as the number of rectangles of I ′ increases and r(α, β) decreases. Roughly speaking, the simple shapes of rectangles of I ′ allows us to construct π 0 with a simple structure. We will denote by π i 0 the set of rectangles packed by π 0 in S i .

Proving steps b) and c)

We now prove that applying steps b) and c) leads to a r(α, β) ratio. In this section, we suppose that we are given a guess v, and a packing π 0 (called the preallocation) of I ′ = L W D ∪ L H that fits in r(α, β)v, and such that rectangles of L W D ⊂ I ′ are already packed in a (1α)-compact way. We consider step b): how to turn π 0 into a (1α)compact packing.

Lemma 4 (Step b)) Let π 0 be the preallocation of

I ′ con- structed in Step a). Let π 1 = π 0 ∩ L W D denote π 0 when keeping only rectangles of L W D . Recall that π 1 is already a (1 -α)-compact packing of rectangles of L W D .
Then, we can complete π 1 into a (1α)-compact packing π 1 of I ′ , such that the height of π 1 is lower or equal to the height of π 0 .

Proof Let us define the LS π0 algorithm that adds rectangles of I ′ \ L W D . Let us consider a single strip S i . Let π i 0 denote π 0 restricted to S i , and π i 1 denote π 1 restricted to S i . Let X = {r 1 , . . . , r p } be the set of preallocated rectangles of I ′ \ L W D that we have to add to S i . We assume that lvl(j) ≤ lvl(j + 1), where lvl(j) is the bottom level of r j in π 0 .

For our considered strip S i , the LS π0 algorithm executes AddAsap(r j , π i 1 ), for 1 ≤ j ≤ p, where AddAsap(r, π i 1 ) adds rectangle r to π 1 (in S i ) at the smallest possible level. Notice first that adding with AddAsap a rectangle r j with w j ≤ α to a (1α)-compact packing creates another (1α)-compact packing. Thus it is clear that π 1 is (1α)compact.

For any 1 ≤ j ≤ p, let ( π i 1 , j) denote the packing in S i just before adding r j with AddAsap, and let (π i 0 , j) denote the packing π i 0 ∩ (L W D ∪ {r 1 , . . . , r j-1 }). Let us prove by induction on j ∈ {1, . . . , p} that u ( π i 1 ,j) (l) ≤ u (π i 0 ,j) (l), for any l ≥ lvl(j). The definition of π 1 implies the property for j = 1 (packings are identical). Let us suppose that the property holds for j, and prove it for j + 1. Let l ≥ lvl(j + 1). The induction property for rank j implies that r j is added by AddAsap at a level lower or equal to lvl(j). Thus, if r j intersects l in ( π i 1 , j + 1), then it also occurs in (π i 0 , j + 1). Thus in this case we have

u ( π i 1 ,j+1) (l) = u ( π i 1 ,j) (l) + w j ≤ u (π i 0 ,j) (l) + w j = u (π i 0 ,j+1) (l) If r j does not intersect l in ( π i 1 , j), then clearly u ( π i 1 ,j+1) (l) = u ( π i 1 ,j) (l) ≤ u (π i 0 ,j) (l) ≤ u (π i 0 ,j+1) (l)
Thus we proved that for any 1 ≤ j ≤ p we have u ( π i 1 ,j) (l) ≤ u (π i 0 ,j) (l) for any l ≥ lvl(j), implying that every r j is added by AddAsap at a level lower or equal to lvl(j). Thus, the height of π 1 is lower or equal to the height of π 0

We now prove in Lemma 5 that after adding rectangles in step c), the height of the packing does not exceed r(α, β)v = ( 11-α + β)v. This explains why the height of the pre-allocation should also be bounded by r(α, β)v.

Lemma 5 (Step c))

Let π 1 be a (1-α)-compact packing of I ′ . Adding to π 1 rectangles of I \ I ′ with a List Scheduling algorithm (LS) leads to a packing π having height lower than max(height(π 1 ), v( 11-α + β)).

Proof The LS algorithm scans all the strips from level 0, and at any level adds any rectangle of I \ I ′ that fits. Notice that the final packing π is (1α)-compact, since we add rectangles r j with w j ≤ α to an (1α)-compact packing.

Let us assume that the height of π is due to a rectangle r j ∈ I \ I ′ that starts at level s. This implies that when packing r j we had l i ≥ s for any strip i (with l i defined as in Definition 3). According to this definition we have u i (l) > 1α for any l ≤ l i . Thus, we have

S(I) > N i=1 l i (1 -α) ≥ N (1 -α)s, implying that s < v 1
1-α , and thus that of height of π is lower or equal to s + max j∈I\I ′ h j ≤ v( 11-α + β).

Thus, we now apply this framework with α = 1 3 and β = 1 2 to get a 2-approximation.

2-approximation

Hardness

As explained before, the 2 + ǫ-approximation in [10] and the 2-approximation we recently proposed in [START_REF] Bougeret | Approximation algorithm for multiple strip packing[END_REF] are rather complexity results than practical algorithms. We aim at constructing a low cost algorithm that could be used in a practical context. Thus, we are considering a restriction of M CSP where the inapproximability bound could be tightened with a fast algorithm, and we consider that all the rectangles have a width lower or equal to 1 /2.

Theorem 6

The M CSP where every rectangle has a width lower (or equal) to 1 2 has no polynomial algorithm with a ratio strictly better than 2, unless P = N P . Proof As in [START_REF] Zhuk | Approximate algorithms to pack rectangles into several strips[END_REF] for the general version, we construct a gap reduction from 2-partition. Let {x 1 , . . . , x n } ⊂ N n and a such that n i=1 = 2a. Without loss of generality, let us assume that for any i, x i < a. In order to only have items with size at least two, we define x ′ i = 2x i for any i, and a ′ = 2a. We construct the following instance I MSP : N = 2 strips, each of size 2a ′ -1. The set of rectangles is {r 1 , . . . , r n , r n+1 , r n+2 }, with w i = x ′ i for 1 ≤ i ≤ n, w n+1 = w n+2 = a ′ -1, and h i = 1 for 1 ≤ i ≤ n + 2. We have w i ≤ 2a ′ -1 2 for any i, as all the x i are strictly lower than a. Notice than any solution of I MSP that packs r n+1 and r n+2 is the same strip has a height of at least 2, as the available width of size 1 in that strip cannot be used by any rectangle.

Obviously, if there is a 2-partition, then Opt(I MSP ) = 1. Otherwise, as r n+1 and r n+2 cannot be packed together, we have Opt(I MSP ) = 2

The previous proof can easily be adapted for any non-trivial restriction on the size of the widest rectangle. Therefore, the fast 2-approximation presented in this section is the best possible result, even for more restricted versions of the MCSP.

Decomposition

We follow the ideas presented in Section 2, and thus we re-use the notion of layer, shelf, bin, and the procedures named CreateLayer and CreateShelf.

Again, we use the dual approximation technique [START_REF] Hochbaum | A polynomial approximation scheme for scheduling on uniform processors: Using the dual approximation approach[END_REF], and we denote by v the guess of the optimal value. Conforming to the dual approximation technique, we will prove that either we pack I with a resulting height lower than 2v, or v < Opt. Then, we will perform a binary search on v to turn the dual approximation algorithm into a classical approximation algorithm. Notice that for the sake of simplicity we did not add the "reject" instructions in the algorithm. Thus we consider in all the proof that v ≥ Opt, and it is implicit that if one of the claimed properties is wrong during the execution, the considered v should be rejected.

Recall that all rectangles have w j ≤ 1 2 . Let us define the following sets:

• let L W D = {r j |w j > 1 /3} be the set of wide rectangles

• let L XH = {r j |h j > 2v /3} be the set of very high rectangles

• let L H = {r j | 2v /3 ≥ h j > v /2}
be the set of high rectangles

• let L B = (L XH ∪ L H ) ∩ L W D
be the set of huge rectangles, and b = Card(L B ).

• let I ′ = L W D ∪ L XH ∪ L H
Notice that we only consider the values v such that

W (L XH ∪ L H ) ≤ N and H(L W D ) ≤ 2N v.
As expected, the set I ′ corresponds in our framework to the set of big rectangles for α = 1 3 and β = 1 2 . The construction of the preallocation π 0 of I ′ is presented from Section 3.3 to 3.5. The final steps to turn π 0 into a 2 3compact packing π 1 and to turn π 1 into the final packing π are quickly described in Section 3.6, as they follow the steps presented in Section 2.2.

We now provide a two phases algorithm that builds the preallocation π 0 of the rectangles of I ′ . Phase 1 (Section 3.3) preallocates rectangles of L W D , and phase 2 (Section 3.5) preallocates rectangles of L H ∪ L XH .

Phase 1

Phase 1 packs the rectangles of L W D by calling for each strip (until L W D is empty) two times CreateLayer(L W D , 2v). Let us denote by Lay 2i-1 and Lay 2i the layers created in strip S i . Let us say that Lay 2i-1 is packed left justified, and Lay 2i is packed right justified. Moreover, each layer is repacked in non increasing order of the widths, such that the narrowest rectangles are packed on the top.

Let N 1 denote the number of strips used in phase 1, and let i 1 denote the index of the last created layer (Lay i1 is of course in S N1 ). Let L 1 H and L 1 XH denote the set of remaining rectangles after phase 1 of L H and L XH , respectively. Thus, for the moment we have

π i 0 = Lay 2i ∪ Lay 2i-1 for all i ≤ N 1 . Lemma 7 If ∃i 0 < i 1 such that H(Lay i0 ) ≤ 3v 2 then it is straightforward to preallocate I ′ . Proof Let i 0 < i 1 such that H(Lay i0 ) ≤ 3v
2 . This implies that we ran out of rectangles of L W D \ (L H ∪ L XH ) while creating layer i 0 . Thus, because of the BF H order there are at least two rectangles of L B in every layer Lay i , for 1 ≤ i < i 1 , implying that the width of high and very high rectangles packed in each of these layers is strictly larger than 2 /3. Thus, W (π i 0 ∩ (L H ∪ L XH )) > 4/3 > 1 for 1 ≤ i < N 1 . Thus, the total width of remaining high and very high rectangles is lower than N -(N 1 -1).

Let us prove that we can pack all the remaining rectangles of I ′ (which are included in (L H ∪ L XH )) in the remaining strips. For each i ∈ [|N 1 + 1, N |] we create two shelves in S i (one at level 0 and one at level v). If there are still some unpacked rectangles, then all the shelves are "full", that is the width of each shelf is larger than 2 /3 (as all the width of any rectangle of L H ∪ L XH is lower than 1/3). Thus, we have

W (π i 0 ∩ (L H ∪ L XH )) > 4 /3 > 1 (for N 1 + 1 ≤ i ≤ N ).
This implies that the total width of remaining rectangles of L H ∪ L XH (including those in strip S N1 ) is now lower than 1. Thus, we can pack all of them in one shelf in S N1 .

From now we assume that H(Lay i ) > 3w 2 for all i < i 1 . This implies that S(

π i 0 ) > v for i < N 1 . Moreover, we have 2N v ≥ H(L W D ) ≥ N1-1 i=1 H(π i 0 ∩ L W D ) > (N 1 - 1)2( 3v /2), implying N 1 < 2 3 N + 1. It remains now to pack L 1 H ∪ L 1 XH . Notice that (L 1 H ∪ L 1 XH ) ∩ L W D = ∅ (we say that (L 1 H ∪ L 1 XH
) contains purely high and very high rectangles).

Packing high and very high rectangles

Preliminaries

Let N 2 = N -N 1 denote the number of free strips after phase 1. Roughly speaking, phase 2 packs shelves of high or very high rectangles in each of the N 2 last strips and merges some high or very high rectangles with the ones packed in strip N 1 (using the M erge procedure).

In this section we present a technique to fill γ empty strips with high or very high rectangles. In the Section 3.5, we use this technique for γ = N 2 (using strips S N1+1 . . . S N ) and an additional merging algorithm (that fills efficiently strip S N1 ) to pack L 1 H ∪ L 1 XH . Let us now introduce the procedure GreedyPack(X, seq). Given an ordered sequence of bins seq, GreedyP ack creates for each empty bin b ∈ seq a shelf of rectangles of X using CreateShelf (X, 1) and packs it into b (an example of a shelf packed in a bin is depicted Figure 3, Page 3). This procedure returns the last bin in which a shelf has been created, or null if no shelf is created. Notice that we will always use sequence of bins that have always width 1, and the same height h b such that max rx∈X h x ≤ h b .

We now define the two sequences of bins seq XH and seq H that will be used by GreedyP ack. Every bin of seq XH (resp. seq H ) will (possibly) contain one shelf of rectangles of L XH (resp. L H ). Notice that in a free strip it is possible to pack two bins of height v (width of bins is always 1), three bins of height 2v /3, or one bin of size v and one bin of size 2v /3. Thus, seq XH is composed of 2γ bins (b 1 , . . . , b 2γ ) of height v, considering that we created two bins of height one in each of the strips S 1 , . . . S γ . More precisely, for all i we locate b 2i-1 and b 2i in S i , with b 2i-x at level v(1x) for x ∈ {0, 1}. The sequence seq H is composed of 3γ bins (b ′ 1 , . . . , b ′ 3γ ) of height 2v /3, considering that we created three bins in each of the strips S γ , . . . S 1 .

It means that for all

i ≥ 1, bins b ′ 3i-2 , b ′ 3i-1 and b ′ 3i are lo- cated in S γ-i+1 , with b ′ 3i-x at level 2xv 3 for x ∈ [|0, 2|
]. This sequences of bins will be used in Lemma 9, and later in phase 2.

Finally, let us define the Add(X, S i last ) procedure that packs the set of rectangles X ⊂ L H \ L W D in S i last . As one can see in Lemma 9, S i last is the last strip where Greedypack created a shelf. Thus, we assume for the moment that S i last may only contain two different shapes of packing, and define the Add procedure accordingly.

In the first case S i last contains a first "full" shelf (full means that the surface of the shelf is at least v /2) of rectangles of L XH at level 0, and a shelf sh of rectangles of L XH packed at level v, right justified. In this case, Add creates a shelf sh 1 using CreateShelf (X, 1 -W (sh)) and preallocate sh 1 at level v, left justified. In the second case (see Figure 5), S i last contains only a shelf of rectangles of L XH packed at level 0, right justified. In this case, Add first moves some rectangles from sh to a new shelf sh until W (sh) ≤ 2 /3. Then, Add packs (right justified) the widest of these two shelves (denoted by sh A ) at level 0, and the other one (denoted by sh B ) at level v. Finally, Add creates two shelves sh 1 and sh 2 using CreateShelf (X, 1 -W (sh A )) and one shelf sh 3 using CreateShelf (X, 1 -W (sh B )). Then, sh i is packed at level 2v(i-1) 3 , left justified. Notice that stacking shelves sh 1 , sh 2 , sh 3 does not exceed 2v.

We end these preliminaries with the following Lemma about the efficiency of Add. Lemma 8 Let X ⊂ L H \ L W D and S i be a strip packed as expected for Add(X, S i ). Let π i 0 denote the rectangles packed in S i before the call Add(X, S i ). If X = ∅ after calling the procedure, then S(π i 0 ∪ X) > v.

Proof Remember that two cases are possible according to what is already packed in S i before the call. Let us first suppose that there is one full shelf (of area strictly larger than v /2) of very high rectangles (at level 0) and another shelf sh of very high rectangles at level v. Then, X = ∅ after the call implies that W (X) > -W (sh), and we have S(

π i 0 ∪ X) > v 2 + W (sh) 2v 3 + W (X) v 2 > v.
Let us now suppose that S i contains only one shelf sh of L XH at level 0. Let sh A , sh B , sh 1 , sh 2 , sh 3 be defined as described in the Add procedure. As W (sh A ) ≤ 2 3 (W (sh b ) ≤ 2 3 is also true), and X ∩ L W D = ∅, sh 1 and sh 2 contain at least one rectangle, implying that W (sh 1 ) and W (sh 2 ) are strictly larger than 1-W (shA) 2 according to Lemma 2. Moreover, X = ∅ after the call implies that after creating sh 1 and sh 2 the total width of remaining rectangles of X was strictly larger than 1 -W (sh B ). Putting this together, we get S(π

i 0 ∪ X) > (1 -W (sh A )) v 2 + (1 - W (sh B )) v 2 + (W (sh A ) + W (sh B )) 2v 3 > v.

Filling γ empty strips with high and very high rectangles

The next lemma shows how to fill γ free strips.

Lemma 9 Let L XH ⊂ L XH \ L W D and L H ⊂ L H \ L W D
be two sets of rectangles that we have to pack. Suppose that we execute the following calls:

1. last = GreedyP ack( L XH , seq XH ) 2. GreedyP ack( L H , seq H ) 3. Add( L H , S i last )
where S i last denotes the strip containing the bin "last".

Then, we get the following properties:

• If L XH = ∅ after 1, then S( L XH ) > (γ + 1 6 )v • Otherwise, if L H = ∅ after 3, then S( L XH ∪ L H ) > γv.
Remark 10 Let X such that X ∩ L W D = ∅ and let sh denote a shelf created by CreateShelf (X, 1), supposing that we did not run out of rectangle while creating the shelf. Then, according to Lemma 2, as at least three rectangles fit we have W (sh) >

i ∈ [|1, i last -1|], S(π i 0 ∩ ( L XH ∪ L H )) > 2 v 2 = v. For all i ∈ [|i last + 1, γ|], S(π i 0 ∩ ( L XH ∪ L H )) > 3 3v 8 > v. According to Lemma 8, L H = ∅ implies S(π i last 0 ∩ L H ) > v. Thus, we get that S( L XH ∪ L H ) > γv.

Phase 2

In phase 1 we preallocated L W D in strips S 1 , . . . , S N1 . Recall that each layer created in phase 1 is sorted with the narrowest rectangles on the top. It remains now to preallocate L 1 XH ∪ L 1 H in S N1 , . . . , S N .

Theorem 11 It is possible to preallocate L 1 XH ∪ L 1 H in S N1 , . . . , S N with a resulting height lower than 2v. Thus, our algorithm is a 2-approximation.

Proof Due to lack of space, we only sketch the proof of the theorem here and refer the reader to [START_REF] Bougeret | Tight approximation for scheduling parallel jobs on identical clusters[END_REF] for the complete case analysis.

An example of the final packing is depicted Figure 6. Phase 1 optimally filled the (N 1 -1) first strips using rectangles of L W D . Phase 2 optimally filled strips N 1 + 1 to i last -1 using two shelves of very high rectangles in each strip, and optimally filled strips i last + 1 to N using three shelves of high rectangles in each strip. Strips N 1 and i last will be carefully filled according to a case by case analysis. According to Section 3.3, we know that the area packed in the N 1 -1 first strips is greater than (N 1 -1)v. The general idea is to pack L 1 XH ∪ L 1 H is to call GreedyP ack (see Lemma 9) on strips S N1+1 , . . . , S N . If L 1

XH and L 1 H are entirely packed by GreedyP ack, then all the rectangles are packed in 2v. Otherwise, we can use Lemma 9 to claim that, S(L 1 XH ∪ L 1 H ) > (N -N 1 )v (or (N -N 1 + 1 6 )v), and the remaining rectangles are packed by carefully studying what was packed in S N1 by phase 1. Thus, we finish the proof using a case distinction according to the shape of the preallocation in S N1 .

Concluding remarks

According to the main steps defined in Section 2.2, the previous 2 phases algorithm that preallocates I ′ is sufficient to get a 2-approximation. Indeed, we simply add rectangles of I \ I ′ using list algorithms defined in Section 2.2.

Let us now sketch the analysis of the running time of the preallocation algorithm. Phase 1 and Phase 2 run in O(n log(n) + N n) as for each strip creating a layer or a shelf can be done in O(n) (by sorting wide rectangles according to their heights before phase 1, and sorting high rectangles according to their widths before phase 2). The list scheduling algorithms of Section 2.2.2 that turn π 0 into the final packing can be implemented in O(n log(n)) by using a global list (i.e. that refers to the N strips) of currently "scheduled" rectangles, instead of scanning level by level and strip by strip (which is in O(N n log(n))).

Finally, taking into account the repetitions due to the binary search on v (where 0 ≤ v ≤ nh max ), the overall algorithm runs in O(log(nh max )n(N + log(n))).

Figure 1 .

 1 Figure 1. Example for 9 jobs and 2 clusters of a feasible solution for M CSP and not feasible for M SP P .

= w 1 + w 2 Figure 3 .

 23 Figure 3. Example of a layer, a shelf, a bin and of the utilization function. sh is packed in b.

Figure 4 .

 4 Figure 4. (1α) compact packing

Figure 5 .

 5 Figure 5. Example of the add procedure.

Figure 6 .

 6 Figure 6. Example of a preallocation.

  .

	Problem Ratio Remarks		Source
	M CSP ,	2ρ	Need solving P ||C max	[10]
	M SP P		with a ratio ρ
	M CSP	5/2	Fast algorithm		[3]
	M CSP ,	2	Need applying the PTAS	[2]
	M SP P		of [1] for a constant (≈
			10 4 ) number of clusters
	M CSP ,	AFPTAS	Additive constant in	[2]
	M SP P		O( 1 ǫ 2 ), and in O(1) for large values of N
	M CSP	3	Fast (and decentralized)	[9]
			algorithm	handling
			clusters of different size
	M CSP	2	Requires max j w j ≤ 1 2

  Proof [Proof oflemma 9] Let us first suppose that L XH = ∅ after step 1. It implies that after creating the first 2γ -1 shelves of width at least 3 /4 (according to Remark 10), the total remaining width of rectangles of L XH was strictly larger than 1. Thus, S( L XH ) > Let us now suppose that L XH = ∅ and L H = ∅ after step 3. Let sh denote the shelf of rectangles of L XH contained in bin last. Remind that i last is the index of the strip containing last. For all

	(γ + 1 6 )v.	3 4 (2γ-1) 2v 3 + 2v 3 =

3

/4 . Moreover, if X ⊂ L XH then S(sh) > v /2, and if X ⊂ L H then S(sh) > 3v /8.

Therefore, when trying to pack a rectangle of width w j at level l of a strip, we only need to check if the occupation at level l is lower than 1w j .

Even if some recent advances in the PTAS design for P ||Cmax allowed to decrease the asymptotic depedencies in 1 ǫ (like 2 O( 1 ǫ 2

log 3 ( 1 ǫ )) )in[START_REF] Jansen | An EPTAS for scheduling jobs on uniform processors: using an MILP relaxation with a constant number of integral variables[END_REF], the running time of these algorithms remain high due to constants.