Adina M Panchea
email: panchea@lix.polytechnique.fr

Alexandre Chapoutot
email: chapoutot@ensta-paristech.fr

David Filliat
email: filliat@ensta-paristech.fr

Extended Reliable Robust Motion Planners *

A new method to plan guaranteed to be safe paths in an uncertain environment, with an uncertain initial and final configuration space, while avoiding static obstacles is presented. First, two improved versions of the previously proposed BoxRRT algorithm are presented: both with a better integration scheme and one of them with the control input selected according to a desired objective, and not randomly, as in the original formulation. Second, a new motion planner, called towards BoxRRT*, based on optimal Rapidly-exploring Random Trees algorithm and using interval analysis is introduced. Finally, each of the described algorithms are evaluated on a numerical example. Results show that our algorithms make it possible to find shorter reliable paths with less iterations.

Introduction

The motion planning problem, as addressed in this paper, consists in finding a path, or a sequence of control policy which drives a mobile robot, with a given dynamics description, from a given initial state region to a given goal region while avoiding collisions with a given set of obstacles.

The development of efficient and intelligent motion planning algorithms used for autonomous navigation, which is at the core of autonomous mobile robotic devices, remains a challenge in particular when trying to guarantee the safety of the vehicle. The guarantee of vehicle's safety implies that the motion planning should take into account uncertainties usually resulting from the approximate initial mobile robot localization, from imperfect embedded sensors or from the approximate models used to describe the behaviour of the robotic devices.

Related work

Many motion planning algorithms have been proposed in the literature. When dealing with complex environments, non-holonomic vehicles or high-dimensional state-space, a popular approach is to use stochastic sampling to discretize the configuration space. Among the existing approaches, we focus on the Rapidly-exploring Random Trees (RRT) ( [START_REF] Lavalle | Rapidly-exploring random trees: a new tool for path planning[END_REF][START_REF] Lavalle | Rapidly-exploring random trees: Progress and prospects[END_REF][START_REF] Kuffner | RRT-connect: An efficient approach to singlequery path planning[END_REF][START_REF] Lavalle | Planning Algorithms[END_REF]) path planning algorithm and its many variants. These algorithms have the advantage of rapidly covering the whole configuration space and of easily integrating complex robot models.

However, RRT lacks in guaranteeing an optimal solution. The (asymptotic) optimality of the solution is provided by the optimal Rapidly-exploring Random Trees (RRT*) first proposed in [START_REF] Karaman | Optimal kinodynamic motion planning using incremental sampling-based methods[END_REF] and recently used, for example to plan trajectories for aerial vehicle ( [START_REF] Pharpatara | Shortest path for aerial vehicles in heterogeneous environment using RRT[END_REF][START_REF] Pharpatara | 3-d trajectory planning of aerial vehicles using RRT*[END_REF]). In [START_REF] Noreen | Optimal path planning using RRT* based approaches: A survey and future directions[END_REF], a survey on motion planning algorithms based on RRT * is given.

Uncertainty in path planning has been considered using several representations such as set-membership ( [START_REF] Page | Robot motion planning for sensor-based control with uncertainties[END_REF][START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF]) or covariance matrices ( [START_REF] Lambert | Safe path planning in an uncertain-configuration space[END_REF][START_REF] Pepy | Safe path planning in an uncertain-configuration space using rrt[END_REF][START_REF] Censi | A bayesian framework for optimal motion planning with uncertainty[END_REF]). While the latter is able to find paths with a collision probability under a given threshold, setmembership approaches can guarantee safe trajectories under a bounded noise assumption.

The localization information provided by imperfect proprioceptive sensors, while represented by Gaussian functions ( [START_REF] Lambert | Safe path planning in an uncertain-configuration space[END_REF][START_REF] Pepy | Safe path planning in an uncertain-configuration space using rrt[END_REF]) can guarantee the safety of the path at a certain confidence threshold. Recently, [START_REF] Pepy | Reliable robust path planner[END_REF] provided the guarantee of a safe path to imperfect proprioceptive sensors, while considering the uncertainties bounded with know bounds. Under the latter uncertainty representation, [START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF] introduced a preliminary conceptual reliable and robust path planner based on RRT principles and solved in an interval analysis ( [START_REF] Moore | Interval Analysis[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF]) framework. The interval analysis framework was previously used ( [START_REF] Jaulin | Path planning using intervals and graphs[END_REF]) along with graph algorithms to find collision-free short paths in a given configuration space.

Contributions

This study presents motion planning algorithms which can guarantee safe paths in an uncertain configuration space, where all approximate initial and final mobile robot localisation are bounded with known bounds. Improving the motion planner (denoted BoxRRT) proposed in [START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF], our first contribution consists in an improved BoxRRT planner which makes use of modern and new tools ( [START_REF] Dit Sandretto | Validated explicit and implicit Runge-Kutta methods[END_REF][START_REF] Dit Sandretto | Validated simulation of differential algebraic equations[END_REF]) for the guaranteed numerical integration employed by the motion planner combined with methods coming from constraint satisfaction problems. The second contribution consists in a second improved BoxRRT algorithm which makes use of the guaranteed numerical integration improvements along with the choice of the control input according to a desired objective and not randomly, as already proposed in the literature. As a third contribution, a preliminary attempt towards a new reliable and robust motion planning algorithm based on RRT * principles is introduced.

This paper is organized as follows. First, Sect. 2 introduces the problem formulation, while Sect. 3 provides two improved versions of BoxRRT algorithm. Next, Sect. 4 describes the new planner based on RRT * . The resulting three proposed reliable and robust path planners are applied to plan paths of a non-holonomic vehicle in Sect. 5, where the simulations results are provided. Finally, some concluding remarks and perspectives are drawn in Section 6.

Problem Statement

This paper considers a mobile robot which has to be driven in a two-dimensional static environment from an initial state to a desired one while avoiding obstacles represented by polygons shapes.

The configuration space S = S free ∪ S obs is therefore composed of two subsets: the free region subset S free = S \ S obs where the mobile robot is allowed to move and the obstacle region subset S obs which the mobile robot needs to avoid. Moreover, uncertainties related to its initial and final position and orientation w.r.t. a frame attached to the environment are considered.

Problem formulation

Consider the differential system which can describe the evolution of a mobile robot system:

ṡ(t) = f(s(t), u(t)) (1)
where s ∈ S is considered to be the measurable state of the system, while u(t) ∈ U is the admissible control input. The exact solution of (1) from the inital condition s 0 is denoted by s(t; s 0 ). From an initial state s 0 which belongs to a known set s 0 ∈ S init ⊂ S free the system needs to reach a given set of goal states S goal ⊂ S free . The problem formulation comes from [START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF].

The purpose of the robust motion planner is to provide a sequence of control inputs u ∈ U ∆t

[u] bounded over intervals of time [K∆t, (K + 1)∆t[, with ∆t > 0 and K ∈ N, which will drive the system to reach S goal while avoiding the non-admissible states S obs whatever the initial state s ∈ S init are. If such a sequence of control input u ∈ U ∆t is proved to drive the system from any initial state s ∈ S init to a final state in S goal then the found robust planned path is reliable.

The formulation of such a robust motion planner for which there exists a sequence of control input u ∈ U ∆t

[u] to drive the system from an uncertain initial state to a set of goal states S goal is as follows:

∃K > 0 and u ∈ U such that ∀s 0 ∈ S init , ∀ s(K∆t; s 0 ) ∈ S goal and ∀t ∈ [0, K∆t], s(t; s 0 ) ∈ S free , [START_REF] Dit Sandretto | Validated simulation of differential algebraic equations[END_REF] with s(t) the solution of (1).

Motion Planner Algorithms

This section recalls all necessary notation deployed in this study regarding interval vectors or boxes ( [START_REF] Jaulin | Applied Interval Analysis[END_REF]) which are being used to represent the environment uncertainties.

Next, two new versions of BoxRRT motion planner: the random control input BoxRRT (rciBoxRRT) motion planner algorithm and the selected control input BoxRRT (sciBoxRRT) motion planner algorithm are introduced. Both algorithms are based on RRT motion planner ( [START_REF] Lavalle | Rapidly-exploring random trees: a new tool for path planning[END_REF][START_REF] Lavalle | Planning Algorithms[END_REF][START_REF] Lavalle | Rapidly-exploring random trees: Progress and prospects[END_REF][START_REF] Kuffner | RRT-connect: An efficient approach to singlequery path planning[END_REF]) which is an incremental-based method with the purpose of efficiently explore all the given configuration space from a given starting configuration.

The idea of BoxRRT is not new, being previously proposed by [START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF][START_REF] Pepy | Reliable robust path planner[END_REF] for the case where the uncertainties related to the configuration space are considered only on the final state and not on the initial configuration state. While in this study the interest is to consider uncertainties in the initial and final configuration state space. Moreover, the improvements made on our new versions of BoxRRT are presented in the followings.

Interval analysis

A scalar (real) interval [x] = [x, x] is a closed and connected subset of R, where x represents the lower bound and x represents the upper bound. Two intervals [u] and [v] are equal if and only if u = v and u = v. An interval vector (or box)

[x] is a subset of R n which is the Cartesian product of scalar intervals [x] = [x 1 ] × [x 2 ] × • • • × [x n ],
where the ith component is the projection of [x] onto the ith axis. The interval hull of a set A is the smallest box which contains A, denoted by Hull(A). The inner approximation of a set A, denoted Int(A), is a box included in A, i.e, Int(A) ⊂ A. The Hausdorff distance [START_REF] Jaulin | Applied Interval Analysis[END_REF][START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]

of two intervals [x 1 ] and [x 2 ] is d([x 1 ], [x 2 ]) = sup{|x 1 -x 2 |, |x 1 -x 2 |} (3)
Validated numerical integration methods are interval counterpart of numerical integration methods. A validated numerical integration of a differential equation, as defined in (1) assuming piece-wise constant input, consists in a discretization of time, such that t 0 • • • t end , and a computation of enclosures of the set of states of the system s 0 , . . . , s end , by the help of a guaranteed integration scheme. In details, a guaranteed integration scheme is made of:

• an integration method Φ( f , s j ,t j , h), starting from an initial value s j at time t j and a finite time horizon h (the step-size), producing an approximation s j+1 at time t j+1 = t j +h, of the exact solution s(t j+1 ; s j ), i.e., s(t j+1 ; s j ) ≈ Φ( f , s j ,t j , h);

• a truncation error function lte Φ (d, s j ,t j , h), such that s(t j+1 ; s j ) = Φ( f , s j ,t j , h) + lte Φ ( f , s j ,t j , h).
Our validated numerical integration method is a two step method starting at time t j and for which i) it computes an enclosure [s j ] of the solution of (1) over the time interval [t j ,t j+1 ] to bound lte Φ (d, s j ,t j , h); ii) it computes a tight enclosure of the solution of (1) for the particular time instant t j+1 . There are many methods for these two steps among Taylor series and Runge-Kutta methods see [START_REF] Nedialkov | validated solutions of initial value problems for ordinary differential equations[END_REF][START_REF] Dit Sandretto | Validated explicit and implicit Runge-Kutta methods[END_REF] and the references therein for more details. As a result, validated numerical integration methods produce two functions depending on time

R : R → IR n t → [s] (4)
with for a given t i , R(t i ) = {s(t i ; s 0 ) :

∀s 0 ∈ [s 0 ]} ⊆ [s], and
R : IR → IR n [t,t] → [s] (5) with R([t,t]) = {s(t; s 0 ) : ∀s 0 ∈ [s 0 ] ∧ ∀t ∈ [t,t]} ⊆ [s].

The rciBoxRRT and sciBoxRRT proposed motion planners

Let's start by introducing the global description which is the same for both algorithms.

Next, the algorithm which gathers improvements regarding to the previously proposed BoxRRT algorithm and which has the same formulation for both our new versions of BoxRRT motion planner is introduced. Finally, each procedure of the algorithm is explained separately. The difference between the proposed algorithms lies in the choice of the control input and it will be explained as follows.

Description: First the given initial configuration [s init ] is added to the exploration tree G (Line 1 Steer procedure: The control input used by this procedure is the same as the one proposed for the sciBoxRRT motion planner, i.e. designed according to a desired behaviour or to a chosen criterion.

Near procedure: In this study, the k-nearest neighbors algorithm is employed to determine the set of vertices nearest to the state [s new ], according to the metric d, defined in (3). At each iteration, S near will contain the closest vertices with the metric d([s new ], G) < r. This means that the vertices contained in S near are searched within the area of a ball of radius r(n) = γ log(n), with γ > ε(1 + 1 dim ) or γ = 2ε as suggested in [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF] where ε is Euler's number, n is the number of vertices in the tree at a given iter- Rewire-Parent procedure: If a better parent is found along with a control input which connects it to [s new ], by the NewParent procedure, the Rewire-Parent procedure will update [s new ] parent and cost value. Even though the tBoxRRT * motion planner is build upon the original RRT * planner principles, the rewire procedure, now denoted Rewire-Parent, searches only potential parents and not potential children, as done in the original RRT * planner. Moreover, the original RRT * planner stops the algorithm when a given number of iterations is reached, while as seen in Algo. 2 our proposed planner stops when a solution is found. Even if the choice of the stop criteria can have an impact in proving a near-optimal solution for the proposed algorithm, the latter benefits of different advantages such as guaranteeing the reliability and robustness of the found solution, if exists.

input : [s init ] ⊂ S free , [s goal ] ⊂ S free , MaxIter ∈ N output: G 1 G.init([s init ]); 2 i ← 0 ; 3 repeat 4 [s rand ] ← random-box-GoalBias (S free ); 5 [s nearest ] ← nearest-neighbor(G, [s rand ]); 6 ([s new ], u) ← steer([s nearest ], [s rand ]); 7 if collision-free-path([s nearest ], [s new ]) then 8 cost([s new ]) ← cost([s nearest ]) + d([s nearest ], [s new ]); 9 S near ← near(G, [s new ]); 10 ([s near ], u) ← NewParent(S near , [s nearest ], [s new ]); 11 G ← Rewire-Parent([s near ], [s new ], G);

Numerical Example

The three proposed motion planners are applied on four different scenarios for which the configuration space size is 0.6m ×0.6m×2πrad. The initial state [s init ] size for each environment is: 0.2m × 0.2m × 0.02rad, 0.3m × 0.3m × 0.02rad, 0.4m × 0.4m × 0.02rad, 0.2m × 0.2m × 0.02rad and has to reach the following final state [s goal ] size: 2m × 2m × πrad, 3m × 3m × πrad, 3m × 3m × πrad, 6m × 6m × πrad.

Each algorithm performs 50 iterations for each proposed scenario on a Intel Core m7-6Y75 CPU at 1.20GHz×4. The used software consists in DynIBEX1 [START_REF] Dit Sandretto | DynIBEX: a differential constraint library for studying dynamical systems[END_REF] which is a library providing operators to deal with constraint satisfaction problems embedding differential equations.

Robot mobile modelling

The considered mobile robot is represented by a simple car model, evolving in a 2D environment. The car moves in a configuration s = (x y θ ) with its position (x y) and orientation θ w.r.t. a frame attached to the environment. The simple car model which involves non-holonomic constraints is as follows: While the control input employed by the rciBoxRRT planner is randomly chosen in the admissible set, the one used by sciBoxRRT and tBoxRRT * planners is designed in two steps, as follows. First, the car is oriented towards the target, case in which the control input is obtained by considering the error between the orientation and direction to the goal equal to zero: atan( y target -y start x target -x start )θ target = 0, and using it to simulate (6) for ∆t. Once the car is oriented towards the target a second control input is designed to move the car straight ahead to the target. The steering angle is equal to 0, while the longitudinal speed is obtained from: v = ẋ sin(θ ) . Finally, ( 6) is simulated for ∆t with this last control input.

ẋ = v cos θ , ẏ = v sin θ , θ = v L tan(δ ), (6)

Simulations

All scenarios are performed with ∆t = 1s, a probability p = 0.33 mentioned in Section 3.1 and the maximum limit of iterations fixed to MaxIter = 20.000. The three proposed motion planners are performed on four different environments denoted env i with i = 1..4 and illustrated on Fig. 3. On the same figure a solution found by each algorithm is represented along with the exhibited total number of vertices and the computation time (CPU [s]). When the complexity of the environment increase, as well the algorithm's performances in terms of CPU, number of vertices and distance for the planned path will increase.

Fig. 1 reports the number of iterations necessary for the convergence of each algorithm. We observe that rciBoxRRT which applies a random control input requires the most iterations for convergence, while sciBoxRRT and tBoxRRT * which use a designed control, presented above, converged after less iterations. Fig. 2 illustrates the mean and standard deviation of computational time, number of vertices and length of the planned path, for all simulations performed by the three planners. For all environments, as Fig. 2 stands for, while comparing the planners two different classifications can be made which is the same for all 4 environments: (a) in terms of CPU time, the order of the planners performance enumerated from the more expensive to the less one is: rciBoxRRT, tBoxRRT * and sciBoxRRT; (b) in terms of number of vertices and length of the planned path, the order of the planners performances given in decreasing order is: rciBoxRRT, sciBoxRRT and tBoxRRT * . These results suggest that the two planners for which a control input is designed have better performances than the one in which a random control input is used. Moreover, it was not a surprise to see that tBoxRRT * is more time consuming than sci-BoxRRT while the first one recalls multiples times the steer procedure (Lines 10-11 in Algo. 2) so that better length path performances to be obtained.

Conclusion and Perspectives

Improved versions of the previously proposed boxRRT algorithm and a new motion planner tBoxRRT * based on RRT * are presented in this paper. All motion planners are able to find reliable and robust paths in an uncertain environment, where the uncertain quantities are assumed to belong to boxes.

If the imperfections on the initial states are too large, the imprecision at each new uncertain state can increase. This issue, which will be the topic of future studies, can be encountered by our proposed planners and by the original BoxRRT planner as well. For this reason, in practical settings the motion planners can be updated from time to time by using observers to estimate the state evolution using informations provided by sensors. This can be very useful in decreasing the large imperfections of these uncertain new states. Also, the use of different control inputs between two states can limite the uncertain new states growth.

In the presented version, the tBoxRRT * has a basic form in which not all the RRT * principles are employed (in particular the stop criteria and the rewire procedure). Applying those principles would improve the planned path length value. The proposed planners can be adapted for cases where the free subspace of the configuration space varies with time, to describe moving obstacles. Moreover, model or/and sensorial uncertainties along with a more complex model of the mobile robot which takes into account its dynamics can be considered. This should form the subject of future studies.

2 π 2 ]

 22 where the control input u = [v δ ] is represented by the longitudinal speed v ∈ [-1 1] and the steering angle δ ∈ [-π . L = 1.5 represents the distance between the front and back axes of the car.

Figure 1 :

 1 Figure 1: Number of iteration for convergence required by each of the three proposed algorithms.

Figure 2 :

 2 Figure 2: The rciBoxRRT, sciBoxRRT and tBoxRRT * path planner found performances for each of the four proposed scenarios.

Figure 3 :

 3 Figure 3: (a) -(d) rciBoxRRT ((a) total vertices 2200 in 28 [s]; (b) total vertices 5880 in 103 [s]; (c) total vertices 3416 in 51 [s]; (d) total vertices 7802 in 141[s]). (e) -(h) sciBoxRRT( (e) total vertices 570 in 11 [s]; (f) total vertices 1149 in 32 [s]; (g) total vertices 278 in 5[s]; (h) total vertices 978 in 26[s]). (i) -(l) tBoxRRT * ( (i) total vertices 156 in 3 [s]; (j) total vertices 1088 in 38 [s]; (k) total vertices 786 in 20 [s]; (l) total vertices 963 in 28[s]).

  ). Then, a state [s rand ] ⊂ S free is randomly chosen by the procedure randombox-GoalBias (Line 4). The nearest-neighbor procedure from Line 5 returns the closest vertex [s near ] to [s rand ] in the tree G, according to a certain metric d. A control input u ∈ [u] is chosen according to a specified criterion or randomly through the select input procedure. Then, in the prediction procedure, (1) is integrated over a fix time interval ∆t with the initial condition [s near ] and a constant control input u (given atLine 6) and will result in a new state [s new ] (Line 7). If it can be proved that all state values along the trajectory between [s near ] and [s new ] lie in S free being a collision free path, then the path between [s near ] and [s new ] is considered reliable and [s new ] is added to G as a new vertex and connected to its parent [s near ] though the G.add-guaranteed-vertex procedure. Otherwise, [s new ] is not added to G. Lines 4 to 11 are repeated until a chosen number of iterations MaxIter is reached or until a path is found meaning [s new ] = [s goal ], or most likely when [s new ] ⊂ [s goal ]. Note that we have [s init ] = Hull(S init ), [s obs ] = Hull(S obs ) and [s goal ] = Int(S goal ) to ensure the soundness of the proposed algorithm. Other techniques towards a more improved random procedure can be thought of, such as the use of the artificial potential field algorithms (APFs) as proposed in[START_REF] Pharpatara | Shortest path for aerial vehicles in heterogeneous environment using RRT[END_REF][START_REF] Pharpatara | 3-d trajectory planning of aerial vehicles using RRT*[END_REF]. Nearest neighbor procedure: Finds the closest vertex to the [s rand ] one according to a chosen metric d. The Hausdorff distance between two intervals as defined in (3) is considered. Prediction procedure: Finds a new state [s new ] while integrating (1) with the selected control input, given by the select input procedure, over an interval of time ∆t. This step is based on validated numerical integration methods as explained in Section 3.1 and using function R(t). Collision free path procedure: If [s init ] and [s goal ] are, respectively, the imperfect initial and final states, one has to show before starting the path planner that both sets of states belong to S free . When it is proved that no collision occurs between any two consecutive vertices of the tree, one proves by induction that the path between [s init ] and [s goal ] is robustly reliable, if it exists. The techniques used in this procedure are based on new tool and functions proposed by[START_REF] Dit Sandretto | Formal Verification of Robotic Behaviors in Presence of Bounded Uncertainties[END_REF], which are capable of testing during the integration procedure if a collision occurred. Therefore this procedure differs from the previously BoxRRT one which uses wrap techniques[START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF][START_REF] Pepy | Reliable robust path planner[END_REF]. More precisely, using the enclosure R(t) of the trajectory of (1), checking that no collision occurs is simply an interval test which checks if R(t) does not intersect [s obs ] for all t. Select input procedure: This procedure is used to find a control input which finds a new state starting from a given initial state. The difference between the two new versions of BoxRRT is made in this procedure: the rciBoxRRT motion planner uses a control input chosen randomly among the set of admissible values u ∈ U, while sciBoxRRT motion planner uses a designed control input according to a desired behaviour or to a chosen criterion (see Section 5 for an example).

	4 Towards BoxRRT* Motion Planner Algorithm

input : [s init ] ⊂ S free , [s goal ] ⊂ S free , ∆t ∈ R + , MaxIter ∈ N output: G 1 G.init([s init ]); 2 i ← 0; 3 repeat 4 [s rand ] ← random-box-GoalBias(S free ); 5 [s near ] ← nearest-neighbor(G, [s rand ]); 6 u ← select input([s rand ], [s near ]); 7 [s new ] ← prediction([s near ], u, ∆t); 8 if collision free path ([s near ], [s new ], u, ∆t) then 9 G.add-guaranteed-vertex([s near ], [s new ], u); 10 return [s new ] 11 return / 0 12 until i + + < MaxIter or ([s new ] = / 0 and [s new ] ⊂ [s goal ]); 13 return G Algorithm 1: BoxRRT motion planner algorithm Random box GoalBias procedure: This procedure, previously proposed in [21], consists in choosing the random state in the final configuration state [s rand ] ⊂ [s goal ] with a probability p > 0. This section introduces a new reliable robust path planner for uncertain environments based on optimal RRTs (RRT * )([8, 24, 23]), which is denoted tBoxRRT * . A general description of the algorithm is presented, followed by the description of the used procedures. Description: As in RRT, the tree G is initialized with the given initial configuration [s init ]. Then, a state [s rand ] ∈ S free is randomly chosen by random-box-GoalBias procedure and its nearest vertex [s nearest ] according to a defined metric d is provided by the nearest-neighbor procedure. Steer procedure designs a control input according to a desired behaviour or according to a specific criterion. Eq. (1) is integrated over a fixed time interval ∆t with the initial condition [s nearest ] and a constant control input u. Then a new state [s new ] is found. If the trajectory between [s nearest ] and [s new ] lie in S free , then the path between [s nearest ] and [s new ] is reliable and [s new ] is a new vertex added to G, with the cost function (cost([s new ])) associated with the distance from [s init ] to [s new ] through its parent [s nearest ]. Next, the near procedure checks if a better parent for [s new ] can be found. Therefore, a list of potential vertices in a neighborhood [s near ] ∈ S near of [s new ] is selected. For each vertex ([s near ]) from the list of potential parents S near is checked if the cost (according to the distance metric as defined in (3)) to arrive in [s new ] through [s near ] is better than cost([s new ]). If it is the case and the path is collision free then the rewire-parent procedure will update [s new ] parent with [s near ] and its cost([s new ]) accordingly. These steps (Lines 4 to 13) are repeated until the algorithm reaches the MaxIter iterations or until [s new ] = [s goal ], or more likely when [s new ] ⊂ [s goal ]. Note that here as in Subsect. 3.2 [s init ] = Hull(S init ), [s obs ] = Hull(S obs ) and [s goal ] = Int(S goal ). Random-box-GoalBias, Nearest-neighbor and Collision-free-path are already described in Section 3.2.

  MaxIter) or ([s new ] = / 0 and [s new ] ⊂ [s goal ]);

	12	return [s new ]
	13	return / 0
	14 until (i + + < 15 return G
		Algorithm 2: Towards BoxRRT

* motion planning algorithm ation and dim represents the dimension of the configuration space. NewParent procedure: When a vertex is added in G its cost is defined as: cost([

s new ]) = cost([s nearest ]) + d([s nearest ], [s new ]),

where cost([s nearest ]) represents the distance from the initial state ([s init ]) to the vertex initial parent [s nearest ]. This procedure verifies if among the vertices in S near a better parent can be found. For each vertex [s near ] ∈ S near it is checked if the total cost to arrive to [s new ], passing through [s near ] is smaller than cost([s new ]). When a better parent for [s new ] is found, the steer procedure is applied from [s near ] to [s new ]. The control input applied can drive [s near ]: (a) to [s new ] state, (b) inside [s new ] state or (c) as close as possible and with as small cost to [s new ].

http://perso.ensta-paristech.fr/∼chapoutot/dynibex/

Acknowledgment

The authors would like to thank Oliver Mullier, Julien Alexandre dit Sandretto, Eric Goubault and Benjamin Martin for useful comments and discussions.

* This work was supported by DGA MRIS.


A Álvarez
email: aalvarez@el.ub.edu

M Bafleur

J-M Dilhac

J Colomer

D Dragomirescu

J Lopez

M Zhu

P Miribel

Self-Powered energy harvester strain sensing device for structural health monitoring

This paper presents an envisaged autonomous strain sensor device, which is dedicated to structural health monitoring applications. The paper introduces the ASIC approach that replaces the discrete approach of some of the main modules.

Introduction

The proposed device is developed in the frame of the SMARTER (Smart Multifunctional Architecture & Technology for Energy-Aware wireless sensor) Project [START_REF] Bafleur | Proc. IEEE International Conference on Green Computing and Communications and IEEE Internet of Things and IEEE Cyber, Physical and Social Computing[END_REF].

The device is self-powered, based on a piezoelectric transducer as an energy-harvesting source, enabling a battery-free solution. The system would be capable to sense the mechanical strain, which is an interesting topic in structural health monitoring applications, particularly in the aerospace field, which is the present field of application, defining the frequencies available for the transducer, harvesting the energy provided by the in-flight vibration of the aircraft's wing [START_REF] Giuliano | Proc. IEEE International Conference on Green Computing and Communications and IEEE Internet of Things and IEEE Cyber, Physical and Social Computing[END_REF].

The nice approach of the device, which defines its Multifunctional concept, is that the same transducer used to harvest energy also allows extracting the mechanical strain. On the one hand, based on a piezoelectric transducer, the mechanical vibrations are harvested to generate electrical energy, in order to start-up and supply the system, but on the other hand, the system is capable to extract the open voltage (VOC) of the piezoelectric harvester, which has a relationship with the mechanical strain [START_REF] Chen | [END_REF]. Then, based on the open circuit voltage concept to adapt the maximum power transferred condition between the piezoelectric mechanical to electrical energy [4,5,6], a detection of the open voltage of the piezoelectric transducer is derived. Based on these approaches, initial fully discrete electronic solutions have been implemented [START_REF] Colomer-Farrarons | Proc. Conference on Design Circuits and Integrated Systems[END_REF][8] as a proof of concept.

The system is based on four main modules, figure 1: a) the energy transducer as a multifunctional device (MFC) based on a piezoelectric generator (PZT); b) the Power Management, Control and Sensing Circuitry Unit, to be implemented in an ASIC solution (ASIC1); c) the Dynamic Storage Module to store the generated energy (SCap), and d) the Communication (RF) and Data Processing Unit (ASIC2). Full electronics that define the Smarter Device would be placed on a flexible PCB.

In this paper we focus the work regarding the design of ASIC1. This implementation has the role to manage the power harvested, the measurements and the management of the other modules. In section 2 is introduced the ASIC1 approach. The initial discrete solutions [START_REF] Colomer-Farrarons | Proc. Conference on Design Circuits and Integrated Systems[END_REF][8] are translated to a microelectronic one, for a more compact solution.

Finally in section 3 are derived some conclusions and future perspectives.

Integrated approach

Thanks to the discrete approaches [START_REF] Colomer-Farrarons | Proc. Conference on Design Circuits and Integrated Systems[END_REF][START_REF] Colomer-Farrarons | Proc. IEEE International Symposiun on Industrial Electronics[END_REF], the contemplate architecture is validated and with the gained experience we proceeded with the ASIC integration. A first ASIC is going to be implemented based on an AMS 0.35 technology. In figure 2 is depicted the ASIC1 at the functional block level and its layout design in figure 3, in an AMS 0.35µm technology, that will be placed in a QFN32 package. To interface ASIC1 and ASIC2, figure 1,a communication protocol has been implemented, working with 8 bits: Two bits of control, and 6-bits of data. This protocol is implemented in [START_REF] Colomer-Farrarons | Proc. IEEE International Symposiun on Industrial Electronics[END_REF] and validated. In t1 data is available on ASIC1. In t2 ASIC activates Start Bit (SB) when data is ready of transmitting. In ta SB is activated during all operation. In t3 ASIC2 is ready to operate. In tb, the bit Transmission Done (TD) is activated during UR-UWB transmission. In t4 TD is grounded when the transmission is completed. In t5 SB is turned off after communication is done. Finally, in t6 the system waits till new data is available on ASIC1.

The system is focused into the extraction of the strain through VOC, and the validation of the architecture that will combine the electronics envisaged in ASIC1 with those of ASIC2. The wireless solution implemented in [START_REF] Colomer-Farrarons | Proc. Conference on Design Circuits and Integrated Systems[END_REF] will be replaced by ASIC2 in the future Smarter prototype, and ASIC1 has a Data Processing Unit, based on a 6-bits SAR-ADC converter, which prepares the data to be sent by the RF module (ASIC2) replacing the TI Microcontroller in [START_REF] Colomer-Farrarons | Proc. IEEE International Symposiun on Industrial Electronics[END_REF]. The strain is measured and related with the open circuit voltage, and the measurement is outputted by 6-parallel bit ADC, which translates the measurement to the Communication. The main blocs are defined by the Power Circuit Unit, which combines linear regulators and the Energy Management module. The system needs a 1.1V regulated voltage source and a 3.3V regulated voltage source. From the AC/DC rectifier two LDOs will generate these voltages. The 1.1V LDO is designed to bias ASIC2, and the 3.3V LDO to bias the rest of the electronics. Then, the Energy Management module has the role to control the Dynamic Storage Element or SCap [START_REF] Mahboubi | Proc. Powermens 2016[END_REF]. Finally, the ACU module and the 6-bit SAR ADC Converter define the Control Circuit and Sensing Circuitry. The Power Circuit Unit is described more in detail. A specific low-power bandgap reference and current reference circuit has been designed; with a current consumption lower than 9µA, nominally and a nominal power consumption of 30µW. It has the capability to work from a lower voltage of 2V till 5V, the maximum operating voltage, with an output voltage of 1.23V with a Temperature Coefficient (TC) lower than 60ppm. The current reference is 105nA with a TC lower than 250ppm, in the temperature range of -40ºC to 125 ºC.

Then, based on the architectures presented in [5], the LDO circuits have been designed. The 3.3V LDO is designed to work with a current up to 1mA, and the 1.1V LDO with a maximum DC current of 10mA. The main parameters for the LDO are introduced. In the case of the 3.3V LDO the Linear Regulation is 68.7•10 -3 (V/V). In the case of the 1.1V LDO is 665•10 -9 (V/V). The ability to maintain the output voltage with variations of the input voltage is reasonable for both designs. The measurement of the capability of the circuitry to keep the regulated output voltage to the desired value against changes of the load, that is, between light and high load conditions (in our case between few µA that define the quiescent operating point of the electronics and the load condition of 1mA for the 3.3V LDO and 10mA for the 1.1V LDO), is the load regulation, which is 81.6•10 -3 (V/A) in the case of the 3.3V LDO and 99•3•10 -6 (V/A) for the 1.1V LDO regulator.

Full Monte Carlo analyses have been carried out for the bandgap reference circuit and LDOs. In the particular case of the Linear Dropout Regulators, it is the 1.1V LDO that has the greater limitation, with an accepted range between 1.05 and 1.18V, with 200 samples. The average regulated voltage is 1.12V with a standard deviation of 34,19 mV (Yield of 95%).

The ACU module, which is quite complex, and the Energy management module, are out of the scope of the present work. In figure 5 is depicted the caption of the simulation, at the extracted ASIC view level, of the performance of the ACU control stage though the GATE signal (VGATE) [7][8].

Conclusions

This paper presents the ASIC approach for the power management, control and sensing for the SMARTER project. Future full test and characterization will be carried out before its implementation in the final prototype on flexible substrate and reported in future publications.

Figure 1 .

 1 Figure 1. Main blocs for the Smarter Device.

Figure 2 .

 2 Figure 2. ASIC1 envisaged version and implementation.

  Figure 3. ASIC1 layout implementation.

Figure 4 .

 4 Figure 4. Communication protocol between ASICS.

Figure 5 .

 5 Figure 5. Caption of the transient simulation, from the ASIC1 extracted view, of the GATE control signal. Generation of the scenario operating phases. From the Top view to the bottom caption: (Top) Start-Up and harvest/measuring phases; VGATE signal that control the phases; Bandgap reference voltage; 3V3 LDO regulated voltage; (bottom) 1V1 LDO regulated voltage.

Acknowledgements

The SMARTER project is supported by European Union under the ERA-Net funding scheme of the FP7 (CHIST-ERA), and in the case of the spaniard partner through MINECO, reference PCIN-2013-069.


Sergey Kovalev
email: skovalev@inseec.com

Xavier Delorme
email: delorme@emse.fr

Alexandre Dolgui
email: alexandre.dolgui@mines-nantes.fr

Ammar Oulamara
email: ammar.oulamara@loria.fr

Minimizing the number of stations and station activation costs for a production line

Keywords: Line balancing, Scheduling, Computational complexity, Integer linear programming, Heuristics

In the problem under study, a paced unidirectional machining line, consisting of a number of stations, has to be configured to produce parts of several types. A given set of operations is required for each part type and the same operation can be required for different part types. Re-assignment of operations, when switching from one part type to another, is not allowed. All operations assigned to the same station are performed simultaneously. The objective is to assign operations to stations in order to minimize the number of stations and the station activation costs, with respect to precedence and zoning constraints. The two objectives are considered in a lexicographic order, the former being the primary objective. Activation costs refer to the costs induced by the energy consumption, equipment maintenance, setup activities or labor requirement which occur whenever a station is used. Computational complexity for various special cases is established. Heuristic algorithms, integer linear programming formulations, and computer experiments are presented. Instances of practical dimension, with 40 to 80 operations, are solved in an hour on a conventional computer.

Introduction and literature review

We study a line balancing problem whose peculiarity is the consideration of the station activation costs in addition to the usual minimization of the number of stations. The line to be configured is intended for a large series production of several types of parts. It is paced, unidirectional and consists of a number of stations. Parts move along the stations one after another and, at each time instant, any station is occupied by at most one part. In the real production case, which motivates our studies, each station is equipped with a single multi-spindle head where rotating tools such as borers, milling cutters, grinding heads and chamfering mills are mounted to perform operations which are required for all part types. The line balancing decision specifies an assignment of the required operations to the stations. When a blank or a partially processed part of a certain type arrives at a station and there are operations to be performed on this station for this part type, the station is activated and the required operations are performed in parallel.

The station activation includes cleaning the working zone, loading, positioning and unloading of the part. It induces a cost associated with energy, equipment maintenance and labor, see Borgia et al. [START_REF] Borgia | STEP-NC compliant approach for setup planning problem on multiple fixture pallets[END_REF]. Over a long production period, the total station activation cost for a given part type is proportional to the quantity of parts of this type to be produced, or its relative quantity among all types, multiplied by the number of stations which have some assigned operations on parts of this type. The costs associated with characteristics of operations other than the type, such as the number of required tools, sum up to a constant and they do not depend on the line balancing decision. Therefore, the station activation costs are assumed to be dependent only on the part type.

For the type of lines which has motivated this study, the multi-spindle heads are very expensive which often leads decision-makers to focus primarily on the minimization of investment costs. The primary criterion for the line design is thus to minimize the number of stations. However, there are often many solutions with the same minimum number of stations. Therefore, we consider a secondary objective, which is to minimize the total station activation cost.

The line takt time can be common for all types or it can be type dependent. If it is type dependent, then, due to the parallel execution of operations, it is equal to the maximum operation time plus the station activation time, among parts of a given type. If it is common for all types, then it is equal to the maximum operation time plus the corresponding station activation time, among parts of all types. Thus, the line takt time depends on the input data of the problem and it does not depend on the line balancing and part sequencing decisions. The studied lines can be used in a mixed-model mode, where parts of different types are produced interchangeably, see Figure 1, or in a multi-model mode, where parts are produced in batches of the same type, see Figure 2.

A formal definition of the problem under study can be given as follows. A paced production line consisting of several stations has to be configured to produce parts of f types. Parts move along the stations in the same direction in a given sequence. Let F = {1, . . . , v, . . . , f } be the set of part types. Each part of type v ∈ F requires each operation of a given set N v to be executed exactly once on the line. Operations of the set N v are called type v operations. The same operation can be performed on parts of different types, thus, different sets N v can contain common operations. Let N := ∪ f v=1 N v = {1, . . . , n} denote the superset of all the required operations, and let T i denote the set of types of the operation i ∈ N , i.e., T i = {v | i ∈ N v , v ∈ F }. Redesign, i.e., re-assignment of operations when switching from the production of one type part to another, is not allowed. Operations of the same type assigned to the same station are performed simultaneously. Examples of simultaneous execution of operations in the mechanical industry are given by Dolgui et al. [START_REF] Dolgui | Optimisation of multi-position machines and transfer lines[END_REF], Belmokhtar et al. [START_REF] Belmokhtar | Integer programming models for logical layout design of modular machining lines[END_REF], Battaia and Dolgui [START_REF] Battaia | Reduction approaches for a generalized line balancing problem[END_REF] and Falkenauer [START_REF] Falkenauer | Line balancing in the real world[END_REF].

Binary precedence relations are given on the superset N . If operation i ∈ N precedes operation j ∈ N , then j cannot be assigned to the station of i or any preceding station. Precedence relations are transitive and irreflexive. They are represented by a directed acyclic graph G = G(N, A), in which there is an arc (i, j) ∈ A if and only if i precedes j. Precedence relations characterize the technological process. They can be defined separately for each part type, but since re-design is not allowed, they must be consistent and can therefore be given on the superset N .

Plural exclusion relations are defined on the superset N . They are represented by a collection E of sets E ⊂ N such that all operations of E cannot be assigned to the same station, but any proper subset of E can be assigned to the same station. If the simultaneous utilization of some tools at the same station is not possible, exclusion relations prevent this from happening. For example, a set E = {1, 2, 3} of three drilling operations cannot be performed simultaneously, because the total weight of the corresponding drilling heads is unacceptable for the same station.

However, operations of any proper subset of E , which may include one or two arbitrary operations from E , can be assigned to the same station. The total weight of one or two drilling heads associated with operations from E is acceptable. Only the installation of all three drilling heads at the same station is not acceptable. Without loss of generality it can be assumed that there are no two sets in E containing one another. For example, if there are two exclusion sets {1, 2, 3} and {1, 3}, then {1, 2, 3} can be removed from the input since the fact that operations 1 and 3 will never be assigned to the same station implies that 1, 2 and 3 will do the same. Exclusion relations can be defined separately for each part type. However, since re-assignment of operations when switching from the production of one type part to another, is not allowed, they must be consistent and can therefore be given on the superset N .

Each operation i has its size s i which can be viewed as the number of tools required simultaneously to perform this operation. The total size of all operations, i.e., the total number of tools assigned to the same station should not exceed the given station capacity r. Since each tool is associated with an elementary operation, each operation i can be viewed as a block of elementary operations to be executed simultaneously on the same station, with s i as the cardinality of this block and N as a set of blocks. Therefore, operations with non-unit sizes are similar to the inclusion sets such that all elementary operations of the same inclusion set must be assigned to the same station. For example, a set I = {4, 5} of two drilling operations can be required to be performed on the same station because the distance between the corresponding holes must be precise. The required precision will be lost if the drilling operations are performed on different stations. Therefore, set I can be viewed as an operation of size two. See Dolgui et al. [START_REF] Dolgui | MIP approach to balancing transfer lines with blocks of parallel operations[END_REF] for an inclusion sets formulation of a transfer line balancing problem.

Exclusion and inclusion relations refer to what is called zoning constraints in line balancing, see Akpinar and Mirac Bayhan [START_REF] Akpinar | A hybrid genetic algorithm for mixed model assembly line balancing problem with parallel workstations and zoning constraints[END_REF] and Falkenauer [START_REF] Falkenauer | Line balancing in the real world[END_REF]. Zoning constraints are called positive if a specific set of operations must be assigned either to the same station or to a specific set of stations. Negative zoning constraints forbid operations of a specific set to be assigned either to the same station or to a specific set of stations. In our case, exclusion relations correspond to negative zoning constraints while inclusion relations represented by operation sizes refer to positive zoning constraints.

A station is activated for part processing if at least one operation on this part will be processed on this station. An activation cost a v is associated with part type v, v = 1, . . . , f . Let x v denote the number of stations activated for any single part of type v, v = 1, . . . , f . A decision has to be made about the total number of stations, k, and an assignment of the operations to the stations 1, . . . , k. The primary criterion is to minimize the total number of stations, and the secondary criterion is to minimize the total activation cost, a 1 x 1 + • • • + a f x f . We denote the problem with the single primary criterion as P (prec, excl, s i | k) and the problem with the two lexicographically ordered criteria as P (prec, excl, s i | lex(k, cost)). Note that the two criteria are not redundant in the sense that minimizing the total activation cost is not sufficient to minimize the number of stations and vice versa. Let (k * , c * ) denote the lexicographical minimum of the objective function values for the problem P (prec, excl, s i | lex(k, cost)). Thus, k * is the global minimum of the number of stations, and c * is the minimum cost, provided that the number of stations is k * .

For an illustration of the effect of station activations, assume that there are three stations, each having capacity 2, and two part types. Any part of type 1 requires unit size operations 1, 2 and 3 and any part of type 2 requires unit size operations 4, 5 and 6. Assume that the other constraints are not effective. Two feasible assignments of the six operations are given in Figure 3. to the total station activation cost, which is equal to 2a 1 + 2a 2 for the assignment 1. Now consider assignment 2. When a part of type 1 is processed, station 1 performs operation 1, station 2 performs operation 2 and station 3 performs operation 3. This contributes 3a 1 to the total station activation cost. When a part of type 2 is processed, station 2 performs operation 4, station 2 performs operation 5 and station 3 performs operation 6. This contributes 3a 2 to the total station activation cost, which is equal to 3a 1 + 3a 2 for the assignment 2. Boysen et al. [START_REF] Boysen | Resequencing of mixed-model assembly lines: Survey and research agenda[END_REF], Battaia and Dolgui [START_REF] Battaia | A taxonomy of line balancing problems and their solution approaches[END_REF] and Sivasankaran and Shahabudeen [START_REF] Sivasankaran | Literature review of assembly line balancing problems[END_REF]. Results on multi-objective line balancing problems can be found in Gamberini et al. [START_REF] Gamberini | A new multi-objective heuristic algorithm for solving the stochastic assembly line re-balancing problem[END_REF], McMullen and Tarasewich [START_REF] Mcmullen | A beam search heuristic method for mixed-model scheduling with setups[END_REF][START_REF] Mcmullen | Multi-objective assembly line balancing via a modified ant colony optimization technique[END_REF], Suwannarongsri and Puangdownreong [START_REF] Suwannarongsri | Multi-objective assembly line balancing via adaptive tabu search method with partial random permutation technique[END_REF], Ramezanian and Ezzatpanah [START_REF] Ramezanian | Modeling and solving multi-objective mixed-model assembly line balancing and worker assignment problem[END_REF], Chica et al. [START_REF] Chica | A multiobjective model and evolutionary algorithms for robust time and space assembly line balancing under uncertain demand[END_REF]. Setup time and setup cost models are well known in scheduling and lot-sizing research, see, for example, Potts and Kovalyov [START_REF] Potts | Scheduling with batching: a review[END_REF] and Dolgui et al. [START_REF] Dolgui | Multi-product lot-sizing and scheduling on unrelated parallel machines[END_REF].

The rest of the paper is organized as follows. In Section 2, we focus on the problem with the sole primary criterion. The computational complexity of this problem, which is NP-hard in the strong sense in the general case, is analyzed. Two greedy randomized heuristics and an integer linear programming (ILP) formulation are suggested. In Section 3, the lexicographic bi-criteria problem is studied. Again, the computational complexity of this problem is analyzed, and two greedy randomized heuristics as well as an ILP formulation are suggested. Computer experiments with the suggested heuristics and ILP formulations are given in Section 4. A solution procedure for the bi-criteria problem, which dichotomically enumerates the number of stations, is proposed in Section 5. The paper concludes with a summary of the results.

Minimizing the number of stations

In this section, the single criterion problem P (prec, excl, s i | k) to minimize the number of stations is studied. The simplest case with unit size operations and no exclusion or precedence relations is easily solvable: k * = n r . In the following subsections, more complex cases are considered.

Given precedence relations, no exclusion relation and unit sizes

The special case in this subsection is denoted as P (prec,

s i = 1 | k).
We will show that it is equivalent to the scheduling problem P m

| p i = 1, prec | C max . Problem P m | p i = 1, prec | C max :
There are n unit-time jobs to be scheduled on m identical parallel machines. Precedence relations are given on the set of jobs such that if job i precedes job j then j cannot be processed in the unit-time processing interval of i or any preceding unit-time interval. No preemption of job processing is allowed and no two jobs can be processed by the same machine simultaneously. The objective is to assign jobs to the unit-time intervals on the machines so that the precedence relations are satisfied and the completion time of the latest job,

C
p i = 1, prec | C max is NP-hard
in the strong sense for a variable number of machines m (Ullman [START_REF] Ullman | NP-complete scheduling problems[END_REF]), and its computational complexity status is open if m is a constant greater than two (Brucker and Knust [12]). Problem (Hu [31]) or a collection of out-trees (Davida and Linton [START_REF] Davida | A new algorithm for the scheduling of tree structured tasks[END_REF]), and it is solvable in O(n 2 ) time, if G is arbitrary and m = 2 (the best known algorithm is by Gabow [START_REF] Gabow | An almost linear algorithm for two processor scheduling[END_REF]).

P m | p i = 1, prec | C max is solvable in O(n) time if the precedence graph, G, is a collection of in-trees
Thus, problem P (prec, s i = 1 | k) is NP-hard in the strong sense for variable station capacity bound r and arbitrary precedence graph G. It is solvable in O(n) time if G is a collection of in-trees or a collection of out-trees, and in O(n 2 ) time if G is arbitrary and r = 2. The chain-like and tree-like precedence relations are typical for machining and assembly operations.

Arbitrary operation sizes and no exclusion or precedence relation

Denote the special case in this subsection as P (s i | k). The problem P (s i | k) can be interpreted as the well-known bin packing problem, in which there are n items of sizes s 1 , . . . , s n to be packed into the minimum number of bins, each having capacity r. The minimum number of stations k * is equal to the minimum number of the required bins.

The bin packing problem is NP-hard in the strong sense. It is solvable in O(n K L ) time if there is at most K distinct item sizes and at most L items can fit in one bin; see, for example, Fernandez de la Vega and Lueker [START_REF] Fernandez De La Vega | Bin Packing can be solved within (1 + ε) in linear time[END_REF]. If r is a constant, then K and L are constants as well, and the bin packing problem is polynomially solvable. From the above discussion we deduce that the problem P (s i | k) is NP-hard in the strong sense. It is polynomially solvable if the number of operations is a constant or if the station capacity r is a constant.

Given exclusion relations, no precedence relation and unit sizes

Denote the special case in this subsection as P (excl, Any Y i -operation excludes any Y q -operation, i.e., {j, l} ∈ E for every Y i -operation j and every Y q -operation l, i = q. The A i -operation and the A q -operation exclude each other if subsets A i and A q share the same element, i.e., if

s i = 1 | k).
A i ∩ A q = φ, i = q.
We will show that for the constructed instance, a solution to X3C exists if and only if there exists a feasible solution to this instance of the problem P (excl,

s i = 1 | k) with at most a -b + 1
stations. It is obvious that the instance construction is pseudo-polynomial.

Part "if ". Assume that there exists a feasible solution to the constructed instance of the problem P (excl, s i = 1 | k) with at most a -b + 1 stations. Observe that there are a -b + 1 mutually exclusive operations: the X-operation and one Y i -operation for each i = 1, . . . , a -b.

Hence, there must be exactly a -b + 1 stations. Without loss of generality assume that an Y ioperation is assigned to station i, i = 1, . . . , a -b, and the X-operation is assigned to the station a -b + 1. Since any Y i -operation excludes any Y q -operation for i = q, and there are a -b + 1 stations, all the Y i -operations must be assigned to station i, i = 1, . . . , a -b. Therefore, the A i -operation, i ∈ {1, . . . , a}, can be assigned to one of the remaining b positions of the station a-b+1 or to the single position of any of the stations 1, . . . , a-b. The structure of such a solution is given in Table 1 where columns represent stations and symbols "•" represent A i -operations.

Y 1 Y 2 • • • • Y a-b X • Y 2 Y 3 • • • Y a-b • Y 1 • Y 3 • • • • • Y 1 Y 2 Y 3 • • • Y a-b • • • • • • • • • • • • • • • • • • • Y 1 Y 2 Y 3 • • • Y a-b •
The b number of A i -operations at the station a -b + 1 must be mutually non-exclusive, which means that no two subsets A i and A q corresponding to these operations must share the same element. We deduce that the subsets A i corresponding to the A i -operations assigned to station a -b + 1 constitute an exact cover of B, as required for the proof of the part "if".

Part "only if ". If there exists a solution S to the problem X3C, then construct a solution to the problem P (excl, s i = 1 | k) with the structure given in Table 1 and assign the A i -operation to the station a -b + 1 if and only if A i ∈ S. Such a solution is feasible and the number of stations is equal to a -b + 1, as required for the proof of the part "only if".

ILP formulation for the general case of P (prec, excl, s i | k)

Let us now formulate this problem as an Integer Linear Program (ILP). In this model, we assume that an upper bound k 0 ≥ k * on the number of stations is given, i.e. there are k 0 stations to be opened or not. In practice, the number k 0 can be obtained by heuristics (see section 2.5). A station is called open if at least one operation is assigned to it.

Introduce 0-1 variables x ij such that x ij = 1 if and only if operation i is assigned to station j, and 0-1 variables y j such that y j = 1 if and only if station j is open. Recall that the directed graph G = G(N, A) represents the precedence constraints.

Calculate the following lower bound LB 1 on the minimum number of stations:

LB 1 = max l max , n i=1 s i r
, where l max is the number of vertices in the longest chain of G.

The problem P (prec, excl, s i | k) can be re-formulated as the following ILP problem, which we denote as M in#.

Problem M in#:

min k 0 j=1 y j , (1)
subject to

k 0 j=1 y j ≥ LB 1 , (2)
j i=1 y i ≥ jy j , j = 2, . . . , k 0 , (3) n i=1
x ij ≥ y j , j = 1, . . . , k 0 , (4)

k 0 j=1 x ij = 1, i = 1, . . . , n, (5)
n i=1 s i x ij ≤ ry j , j = 1, . . . , k 0 , (6)
i∈E x ij ≤ |E | -1, ∀E ∈ E, j = 1, . . . , k 0 , (7)
h j=1 x bj + k 0 j=h x aj ≤ 1, ∀(a, b) ∈ A, h = 1, . . . , k 0 , (8)
x ij , y j ∈ {0, 1}, i = 1, . . . , n, j = 1, . . . , k 0 .

In the problem M in#, there are k 0 + nk 0 variables and n + k 0 (3 + |E| + |A|) constraints without the 0-1 constraints. Barring the exclusion constraints [START_REF] Borisovsky | Algorithms and implementation of a set partitioning approach for modular machining line design[END_REF], this problem is similar to the simple assembly line balancing problem, SALBP-1, see, for example, Boysen and Fliedner [START_REF] Boysen | A versatile algorithm for assembly line balancing[END_REF].

Constraint [START_REF] Andres | Balancing and scheduling tasks in assembly lines with sequence-dependent setup times[END_REF] incorporates the lower bound LB 1 . Constraints (3) ensure that the stations are opened consecutively in the increasing order of their indices and (4) that at least one operation is assigned to an open station. Note that constraints (2)-( 4) are not necessary for the model to be correct but they accelerate the solution process when using an ILP solver. Constraints [START_REF] Belmokhtar | Integer programming models for logical layout design of modular machining lines[END_REF] ensure that each operation is assigned to a station. Constraints [START_REF] Borgia | STEP-NC compliant approach for setup planning problem on multiple fixture pallets[END_REF] guarantee that the capacity of each station is not exceeded and that any operation can be assigned only to an open station.

Constraints [START_REF] Borisovsky | Algorithms and implementation of a set partitioning approach for modular machining line design[END_REF] prevent operations of the same exclusion set to be all assigned to the same station.

Constraints (8) together with (5) guarantee that the station index of a is strictly smaller than that of b for (a, b) ∈ A, which is the requirement of the precedence constraints.

The minimum number of stations can be calculated as k * = k 0 j=1 y * j where y * j are optimal y j values in the problem M in#.

2.5 Heuristics for the general case of P (prec, excl, s i | k)

As indicated, the problem considered in this paper is NP-hard. In addition, the ILP proposed requires an upper bound on the number of stations. Thus, in this section we suggest two randomized greedy heuristics to solve the problem P (prec, excl, s i | k) approximately. Let the stations be numbered t = 1, 2, . . . Heuristic GreedyN umber1 for the problem P (prec, excl, s i | k):

Step 1 (Elimination of exclusion relations) If there is no exclusion set (E = φ), then set t = 1 and go to Step 2. Otherwise, select a set E ∈ E. Randomly select a pair of operations in E and introduce an arbitrary precedence relation between them. Note that if this precedence relation is satisfied, then the exclusion relation given by E is satisfied as well. Update the precedence graph G = G(N, A). Update E by removing E and any other exclusion set which contains a pair of operations with a directed path between them in the new graph G. Repeat Step 1.

Step 2 (Operation assignment by solving Subset Sum problem) Calculate the set N + of vertices of the precedence graph G which have no predecessor. Assign to station t operations of a subset X ∈ N + , which is selected by solving the following Subset Sum problem: max

X∈N + i∈X s i , subject to i∈X s i ≤ r.
Update graph G by removing vertices of the set X and their outgoing arcs. If G is empty, then stop: a feasible solution is constructed. Otherwise, re-set t := t + 1 and repeat Step 2.

Let us illustrate heuristic GreedyN umber1 on the example with 5 operations, whose sizes are Iterative solving of the Subset Sum problem is proven to be an efficient solution approach for the bin packing problems; see, for example, Caprara and Pferschy [START_REF] Caprara | Modified subset sum heuristics for bin packing[END_REF] and Haouari and Serairi [START_REF] Haouari | Heuristics for the variable sized bin-packing problem[END_REF]. To solve the Subset Sum problem, we use the standard dynamic programming algorithm, which requires O(|N + |r) time.

s 1 = 1, s 2 = 2, s 3 = 3,
Heuristic GreedyN umber1 can be independently run several times to increase the probability to obtain a feasible solution with a small number of stations.

Heuristic GreedyN umber2 differs from GreedyN umber1 in that in Step 2, instead of solving the Subset Sum problem, the subset X ∈ N + such that i∈X s i ≤ r is generated by selecting its elements by random.

Minimizing the total activation cost

In this section, we establish computational complexity and present heuristics and an ILP formulation for the problem P (prec, excl, s i | lex(k, cost)).

Computational complexity

Note that any special case of the problem P (prec, excl, s i | lex(k, cost)) is not easier than the same special case of the problem P (prec, excl, s i | k). Therefore, NP-hardness results for the problem P (prec, excl, s i | k) apply for the problem P (prec, excl, s i | lex(k, cost)).

Kovalev et al. [START_REF] Kovalev | Line configuration to minimize setup costs[END_REF] presented a solution algorithm for the problem with unit size operations and no precedence or exclusion relation, which we denote as P (s i = 1 | lex(k, cost)). The run time of this algorithm depends solely on the number f of part types. Therefore, the problem

P (s i = 1 | lex(k, cost)
) is polynomially solvable for a constant number of part types. We now show that the problem P (s i = 1 | lex(k, cost)) is NP-hard in the strong sense if f is variable (a part of the problem instance) even if all the cost coefficients are unit.

Theorem 2 The problem P (s

i = 1 | lex(k, cost)) is NP-hard in the strong sense even if a v = 1, v = 1, . . . , f .
Proof: We use a reduction from the strongly NP-complete problem 3-Partition; see Garey and Johnson [START_REF] Garey | Computers and intractability. A guide to the theory of NP-completeness[END_REF]: Given 3q + 1 positive integer numbers b 1 , . . . , b 3q and B such that B/4 < b i < B/2, i = 1, . . . , 3q, and 3q i=1 b i = qB, is there a partition of the set {1, . . . , 3q} into q disjoint sets X 1 , . . . , X q such that i∈X l b i = B for l = 1, . . . , q?

Given an instance of 3-Partition, construct the following instance of the problem P (s i = 1 | lex(k, cost)). Define the number of operations n = qB and the station capacity r = B. There are f = 3q types of operations, each operation has a single type and a unit size, and there are b i operations of type i, i = 1, . . . , 3q. Thus, there are qB operations in total. We show that for the constructed instance, a solution to 3-Partition exists if and only if there exists a feasible solution to the constructed instance of the problem P (s i = 1 | lex(k, cost)) with q stations and total activation cost equal to 3q. It is obvious that the construction is pseudo-polynomial.

Part "if ". Assume that there exists a feasible solution to the constructed instance of the problem P (s i = 1 | lex(k, cost)) with at most q stations and at most 3q activation cost. It is easy to see that the minimum number of stations is equal to q. Since there are 3q part types, operations of the same type must be assigned to the same station. Let X h denote the set of part types assigned to station h. We must have i∈X h b i ≤ r = B, h = 1, . . . , q, which, together with the equality q h=1 i∈X h b i = 3q i=1 b i = qB implies i∈X h b i = B, h = 1, . . . , q, as required for the proof of the part "if".

Part "only if ". If there exists a solution X 1 , . . . , X q to the problem 3-Partition, then construct a solution to the problem P (s i = 1 | lex(k, cost)) in which all the operations of the types from the set X h are assigned to the station h, h = 1, . . . , q. This solution has q stations and the total cost of 3q, as required for the proof of the part "only if".

ILP formulation

We now present an ILP formulation for the general case of the problem P (prec, excl, s i | lex(k, cost)). Assume that the minimum number of stations k * has been found. Introduce 0-1 variables z vj such that z vj = 1 if and only if at least one operation of type v is assigned to station j. As before for the problem M in#, x ij are 0-1 variables such that x ij = 1 if and only if operation i is assigned to station j.

Calculate the following lower bound LB 2 on the minimum total cost:

LB 2 = f v=1 a v max l (v) max , i∈Nv s i r
, where l

(v)
max is the maximum number of type v op-erations in any chain of G. The problem P (prec, excl, s i | lex(k, cost)) can be re-formulated as the following ILP problem, which is denoted as M inAct.

Problem M inAct:

min k * j=1 f v=1 a v z vj , (10)
subject to

k * j=1 f v=1 a v z vj ≥ LB 2 , ( 11
) i∈Nv x ij ≥ z vj , v = 1, . . . , f, j = 1, . . . , k * , ( 12
) n i=1 s i x ij ≤ r, j = 1, . . . , k * , (13)
i∈Nv

x ij ≤ |N v |z vj , v = 1, . . . , f, j = 1, . . . , k * , (14) k
* j=1 x ij = 1, i = 1, . . . , n, (15)
i∈E x ij ≤ |E | -1, ∀E ∈ E, j = 1, . . . , k * , (16) h j=1
x bj +

k * j=h x aj ≤ 1, ∀(a, b) ∈ A, h = 1, . . . , k * , (17)
x ij , z vj ∈ {0, 1}, i = 1, . . . , n, j = 1, . . . , k * , v = 1, . . . , f.

In the problem M inAct, there are k * (f + n) variables and 1 + n + k * (1 + 2f + |E| + |A|) constraints without the 0-1 constraints.

Constraint [START_REF] Bukchin | A branch-and-bound based solution approach for the mixed-model assembly line-balancing problem for minimizing stations and task duplication costs[END_REF] incorporates the lower bound LB 2 . Constraints (12), though redundant because z vj will be set to zero in an optimal solution if x ij = 0 for all i ∈ N v , are included to be used by a solver in the solution process. Constraints (13) address the station capacities. Constraints [START_REF] Chica | A multiobjective model and evolutionary algorithms for robust time and space assembly line balancing under uncertain demand[END_REF] guarantee that any operation of a certain type can be assigned to a station only if this station is open for this type. Constraints ( 15) -( 17) have the same meaning as in the problem M in#.

Heuristics

For the problem P (prec, excl, s i | lex(k, cost)), we also suggest a randomized greedy heuristic, denoted as GreedyCost, to provide an upper bound. Let c i denote the total activation cost of the operation i ∈ N : c i = v∈T i a v . The only difference between the heuristics GreedyCost and GreedyN umber1 is that in GreedyCost, instead of solving the Subset Sum problem, the following Knapsack problem is solved:

max X∈N + i∈X c i , subject to i∈X s i ≤ r.
An optimal solution of this problem assigns operations of maximum total cost to the same station.

Again, we use the standard dynamic programming algorithm to solve the Knapsack problem.

It requires O(|N + |r) time. Heuristic GreedyCost can be independently run several times to increase the probability to obtain a feasible solution with a low cost.

Note that heuristic GreedyCost can also be used to solve the problem P (prec, excl,

s i | k)
approximately, and heuristics GreedyN umber1 and GreedyN umber2 can be used to solve the problem P (prec, excl, s i | lex(k, cost)) approximately.

Computer experiments

We used ILOG CPLEX version 12. Each operation i ∈ N is randomly assigned its size s i ∈ {1, 2, 3} and a set of types

T i ∈ {1, 2, 3}, {1, 2}, {1, 3}, {2, 3}, {1}, {2}, {3}
with uniform probability distribution. The exclusion sets are defined as follows. For each pair (a, b) of operations such that they have at least one common type, it is decided with probability 1/3 that a excludes b. Thus, all the exclusion sets are of cardinality two.

Each family is subdivided into four series n S , n M , n L and n X with small, medium, large and very large density of precedence constraints, respectively. For each series n J , J ∈ {S, M, L, X}, and each pair (a, b) of operations such that a < b and a and b have at least one common type, it is decided with probability p J that a precedes b. The corresponding probabilities are p S = 0.2, p M = 0.4, p L = 0.6 and p X = 0.8. Since each arc in the obtained precedence graph always goes from a vertex with smaller index to a vertex with larger index, this graph is acyclic.

Each series consists of 20 instances. Their average numerical characteristics are given in Table 2. There, %P rec is the density of the precedence constraints defined as The idea of these experiments is in line with similar studies on the transfer line balancing problems, see, e.g., Belmokhtar et al. [START_REF] Belmokhtar | Integer programming models for logical layout design of modular machining lines[END_REF], Borisovsky et al. [START_REF] Borisovsky | Algorithms and implementation of a set partitioning approach for modular machining line design[END_REF] and Dolgui et al. [START_REF] Dolgui | Minimizing setup costs in a transfer line design problem with sequential operation processing[END_REF]. There are no benchmark instances, however, because problems M in# and M inAct have never been studied before.

|A| max # arcs × 100%= 2|A| n 2 -n × 100%.
The following procedure is used to solve any generated instance.

Solution procedure for the problem P (prec, excl, s i | lex(k, cost)):

Step 1 Apply heuristics GreedyN umber1, GreedyN umber2 and GreedyCost 1000 number of times each. Determine the minimum number of stations, k 0 , delivered by these heuristics.

Step 2 Apply CPLEX to solve the problem M in# with upper bound k 0 on the number of stations. Determine the minimum number of stations k * .

Step 3 Apply CPLEX to solve the problem M inAct with the number of stations k * .

Step 4 Select the best solution for the problem P (prec, excl, s i | lex(k, cost)) among the solutions obtained in Steps 1, 2 and 3.

The possible outputs of each CPLEX run for an instance I are: 1) an optimal solution with value Opt(I) is found, 2) an approximate solution with value Appr(I) with a known relative error %Gap(I) and a lower bound LB(I) are obtained, 3) no feasible solution is found when the time limit is achieved. The relative error %Gap(I) is defined as Appr(I)-LB(I) LB(I) × 100%. In our experiments case 3) never happened.

Table 3 shows the number of instances of problems M in# and M inAct solved by CPLEX to optimality (row "Opt") and approximately (row "Appr"), for each series. Furthermore, the average (row "Avg.Gap") and the worst (row "W orst Gap") relative errors are listed for M in# problem. The average relative error Avg.Gap is defined as

1 |n J |
I∈n J %Gap(I) and the worst relative error W orst Gap as max I∈n J {%Gap(I)} for each series n J , J ∈ {S, M, L, X}. For the problem M in#, 237 instances out of 240 were solved optimally and the remaining 3 instances were solved approximately with an acceptable relative error. As for the problem M inAct, all instances were solved to optimality for all series. Table 3 also provides average (row "Avg.T ime") and maximal (row "M ax.T ime") computational times of CPLEX application to problems M in# and M inAct for all series. Performance of the heuristics GreedyNumber1, GreedyNumber2 and GreedyCost is presented in Table 4. More precisely, the row GreedyNumber AvgGap (GreedyNumber Worst.Gap) provides the average gap (worst gap) of the best results found on each instance with heuristics GreedyNumber1, GreedyNumber2. The section M in# of Table 4 contains the best results of the two heuristics, GreedyNumber1 and GreedyNumber2, with respect to the number of stations. The section M inAct of Table 4 shows the results of the heuristics GreedyCost with respect to the total activation cost.

Two conclusions can be drawn from these results. First, for almost all instances where optimal solutions were found by CPLEX, values of heuristic solutions range between 10% and 40% of the optimum. Second, there is a high correlation between density of precedence constraints and heuristic solution quality. The higher the density, the better is the solution. It is logical since the higher density means a fewer number of alternate assignments to evaluate. We deduce that CPLEX can be used to solve instances of practical dimensions in acceptable time. Heuristics GreedyNumber1 and GreedyNumber2 can be used to obtain an initial number of stations k 0 . All the heuristics can be used to obtain an approximate solution if the solution time is important. This is the case when there is a highly qualified and well paid expert whose job is to evaluate variants of the optimal line design for different types of equipment. Faster solutions save working time for the expert.

Note that earlier methods developed for similar problems observed in the line manufacturing industry were able to solve instances with up to 130 unit-sized operations, see, for example, Borisovsky et al. [START_REF] Borisovsky | Algorithms and implementation of a set partitioning approach for modular machining line design[END_REF]. Since the average size of operations in the experiments of the current paper is two, the largest instances solved in this paper include about 160 unit size operations.

Dichotomic procedure

In this section, we show how to take advantage of the difference in computational times observed between both problems. As shown in the previous section, when CPLEX could not find the optimal solution for problem M in#, it provided a lower bound on the minimum number of stations and an approximate solution, which delivers an upper bound on this value. Therefore, an optimal solution for the problem P (prec, excl, s i | lex(k, cost)) can be found by solving the As shown in these three examples, this dichotomic procedure can in some cases provide better solutions. However, the resolution of the problem M inAct appears much more difficult when CPLEX cannot find a feasible solution. Overall, the procedures described in this paper have permitted to solve optimally 238 instances out of 240, and the 2 remaining instances were solved approximately with a gap of 1 and 2 workstations respectively (and the optimal total cost for the corresponding number of stations).

Conclusions

In this paper, we studied the line balancing problem P (prec, excl, s i | lex(k, cost)) in which a line has to be configured to produce parts of different types. Each part of a specific type moves between stations in the same direction and a station is activated if at least one operation is assigned to it. Operations on the same part assigned to the same station are launched simultaneously and executed in parallel. The primary objective is to minimize the number of stations and the secondary objective is to minimize the total activation cost. Several results on computational complexity have been presented for this problem and are summarized in Table 8.

In order to solve the general case of the problem P (prec, excl, s i | lex(k, cost)), we suggested strongly NP-hard a four-stage solution approach, which combines heuristic algorithms with integer linear programming. Computer experiments show that this approach is able to solve instances of practical sizes.

A dichotomic procedure has also been presented and has permitted to obtain better solutions in the few cases where the four-stage solution approach was not able to solve the problem to optimality.

It is also interesting to note that the proposed dichotomic procedure can be used for solving the problem P (prec, excl, s i | lex(k, cost)) without solving the problem M in#. For example, LB 1 or the optimal value of the continuous relaxation of M in# can be used as a lower bound, and the best solution provided by the heuristics can be used as an upper bound. In this case, the dichotomic procedure can possibly require the resolution of a few more M inAct problems, but these additional resolutions can consume less time. Considering the results obtained for the examples presented in Section 5, it appears that the quality of the lower bound is more critical for the computational time of the dichotomic procedure than the quality of the upper bound. As a consequence, besides a search for more efficient heuristics, improvement of the lower bound for the problem M in# seems to be a more promising way to solve problems of larger dimensions.

Finally, the aforementioned dichotomic procedure could be easily extended to deal with a multiobjective version P (prec, excl, s i | k, cost) of this problem. In this case, the procedure would allow to generate a minimal complete set of efficient solutions in the sense of Pareto.

Figure 1 :

 1 Figure 1: Mixed-model mode. No operation for part type 3 is assigned to station 2 and this station is not activated when part of type 3 is processed

Figure 2 :

 2 Figure 2: Multi-model mode. Takt time is adjusted for each part type. When takt time decreases, the production of parts of the next type -type 2 in this example -starts before the change of the takt time

Figure 3 :

 3 Figure 3: Two assignments of operations

  max , is minimized. Denote the minimum C max value as C * max . For problem P (prec, s i = 1 | k), let us interpret operations as jobs, the station capacity bound r as the number of machines m, and stations 1, 2, . . . as unit-time intervals [0, 1], [1, 2], . . . in the problem P m | p i = 1, prec | C max . It can be seen that any feasible solution of problem P (prec, s i = 1 | k) with the number of stations k can be represented as a feasible solution of the corresponding scheduling problem with the objective value C max = k and vice versa. Therefore, the two problems are equivalent and k * = C * max . If precedence relations are arbitrary, then problem P m |

  It can be reduced to solving n parallel machine scheduling problems P m | d j = r | •, m = 1, . . . , n. In the problem P m | d j = r | •, there are n non-preemptive jobs with processing times s 1 , . . . , s n to be scheduled on m identical parallel machines, and the question is whether all the jobs can be completed by the common deadline r. We have k * = m * where m * is the minimum m for which the question in the problem P m | d j = r | • has the affirmative answer. The problem P m | d j = r | • can be solved in O(nmr m-1 ) time by a dynamic programming algorithm. Hence, the problem P (s i | k) can be solved in O(n 2 r n-1 log 2 n) time by a bisection search over the range 1, . . . , n of values m.

Theorem 1

 1 The problem P (excl, s i = 1 | k) is NP-hard in the strong sense if |E | = 2 for each E ∈ E. Proof: We use a reduction from the problem Exact Cover by 3-Sets (X3C); see Garey and Johnson [27]: Given a family A = {A 1 , . . . , A a } of 3-element subsets of the set B = {1, . . . , 3b}, does A contain an exact cover of B, i.e., a subfamily S ⊆ A such that each i ∈ B belongs to exactly one 3-element set in S? Assume a > b because otherwise there exists a trivial solution. Given an instance of X3C, construct the following instance of the decision version of the problem P (excl, s i = 1 | k). Define the number of operations n = (b + 1)(a -b + 1) and the station capacity r = b + 1. The set of operations includes a single X-operation, b number of the same Y i -operations, i = 1, . . . , a-b, and A i -operations associated with the subsets A i , i = 1, . . . , a.All operations are of the same type. The exclusion relations are given as follows. The Xoperation excludes any Y i -operation, i.e., {X, j} ∈ E for every Y i -operation j, i = 1, . . . , a -b.

s 4 = 2 and s 5 = 1 .

 51 Further, station capacity is r = 4, exclusion relations are given by E = {E 1 , E 2 }, where E 1 = {1, 3} and E 2 = {1, 4}, and precedence relations are given by arcs[START_REF] Akpinar | A hybrid genetic algorithm for mixed model assembly line balancing problem with parallel workstations and zoning constraints[END_REF][START_REF] Andres | Balancing and scheduling tasks in assembly lines with sequence-dependent setup times[END_REF],[START_REF] Akpinar | A hybrid genetic algorithm for mixed model assembly line balancing problem with parallel workstations and zoning constraints[END_REF][START_REF] Battaia | Reduction approaches for a generalized line balancing problem[END_REF],[START_REF] Battaia | Reduction approaches for a generalized line balancing problem[END_REF][START_REF] Battaia | A taxonomy of line balancing problems and their solution approaches[END_REF],[START_REF] Battaia | Reduction approaches for a generalized line balancing problem[END_REF][START_REF] Belmokhtar | Integer programming models for logical layout design of modular machining lines[END_REF] and[START_REF] Battaia | A taxonomy of line balancing problems and their solution approaches[END_REF][START_REF] Belmokhtar | Integer programming models for logical layout design of modular machining lines[END_REF] in graph G with 5 vertices. According to the heuristic, exclusion relation {1, 3} can be eliminated by introducing precedence relation[START_REF] Akpinar | A hybrid genetic algorithm for mixed model assembly line balancing problem with parallel workstations and zoning constraints[END_REF][START_REF] Battaia | Reduction approaches for a generalized line balancing problem[END_REF].Since directed path (1, 3, 4) appeared, exclusion relation {1, 4} can be eliminated as well. Then, t = 1, and operations are assigned to station 1. Since set N + = {1} and s 1 ≤ r, operation 1 is assigned to station 1. It is removed from G, and t := 2. Operations are assigned to station 2. Now,N + = {2, 3} and s 2 + s 3 > r. Solution of the Subset Sum problem is X = {3}.Operation 3 is assigned to station 2. It is removed from G, and t := 3. Operations are assigned to station 3. Now, N + = {2, 4} and s 2 + s 4 ≤ r. Solution of the Subset Sum problem is X = {2, 4}. Operations 2 and 4 is assigned to station 3. Only operation 5 with s 5 ≤ r remains unassigned. It is assigned to station 5.

  4 to solve the problems M in# and M inAct. The experiments were run on PC with Intel(R) Core(TM) Quad 2.83 Mhz processor and 8 Gb of RAM under MS Windows 7 Professional (64 bit). The data generation algorithm for CPLEX and the heuristics were coded in C++. In the experiments, the time limit for solving each instance of M in# or M inAct was set to one hour. Three families of instances were randomly generated. A family is associated with the number of operations n ∈ {40, 60, 80}. For all instances, there are three part types, the activation costs are a 1 = 3, a 2 = 2 and a 3 = 1 and the station capacity is equal to r = 10. Three part types are chosen because two or three modifications of the same product are typical in the application we are concerned with. Values of activation costs are determined based on a possible proportion 3 : 2 : 1 between the quantities of parts of three types in their mass production. Station capacity 10 is one of possible station capacities in the application.

  s i = 1 | k), arbitrary r strongly NP-hard P (prec, s i = 1 | k), constant r open P (in -tree, s i = 1 | k) O(n) P (out -tree, s i = 1 | k) O(n) P (prec, s i = 1 | k), r = 2 O(n 2 ) P (s i | k), arbitrary r strongly NP-hard P (s i | k), constant r O(n const ) P (excl, s i = 1 | k)strongly NP-hard P (s i = 1 | lex(k, cost))

  Kovalev et al.[START_REF] Kovalev | Line configuration to minimize setup costs[END_REF] studied a special case of the problem P (prec, excl, s i | lex(k, cost)) with unit size operations and no precedence or exclusion relation, which we denote as P (s i = 1 |

	lex(k, cost)). They suggested a solution procedure with running time depending solely on the
	number of part types f . This procedure combines the enumeration of different assignments of part
	types with the solution of a system of linear inequalities. Dolgui et al. [17] investigated a problem
	similar to P (prec, excl, s i | lex(k, cost)), in which operations of the same type assigned to the same
	station are performed sequentially and there are no exclusion relations. The takt time constraints
	were explicitly present in the formulation of this problem because of the sequential execution of
	operations at the same station. NP-hardness proofs, polynomial time algorithms for special cases,

heuristics and integer linear programming formulations for the general problem were suggested in

[START_REF] Dolgui | Minimizing setup costs in a transfer line design problem with sequential operation processing[END_REF]

. These results cannot be applied to the problem P (prec, excl, s i | lex(k, cost)) because the difference between sequential and simultaneous processing makes the solution structures and feasible domains of these problems substantially different. The fact that the type of processing drastically changes the mathematical model and affects solution methods is widely observed in the line balancing problems.

Problem P (prec, excl, s i | k) is related to multi-product line balancing problems, which are studied by Sawik

[START_REF] Sawik | Integer programming models for the design and balancing of flexible assembly systems[END_REF]

, Kabir and Tabucanon

[START_REF] Kabir | Batch-model assembly line balancing: A multiattribute decision making approach[END_REF]

, Bukchin and Rabinowitch

[START_REF] Bukchin | A branch-and-bound based solution approach for the mixed-model assembly line-balancing problem for minimizing stations and task duplication costs[END_REF]

, Vilarinho and Simaria

[START_REF] Vilarinho | ANTBAL: An ant colony optimization algorithm for balancing mixed-model assembly lines with parallel workstations[END_REF]

, Andres et al.

[START_REF] Andres | Balancing and scheduling tasks in assembly lines with sequence-dependent setup times[END_REF]

, Yuan et al.

[START_REF] Yuan | An effective hybrid honey bee mating optimization algorithm for balancing mixed-model two-sided assembly lines[END_REF]

, among others. Surveys of line balancing research are given by Erel and Sarin

[START_REF] Erel | A survey of the assembly line balancing procedures[END_REF]

, Rekiek et al.

[START_REF] Rekiek | State of art of assembly lines design optimization[END_REF]

, Boysen et al.

[START_REF] Boysen | Sequencing mixed-model assembly lines: Survey, classification and model critique[END_REF]

, Dolgui and Proth

[START_REF] Dolgui | Supply chain engineering: useful methods and techniques[END_REF]

,

Table 1 :

 1 Structure of a feasible solution with a -b + 1 stations.

Table 2 :

 2 Average numerical characteristics

	Series n J	40 S	40 M	40 L	40 X	60 S	60 M	60 L	60 X	80 S	80 M	80 L	80 X
	%P rec	9.35% 19.25% 29.74% 40.25% 9.88% 19.74% 29.74% 41.24% 10.09% 19.93% 30.09% 39.43%
	n i=1 s i	78	79	79	80	117	120	120	120	160	160	158	159
	|E|	189	186	202	199	444	451	443	458	797	790	791	786

Table 3 :

 3 CPEX: solution quality and run time

	Series n J	40 S	40 M	40 L 40 X	60 S	60 M	60 L	60 X	80 S	80 M	80 L	80 X
								M in#
	Opt	20	20	20	20	20	20	20	20	18	19	20	20
	Appr	0	0	0	0	0	0	0	0	2	1	0	0
	Avg.Gap	-	-	-	-	-	-	-	-	11.03	16.67	-	-
	W orst Gap	-	-	-	-	-	-	-	-	13.1	16.67	-	-
	Avg.T ime (sec)	5.99 4.91 3.06 2.08 151.80 94.58 27.20 26.61 1416.60 1144.04 718.22 249.84
	M ax.T ime (sec) 18.11 13.00 12.18 7.02 312.19 280.97 68.95 173.29	3600	3600	3137.18 1312.09
								M inAct
	Opt	20	20	20	20	20	20	20	20	20	20	20	20
	Appr	0	0	0	0	0	0	0	0	0	0	0	0
	Avg.T ime (sec)	0.84 0.30 0.37 0.82	5.08	1.41	1.96	5.83	145.56	65.45	9.43	38.86
	M ax.T ime (sec) 2.31 0.86 0.64 2.11 17.96	3.26	3.9	28.77 1037.45 1182.77 18.63	309.05

Table 4 :

 4 Heuristics GreedyCost and GreedyNumber : solution quality and run time GreedyN umber Avg.Gap (%) 42.63 21.42 20.98 11.09 32.99 23.26 19.33 11.32 31.79 22.47 16.29 09.74 GreedyN umber W orst Gap (%) 66.67 38.53 37.57 25.52 46.15 39.83 32.78 20.90 50.12 39.35 31.67 24.16General observations derived from the computer experiments are the following. First of all, CPLEX solves large scale instances of the problem P (prec, excl, s i | lex(k, cost)) in a short time.On average, problems are solved within 25 minutes. Secondly, problem M in# needs more computation effort than MinAct. It may be related to the high values of k 0 provided by the heuristics GreedyNumber. The third observation concerns the density of precedence constraints. More computation time is usually needed to solve instances with smaller density of precedence constraints.

	Series n J	40 S	40 M	40 L	40 X	60 S	60 M	60 L	60 X	80 S	80 M	80 L	80 X
							M in#
	Avg.T ime (sec)	95.8 95.1 94.9 94.5 142.7 142.2 141.7 141.1 189.9 189.5 188.2 187.6
							M
	GreedyCost Avg. Gap (%)	32.73 26.13 17.86 9.64 26.49 21.78 16.78 11.32 29.74 24.48 22.28 9.31
	GreedyCost W orst Gap (%)	54.60 41.18 34.05 17.78 42.96 35.86 29.69 23.14 53.87 41.33 32.35 23.60
	Avg.T ime (sec)	125.6 125.2 124.7 124.3 202.8 202.3 201.9 201.4 342.1 341.6 341.1 340.8

Table 5 :

 5 Illustration of dichotomic procedure on the first instance unsolved for the series 80 S Employ dichotomic procedure with LB = 18 and U B = 19:Apply CPLEX to solve M inAct with 19 stations Output: optimal solution for 19 stations in 318 sec Employ dichotomic procedure with LB = 18 and U B = 18:Apply CPLEX to solve the problem M inAct with 18 stations Output: optimal solution for 18 stations in 389 sec Output: optimal solution in 707 sec

Table 6 :

 6 Illustration of dichotomic procedure on the second instance unsolved for the series 80 S Employ dichotomic procedure with LB = 18 and U B = 20: Apply CPLEX to solve M inAct 19 stations Output: optimal solution is for 19 stations found in 679 sec Employ dichotomic procedure with LB = 18 and U B = 18: Apply CPLEX to solve M inAct with 18 stations Output: no feasible solution for 18 stations within the time limit Output: Approximate solution with 19 stations in 4,279 sec A gap of 1 station (5.56%) For 19 stations, solution with minimum total cost

Table 7 :

 7 Illustration of dichotomic procedure on the instance unsolved for the series 80 M Employ dichotomic procedure with LB = 26 and U B = 30: Apply CPLEX to solve M inAct with 28 stations Output: optimal solution for 28 stations is found in 1,233 sec Employ dichotomic procedure with LB = 26 and U B = 27: Apply CPLEX to solve M inAct with 27 stations Output: no feasible solution for 27 stations within the time limit Employ dichotomic procedure with LB = 26 and U B = 26: Apply CPLEX to solve M inAct with 26 stations Output: no feasible solution for 26 stations within the time limit Output: Approximate solution with 28 stations in 8,433 sec A gap of 2 stations (7.69%) For 28 stations, solution with minimum total cost

Table 8 :

 8 Computational complexity ProblemComplexity P (prec,

problem M inAct for each number of stations between these two bounds. In order to solve problem P (prec, excl, s i | lex(k, cost)) more efficiently in this case, we suggest the following dichotomic procedure.

Dichotomic solution procedure for the problem P (prec, excl, s i | lex(k, cost)) with a lower bound LB and an upper bound U B on the minimum number of stations:

Step 1 Apply CPLEX to solve the problem M inAct with the number of stations LB+U B 2 .

Step 2 Consider the following cases. Step 3 Among the solutions obtained in Step 1 and 2, select the best solution for the problem P (prec, excl, s i | lex(k, cost)), which has the minimum cost among the solutions with the minimum number of stations.

The possible outputs of this procedure for an instance I of the problem P (prec, excl, s i | lex(k, cost)) are: 1) an optimal solution is found, 2) an approximate solution with a known relative error on the number of stations, or on the total cost, or on both criteria is obtained, 3)

no new feasible solution is found when the time limit is achieved. In the latter case, however, the relative error can possibly be reduced if CPLEX proves that no feasible solution exists for some of the numbers of stations considered.

Let us now illustrate this procedure on the three instances that remained unsolved in the previous section (see Tables 5,[START_REF] Borgia | STEP-NC compliant approach for setup planning problem on multiple fixture pallets[END_REF][START_REF] Borisovsky | Algorithms and implementation of a set partitioning approach for modular machining line design[END_REF].