Laurent Houssin
email: houssin@laas.fr

Christian Artigues
email: artigues@laas.fr

Erwan Corbel
email: erwan.corbel@thalesaleniaspace.com

Frequency allocation problem in a SDMA satellite communication system

Keywords: SDMA system, radio resource management, frequency assignment problem, integer programming, greedy algorithms

come

Introduction

Satellite communication systems move towards greater capacity, higher flexibility (with respect to the position of the users) and better service to the end-user. SDMA (Spatial Division Multiple Access) appears to be one way to achieve these requirements at the same time [START_REF] Liberti | Smart Antennas for Wireless Communications[END_REF]. SDMA is a principle of radio resource sharing that relies on the division of the space dimension into separated communication channels. It can be used with common Frequency Division Multiple Access (FDMA), Time Division Multiple Access (TDMA) or Code Division Multiple Access (CDMA) techniques. Main future terrestrial communication standards (such that WIMAX, 3GPP, LTE) implement SDMA.

Today SDMA is currently used by IRIDIUM system in L-band, a constellation of 66 Low-Earth Orbit satellites, thanks to time beam-switching. SDMA is also foreseen as a key enabling technique to increase the capacity of future two-way satellite communications systems in low-frequency bands (typically lower than 5-6 GHz) through the interference mitigation and high frequency reuse [START_REF] Corbel | Hybride satellite and terrestrial mobile network for 4g: Candidate architecture and space segment dimensioning[END_REF]. It is also expected to play an important role in future systems devoted to Public Protection and Disaster Relief (PPDR) and Global Monitoring for Environment and Security (GMES) missions which require on-demand beam-forming [START_REF] Gayrard | The role of satellites in broadband wireless access for public protection, disaster relief and gmes missions[END_REF].

SDMA basically relies on adaptive and dynamic beam-forming associated to a clever algorithm in charge of resource allocation. The satellite beamformer optimizes the antenna diagram with respect to the positions of the users in order to maximize the gain while mitigating interferences. The resource allocation algorithm carefully designs a frequency plan that

• prevents or limits interferences between users,

• tailors the allocated bandwidth to the user need in order to save the spectrum.

The spatial filtering capability of adaptive antenna can be used to implement SDMA in communication systems where users share the same channels. Moreover, the resource management by spatial filtering of users give the opportunity to massively reuse radio channels. It makes the channel allocation strategy play an important role in the system performance. This class of problem is well known as Frequency Assignment Problem (FAP). The satellite telecommunication system that we study in this paper aims at establishing bi-directional communications involving a gateway connected to a terrestrial network and some user terminals located in a service area. This paper studies static frequency allocation problems in this system and two scenarios are considered: one based on SDMA and the other based on "traditional" spot coverage. We propose original Integer Linear Programming (ILP) formulations and greedy allocation algorithms for these problems. The difficulty for solving the FAP is increased by considering cumulative interference constraints. We then compare the performance of the two scenarios.

The remainder of the paper is organized as follow: Section 2 is dedicated to the telecommunication system and the description of the scenarios. In Section 3, we present a cumulative formulation of the FAP interference constraints. In Section 4, ILP formulations and greedy algorithms proposed for both scenarios are described. Section 5 presents the results obtained by the different algorithms on the scenarios. Concluding remarks are drawn in Section 6.

Telecommunication system and scenarios for frequency allocation

System description

The service area is a rectangular grid where users position is uniformly distributed and where beams are directed. The grid size is u = [-0.043980, 0.048870] and v = [-0.021152, 0.012702] (Cartesian coordinates are considered although the same study can be done with spherical coordinates) with a step of 5 × 10 -4 . The satellite orthogonal projection onto the service area defines point (0, 0) of the grid.

Beams have two particular characteristics which are the direction that we consider through the position of the beam center in the service area, and the radiation pattern. An analytic representation of the radiation pattern enables to compute the directive gain of the antennas for a considered direction. This description is such that

G Sat (u, v, u 0 , v 0 ) = G 1 × G 2 (u, v, u 0 , v 0 ) × G 3 (u, v, u 0 , v 0 ) (1)
with

G 1 = η πD λ 2 , G 2 (u, v, u 0 , v 0 ) =   2J 1 πD λ (u -u 0 ) 2 + (v -v 0 ) 2 πD λ (u -u 0 ) 2 + (v -v 0 ) 2   2 and G 3 (u, v) = 2J 1 πd λ √ u 2 + v 2 πd λ √ u 2 + v 2 2 ,
where J 1 (x) is the Bessel functions of the first kind. A singular property of Bessel functions of the first kind is that The equation [START_REF][END_REF] was provided by the satellite company Thales Alenia Space. To the best of our knowledge, this equation has a particular shape and is not known in the literature. The term G 1 corresponds to the maximum gain antenna whereas G 2 depends on the distance between the user and the beam center, i.e.

J 1 (x) x = 1/2 if x = 0. It leads to G 2 (u, v, u 0 , v 0 ) = 1 if (u, v) = (u0, v0) (G 2 (u, v, u 0 , v 0 ) < 1 otherwise) and G 3 (u, v) = 1 if (u, v) = (0, 0) (G 3 (u, v) < 1 otherwise) .
(u -u 0 ) 2 + (v -v 0 ) 2 and G 3 depends on user position related to the satellite, i.e. 1 and by a 2D representation in fig. 2 where the left pattern is the pattern of beam centered in (0.03,-0.03) and the right one is centered in (0,0), i.e directly under the satellite. It emphasizes that the more the beam center is closed to the satellite, the more the gain is high. We can also remark in fig. 2 that closed to the sharp summit, small waves exist and will disturb users that share the same frequency. A minimum quality for a communication between a user and the gateway is required. It corresponds to C

√ u 2 + v 2 . The product G 2 (u, v, u 0 , v 0 )×G 3 (u, v, u 0 , v 0 ) is illustrated by a 3D representation in fig.
N RsM odCod that depends on the modulation scheme and the code scheme (the digital communication scheme considered here is a phase key shifting). More precisely, a user will be served if

C N + I ≥ C N RsM odCod (2)
where C N +I is the user Signal-to-Interference-plus-Noise Ratio (SINR) which determines quantitatively the signal quality.

The link budget enables to compute the SINR :

C N +I -1 = C N -1 F eeder + C I -1 F eeder + C IM -1 + C N -1 U ser + C I -1 U ser , (3)
in which

C N U ser = K 1 × G Sat (User Beam → User) K 2 and C I U ser = G Sat (User Beam → User) j∈Interf (G Sat (User Beam → User(j))
.

Terms K 1 and K 2 involve technical parameters, such as the atmospheric loss, the antenna temperature and the Equivalent Isotropically Radiated Power (EIRP) which are considered as constant. Consequently, C N U ser and C I U ser only depends on user position and beam center position. Indeed, G Sat (User Beam → User) describes the gain for a user terminal and its beam (see equation ( 1)). The set Interf is the set of users sharing the same channel. Moreover, C N f eeder and C I f eeder which are gateway characteristics, and C IM that is a satellite characteristic are also constant. Remark 1. For simplification reasons, considerations on signal processing are not tackled. However, for the different ratios C I , C N and C N +I , the variables C, I and N mean respectively "carrier", "interferences" and "noise". Specifications of the system are the following.

• Beams are only adaptive in direction and not in shaping (although it is technically feasible).

• For the first scenario (involving SDMA), each user has a beam directly centered on him. It leads to consider as many beams as users in the system.

• We only focus on the case of users-to-feeder (gateway) link.

Scenarios description

In the users-to-feeder link, interferences can occur in the link budget when several users share the same frequency. We illustrate this phenomena in fig. 3 where black diamonds are users sharing the same frequency. It shows that users u 3 and u 2 are interferers for user u 1 . In the first scenario, named "scenario 1", we consider that each user has a beam directly centered on him (which is possible thanks to SDMA). The number of available channels is 8. In this case,

G 2 (u, v, u 0 , v 0 ) is always equal to the maximum value (that is G 2 (u, v, u 0 , v 0 ) = 1) since (u, v) = (u 0 , v 0 ).
Contrary to scenario 1, the second scenario involves fixed beams, however the frequency assignment is variable according to the demand (the channel number can be adjusted for a beam). The service area is composed of 40 fixed beams which form a spot-based coverage. We also have 8 available channels. Consequently, (since no channel can be used more than once in a spot) we can not serve more than 8 users in a spot.

For both scenario, the FAP considered in this paper consists in finding and interference-free frequency allocation to the users maximizing the number of served users. This problem is static: the set of users is known in advance and we do not take dynamic arrivals and departures of users into account.

A cumulative interference representation

In a FAP where channels are limited, results depend on the ability of the system to allocate the same channel to several users.

In this section, we show that it is possible to obtain a cumulative representation of interferences from equations ( 2) and (3).

Involving the link budget, constraint (2) for a user i, becomes 1

A + B i + P j∈Interf (i) G Sat (U ser Beam(i)→U ser(j)) G Sat (U ser Beam(i)→U ser(i))
≥ D

where

A = C N -1 F eeder + C I -1 F eeder + C IM -1 , B i = C N -1 U ser i = K 2 K 1 × G Sat (User Beam(i) → User(i)) and D = (C/N) RsM odCod .
Previous inequation leads to

1 ≥ AD + DB i + P j∈Interf (i) DG Sat (U ser Beam(i)→U ser(j)) G Sat (U ser Beam(i)→U ser(i))
and finally, the cumulative representation follows

j∈Interf (i) DG Sat (User Beam(i) → User(j)) ≤ G Sat (User Beam(i) → User(i)) (1 -AD -B i D).
For scenario 1, we write last constraint

j∈Interf (i) δ ij ≤ α i , (4)
where

α i = G Sat (User Beam(i) → User(i)) (1 -AD -B i D)
and δ ij = DG Sat (User Beam(i) → User(j)) .

Let the system with n users, α is a n-row vector and δ is a n × n matrix. We can deduce that α i represents the maximum level of interferences that the user i can support. In this way, δ ij describes the interference level of user j if users j and i share the same frequency. The cumulative representation is motivated by a linear representation (although constraint (2) was not linear) that enable the use of integer linear programming (see §4.2).

Concerning scenario 2, we denote the interference inequality by

j∈Interf (i) γ ij ≤ β i . (5)
It is worth emphasizing that α i and β i are different since the center of the beam related to user i does not correspond to user i coordinates in scenario 2 whereas in scenario 1, it does. The same remark can be done for δ ij and γ ij .

Modeling and solving scenarios 1 and 2 FAP

FAP literature overview

Most approaches dealing with interference minimization FAP consider binary interference constraints, i.e. involving only two users. Because of the strong links between graph coloring and frequency allocation with binary interference constraints, most methods found in the literature are inspired by coloring algorithms. We also know unfortunately the graph coloring problems, and consequently the FAP, are NP-hard. Among the proposed methods, the constructive (greedy) algorithms are widely used since they are simple, fast and also able to solve dynamic FAP. In this category, we find the generalization of DSATUR procedure [START_REF] Brélaz | New methods to color the vertices of a graph[END_REF]. Other more sophisticated algorithms, such as local search, metaheuristics, ILP and constraint programming approaches, are frequently encountered, see [START_REF] Resende | Handbook of Optimization in Telecommunications[END_REF] [START_REF] Murphey | Handbook of Combinatorial Optimization: Supplement Volume A[END_REF] and [START_REF] Aardal | Models and solution techniques for the frequency assignment problem[END_REF] for a comprehensive review of the state-of-the-art methods for the FAP with binary interference constraints.

One of the difficulties appearing in the telecommunication system considered in this study (for both scenarios 1 and 2) lies in the explicit consideration of non-binary interference constraints. In terms of graph coloring, deciding whether a given coloring is feasible or not cannot be made anymore by checking pairwise user color assignments. Instead, for a given user, the cumulative interferences of the users assigned to the same color has to be computed. Then, the coloring is feasible if this cumulative interference remains under a user-dependent threshold (see Section 3). In the literature, only a few approaches take account explicitly of such interferences for frequency assignment [START_REF] Dunkin | Towards high order constraint representations for the frequency assignment problem[END_REF]; [START_REF] Mannino | An enumerative algorithm for the frequency assignment problem[END_REF]; [START_REF] Palpant | Models and methods for frequency allocation with cumulative interference constraints[END_REF]. This study is partly based on integer linear programming formulations proposed in [START_REF] Palpant | Models and methods for frequency allocation with cumulative interference constraints[END_REF].

Integer linear programming formulations (ILP)

Taking account of hypothesis and simplifications presented in Section 2, FAP corresponding to scenarios 1 and 2 is at first sight similar to coloring problems and thus formalized as the corresponding combinatorial optimization problems. Each user has to be assigned a color, representing the allocated carrier. However there are some significant differences. There cumulative interferences and the number of color is fixed while the number of assigned users must be maximized. This correspond to the MaxFAP problem in the classification described in [START_REF] Aardal | Models and solution techniques for the frequency assignment problem[END_REF].

SDMA FAP (scenario 1)

For the SDMA FAP, (scenario 1), the following data are considered. n denotes the number of users. U = {1, . . . , n} is the set of users. C is the number of colors (channels). α i denotes the interference threshold for user i. δ ij is the interference component from user j on user i, if i and j are assigned the same color.

Binary decision variables x ic are defined for i ∈ {1, . . . , n} and c ∈ {1, . . . , C} with n, the number of users and C the number of available colors.

x ic = 1 if color c is allocated to user i and x ic = 0 otherwise. The problem can be represented by the following ILP:

max n i=1 C c=1 x ic (6) C c=1 x ic ≤ 1 i = 1, . . . , n (7)
n j=1 δ i,j x jc ≤ α i + M i (1 -x ic ) i = 1, . . . , n c = 1, . . . , C (8)
x ic ∈ {0, 1} i = 1, . . . , n c = 1, . . . , C (9)
Objective ( 6) consists in maximizing the number of accepted users. Constraints ( 7) state that at most one color has to be selected for each user. Constraints (8) are the cumulative interference constraints. They represent, in case color c is allocated to user i, the respect of the threshold for user i taking account of users that are assigned color c, i.e. possible interferers. Constant M i has to be large enough to withdraw the constraint if i is not assigned color c (x ic = 0). More precisely, we set

M i = n j=1 δ ij -α i .

Fixed beam FAP (scenario 2)

For the fixed beam FAP (scenario 2), the data are similar with additional features concerning the spots (a spot designing the area covered by a given beam): m denotes the number of spots. S = {1, . . . , m} is the set of spots. U s is set of users covered by spot s ∈ S. A user i belong to the set U s corresponding to his closest beam. β i denotes the interference threshold for user i. γ ij is the interference component from user j on user i if i and j are assigned the same color.

We define the fixed beam FAP as a combinatorial optimization problem resembling the SDMA FAP preventing two users covered by the same spot from being assigned the same non-zero color. We obtain the following ILP:

max n i=1 C c=1 x ic (10) C c=1 x ic ≤ 1 i = 1, . . . , n (11)
i∈Us

x ic ≤ 1 s = 1, . . . , m c = 1, . . . , C (12)
n j=1 γ i,j x jc ≤ β i + N i (1 -x ic ) i = 1, . . . , n c = 1, . . . , C (13)
x ic ∈ {0, 1} i = 1, . . . , n c = 1, . . . , C (14)
Scenario 2 ILP differs with the scenario 1 ILP via constraints ( 12) and ( 13). Constraints ( 12) prevent any color from being allocated more than once in a given spot. Hence a maximum number of C users may be served in the same spot. In constraints ( 13), values β i and γ ij differ from α i and δ ij since beams are not centered on the users (see Section 3). As for scenario 1, N i has to be large enough to ensure the constraint is verified when color c is not allocated to user i. For this purpose, we set N i = n j=1 γ ij -β i . Both above-defined ILP can be solved by an integer linear programming solver, via branch and bound.

Greedy algorithms

Solving the ILP formulations provides optimal solutions only for small problems. For large-sized problems it is necessary to use a heuristic. We propose greedy algorithms to solve scenarios 1 and 2 FAP. For both scenarios, the principle of the greedy algorithms is, first, to consider the users sequentially according to a given criterion named the user priority rule. Second, either the selected user is assigned a frequency or rejected according to a second criterion, the frequency priority rule. Let Q denote the set of users that have not been assigned a color yet. Initially, we have Q = U. At each step of the greedy algorithm, a user i is removed from Q and is either rejected or assigned a color. In the following sections we described the variants proposed for each scenarios.

SDMA FAP (Scenario 1)

For scenario 1, the principle of the greedy algorithm is summarized in Algorithm 1, where F j denote the frequency (color) allocated to user j if 1 ≤ F j ≤ C and F j = 0 indicates user j is rejected.

Algorithm 1 Greedy Algorithm Require: n, C, α, δ F i ← 0, ∀i = 1, . . . , n for q = 1 to n do i ← SelectUser(n, C, α, δ, F ) F i ← SelectColor(i, n, C, α, δ, F ) end for Return F
For both the user and priority rules, we may use the frequency margin, where the margin M(i, c) of a user i ∈ Q for a color c is given by M(i, c) = α i -j∈U \Q∪{i},F j =c δ ij . Namely, this margin corresponds to the positive or negative slack of the cumulative interference constraint for user i terminal if it is assigned color c.

As a preliminary result, we observed that the user priority rule aiming at selecting first the most constrained users in terms of available colors, as for the well-known DSATUR algorithm for standard graph coloring problems gives bad results for this problem. In fact we were brought to consider a kind of hybrid reverse DSATUR rule, aiming a assigning users with a large number of available colors, alternating this rule with the selection of the user having maximum interferences with the previously assigned user. More precisely, we tested the two following user priority rules.

• Lexicographic. The user with the smallest number is selected.

• Hybrid. The proposed rule selects the user having the largest number of available colors, a color c being available for user i ∈ Q if M(i, c) ≥ 0 and if for all users j ∈ U \ Q having already been assigned color c, M(j, c) ≥ 0. In case of a tie, we select the user having the largest total margin for all its available colors. Let i denote the user selected with this rule. For the next iteration, we select the user having maximum inteference with i, i.e. the user j maximizing δ i,j + δ j,i and we alternate the two rules.

For the frequency selection, we propose to use two priority rules:

• Lexicographic. The smallest available frequency is selected.

• Most used. The most used available frequency is selected. In case of tie we select the color c that maximizes the sum of margins M(j, c) for all users j ∈ Q.

The proposed greedy algorithms run in O(n 2 C) time.

Fixed beam FAP (Scenario 2)

For scenario 2, the above-described user and frequency priority rules can be applied with the notable difference that

• values β i and γ i,j are used instead of α i and δ i,j , respectively.

• color availability for a user i ∈ Q in a spot s has to consider the other users in U s . As soon as a user in s has been assigned a color c, then c becomes unavailable for the other users in U s .

Computational experiments and simulations

The ILP formulations have been solved using IBM/ILOG CPLEX 12.1 CPLEX (2009) and the greedy algorithms have been coded in C++. We tested the proposed algorithms with C = 8, increasing stepwise the numbers of users by 20 from 20 to 200 users (which corresponds to reuse rates from 20/8=2.5 to 200/8=25), yielding 1000 instances. With regard to scenario 2, we consider 40 fixed beams (m = 40). For each number of users, a set of 100 FAP data instances was obtained by randomly generating the user positions on the service area (uniform distribution). The results were obtained on a Intel Core 2 Duo processor with 2,93GHz and 1.9 GB RAM. The CPU time for the ILP resolution has been limited to 60 seconds afterwhich the best obtained integer solution is returned. The CPU times for the greedy algorithms were negligible.

We first present a comparison of the greedy algorithms. We selected only the algorithms performing better than the simple lexicographic rule (for user and frequency). Table 1 reports the average number of accepted users over the 1000 instances for the three best algorithms.

Fig. 4 displays, for each scenario/algorithm/number of users, the average number of accepted users in the computed frequency allocation plans. For the greedy algorithms we selected for each instance the best value obtained by the three best algorithms displayed in Table 1.

Analysis

Table 1 shows that the results of the greedy algorithm is very close. The hybrid method we propose for the user selection rule gives the best results with the most used frequency rule. However it was difficult to give better results than the simple lexicographic rules.

Fig. 4 clearly shows that the best results in terms of quality of service using the proposed algorithms are obtained for scenario 1 (SDMA-based system). ILP-based algorithms obtain significantly better results than the simple greedy algorithms for both scenarios. However, the greedy algorithm for scenario 1 performs better than the ILP-based algorithm for scenario 2.

Fig. 5 shows as expected that the gap increases as the number of users increases. The gap is also more important for scenario 2 than for scenario 1. This shows that the optimization problem of scenario 1 is easier to solve than the optimization problem of scenario 2, at least for the tested methods. This is confirmed by the ILP results in optimality terms. Indeed, Table 2 shows that the number of proved optimal solutions falls from 100% to 82% then 6% when the number of users increase from 60 to 100 for scenario 2 while this number stays above 90% of up to 120 users for scenario 3. Forbidding identical colors in the same spot makes the search of an optimal solution harder. Furthermore, crossing Figure 4 and Table 2, we see that for instances where 100% of optima are proved for both scenarios, i.e. for instances having up to 60 users, no difference can be observed between the scenarios since all users are accepted. To rigorously compare the performance of the scenarios, we gather the instances with n ≥ 80 and for which the optimum number of accepted users was proved by ILP for both scenarios. We obtain 63 instances of known optima in both cases. Among these 63 instances, Greedy algorithm for scenario 1 obtains solutions strictly better that ILP for scenario 2 for 24 instances. It obtains only one strictly worse solution. Hence the simple greedy algorithm for scenario 1 outperforms the optimal solutions for scenario 2.

This numerical example with real parameters for the simulation shows qualitatively the benefits of scenario 1 even with a simple greedy algorithm. We can also note that beyond 80 users, which correspond to a reuse rate of 10, differences between algorithms increase.

Conclusion

In this paper we have developed integer linear programming formulations and greedy algorithms for frequency assignment problems involving cumulative interferences. The cumulative approach enables to take into account the non linear characteristics of interferences. Combining the SDMA system and the cumulative approach, we proved in section 6 the efficiency of scenario 1 with the two algorithms. Even better results and system optimizations using SDMA could be obtained by allowing shifts of the beam centers around the users and adjusting the EIRP parameters, yielding as a counterpart harder FAP. These features, together with consideration of dynamic aspects, constitutes the basis for further research.

  We use the following notations: u, v Cartesian coordinates of the user terminal; u 0 , v 0 Cartesian coordinates of the beam center; η antenna gain; D antenna diameter; d primary source diameter; λ wavelength.

Figure 1 Figure 2 :

 12 Figure 1: 3D representation of beam patterns

Figure 3 :

 3 Figure 3: Users-to-feeder link.

Figure 4 :

 4 Figure 4: Average numbers of accepted users obtained by the proposed algorithms for each number of users

Figure 5 :

 5 Figure 5: Average gap obtained by the proposed algorithms for each number of users

Table 1 :

 1 Average number of accepted users over 1000 instances

	Scenario 1 Scenario 2

Table 2 :

 2 Number of optima proved by ILP


Laurent Houssin
email: laurent.houssin@laas.fr

Cyclic jobshop problem and (max,plus) algebra

Keywords: Cyclic scheduling, (max, plus) algebra

AbstractIn this paper, we focus on the cyclic job-shop problem. This problem consists in determining the order of a set of generic tasks on machines in order to minimize the cycle time of the sequence. We propose an exact method to solve this problem. For each solution, a linear max-plus model (possibly non causal) is obtained. To evaluate the performance of a considered schedule, we build the causal max-plus representation and compute the eigenvalue of the evolution matrix. A branch and bound procedure is presented.

INTRODUCTION

The main idea of cyclic (or periodical) scheduling is to perform a set of generic tasks infinitely often (or with a infinite time horizon). The schedule has no end and the production is specified by cycle periodically repeated one behind the other. This kind of problems arises in many contexts such as robotics [START_REF] Kats | A strongly polynomial algorithm for no-wait cyclic robotic flowshop scheduling[END_REF]), manufacturing systems [START_REF] Pinedo | Planning and Scheduling in Manufacturing and Services[END_REF], [START_REF] Hillion | Performance evaluation of job-shop systems using timed event graphs[END_REF]) or multiprocessor computing (Hanen and Munier [1995b]). Up to now, periodic scheduling problems have been studied from several points of view. Most of them consider graph based approaches, integer linear programming, Petri nets. Although these approaches are different, these tools usually leads to the same results. An overview of the cyclic scheduling problem is available in Hanen and Munier [1995a] and [START_REF] Brucker | A general model for cyclic machine scheduling problems[END_REF].

In this paper, we tackle the problem of cyclic jobshop scheduling which is proven to be NP-hard by [START_REF] Hanen | Study of a np-hard cyclic scheduling problem: The recurrent job-shop[END_REF]. More precisely, we propose a max-plus linear model (possibly non causal) for each schedule. This paper presents a method to obtain a simple (and causal) recurrence of the form x(k) = Ax(k -1) for the date of each task. The cycle time of the schedule is given by computing the eigenvalue of the evolution matrix. More precisely, decision variables correspond to an order of tasks on each machine (or dedicated processor). For each assignement of these decison variables, we associate a max-plus linear model. A branch and bound procedure is then presented to obtain the schedule characterized by the minimum cycle time. The remainder of this paper is as follows. In section 2, we introduce algebraic tools used throughout the paper and (max, +)-linear systems. In section 3, the cyclic scheduling problem problem is stated. We first describe the basic scheduling problem and then we propose an optimal solution for the cyclic scheduling problem with resource constraints. A simple numerical example is also developed during all of the paper.

MAX-PLUS MODEL FOR DISCRETE EVENTS SYSTEMS

Max-plus algebra enables to model discrete events systems that involve synchronization and delay phenomena. Applications of this theory have essentially concerned manufacturing systems [START_REF] Menguy | Just in time control of linear systems in dioid: cases of an update of reference input and uncontrollable input[END_REF], [START_REF] Lahaye | Models combination in (max,+) algebra for the implementation of a simulation and analysis software[END_REF], communication networks Boudec and P. [2001] and transportation networks [START_REF] Houssin | Just-in-time control of constrained (max, +)-linear systems[END_REF], de [START_REF] Vries | On max-algebraic models for transportation networks[END_REF]. Few works on scheduling can be found however we can mention [START_REF] Gaubert | Modeling and analysis of timed petri nets using heap of pieces[END_REF] and more recently [START_REF] Bouquard | Application of an optimization problem in max-plus algebra to scheduling problems[END_REF].

Max-plus algebra and discrete event systems

We first recall some algebraic tools. An exhaustive presentation of this theory can be found in [START_REF] Baccelli | Synchronization and Linearity[END_REF]. The (max, +) semiring is the set R∪{-∞} endowed with the max operator, written a ⊕ b = max(a, b), and the usual sum written a ⊗ b = a + b. The sum (resp. product) admits a neutral element denoted ε = -∞ (resp. e = 0), it leads to a ⊕ ε = a and a ⊗ e = a. For matrices, additions and products give (

A ⊕ B) ij = A ij ⊕ B ij and (A ⊗ B) ij = n k=1 A ik ⊗ B kj .
The behavior of such systems can be represented by some discrete functions called dater functions. More precisely, a discrete variable x(•) is associated to an event labeled x. This variable represents the occurring dates of event x. The numbering conventionally begins at 0: x(0) corresponds to the date of the first occurrence of x. In this way, autonomous DEDS can be modeled by a linear state representation

x(k) = Ax(k -1), (1) where x is a vector and A is the evolution matrix. Moreover, to each square matrix A, a graph G = (N, E) can be associated. Indeed, there is no arc (i, j) if A ji = ε, and (i, j) is labelled with A ji otherwise. In the same way, a dedicated (max, +) matrix is associated to each graph.

We now recall the classical spectral theorem of an irreducible square matrix. Theorem 1. Any irreducible matrix A of size n × n with entries in R ∪ {-∞} has a unique eigenvalue which is the maximal cycle mean of the graph associated to A.

Several algorithms for computing the eigenvalue are possible. We can mention the Karp's algorithm in O(|n| 3 ) [START_REF] Karp | A characterization of the minimum cycle mean in a digraph[END_REF], [START_REF] Dasdan | Faster maximum and minimum mean cycle algorithms for systemperformance analysis[END_REF]) and the Howard's algorithm [START_REF] Cochet-Terrasson | Numerical computation of spectral elements in max-plus algebra[END_REF]) of unproved complexity but showing excellent performance.

Max-plus model for non causal systems

We now focus on non causal (max, +)-systems and we establish conditions on the existence of a causal representation of these systems.

The general form of an autonomous (max, +)-system can be written as

x(k) = A q x(k -q), q ∈ Z.
By means of an extension of the state vector, last equation leads to this compact form

x(k) = A 1 x(k -1) ⊕ A 0 x(k) ⊕ A -1 x(k + 1).
(2) We are confronted to a fixed point equation (see [Baccelli et al., 1992, Th. 4.70]) and the least solution (that corresponds to the earliest functionning) is given by

x(k) = A ′ 1 x(k -1) ⊕ A ′ -1 x(k + 1), (3) with A ′ 1 = A * 0 A 1 and A ′ -1 = A * 0 A -1 considering X * = i∈N X i .
Proposition 1. If the matrix A -1 contains only one element i, j such that A -1ij = ε and i = j, the causal representation of the system is given by

x(k) = (A ′ -1 A ′ 1 ) * A ′ 1 x(k -1). (4)
Proof : First, we state that if A ′ 1 has also only one element different from ε and if it is not a digonal element, it involves that (A ′ -1 ) 2 = ε. If it is a diagonal element, the system can not admit a causal representation since it leads to

x i (k) A ii x i (k + 1). From (3), we deduce x(k) = A ′ 1 x(k -1) ⊕ A ′ -1 x(k + 1) = A ′ 1 x(k -1) ⊕ A ′ -1 (A ′ 1 x(k) ⊕ A ′ -1 x(k + 2)) = A ′ 1 x(k -1) ⊕ A ′ -1 A ′ 1 x(k) ⊕ (A ′ -1 ) 2 x(k + 2)) Last equality leads to x(k) = A ′ 1 x(k -1) ⊕ A ′ -1 A ′ 1 x(k), since (A ′ -1 ) 2 = ε.
In accordance with the earliest functionning rule (an event occurs as soon as possible), we select the least solution of last equality which is given by

x(k) = (A ′ -1 A ′ 1 ) * A ′ 1 x(k -1).
2 Remark 1. For the meantime, we didn't prove neatly the general case, i.e. no restrictions on the number of elements of A -1 different from ε.

CYCLIC SCHEDULING PROBLEMS

This section is devoted to cyclic scheduling problems. A schedule is called periodic (or cyclic) with cycle time ω if

x i (k) = ω + x i (k -1) = k × ω + x i (0)
where x i (k) denotes the k-th occurrence of task i. We describe these problems from a (max, +) point of view. In (max, +)-algebra the last equation becomes

x i (k) = ω ⊗ x i (k -1) = ω k ⊗ x i (0).
We first recall the basic periodic scheduling problem and then we tackle the problem of the jobshop cyclic scheduling problem.

The basic cyclic scheduling problem

We recall the basic cyclic scheduling problem. In this framework, there is no resource constraints since each task is performed without preemption on a dedicated machine.

In this problem we have:

• operations (or tasks) i = 1, . . . , n and the associated processing times p 1 , . . . , p n , • generalized precedence constraints (or conjunctive constraints) between the k-th occurrence of task i and the k + H ji -th occurrence of task j given by

x i (k + H ij ) ≥ p j x j (k).
The event shift function H is called the height. We denote T the set of operations. A uniform graph G u (T, U ) can be associated to the generalized precedence constraints in which a node is a task and an arc corresponds to a constraints. Each arc (i, j) ∈ U is supplied by two values p i and H ji .

A cyclic scheduling problem is said to be consistent, if and only if every circuit of G u has a positive height (Hanen and Munier [1995a], [START_REF] Brucker | A general model for cyclic machine scheduling problems[END_REF]). Let us define a matrix

H ′ such that H ′ ij = -H ij if (i, j) ∈ U and H ′ ij = ε, otherwise. Proposition 2.
The basic scheduling problem is consistent if and only if the largest eigenvalue of H ′ is negative.

Proof : It also means that the problem is consistent if and only if the minimum cycle mean of G u is positive. We can build a conjugate graph G ′ u (T, U ) in which every (i, j) ∈ U is labeled with -H ji . Consequently, the problem is consistent if and only if the maximum cycle mean of G ′ u is negative. Then, the matrix associated to G ′ u is H ′ . Hence, the consistency of the problem can easily checked by means of Karp's algorithm or Howard algorithm.

2

Of course, the k +1-occurrence of a task i can not start before the end of the k-occurrence of a task i. Consequently, we get the following non-reentrance constraints

x i (k + 1) ≥ p i x i (k) ∀i ∈ T
Such systems, if it is consistent, can easily be expressed as a causal (max, +)-linear system and the average cycle time is given by the largest (max, +)-eigenvalue.

Example 1. Let us consider the following basic scheduling problem task 1 2 3 4 5 processing time 2 3 1 2 2 and the following uniform constraints :

x

1 (k) ≥ 2x 4 (k -1) x 2 (k) ≥ 2x 1 (k) x 3 (k) ≥ 3x 2 (k) x 3 (k) ≥ 2x 5 (k -1) x 4 (k) ≥ 3x 2 (k) x 4 (k) ≥ 2x 5 (k -1) x 5 (k) ≥ 2x 1 (k + 1).
We first check the consistency of the problem. There is a solution since the eigenvalue of the matrix H ′ (in which non-reentrance constraints are included) associated to the constraints is equal to -1 3 . The (max, +) formulation leads to a non causal system of the form (2). Nonetheless, a causal representation of the earliest schedule can be obtained as shown in §2.2. The new associated model is

x(k) =      2 ε ε 2 ε 4 3 ε 4 ε 7 6 1 7 2 7 6 ε 7 2 11 10 ε 11 6      x(k -1).
The eigenvalue of the evolution matrix leads to 7. It corresponds to the optimal cycle time of this basic scheduling problem. The solution of the schedule is depicted in figure 1.

Figure 1. Optimal solution of the basic cyclic scheduling problem.

Remark 2. Notice that the number of work-in-process is not limited to one. For instance, the solution of example 1 presents two work-in-process for t ∈ [7, 11].

The jobshop cyclic scheduling problem

In a jobshop cyclic scheduling problem, tasks are mapped onto m machines. Naturally, we have m < n. It means that the number of machine is lower than the number of tasks (n = |T |). A presentation of this problem is proposed by Hanen and Munier [1995a].

We denote M (i) ∈ {1, . . . , m} the machine dedicated to the task i. Occurrences of operations to be processed on the same machine can not overlap. Indeed, new constraints have to be took in account to consider this allocation. Let us consider two tasks i, j performed on the same machine such that the next occurrence of j to be performed after the l-th task i is the k-th occurrence. Therefore, we can consider x j (k) ≥ p i + x i (l). Since the schedule is periodic, last equation is equivalent to x j (0) + kω ≥ p i + x i (0) + lω.

It is also equivalent to x j (0) ≥ p i + x i (0) + (lk)ω and we can state

x j (s) ≥ p i + x i (s + l -k) ∀s ∈ N. (5)
Moreover, the next occurrence of i to be performed after the k-th task j is the l + 1-th occurrence. So, we have x i (l + 1) ≥ p j + x j (k). In the same manner of (5), we obtain

x i (s) ≥ p j + x j (s + k -l -1) ∀s ∈ N. (6)
Constraints ( 5) and ( 6) can be rewritten as follows : for all different tasks i, j ∈ T such that M (i) = M (j), we have

∃K ij ∈ Z x i (s) ≥ p j + x j (s -K ij ) ∀s ∈ N (7)
and K ij + K ji = 1. The event shifts K ij and K ji can be interpreted as the ordering on machine M (i). For instance, the selection K ij = 0 and K ji = 1 means that the first occurrence if task i starts after the first occurrence of task j and, since the schedule is periodic, we can state that the k-occurrence if task i starts after the k-occurrence of task j. A solution of the jobshop cyclic scheduling problem is a complete selection of all K ij . These selection constraints, also called disjunctive constraints, defines a selection graph G s (T, S). This graph is non connected and there is one non connected component per machine.

Thus, the jobshop cyclic scheduling problem leads to two graphs G s (T, S) corresponding to selection constraints and G u (T, U ) that characterizes the conjunctive constraints as in the basic cyclic problem. The total graph, denoted G t (T, U ∪S) of the problem can be obtain through merging G s (T, S) and G u (T, U ). If (i, j) ∈ U and (i, j) ∈ S, arc (i, j) is supplied by p i and min(H ij , K ij ) (by selecting the minimum value of H ij and K ij , both constraints are met).

As for the basic cyclic scheduling problem, we introduce the matrix After a complete selection, i.e. a choice for each K ij , the functionning of the system is described by an equation of the form (2). Therefore, a causal representation of the form (4) can be obtained. As in the basic cyclic scheduling problem, the optimal period of the schedule is given by the largest eigenvalue of the system. Therefore the problem can be summarized as follows : find a complete selection that minimizes the largest eigenvalue of the evolution matrix. Let us denote A u the matrix of conjunctive constraints such that x(k) ≥ A u x(k -1) and A s the evolution matrix of selection constraints such that x(k) ≥ A s x(k-1).

H ′ t defined such that H ′ tij = min(H ij , K ij ) if (i, j) ∈ U ∪ S
We now introduce the following property that bounds the optimal period of this problem. Proposition 4. The largest eigenvalue of A u and the largest eigenvalue of A s are lower bounds of the optimal cycle time of the periodic jobshop problem.

Proof : Considering the conjunctive constraints, this lower bound is the optimal solution of the basic cyclic scheduling problem, i.e. without resource constraints. The optimal schedule of a complete selection is given by the largest eigenvalue of matrix A = A u ⊕ A s . Hence, we have

x(k) ≥ A u x(k -1) and consequently x(k) ≥ λ(A u )x(k -1)
where λ(A u ) denotes the largest eigenvalue of A u . The same reasoning holds for A s . However, the computation of λ(A s ) is not so easy because A s is not defined while a complete selection has not been realized. We propose a simple way to compute it. Let us consider the set of tasks {i 1 , i 2 , . . . , i l } performed on a machine. It corresponds to a strongly connected part of A s . The functionning of a sequence on this machine is obviously limited by the overall time of processing of all the tasks {i 1 , i 2 , . . . , i l }, so the cycle time of this machine cannot exceed the sum p i1 +p i2 +. . .+p i l that is the eigenvalue of this component. Thus, the largest eigenvalue of A s is given by maximum of the overall processing times of each machine. It also could be proven that every circuit of A s is of height 1 and the sum of all processing times is an eigenvalue of the component.

2

Example 3. Considering the example 2, we get λ(A u ) = 7 (already computed in example 1). For machine 1 (resp. 2), the overall processing time is 6 (resp.4). So we can state that 7 is a lower bound of the optimal cycle time for this cyclic jobshop problem.

As seen previously, every K ij ∈ Z, so we now propose some bounds on this event shift.

If the basic scheduling problem (i.e. without constraints) of the jobshop scheduling problem is consistent, then H ′ * can be computed (the computation of H ′ * converges since the largest eigenvalue is negative, see [START_REF] Baccelli | Synchronization and Linearity[END_REF]). Therefore, (-H ′ * ) ij , if it is not equal to ε, is the minimum shift event between task i and task j induced by the uniform constraints. More precisely, we have

x i (k) ≥ P ij + x j (k --H ′ * ij ) and so x i (k) ≥ P ij + x j (k + H ′ * ij ), (8)
in which P ij is a sum of several processing times (at least p j ).

Let us consider two tasks i and j such that M (i) = M (j) and suppose that (-H ′ * ) ij = ε, consequently we get (8) and since

x j (k) ≥ p i + x i (k -K ji ), we can state x i (k -H ′ * ij ) ≥ P ij + x j (k), ⇒ x i (k -H ′ * ij ) ≥ P ij + p i + x i (k -K ji ), ⇒ x i (k -H ′ * ij ) > x i (k -K ji ).
As x(k) is monotone, it requires K ji > H ′ * ij and for integrity reason we obtain the following bound

K ji ≥ H ′ * ij + 1. Consequently, we also have a bound for K ij since K ij + K ji = 1.
Example 4. Let us consider tasks 1 and 5 of example 2. We have H ′ * 15 = -2. It leads to K 51 ≥ -1 and K 15 ≤ 2. Moreover H ′ * 51 = 1, and we get K 15 ≥ 2 and K 51 ≤ -1. We can conclude that K 51 = -1 and K 15 = 2. Generally, we only obtain upper bounds and lower bounds for the K ij .

A branch and bound procedure for the jobshop cyclic scheduling problem

In previous section, we have seen that the definition domain of the variables K ij can be reduced by using the uniform constraints. Besides, a lower bound of the optimal cycle time has been computed. A branch and bound procedure (not detailed here) can be realized to solve this kind of problem. The branching is realized on variables K ij . Since, there is no need to branch K ji if K ij is already branched, the depth of the exploration tree is equal to the half of the number of selection variables. For each node, the consistency is checked. If the solution is consistent, we build the causal (max, +)-representation of the system and compute the largest eigenvalue of the evolution matrix. Example 5. For our example, we obtain the following results : K 14 = 1, K 41 = 0, K 15 = 2, K 51 = -1, K 32 = 0, K 23 = 1, K 45 = 1 and K 54 = 0. The lower bounds is reached and we obtain an optimal cycle time of 7. The corresponding schedule can be found in figure 3.

CONCLUDING REMARKS

We presented a general (max, +)-linear model for cyclic scheduling problems. This class of systems leads frequently to a non causal representation in (max, +)-algebra. In light of those observations, we propose a way to obtain a causal representation with recurrence of order 1. However, this method involves strong constraints on the system structure. We first apply it to the basic scheduling problem and then we try to formalize the cyclic jobshop problem with (max, +)-algebra. Bounds are proposed, both for the optimal cycle time and for the selection variables K ij . Future works should refine the method and compare it to mixed integer linear programs of [START_REF] Hanen | Study of a np-hard cyclic scheduling problem: The recurrent job-shop[END_REF] and [START_REF] Brucker | A general model for cyclic machine scheduling problems[END_REF] for larger instances of the problem. We also plan to solve the general case of proposition 1.

One may note that this kind of approach suggests that some performance evaluation could be performed (robustness to delay for instance) and it is not necessary relevant to schedule the minimum cycle time if the solution is not robust. It might be interesting, in some cases, to minimize the deviation of the cycle time when delays occurrs. It is worthwhile to study this connection in more details.

Figure 2 .

 2 Figure 2. Total graph of the cyclic jobshop problem.

Figure 3 .

 3 Figure 3. Optimal solution of the jobshop cyclic scheduling problem.

  and H ′ tij = ε otherwise. Nonetheless not all values for each K ij are possible and we now propose a consistency property. Proposition 3. A selection is consistent if and only if the largest eigenvalue of H ′ t is negative.

	Proof : As for proposition 2, we can build the graph G ′ t and the problem is consistent if and only if the maximum
	mean cycle of this graph is negative.
	2
	Example 2. Let us consider the cyclic scheduling problem
	of example 1. We now consider that tasks 1,4 and 5 are
	performed by machine 1 and tasks 2 and 3 are performed
	by machine 2.

The selection graph (resp. uniform graph) of this problem is depicted in dashed lines (resp. solid line) of figure 2. For more clarity, we don't report in figure

2

arcs from non reentrance constraints). The merge of these two graphs leads to the total graph.