Olivier Guyon
email: guyon@emse.fr

Pierre Lemaire
email: pierre.lemaire@grenoble-inp.fr

Éric Pinson
email: pinson@uco.fr

David Rivreau
email: rivreau@uco.fr

Solving an integrated job-shop problem with human resource constraints

Keywords: Employee Timetabling Problem, Job-shop, Job-shop with resource availability constraints, Probing, Cut generation, Branch and Bound

We propose two exact methods to solve an integrated employee-timetable and job-shop-scheduling problem. The problem is to find a minimum cost employee-timetable, where employees have different competences and work during shifts, so that the production, that corresponds to a job-shop with resource availability constraints, can be achieved. We introduce two new exact procedures: 1) a decomposition and cut generation approach and 2) a hybridization of a cut generation process with a branch and bound strategy. We also propose initial cuts that strongly improve these methods as well as a standard MIP approach. The computational performances of those methods on benchmark instances are compared to that of other methods from the literature.

Introduction

The purpose of a manufacturing factory is merely to produce some goods to meet some demand. Due to limited resources, an optimal production plan is usually hard to compute. The complexity lies mainly in two intertwined aspects: 1) a production schedule, that is an allocation of human and material resources to the different tasks (or jobs) that have to be processed, and 2) an employee timetable, that ensures that the human resources required by the production schedule are actually met.

Except for special cases, each aspect is per se a difficult problem, as shows the huge literature on both scheduling problems (e.g., see Pinedo [Pin04] and Leung [Leu04]) and timetabling problems (e.g., see Ernst et al. [EJKS04] and Soumis et al. [SPR05] for states of the art). As a consequence, the resulting integrated problem has long been considered as too complex to be solved in practice and it is usually decomposed into an assignment part and a scheduling part, resulting in sub-optimal solutions.

Because of economical and financial pressure, higher performances must be sought for, through elaborate and adequate optimization procedures. That is why some recent efforts have been made to actually tackle the integrated problem, and this paper is one of them (indeed, this is the continuation of a previous work on a simpler problem [GLPR10] to try out the performances of cut generation processes for integrated problems).

In this paper, we propose two exact methods to solve an integrated employee-timetable and job-shob-scheduling problem. First, we provide a precise description of the problem we intend to solve, together with an integer linear formulation (Section 1), and we stress its links with the existing works (Section 2). Then, we introduce our two new procedures: a decomposition and cut generation approach (Section 3), and a hybridization of a cut generation process with a branch and bound strategy (Section 4). The computational performances of those methods on benchmark instances are compared to that of other methods from the literature (Section 5). Some conclusions are finally drawn in Section 6.

A model integrating an employee timetabling and a job-shop scheduling problems

Problem description

Our purpose is to solve two decision levels which interact in a global optimization process: a timetabling and a scheduling problem.

The job-shop problem consists in processing a set J of n jobs on a set K of m different machines. Each job i is made of a sequence of operations O i1 , O i2 , . . . , O im which have to be scheduled according to a given order. Each operation O i j has a processing time p i j ∈ N. It can be processed exactly by one of the m available machines. This performing machine is denoted by m i j . For the sake of clarity, we denote by ρ ik the processing time of job i on machine k. It is not allowed to preempt operations.

An operation can only be processed on machine k if an operator, able to use this resource, is available. We thus introduce a set E of µ operators where each operator e masters a subset K e of machines. We assume that the operators work on a three-shift system. Consequently, the time horizon H we use is divided into a set S of σ consecutive and identical shifts s 0 , s 1 , . . . , s σ -1 . Each shift corresponds to a time period and has a fixed duration π (thus: H = σ • π). Each employee e is assumed to be available on a subset of shifts S e . The cost of its assignment to machine k during shift s is denoted by c eks . Furthermore, regulation constraints impose that an operator must not work more than one shift among each triplet of consecutive shifts.

The objective of the problem is to find a minimum cost assignment of operators to both machines and shifts in such a way that a feasible operating sequence exists for each machine for a target makespan C max ≤ H (maximum of the completion times of operations).

A small illustrative instance is provided by Example 1 (shown in Figure 1). It consists of 2 jobs, 2 machines, 3 employees and 3 shifts. The shift duration is π = 8 hours (H = 24). We want to schedule each operation before the end of the day, i.e., C max = 24.

In Figure 1.a), the job data (machines and durations) are provided. Employee data are displayed in Figure 1.b). The costs to assign an employee to this machine are displayed for all shifts for which the employee is assumed to be available.

Figure 1.c) illustrates an optimal solution (with a total cost of 9) for Example 1. In this solution, a production plan with a makespan of 20, lower than C max , is obtained. Employee e 1 (respectively e 2 ) is assigned to machine m 2 (resp. m 1 ) during s 0 in order to perform operations O 11 and O 22 (resp. O 21 ). During the last shift s 2 , operation O 12 is processed on m 1 by e 3 .

i mij pij 1 2 1 3 4 2 1 2 2 5 O11 16 0 8 24 = Cmax m1 O21 O12 O22 e2 (4) e1 (2) e3 (3) s0 s1 s2 m2 c) An optimal solution a) Job data b) Employee data e Se Ke c eks 1 0 1 2 2 2 4 7 2 0 2 1 4 -3 3 1 2 1 -5 3 2 -6 8
Fig. 1 Example 1

An integer linear programming model

In this section, we propose an integer linear programming model [P] associated with the problem at hand. For the sake of clarity, this mathematical formulation is split into four parts: the objective-function, the employee timetabling sub-problem, the job-shop sub-problem, and the coupling constraints. Beforehand, we define the release date, or head (respectively the latency duration, or tail), denoted r ik (respectively d ik ), associated with the processing of job i on machine k. These dates are simply induced by the conjunctive constraints associated with the job sequences. They can be computed in O(n • m) time using the following recurrence relations:

r ik = 0 i = 1, . . . , n k = m i1 r il = r ik + ρ ik i = 1, . . . , n j = 1, . . . , m -1 k = m i j l = m i( j+1)
(1)

d ik = C max i = 1, . . . , n k = m im d ik = d il -ρ il i = 1, . . . , n j = 1, . . . , m -1 k = m i j l = m i( j+1)
(2)

Let us now define the two groups of binary decision variables involved in the model:

-x eks = 1 if and only if operator e is assigned to machine k during shift s, 0 otherwise -y ikt = 1 if and only if job i is processed on machine k at time t, 0 otherwise

Objective-function

[P] minΘ = ∑ e∈E ∑ k∈K e ∑ s∈S e c eks • x eks (3) Employee timetabling specific constraints ∑ k / ∈K e σ -1 ∑ s=0 x eks = 0 e = 1, . . . , µ (4)
∑

k∈K e ∑ s / ∈S e x eks = 0 e = 1, . . . , µ (5)
∑ k∈K e (x eks + x ek(s+1) + x ek(s+2) ) ≤1 e = 1, . . . , µ s = 0, . . . , σ -3 (6)

x eks ∈ {0, 1} e = 1, . . . , µ k = 1, . . . , m s = 0, . . . , σ -1 (7)
Job-shop specific constraints

d ik -ρ ik ∑ t=0 t • y ikt + ρ ik ≤ C max i = 1, . . . , n k = m im (8) d ik -ρ ik ∑ t=r ik y ikt = 1 i = 1, . . . , n k = 1, . . . , m (9) r ik ∑ t=0 y ikt + C max ∑ t=d ik -ρ ik +1 y ikt = 0 i = 1, . . . , n k = 1, . . . , m (10) t ∑ u=r ik +ρ ik y ilu - t-ρ ik ∑ u=r ik y iku ≤ 0 i = 1, . . . , n j = 1, . . . , m -1 k = m i j l = m i( j+1) t = r ik + ρ ik , . . . , d il -ρ il (11) n ∑ i=1 min(d ik -ρ ik ,t) ∑ u=max(r ik ,t-ρ ik +1) y iku ≤ 1 k = 1, . . . , m t = 0, . . . ,C max (12) y ikt ∈ {0, 1} i = 1, . . . , n k = 1, . . . , m t = 0, . . . ,C max (13) Coupling constraints ∑ e∈E x eks - n ∑ i=1 min(d ik -ρ ik ,t) ∑ u=max(r ik ,t-ρ ik +1) y iku ≥ 0 k = 1, . . . , m t = 0, . . . ,C max s = ⌊t/π⌋ (14)
In this formulation, assignment variables x eks are fixed to 0 if employee e does not master resource k (4) or is not available during shift s (5) 1 . Constraints (6) formalize the regulation rules stating that an employee must not work more than one shift among each triplet of consecutive shifts. Constraints (8) ensure that all jobs are completed prior to the target makespan C max . Each operation has to be processed within its time window ( 9)-( 10) and cannot start prior to the completion of its job predecessor (11). Furthermore, at most one operation can be processed on a given machine at each instant (12). Lastly, coupling constraints ( 14) ensure that an operation on a machine can be processed only if a qualified employee is assigned to the machine at hand. This model is at the core of this work: both procedures of Sections 3 and 4 aim at solving it, using different decomposition strategies.

For the scheduling stage, one may remark that operators may be unavailable and/or not assigned to some machines for some shifts over the time horizon. Therefore, we have to take into account for each machine a collection ϒ of time periods on which it cannot process any operation. Obviously, these periods depend on an instantiation x of decision variables x eks , since x corresponds to a particular assignment of operators to both machines and shifts. This leads to the so called job-shop scheduling problem with availability constraints.

However the problem can be reduced to a classical job-shop scheduling problem by creating additional fictional jobs, as stated in the following proposition:

Proposition 1 Let us consider a time slot υ = α υ , β υ corresponding to an unavailability period on machine k. To take into account this fixed inactivity period is equivalent to consider the additional job i F composed of m operations with execution domain 0,C max and a null processing time except one of them (without loss of generality the first one) with execution domain α υ , β υ and a duration equals to

ρ i F k = (β υ -α υ ).
Proof Operations of i F with an execution domain 0,C max and a null processing time can be processed at any moment. Scheduling i F thus only depends on the scheduling of the only operation with a non-null processing duration, i.e., the one which has to be processed on k during υ. Let us denote by O k i F this specific operation. Because all other operations of i F have a null duration, O k i F can be placed at any rank of the sequence of operations of i F . O k i F has a duration (β υ -α υ ) and a processing domain α υ , β υ . Scheduling O k i F therefore implies that machine k is fully dedicated to O k i F during υ. No real job can indeed be processed on k during υ. Machine k is in fact not used during υ. This complete the proof.

In the following, we will denote by J F the set of n F fictional jobs associated with the fixed inactivity periods of ϒ , with n F = card {ϒ }. Given an assignment vector x, the schedul- ing component associated with our problem consists in checking whether a feasible operating sequence associated with both real and fictional jobs (J ∪ J F ) exists for each machine and for the target makespan C max .

Existing works

As it has already been mentioned, the integrated approach is usually considered as too complex to be solved in practice and, as a consequence, much of the literature deals either with production scheduling problems (e.g., see [Pin04,Leu04] for states of the art) or with timetabling problems (e.g., see [EJKS04,SPR05]). In what follows, the main works dealing with integrated problems and/or decomposition procedures are presented in Subsection 2.1. Subsection 2.2 is dedicated to solution procedures for job-shop problems.

Integrated problems and decomposition procedures

Only few attempts exist that cope with an integrated solution method. To the best of our knowledge, the existing works either deal with different production systems and/or solve it using fundamentally different approaches.

In Daniels and Mazzola [DM94], flow-shop scheduling contexts are considered, and the authors assume that each operation must be processed by its machine parallel to a set of operators devoted to this operation during its entire processing time. Moreover, each employee can only be assigned to a subset of operations (because of skills considerations) and no shift partitioning of the time horizon is considered. Consequently, operators can be simply considered as additional parallel machines. They propose an enumerative approach of brand and bound type for solving this problem.

In the context of project scheduling, Alfares and Bailey [AB97] propose an integer linear programming model and dynamic programming based heuristics. Their problem is rather simple since the only constraints are precedences among operations, and that the volume of human resources on any given day is big enough for the operations of that day.

More recently, Chen [Che04] considers a parallel machine scheduling problem that involves job scheduling coupled to a resource allocation component; the processing times are inversely related to the amount of resources allocated and the objective is to minimize the total cost associated with both scheduling and resource allocation. A column generation based branch and bound is proposed for solving this problem.

Several other authors use decomposition procedures. Such procedures originate from Benders decomposition [Ben62]. In its basic version, it is a generic method for solving problems that contain groups of variables of different natures (e.g., MIPs with integer and continuous variables). The key idea is to assign some trial values to a group of variables and to find the best solution consistent with this particular instanciation. The major underlying idea relies on dual inference with respect to Lagrangian relaxation and Dantzig-Wolfe decomposition. Thus, the key element of this approach is the derivation of Benders cuts that exclude superfluous solutions. Classical Benders cuts are formulated by solving the dual of the subproblem obtained when the trial values are fixed, by exploiting a "nice" substructure in the global problem. Such an approach is commonly used for problems involving strategic as well as operational decisions.

During the last two decades, extensions of this powerful approach have been proposed: in particular, hybridization with constraint programming techniques or logic based Benders decompositions. One of the main pioneering works exploiting these ideas is due to Hooker [HO03,Hoo05,Hoo07] jointly with some other contributors. For instance, in [Hoo05,Hoo07], an hybrid method combining constraint programming techniques, mixed integer linear programming and logic based Benders decompositions is designed for solving a "planning and scheduling problem" in which tasks are to be assigned to facilities with some consumption considerations.

Similar solution procedures have been used by Artigues et al. [AGRV09] to solve a model linking a job-shop scheduling stage to an employee timetabling component. To the best of our knowledge, their model is the previously studied model that is the closest to the one developed in this paper. More precisely, a set A of additional activities to be completed by a set E of operators is introduced. The subset of activities each operator is able to perform is known. Each operation associated with the job-shop problem is then assumed to require, during its processing, a predefined number of these employees for each activity a ∈ A. As in our case, the time horizon is assumed to be partitioned in identical shifts. Clearly, the processing of a given operation O i j during shift s induces a fixed requirement δ as of operators devoted to the execution of each activity a on this particular time slot. Notice that this model allows an operator to work simultaneously on several machines. Moreover, an employee shift can cover more than one scheduling period. As a consequence, the authors allow the aggregation of activity demands for the operations processed during a given shift.

A lexicographic optimization problem where the makespan minimization is the primary objective whereas the employee cost minimization is the secondary objective is considered. The problem tackled in our study corresponds to the second stage of the problem described in [AGRV09].

As for the solution procedure, Artigues et al., unlike we, heavily rely on constraint programming. Indeed, two constraint programming based formulations, in association with linear programming relaxations, are proposed. The first (direct) constraint programming formulation models both job-shop scheduling and employee timetabling components by means of global constraints. A linear programming relaxation based on the employee timetabling part is used for search tree pruning and reduced cost reduction is proposed. The second constraint programming formulation relies on a decomposition involving the domains of the starting time variables, and leading to the definition of demand intervals associated with each activity a ∈ A. A new linear programming relaxation, including the scheduling component, is introduced by means of additional binary variables. Energetic reasoning feasibility checking techniques can thus be exploited. Both constraint programming formulations are solved using a backtrack search strategy.

In [GLPR10], Guyon et al. investigate the integration of an employee timetabling and a production scheduling problems. At the first level, they manage a classical employee timetabling problem whereas at the second level, they aim at supplying a feasible production schedule for a set of interruptible tasks with qualification requirements and time-windows. They propose two exact methods to solve the resulting problem. The former is based on a Benders decomposition while the latter relies on a specific decomposition and a cut generation process. Although this second approach has some basic similarities with the first one presented in this paper, most of the key components (cut structure, master and slave problems, ...) differ. In particular, the scheduling stage in [GLPR10] (scheduling on parallel machines with preemption) can be reduced to a polynomially solvable flow problem, while in the present work the scheduling stage (a job-shop) is NP-complete and notoriously hard to solve in practice.

Solution procedures for job-shop problems

As it is mentioned above, the scheduling stage is a job-shop problem, which is an NPcomplete problem. Hence an important issue is how to solve it efficiently. For this particular sub-problem, we rely on the existing litterature.

More precisely, we want to solve a job-shop with availability constraints. This is NPhard since job-shop without unavailability periods is already strongly NP-hard [RK76,GJ79]. Due to their evident utility for applications, scheduling problems with availability constraints are specifically adressed in the scheduling literature [BFHS88, DM94, Lee96, SS98, Sch00, Agg04a, GH05, WSC05, LS08, MCZ10].

In particular, Mauguière et al. [MBB05] show that 5 categories of the job-shop scheduling problem with availability constraints have to be distinguished. The type of problems we consider in this paper corresponds to the class denoted JP2 which refers to the jobshop problem with non-crossable unavailability periods and non-resumable operations. In [Agg02], Aggoune defines a method (based on the search of a shortest path problem) to solve up to optimality problems of JP2 with only two jobs. He also proposes an exact branch and bound method (that exploits the classical disjunctive graph representation of the jobshop problem [RS64]) for any problem of JP2 [Agg04b]. Other effective attempts [MBB05,Zri05] to solve the job-shop scheduling problem with availability constraints clearly show the growing interest of researchers for it.

As said before, the scheduling part of our problem can be seen as a classical job-shop problem with fictional jobs. Since the concern of our paper was not to investigate that specific sub-problem but mainly to deal with the integrated problem, we chose to use the dedicated job-shop solver of [Riv99] as a black box. In brief, this solver is a classical branch and bound process where the key point is the use of the immediate selections on disjunctions [Car75], the immediate selections on ascendant/descendant sets [CP89, CP90] and the shaving procedure [CP94, MS96] as elimination rules. A survey on these elimination rules can be found in [CPPR04]. Since the aim of the solver is either to find a solution within a given C max or to prove that there is no such solution in a given amount of CPU time, in each node we try to build a feasible solution using serial scheduling based on the current deadlines of operations. If this procedure fails, we relax the precedence constraints on jobs and we build a feasible solution on each machine (using [Car82]). These solutions provide us with earliest starting times and latest starting times for each operation. Clearly, most of the time, these dates do not fulfill the precedence relations constraints of jobs (otherwise we have a feasible solution and we are done). So in order to identify the greatest violation, the earliest starting times and latest starting times are amended through the precedence relations of each job by a forward and backward process. This gives us an evaluation of violation for each job: we select the job with the lowest compliance level and then we develop two nodes by halving the time window of one of its operations.

From a functional point of view, the solver is given a CPU time limit, so it can fail to prove that there is a solution or not (job-shop is not a particularly easy scheduling problem). We can also select to use the shaving procedure or not: indeed, if for hard problems this elimination procedure is essential, for easy problem it induces a CPU over-consumption that may be sometimes avoided.

A decomposition and cut generation approach

Problem decomposition

Due to the intrinsic two-level decision structure of [P], it seems quite natural to investigate decomposition methods based on the splitting of [P] into two interacting sub-problems:

-A master problem corresponding to the employee timetabling component of [P]: [ET P] -A slave or satellite problem corresponding to the scheduling component of [P]: [JobShop] This approach clearly relies on the relaxation of the coupling constraints (14).

Master problem

The master problem [ET P] is an employee timetabling problem that attempts to find a minimal cost assignment of operators to both machines and shifts, where coupling constraints binding operators and machines are relaxed. This problem, relying only on the decision variables x eks , is stated as follows: 5), ( 6), (7) where Cut denotes a set of valid inequalities that are iteratively added to the model. They invalidate solutions of [ET P] which do not lead to a feasible solution for the scheduling component of [P].

[ET P] : minΘ = ∑ e∈E ∑ k∈K e ∑ s∈S e c eks • x eks Cut (4), (

Slave problem

Let us denote by x an optimal solution to the current master problem [ET P]. Clearly, x is a particular assignment of operators to both machines and shifts. Moreover, x is a feasible solution for the global problem [P] if it satisfies the scheduling component [JobShop]. To perform this test, we first define aggregate information based on x in the following way :

zks = min( ∑ e∈E xeks , 1) k = 1, . . . , m s = 0, . . . , σ -1 (15)
Clearly, zks = 1 if at least one operator is assigned to machine k during shift s, 0 otherwise (the machine k is not used). Each machine requires an employee for its use. Therefore, machine k can be used during shift s if and only if zks = 1. Considering those entries, checking whether or not [JobShop] is feasible for a given assignment x can be performed by solving the following decision problem:

[JobShop(z)] : Is it possible to find an assignment of decision variables satisfying:

n ∑ i=1 min(C max ,( ῡ+1)•π) ∑ t=max(r ik , ῡ•π-ρ ik ) y ikt = 0 k = 1, . . . , m ῡ ∈ Sk ( 16
)
(8), ( 9), ( 10), ( 11), ( 12), ( 13)

where Sk = {s ∈ S|z ks = 0} denotes the set of time slots during which machine k cannot be used with respect to employee resources defined by z.

[JobShop(z)] is a job-shop scheduling problem with fixed inactivity periods. It can be reduced to a classical job-shop scheduling problem (see Proposition 1, Section 1.2). The only difference with the job-shop specific constraints presented in Section 1.2 relies on constraints (16). These constraints prevent, with respect to z, the use of any machine during a shift for which no operator has been assigned.

Cut generation process

Let us denote by z the vector determined, according to (15), by a feasible solution x of [ET P]. If [JobShop(z)] is unfeasible, it implies that jobs cannot be processed with the available resources defined by z. Consequently, [JobShop(z)] fails when x does not lead to a feasible solution for problem [P]. x must therefore be discarded from the set of feasible solutions.

Trivially, z (and thus x) is unfeasible because at least one machine among those who are temporarily not used should be available for executing a job. A valid cut translating this simple finding in mathematical terms is:

∑ e∈E ∑ k∈K e ∑ s∈S e βks • x eks ≥ 1 ( 17
)
where βks = 1 if zks = 0, 0 otherwise. This inequality should be added to the pool Cut of the master problem [ET P]. Since there is a finite number of eligible assignments of employees to both machines and shifts, the process stops in a finite number of steps.

Initializing the pools of valid inequalities

Introduction

In order to speed up the convergence of the process in such cut generation approaches, it is often interesting to initialize the pool of cuts with a subset of inequalities allowing to (partially) recover the feasibility when solving the master problem. In this section, we propose three different valid inequalities leading to such an initialization of the set Cut of [ET P].

Probing cuts

Basically, probing refers to an elimination technique used for enhancing MIP solvers [Sav94]. The underlying idea is to fix a decision variable to one of its bounds, and to measure the consequence in terms of logical implications. If fixing this variable leads to an unfeasibility, it clearly implies that no solution with the current variable value exists. In some cases, some variables (namely binary variables) can thus be definitively fixed.

In our approach, this idea is exploited to identify shifts during which some machines must be used. For this purpose, we specify for each machine k the list ϒ k of shifts during which k cannot be used. The exhaustive list of shifts where there is a machine which cannot be used is denoted by ϒ = ∪ m k=1 ϒ k . Notice that according to proposition 1 (see Section 1.2), each inactivity shift s on machine k can be modeled by an additional fictional job. Clearly, modeling this inactivity by a fictional job is equivalent to set zks = 0. If the resulting job-shop scheduling problem (with fixed inactivity periods) is unfeasible, then at least one inactivity pair ( k, s) involved in ϒ must be removed, i.e., the associated machine k must be used during shift s. That leads to the following so-called probing cut:

∑ e∈E ∑ k∈K e ∑ s∈S e α ks • x eks ≥ 1 ( 18
)
where

α ks = 1 if s ∈ ϒ k , 0 otherwise.
In our experiments, and for obvious computational considerations, we focused only on configurations (machine k, shift s) of ϒ corresponding to singletons (e.g. (k 1 , s 3 )) or pairs of couples (e.g., (k 1 , s 3 ) ∧ (k 2 , s 0 )). Notice that if a probing cut is generated for a given singleton, this implies that machine k has to be used during shift s. In the following, such singleton will be refered to as fixed by probing.

Minimal number of working shifts per machine

The valid initial inequalities presented in this section are based on the computing of a lower bound LB k for each machine k representing the minimum number of worked shifts which are required to ensure the feasibility of [P]. The corresponding cut simply makes sure that machine k will be used during at least LB k shifts:

∑ e∈E|k∈K e ∑ s∈S e x eks ≥ LB k k = 1, . . . , m (19)
LB k relies on the computing of an upper bound σ ′ ≤ σ on the maximal number of shifts during which machine k can be used. By symmetry, LB k can hence be defined as follows:

LB k = (σ -σ ′ ). It is computed by Algorithm 1. Algorithm 1 Computing of LB k σ ′ ← max(LB 1 k , LB 2 k ) unfeasible ← false repeat
O F ← σ ′ fictional operations which have to be processed on k for a duration of π [JobShop F ] ← job-shop problem that aims to schedule each job i ∈ J and each operation of

O F if [JobShop F ] is unfeasible then unfeasible ← true else σ ′ ← σ ′ + 1 end if until (unfeasible) ∨ (σ ′ > σ ) LB k ← σ -σ ′
At the beginning of Algorithm 1, we use two auxiliary lower bounds: LB 1 k and LB 2 k . Indeed, computing LB k from scratch can be time-consuming and those auxiliary bounds greatly improve the performances. They are defined as follows:

-LB 1 k is the number of pairs (machine k, shift s) fixed by probing (see 3.5.2):

LB 1 k = card{s ∈ S|(k, s) is fixed by probing} -LB 2
k is based on the sum of processing times of operations which have to be scheduled on machine k:

LB 2 k =
∑ j∈J ρ ik π

Preventing multiple assignments

The two groups of cuts defined above do not prevent from generating solutions assigning multiple operators to the same pair (machine, shift). Such solutions are clearly sub-optimal since assignment costs c eks are strictly positive and a machine needs only one operator to be handled. This third collection of initial cuts ensures that at most one operator is assigned to any pair (machine, shift) :

∑ e∈E|(k∈K e )∧(s∈S e )

x eks ≤ 1 k = 1, . . . , m s = 0, . . . , σ -1 (20)

Even if they could seem redundant, according to our experiments, these additional cuts are quite efficient in solving the master problem [ET P].

Algorithm

The overall process underlying the cut generation approach presented in this section is summarized in Algorithm 2: Algorithm 2 Exact cut generation process

LB ← 0 repeat x ← optimum of [ET P], with a cost: Θ x if x is defined then LB ← Θ x ∀(k, s) ∈ {K × S} zks ← min(∑ e∈E xeks , 1) feasible ← solve[JobShop(z)] if ¬ feasible then add cut (17) to the pool Cut of [ET P] end if end if until (feasible) ∨ ( x is not defined) return LB

An hybridization of a cut generation process with a branch and bound strategy

Introduction

In this section, we present a second exact solution strategy. Two classic approaches for tackling hard combinatorial optimization problems are hybridized: branch and bound and cut generation process. Like in the decomposition and cut generation process (see Section 3), we also decompose the problem [P] into two independent sub-problems. The method's main particularity relies on the fact that the global process control is ensured by a vector z of indicatrices which decides whether a machine k is available for processing jobs during a given shift or not. z is instantiated by an enumerative process.

The underlying global process is explained in Section 4.2. The main components associated with the branch and bound approach are detailed in sub-Sections 4.2.1 to 4.2.4, while implementation issues and the overall algorithm is presented in sub-Section 4.2.5.

The global process

Let us define, for any couple (k, s) of machine/shift, an indicatrice representing the availability of machine k during shift s (implying that at least one operator is assigned to k during s):

∀(k, s) ∈ {K × S} zks =    -1 if machine k is not constrained during shift s 0 if machine k must not be used during shift s +1 if machine k must be used during shift s
In this alternative approach, the vector z is exploited to guide the global process, and its instantiation is performed by means of an enumerative approach of branch and bound type, as described in Sections 4.2.1 to 4.2.4. Thus, each node of the search tree corresponds to a partial availability planning on the set of machines M. A given vector z leads to two sub-problems [ET P(z)] and [JobShop(z)] derived from the decomposition principle exploited previously. The employee timetabling sub-problem [ET P(z)] can be formulated as follows:

[ET P(z)] : min Θ z = ∑ e∈E ∑ k∈K e ∑ s∈S e c eks • x eks ∑ e∈E x eks = zks ∀(k, s) ∈ {{M × S} |z ks = -1} (21)

Cut

(4), ( 5), ( 6), ( 7)

where Cut denotes a set of valid inequalities (17), already defined in Section 3, that are iteratively added to the model.

Clearly, [ET P(z)] aims at finding a minimal cost assignment of operators to both machines and shifts, the solution space being limited to assignments x consistent with constraints ( 21). [JobShop(z)] is exactly the same as in Section 3.3.

In our approach, [JobShop(z)] is first solved. If it is feasible, the linear relaxation of [ET P(z)] is then solved up to optimality thanks to a LP solver. The corresponding optimal solution clearly gives a lower bound LB(z) associated with the current node (z).

Let us denote by UB the value of the incumbent solution (best solution found so far). If either 1) LB(z) > UB, 2) the LP-relaxation of [ET P(z)] has no feasible solution or 3) [JobShop(z)] is unfeasible, then the current distribution z is not consistent. The associated node in the enumeration search tree is discarded and a backtracking occurs.

Otherwise, the collection of implication rules detailed in Section 4.2.3 is performed, as well as a consistency test consisting in forcing any non constrained pair (machine k, shift s) (z ks = -1) to be inactive (z ks = 0), and solving the resulting sub-problems [ET P( z Section 4.2.5 provides the detailed implementation issues related to this procedure.

Branching strategy

The enumerative strategy we use is a binary branching scheme. To expand the current active node, a pair (machine k, shift s) is selected by the means of a branching heuristic. Two child nodes are then created, corresponding respectively to zks = 0 (machine k must not be used during shift s) and zks = 1 (machine k must be used during shift s).

In our implementation, the node associated to zks = 0 is systematically explored first. The reason for this choice is that prohibiting machine activity leads to a more constrained problem, and thus to a much more reduced search tree size.

The branching heuristic we use selects the pair of machine/shift ( k, s) defined as follows:

k is the machine k ∈ K with the minimal difference between the number of available shifts and the lower bound LB k of the number of worked shifts on k (see 3.5.3):

k = argmax k∈K|∃s∈S,z ks =-1 (card {(k, s) ∈ {K × S} |z ks ∈ {-1, 1}} -LB k )
-s is the most tardive shift such that no decision has been taken on pair ( k, s):

s = argmax (s ∈ S|z¯k s = -1)
Several other branching heuristics have been investigated. The one presented above has been the most efficient in the experiments.

Lower bound

As mentioned before (see 4.2), [ET P(z)] aims at finding a minimal cost assignment of operators to both machines and shifts consistent with the current distribution z. Clearly, an optimal solution to problem [ET P(z)] gives a valid lower bound associated to the related node in the search tree. For obvious computational reasons, we only solve the linear relaxation of [ET P(z)] for each node of the search tree.

Implication rules

In order to reduce the search space and thus to speed up our process, two implication rules have been defined.

Probing based implication rule

The first implication rule is based on the probing strategy detailed in Section 3.5.2. It aims at checking if successive decisions about the absence of work for some pairs (machine, shift) are not too restrictive. Such restrictions could indeed prevent any descendant of the branching node from leading to a feasible complete solution. For this purpose, each non constrained pair (machine, shift) in the current distribution z (i.e., zks = -1) is in turn forced to inactivity (z ks = 0). The resulting sub-problem [JobShop(z)] is solved, and if unfeasibility is detected, we can set zks to 1, i.e., any subsequent distribution in the sub tree z satisfies that machine k is used.

Implication rule based on the minimal number of working shifts per machine

This implication rule ensures that, for each machine k ∈ K, the lower bound LB k on the minimal number of working shifts (see 3.5.3) is verified. This logical rule can be expressed as follows:

(card {s ∈ S|z ks = 0} = LB k ) ⇒ (z ks = 1 ∀s ∈ {S|z ks = -1})

Consistency checking

As mentioned in Section 4.2, we perform a consistency test for each non pruned node of the search tree. This test aims at checking if the partial distribution z associated with the current branching node is already a complete feasible solution for the problem.

To do so, we first get a complete instantiation of z by setting zks = 0 for any pair (machine k, shift s) that the current partial distribution z does not constraint (i.e., zks = -1). For the resulting complete distribution z, we check the feasibility of both [ET P(z)] and [JobShop(z)]. If these two sub-problems are feasible, the associated optimal solutions yield a feasible solution for the global problem [P]. They can hence be compared to the best incumbent solution. Otherwise, adding the cut (17) to the pool Cut of valid inequalities in problem [ET P] prevents us from generating such a non consistent distribution by enumeration.

Algorithm

The hybrid branch and bound / cut generation method described above is summarized in Algorithms 3 to 8. Algorithm 3 is the main algorithm; it refers to algorithms 4 to 8 described hereafter. In our implementation, we used the two following specific data structures:

-P: a typical stack which is used to store the nodes of the search tree in a Last In, First Out (LIFO) order -γ: a node of the search tree. Each node γ is characterized by 5 attributes:

-k ∈ M: the machine of the selected branching variable used to construct γ -s ∈ S: the shift of the selected branching variable used to construct γ -explored ∈ {true, false}. If true, the node has already been explored and its two children have been (if necessary) generated. If false, it is the first time the node is met. -value ∈ {-1, 0, 1}: fixed value for the branching variable used to construct γ -decisions: list of the decisions (due to the branching variable selection and the implications rules) related to node γ. By construction, each descendant node of γ respects all these decisions.

Computational experiments

In this section, we first describe the instances used as a test bed and the methods used as comparison (section 5.1). We then briefly evaluate the interest of the initial inequalities (section 5.2). Finally, we investigate the experimental behavior of the methods proposed in this paper (section 5.3).

Test bed and methods

To test the methods proposed above, we used the generator of instances of [AGRV09]. We have defined 6 categories of instances of increasing size. For each category, two sets of 25 "feasible" (respectively "unfeasible") instances have been generated. The parameters used as inputs by the instance generator are displayed in Table 1.

In Table 1, column µ refers to the total number of employees, of whom µ extra are extra employees. Those extra employees ensure that the production requirements can always be mathematically met (as in [AGRV09]). However we make a distinction between the two types of instances. In the case of "feasible" instances, the production requirements can actually be met in pratice: virtually no extra employees are needed and, even though more expensive, they do not incur a prohibitive cost (they are typically outsourced personel or regular staff working overtime). In the case of "unfeasible" instances, it is not possible to meet production requirements without extra employees, and those have a prohibitive cost, making the instances unfeasible in practice (an extra employee has a cost of 1000, that dominates the sum of all other costs; hence, if a solution has a cost of e.g., 8028, that means that 8 extra employees are needed).

In Table 1, Θ * is the best known solution for [P]. It is the optimum for 57.7% of all cases; otherwise it corresponds to the best solution found by either the hybrid branch and bound / cut generation method described in Section 4 or a MIP solver applied to the formalisation [P] (within a CPU time limit of 2 hours for each method). For the 15 compatible instances picked out from [AGRV09], C max is the optimal makespan of the job-shop problem. For the 285 other generated instances, C max is computed as follows: 1) a heuristic makespan C heur max to the job-shop problem is found thanks to our job-shop solver (cf. Section 2.2), 2) the number δ of shifts is fixed to C heur max /π and 3) C max is fixed to (δ • π). We try to solve every instance with the two exact methods proposed above:

-Cut: the decomposition and cut generation approach described in Section 3 ; -HyBB: the hybrid branch and bound / cut generation method described in Section 4.

As a basis for comparison, we also try to solve each instance with :

-MIP: a MIP solver applied to the formalisation [P].

-AGRV09: a slightly modified version of the method proposed in [AGRV09]. Indeed, Artigues et al. solve a lexicographic optimization problem where the makespan minimization is the primary objective whereas the employee cost minimization is the secondary objective. We thus define AGRV09 as the second stage, the minimal makespan value being directly fixed in a preprocessing stage. -Heur: a heuristic solution. A heuristic solution to the job-shop problem is found thanks to our job-shop solver (cf. Section 2.2) and every pair (machine, shift) that is worked is fixed. The employees are then assigned optimaly using a MIP solver. (Note: a simpler heuristic, where employees are assigned greedily instead of optimaly, has also been tried out; it is barely faster than Heur and produces drastically worse solutions, and thus it has been discarded.)

CPU time has been limited to 5 minutes for the two categories of small instances ejs6 × 4 × 25 and ejs6 × 6 × 25. Is has been limited to 10 minutes for the four categories of big instances ejs8 × 8 × 40, ejs8 × 8 × 50, ejs10 × 10 × 40 and ejs10 × 10 × 50.

In the remaining (if it is not specified), the initial inequalities of Section 3.5 are added to the three methods Cut, HysBB and MIP. We did not experiment the interest of the initial cuts for AGRV09.

Initial inequalities

To evaluate the initial inequalities described in Section 3.5, we first compare in Table 2 For each type and each category of instances, the LP-relaxation is much better with the initial inequalities: the average ratio LP/Θ * is indeed 17% higher when the initial inequal-ities are used, and the extra computing time remains moderate (24s in average and about 128s for the worst case).

To further chalenge the initial inequalities, we compare the rate of success of the three exact methods (Cut, HyBB, MIP), depending on whether the initial inequalities are used (Table 3). It is very clear from these results that the initial inequalities proposed in Section 3.5 strongly improve all three methods. The improvements is particularly important for Cut and HyBB. Indeed, these two methods cannot solve instances larger than 6 jobs without initial inequalities. As a matter of fact, the initial cuts are an integrated part of the two exact methods Cut and HyBB proposed in this article and they are nicely complementary to the other components of those methods.

As stated before, the initial inequalities are added to each method described in the remaining (except for AGRV09).

Exact methods

In this section, we compare the resultats of the two exact methods we have proposed in this paper (Cut and HyBB) with other approaches (MIP, AGRV09 and Heur).

Table 4 reports the performances of each method. Column success gives the percentage of instances solved to optimality within the time limit. Column time displays the average computation time. Column gap gives the average deviation to the best known solution Θ * 2 . The statistical significance of those results is verified by an exhaustive computation of Kruskal-Wallis tests for every pair of methods and every measure; when relevant, p-values are provided in the following discussion.

Table 5 compares the successes of the different methods. This table can be read as a matrix a i j where a i j is the number of instances the method i solves to optimality within the time limit whereas the method j does not. a ii gives the number of successes of the method i.

Furthermore, Table 6 provides a direct comparison between the two best performing methods: HyBB and MIP. In this table, success is the number of instances solved to optimality by both methods; time(S) is the CPU time in case of success (that is for the instances solved to optimality by both approaches); and gap(F) is the average deviation to the best known solution in case of failure (that is when at least one method fails). From all those experiments, the first general conclusion is that HyBB is the best performing methods on any kind of instances. Then come MIP and Cut, then AGRV09 and the poorest method is Heur. Typically, HyBB solves to optimality more instances than the other methods (but MIP on small instances) and besides it solves them faster when it succeeds or provides tighter bounds when it fails.

From Table 4, it is clear that HyBB achieves better success and gap than Heur, and this is statistically significant (p-values < 10 -7 ). HyBB is also better than AGRV09 (all p-values < 10 -3 ). Compared to Cut, HyBB is more successful and provides tigher gap (p-values < 0.05 for success, and < 10 -7 for gap); it also seems to be faster, but this is not statistically significant (p-value> 0.6 on feasible instances). Furthermore, none of Cut, AGRV09 or Heur solves to optimality an instance that HyBB does not (see Table 5).

As a consequence, HyBB is only challenged by MIP for the title of best method. From Table 4, it is not clear whether HyBB has better success than MIP; indeed, the two methods cannot be distinguished on this criteria (p-value> 0.8), nor on the computation time. The dominance of HyBB is however quite clear for the gap (p-values < 0.05). Table 6 allows for a more detailed analysis. In case of success, HyBB is quicker (however, the statistical relevance is not so good). In case of failure, HyBB provides a significant better gap (p-values < 0.003).

All in all, HyBB performs generaly better than MIP, with exceptions. It is not dominated for the success, nor for the computation times, and achieves better gaps. In particular, on hard instances that none of the methods solve, HyBB is clearly better: the gap is significantly lower, and does not explode when the size of the problem increases; moreover it must be noted that HyBB provides feasible solutions for all 300 instances, while MIP failed for 74 of them. Nevertheless, on some instances, MIP may outperform HyBB (e.g., for small instances).

Similar detailed analysis confirm 1) that HyBB also outperforms Cut, AGRV09 and Heur (all results are statistically significant); 2) that MIP outperforms Cut (except for big instances), AGRV09 and Heur; 3) that Cut outperforms AGRV09 and Heur.

Note that Heur, which is never very good, is very consistent in providing very quickly feasible solutions with a gap of about 30%. Such solutions are not good for small instances, but become competitive for the big ones. Actually, Heur outperforms every method but HyBB on the most challenging instances. This is essentially due to the fact that, within the given time, methods Cut, MIP and AGRV09 may fail to find a feasible solution for such instances.

One should notice here that the comparison with [AGRV09] is not totally fair. Indeed, [AGRV09] solves more general problems than the methods described in this article. We remind the reader that [AGRV09] tackles a problem where each job can require more than one employee to be processed. It also addresses a problem with two hierarchical objectives: 1) minimizing the makespan and 2) minimizing the employee cost. Furthermore, AGRV09 is not competitive at all on unfeasible instances because of the very high costs of the extra employees. The method (dichotomy) has not been fitted for such test problems.

To conclude with these experiments, we point out that really large instances of the problem tackled in this paper are intractable with the current methods. Indeed one can observe that none of the methods experimented here can solve instances with more than 10 jobs and 10 machines; starting from 8 jobs and 8 machines, all methods have serious difficulties. This is not a surprise since the job-shop sub-problem itself is awfully difficult to solve, even with only 10 jobs and 10 machines.

Conclusion

In this paper we have proposed two different procedures to solve an integrated employeetimetable and job-shob-scheduling problem, and we have evaluated them through computational experiments. In the process, we have outlined initial cuts that significantly help the finding of a solution for our own methods, but also for a standard MIP approach.

Our first method, the decomposition and cut generation procedure (Cut), is rather disappointing. However, its mixed results can easily be explained: the generated cuts do not contain enough global information and tend to eliminate only the current solution. Stronger cuts would be necessary; however the design of such cuts appears to be quite difficult and unpredictable perspectives (a tighter collaboration with the job-shop solver) would probably be the key. This phenomenon illustrates the fact that integrated problems require ad-hoc, well-tuned procedures that are capable of using the very particularities of the problem.

The second method, the hybridization of a cut generation process with a branch and bound strategy (HyBB) provides very good results and greatly improves over the performances of the best methods of the literature.

As for future research directions, the relevance of the hybridization of a cut generation process with a branch and bound strategy should be proved on other integrated problems, e.g., that incorporate different scheduling and/or timetabling constraints. Furthermore, variants of this hybridization should be tried out to provide a fast initial convergence that would allow to quickly get good feasible solutions.

  )] and [JobShop(z)]. If the solution ȳ of [JobShop(z)] is once again feasible, then the optimal solution ( x) of [ET P(z)], if defined, yields a feasible solution ( x, ȳ) of value Θ z to the global problem [P]. If Θ z < UB, then an improved solution to [P] has been exhibited: it becomes the new incumbent solution, and UB is updated. Otherwise, a cut (17) is added to the pool Cut in problem [ET P]. Notice that inequalities (17) are valid for any node of the search tree.

Table 1

 1 Parameters of the 300 instances (25 instances per category)

	type	category	n	m	µ (incl. µ extra )	µ extra	# max skills per employee	π	C max	Θ *
		ejs6 × 4 × 25	6	4	25	10	2	8	∈ [37, 56]	∈ [17, 33]
	feasible	ejs6 × 6 × 25 ejs8 × 8 × 40 ejs8 × 8 × 50 ejs10 × 10 × 40	6 8 8 10	6 8 8 10	25 40 50 40	10 20 20 20	4 4 4 4	10 8 10 8	∈ [48, 70] ∈ [44, 64] ∈ [60, 80] ∈ [72, 80]	∈ [24, 73] ∈ [70, 106] ∈ [40, 64] ∈ [174, 286]
		ejs10 × 10 × 50	10	10	50	20	4	10	∈ [80, 100]	∈ [81, 133]
		ejs6 × 4 × 25	6	4	25	21	2	8	∈ [40, 56]	∈ [8009, 12013]
	unfeasible	ejs6 × 6 × 25 ejs8 × 8 × 40 ejs8 × 8 × 50 ejs10 × 10 × 40	6 8 8 10	6 8 8 10	25 40 50 40	19 32 42 30	4 4 4 4	10 8 10 8	∈ [50, 70] ∈ [64, 64] ∈ [60, 90] ∈ [72, 80]	∈ [8028, 18015] ∈ [22040, 32034] ∈ [19040, 28031] ∈ [39056, 50055]
		ejs10 × 10 × 50	10	10	50	40	4	10	∈ [80, 100]	∈ [35056, 47068]

Table 2

 2 the results of the LP-relaxation of [P] with and without the initial inequalities. In this table, column LP/Θ * gives the deviation between the optimum of the LP-relaxation of [P] (LP) and the best known solution for [P] (Θ * ). Impact of the initial inequalities on the LP-relaxation. Instances are not distinguished by type, since the results are almost identical in both cases.

	category	linear relaxation (LP) of [P] LP/Θ * time	linear relaxation (LP) of [P] with initial inequalities # initial preprocess LP total LP/Θ * inequalities time time time
	ejs6 × 4 × 25	83.3%	0.3s	98.1%	48.3	0.2s	0.2s	0.4s
	ejs6 × 6 × 25	66.9%	1.0s	87.8%	69.3	1.1s	1.0s	2.1s
	ejs8 × 8 × 40	65.8%	4.4s	87.7%	127.8	5.7s	3.9s	9.6s
	ejs8 × 8 × 50	70.3%	8.7s	85.9%	123.4	10.3s	9.9s	20.2s
	ejs10 × 10 × 40	62.4%	21.5s	82.3%	211.6	49.1s	22.3s	71.4s
	ejs10 × 10 × 50	61.9%	44.4s	70.6%	194.5	74.7s	52.3s	127.0s
	Total	68.4%	13.4s	85.4%	129.1	23.5s	14.9s	38.4s

Table 3

 3 Impact of the initial inequalities on the exact methods. Instances are not distinguished by type, since the results are almost identical in both cases.

Table 6

 6 Comparison of HyBB and MIP in case of success or failure (average values)

	2 gap = min 1, |value-Θ * | Θ *	if a solution value has been found, 1 otherwise

In our implementation, those variables are not created. Variables x eks are thus created if and only if the employee e can work on machine k (k ∈ K e ) and is available on shift s (s ∈ S e )

Acknowledgements

The authors wish to thank Christian Artigues for kindly giving us access to his code and test problems. Thanks are also due to the anonymous referees who provided useful comments and suggestions that led to improvements in the paper.

A Computing environment

All experiments have been carried out on a standard PC (Intel(R) Core(TM) i3 CPU M 370 @ 2.40GHz, 2.39 GHz, 3.42 Go RAM) running MS Windows XP. Our own procedures are implemented in Java. The software of [AGRV09] is the original one, kindly provided by the authors, and is written in C++. Mixed integer and/or continuous linear programs have been solved with Ilog Cplex 12.2; for constraint based scheduling, we used Ilog Solver 6.7 and Scheduler 6.7.


B Lebental

P Chainais

P Chenevier

E Delevoye

J-M Fabbri

S Nicoletti

P Renaux

A Ghis

Carbon nanotubes based ultrasonic microtransducers for durability monitoring

Keywords: 07.07.Df, 68.37.-d, 85.35.Kt, 85.85.+j

Structural health monitoring of porous materials such as concrete is becoming a major stake in our resource-limited economy, as it conditions durable exploitation of existing facilities. Durability in porous materials depends on nanoscale features which need to be monitored in-situ with nanometric resolution. To address this problem, we put forward an approach based on the development of a new nanosensor, namely a capacitive micrometric ultrasonic transducer whose vibrating membrane is made of aligned single-walled carbon nanotubes (SWNT). Such sensors are meant to be embedded in large number within a porous material in order to provide with information on its durability by monitoring i n -situ neighboring individual micropores. In the present paper, we report on the feasibility of the key building block of the proposed sensor: we have fabricated well-aligned, ultra-thin, dense SWNT membranes that show above-nanometer amplitudes of vibration in a large range of frequencies spanning from 100 kHz to 5 MHz.

Introduction

Because cement manufacturing represents 7% of the world-wide industrial energy consumption and 4% of the worldwide industrial CO 2 emission, reducing its ecological footprint has become an environmental priority [1]. By improving structural efficiency and durability, nanotechnologies open the path toward a greener construction industry [START_REF] Van Broekhuizen | Nano-products in the European Construction Industry -State of the Art[END_REF]. Two domains appear most promising: nanocomposites construction materials [3], which have shown both high performance [4] and new functionalities [1], and nanosensors, which are used for structural health monitoring, notably for the measurement of macroscopic structural properties like damage, strain or corrosion [5]. The in-situ monitoring of micro and nano-scale features by a large number of embedded nanosensors is not yet addressed, although the random micro-and nanoporosities of cement-based structures play a most significant role in their durability issues [6]. To fill this gap, a relevant path consists in probing individual capillary pores (10 nm to 10 µm large water-or air-filled pores) of cementitious materials by using ultrasonic waves generated by embedded Micromachined Ultrasonic Transducers (MUT) (figure 1a) [7]. Such approach fosters the use of MUT devices with an active surface in the 1 µm² range. As state-of-the-art MUT devices do not meet these surface requirements [8], we designed a new type of capacitive MUT (cMUT). In a cMUT device, ultrasonic transduction is achieved by using the vibrations of a conducting membrane suspended above an actuation electrode. For our application, the span of the membrane had to be lowered by about one order of magnitude compared to existing devices, down to a micrometric size. To maintain significant amplitudes of vibration [START_REF] Timoshenko | Theory of elasticity[END_REF] and thus avoid sensitivity loss, the thickness of the membrane also had to be reduced to a nanometer size. With conventional materials used for microelectronics, the fabrication of such a suspended membrane with micrometric size and nanometric thickness for ultrasonic transduction purposes presented a real challenge. We overcame it by using single-walled carbon nanotubes (SWNTs) instead as membrane material: in the new cMUT device discussed here, the ultra-thin membrane consists of a dense, fewlayered assembly of well-aligned SWNTs [START_REF] Lebental | [END_REF] (figure 1b). This paper presents some novel results we obtained while realizing such SWNT-based cMUT devices: after a brief description of the fabrication process, we present a quantitative analysis of the significant morphological features of the membranes as observed on SEM and AFM images. We also study the device vibrational behavior by laser vibrometry.

Figure 1. Instrumentation principle: (a) schematic representation of a cubic sample of cement gel (matrix in gray, air-filled pores in white, water-filled pores in blue) instrumented with a large number of randomly embedded ultrasonic microtransducers (black squares); (b) ultrasonic microtransducer (black rectangle) probing a water-filled pore in pulse-echo mode (emitted waves in red, reflected waves in green); c) Schematic representation the microtransducer: (top) single-walled carbon nanotubes are aligned to form the dense membrane of a capacitive ultrasonic microtransducer; (down) due to the electrostatic interaction between the gate electrode and the metallic membrane, the membrane vibrates upon application of a varying gate voltage.

Device fabrication

The devices are fabricated using SWNTs from Carbon Solution Inc synthesized by electric arc from Ni/Y catalysts (the AP-SWNT product). SWNT are first purified by acidic treatment following the process described in [11] and then dispersed in N-methyl-2-pyrrolidone at low concentration using an ultrasonic bath (typically 70 W for 60 min). The still partially bundelized [12] nanotubes are deposited by dielectrophoresis [13]: an AC electric field (typically 5 V/µm at 5 MHz during 5 min) is applied between two parallel, 1 mm long, metallic (Pt or Ti/Au) electrodes deposited on a SiO 2 layer and separated by 1 µm or 2 µm; polarized by the electric field, the nanotubes are attracted to the electrodes and settle perpendicularly to them (figure 2a). After deposition, the SWNT are anchored to the electrodes by 50 nm thin, 20 µm large ribbons of Pt. These ribbons are patterned by e-beam lithography of PMMA 4% resin, evaporation of 50nm of Pt at 0.1 nm/s and acetone lift-off of the resin. Finally the membrane formed by the anchored SWNTs is released by removing the SiO2 layer either by dry etching (typically a focalized CHF3 plasma applied for 20 s at 40 W in a RIE reactor NEXTRAL NE 90A) or by wet etching (typically HF 10 % during 5 min) (figure 2b). We have no evidence of SWNT degradation after etching. However we observed that a longer exposure (more than a few minutes) of the SWNT to the CHF3 plasma will result in their complete etching. The material that glues the bundles together is interpreted as a thin PMMA layer remaining from the patterning process of the anchoring layer. Note that in both SEM pictures, the sample is tilted by 45° to the electronic beam.

Highly aligned and dense SWNT membranes

Optimization of the dielectrophoresis process

Since the application is very demanding regarding to the membrane properties, the capabilities of the dielectrophoresis deposition process in terms of SNWT alignment and density have been extensively studied: over eighty depositions were realized and the resulting membranes were then imaged by SEM, until a proper control over the membrane features was achieved. Several images (3 to 15) were needed for a proper characterization of each deposition, as the deposition area ranges from 1,000 to 2,000 µm², whereas each SEM image corresponds to a surface of only 1 to 2 µm². The whole set of SEM images (660 images) was taken in a SEM Hitachi S4000 or 4100 under 30 kV acceleration voltage and 12 µA emission current on samples tilted by 45° to the electronic beam with a magnification chosen between 30,000 and 60,000. Considering the large quantity and the complexity of the images, the often used "by hand" method [14] was not sufficient for quantitative analysis and objective comparisons. An automated image processing technique providing statistical quantities related to nanotubes orientation and density appeared highly desirable to validate the membrane fabrication process. It had to be able to process large sets of images with significant diversity.

SEM images processing technique

We implemented such a processing sequence (figure 3a-d). While requiring very little user supplied information, our procedure extracts quantitative characteristics of an assembly of filamentary objects. The initial images (figure 3a) are first denoised using the software GREYCStoration (based on an anisotropic diffusion equation) [15]. Resulting images are analyzed using the software MATLAB, and more specifically its IMAGE and DIFFC [16] toolboxes. Two main issues have to be solved, namely the segmentation of the filaments skeleton and the estimate of their local orientation.

To begin with, the images are convoluted by short lines with every different possible orientation. Because nanotubes (or nanotubes bundles) are linear objects, this operation enhances the nanotubes, enabling the extraction of the pixels of the image that belong to nanotubes (figure 3b). The skeleton of the nanotubes is obtained by morphological thinning of each nanotube down to a line (figure 3c). The ratio of the skeleton total length to the gap width yields the number of nanotubes (or bundles) in the image, from which we derive the nanotubes linear density d l (number of nanotubes by unit length of deposition). The orientation of each nanotube is derived at the position of each pixel of its skeleton by using local structure tensors (figure 3d). We thus obtain the skeleton angular distribution p(α), defined by the percentage of the nanotubes total skeleton length whose orientation is comprised between α and α+180/256°, for -90°<a<90°-180/256°. From p(α), we determine the order parameter of the deposition S=<cos²(α)>=∫p(α)cos²(α)dα, similar to the order parameters used in [17]. Note that the image tilt is taken into account in the analysis.

Results of the SEM image analysis

This image processing method enables the quantitative comparison of various depositions and thus the optimization of the deposition process. It makes the whole device fabrication process more efficient, as it enables the objective selection of only the samples with the best-aligned and densest nanotubes to be post-processed (anchoring and suspension steps). The result of the image processing applied to the whole set of images is presented figure 4: the average order parameter (0.5 in absence of any preferential order, 1 in case of perfect ordering) is equal to 0.74. Over single images, best alignment was found to be as high as 0.93: more than 80 % of the nanotubes are deposited with a tilt of less than 20° to the horizontal axis. These values indicate the high quality of the SWNTs alignment [18]. In the literature, the linear density of an assembly of aligned nanotubes (bundles) is found to range from 2 bundles by micrometer (by in-place growth [19]) to 60 bundles (by Langmuir-Blodgett [20]) by micrometer. Here, the density of nanotubes reaches the intermediate value of 19 bundles by micrometer. Deposited layers appear to be made up of a single layer of SWNTs bundles (fig. 2a; fig. 3a-d), probably as a result from the screening effect commonly seen in SWNTs dielectrophoresis [21].

Ultra-thin SWNT membranes

As membrane thickness is critical to the application, we evaluated before suspension the thickness of three chosen samples by AFM profilometry in tapping mode [22]. We used a Multimode AFM equipped with a Nanoscope IIIa electronics from Bruker AXS with Aluminium-coated silicon tips (TAP300 Al-G, Budget Sensors). The typical tip radius is lower than 10 nm. 512x512 pixels images of the sample height were acquired in tapping mode at the typical scan rate of 0.5 Hz (2 seconds by line). Pixel by pixel error on the height measurement is lower than 0.5 nm. We defined the sample thickness as the average of the height distribution, the substrate height being chosen as reference height. The convolution of the shape of the tip with small-size features (small compared to the tip radius) leads to an overestimation of the thickness (broadening effect). To correct this effect, images were first processed before the thickness is computed. Among the three samples, two showed an average thickness of 10.7 nm (figure 5a). The third one showed an average thickness of 1.4 ± 0.4 nm over a surface of about 0.06 µm² (figure 5b): hence this sample appears to be made of a monolayer of individual SWNTs. It is an exceptional feature among reported devices made of deposited SWNTs [20].

Above-nanometer amplitudes of SWNT vibrations

Laser vibrometry

To test the suitability of the fabricated devices for acoustic wave generation, we studied their vibratory behaviour. SWNT vibrations are often detected by electrical methods [23], where the amplitudes of vibration are determined indirectly by measuring some electrical parameters and feeding them to theoretical models. To measure directly the vibration amplitudes, an alternative method consists in using AFM vibration imaging [24]. Here, we chose to implement scanning laser Doppler vibrometry. The amplitudes of the membrane out-of-plane vibrations are measured by a scanning laser (He-Ne) Doppler vibrometer Polytec MSA-400. Points are typically acquired on a grid with controllable spacing (from 250 nm to 1 µm in both directions). The instantaneous vertical speed of the membrane at each grid point is measured during 100 periods at a rate of 10 acquisitions by period. Then Polytec's dedicated PSV software derives the frequency spectrum of vibration and calculates the RMS amplitude of vibration (quadratic mean of the position over the measurement duration) at each grid point. Until the present work, this technique had not been used with SWNT, as the laser beam is much larger than the SWNT diameter. Here, its use is enabled by the assembly of the SWNTs into a membrane with a density high enough to reflect the laser beam.

Above nanometer amplitude of vibration

For the measurements, we considered a membrane suspended over a 2 µm wide trench and placed in air at ambient pressure and room temperature. The membrane was actuated electrostatically [23] by applying an AC voltage of typically 8V peak-to-peak and an additional DC bias of typically 1V to the actuating electrode underlying the 300 nm deep trench. The membrane itself was grounded via the top electrodes. Under these conditions, we measured the vibrations of the membrane over a large range of frequencies spanning from 100 kHz to 5 MHz. Figure 6 (left) shows an example of vibration profiles obtained at 100 kHz. The maximum peak-to-peak amplitude of vibration measured along the membrane was 5.6 nm at 3.7 MHz (figure 6

SWNT-based cMUT devices for durability monitoring in multiphasic porous materials

According to an advanced numerical microfluidic model of the proposed cMUT devices [START_REF] Lebental | [END_REF], their low-frequency amplitudes of vibration are expected to be four orders of magnitude higher in air than in a liquid environment. Their first resonance frequency in water should be in inverse proportion to the pore size. Hence, if the proposed sensors were to be embedded in large number in a multiphasic porous material (concrete, bones or oil reservoirs), they could differentiate air-filled micropores from other liquidfilled pores and determine the size of water-filled pores. Such in-situ data could then be effectively exploited to monitor hydration, water-related degradations and durability [27] with an unprecedented spatial resolution.

Conclusion

As a summary, we fabricated particularly well-aligned, ultra-thin, dense suspended membranes of SWNT. To our knowledge, some of our membranes feature the highest order parameter and lowest thickness ever reported for dielectrophoresis deposition of SWNT. By using scanning laser Doppler vibrometry, we showed that these membranes can feature above-nanometer amplitude of vibrations. It appears to be the first report on laser vibrometry being applied to measure the amplitude of vibrations of SWNTs. The large amplitudes of vibrations we measured indicate the feasibility of SWNT-based cMUT devices. They are the first type of nanodevices whose use is realistically envisioned for in-situ durability monitoring of porous materials at the micro-scale.

Figure 2 .

 2 Figure 2. a) SWNT deposition: (top) schematic view of the substrate and of the deposited SWNTs; the electrode material can be either Ti/Au, as displayed here, or Pt. (bottom) SEM picture of a membrane of well-aligned, thin SWNT bundles deposited by dielectrophoresis between Ti/Au electrodes on a SiO 2 layer. b) Suspension of the SWNT membrane: (top) schematic view of the device featuring suspended SWNTs anchored on the electrodes by a thin metallic ribbon; (bottom) SEM picture of a membrane of SWNT bundles deposited between Ti/Au electrodes and suspended by RIE etching of the SiO 2 layer.The material that glues the bundles together is interpreted as a thin PMMA layer remaining from the patterning process of the anchoring layer. Note that in both SEM pictures, the sample is tilted by 45° to the electronic beam.

Figure 3 :

 3 Figure 3: application of the image processing technique on a SWNT deposition of average quality (linear density 7 SWNT/µm and order parameter 0.74): a) SEM picture of the considered SWNT membrane (Magnification: 24,900; 45° tilt to the electronic beam). b) Extraction of the nanotubes pixels (in white) from the background (in black) by anisotropic contrast enhancement. c) Extraction of the skeleton by morphological thinning. d) Determination of the nanotubes orientation by gradient computation (colour code: from red (-90°) to light green (0°, i.e. alignment) to blue (+90°)).

Figure 4 :

 4 Figure 4: Result of the image processing applied on 660 SEM images taken on 87 different depositions: for each image, the graph presents the order parameter <cos²(α)> of the nanotubes angular distribution versus linear density (black dots); the pink square marks the position of the average of the order parameter and linear density distributions. The pink errors bars feature the standard deviation of these two distributions.

Figure 5 :

 5 Figure 5: Height distribution of two SWNTs membranes obtained from AFM images: a) 10.7 nm and b) 1.4 nm thick membranes. The height of the SiO2 layer is the reference ("zero") height. The height distribution is expressed in percentage of the pixels number. The green line marks the average of the height distribution, which defines the sample thickness.

  right). Such an amplitude of vibration, about as high as 4 times the diameter of the SWNTs composing the membrane (1.4 nm), makes t h e d e v i c e c o m p a t i b l e w i t h low-frequency (below 20 MHz) [25] vibrations detection system.

Figure 6 :

 6 Figure 6: Demonstration of membrane vibrations using scanning laser vibrometry. (Left top) top view and (left bottom) side view of the RMS amplitude of vibration of a 2 µm wide membrane at 100 kHz under a bias voltage V DC =1 V and an alternating voltage of V AC =8 Vpp. The membrane, delimited by the two horizontal dashed lines, vibrates with different amplitudes along its length, depending on its local properties such as thickness, bundles crossing and trapped impurities. The superimposed white grid shows the effective points of measurement by the scanning vibrometer. The RMS amplitude of noise is 30 pm (RMS amplitude of vibration measured on the motionless anchoring layer). The maximal RMS amplitude of vibration (measured approximately at the center of the dashed circle) is 200 pm. (Right) Focus on the membrane location with largest amplitude of vibrations (located within the red circle): Frequency spectrum of vibration measured under a bias voltage V DC =1 V and an AC voltage V AC =8 Vpp applied at 3.7 MHz. As a bias voltage is applied, the vibration spectrum is made of two Dirac peaks signaled by the red lozenges, one at 3.7 MHz, the second at twice this frequency, 7.4 MHz. RMS and peak-to-peak amplitudes of vibration are respectively 1.7 nm and 5.6 nm. Noise level is below 65 pm.

The first author performed the work while at CEA-LETI and IFSTTAR.

The second author performed the work while at LIMOS UMR CNRS 6158, Clermont Université, France.

Acknowledgments

The authors thank J. C. Gabriel for helpful discussion. This work was in part achieved at the PTA (Up-line Technological Platform) facility in Grenoble (http://www.pta-grenoble.com/).