Claude Yugma
email: claude.yugma@emse.fr

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

Christian Artigues

Alexandre Derreumaux

Olivier Sibille

A Batching and Scheduling Algorithm for the Diffusion Area in Semiconductor Manufacturing

Keywords: Batch, Scheduling, Disjunctive graph, Local search, Simulated Annealing, Wafer fabrication

This paper proposes an efficient heuristic algorithm for solving a complex batching and scheduling problem in a diffusion area of a semiconductor plant. Diffusion is frequently bottleneck in the plant and also one of the most complex areas in terms of number of machines, constraints to satisfy and the large number of lots to manage.

The purpose of this study is to investigate an approach to group lots in batches and to schedule these batches on machines. The problem is modeled through a disjunctive graph formulation. A constructive algorithm is proposed and improvement procedures based on iterative sampling and Simulated Annealing are developed. Computational experiments, carried out on actual industrial problem instances, show the ability of the iterative sampling to significantly improve the initial solution. The Simulated Annealing enhances the results of the iterative sampling. The constructive algorithm has been embedded in a software and is currently being used.

Introduction

Semiconductor wafer fabrication can be described as a multistage process with re-entrant flows. The processing is done layer by layer. Each layer requires several steps of processing such as chemical-mechanical polishing, diffusion, film deposition, photolithography, implant (doping) and etching. For each of the product types, and depending on the technology, a wafer goes through more than 400 process steps over a period of a few weeks. Wafer fabrication planning and scheduling is a complex task due to the large number of products and machines involved. It is further complicated by additional constraints such as re-entrant flow of operations Kumar [START_REF] Kumar | Scheduling semiconductor manufacturing plants[END_REF], setup issues, preventive maintenances, random machine breakdowns, Ovacik and Uzsoy [START_REF] Ovacik | Decomposition methods for complex factory scheduling problems[END_REF], Sze [START_REF] Sze | Semiconductor devices: Physics and technology[END_REF], etc. The importance of scheduling on the performance of semiconductor wafer fabrication facilities (fabs) is known for many years see Wein [START_REF] Wein | Scheduling semiconductor wafer fabrication[END_REF] and Sze [START_REF] Varadarajan | A survey of dispatching rules for operational control in wafer fabrication[END_REF].

In this paper, we focus on an important part of the manufacturing process. Among the complex operations involved in the fabrication of a wafer, the diffusion phase is of critical importance since the batching decisions that are involved may affect the performance of the entire wafer fab Ibrahim et al. [ICN + 03] and Mönch and Habenicht [START_REF] Mönch | Simulation-based assessment of batching heuristics in semiconductor manufacturing[END_REF]. The processing time of the operations in the diffusion area can be large (10h) compared to other operations in the fab. Mehta and Uzsoy [START_REF] Mehta | Minimizing total tardiness on a batch processing machine with incompatible job families[END_REF] state that optimizing batching operations results in good performance measures of the whole production process. Lots regularly arrive in the diffusion area and the diffusion phase is primarily used to alter the type and level of conductivity of semiconductor materials. The purpose of this article is to develop efficient methods to partition lots in batches and to schedule batches on machines in the diffusion area of a semiconductor plant while taking into account numerous constraints and optimizing three main production criteria: maximize the number of operations (moves), maximize the batch sizes and minimizing the total tardiness.

The remaining sections of this paper are organized as follows. In Section 2, we provide some background on existing related batch and scheduling problems. In Section 3, the problem is stated. We present in Section 4, a disjunctive graph representation of the problem. The method for computing an initial solution and improvement procedures based on iterative sampling and Simulated Annealing are described in Section 5. Experimental results on real problem instances and comparison with a literature algorithm are given and discussed in Section 6. Section 7 concludes the paper with recommendations for further research.

Previous related work

The operations of batching and scheduling jobs is a common practice in a manufacturing system especially in semiconductor manufacturing systems, see Mathirajan et al. [START_REF] Mathirajan | A literature review, classification and simple metaanalysis on scheduling of batch processors in semiconductor manufacturing[END_REF]. For examples reasons for batching are the reduction of setups, the ability of machines to process several jobs simultaneously, etc. With the classification given by authors by Mathirajan et al. [START_REF] Mathirajan | A literature review, classification and simple metaanalysis on scheduling of batch processors in semiconductor manufacturing[END_REF], we notice that there is not much literature on scheduling batch problems taking into account the re-entrance features of the system as for example we can name Cigolini et al. [START_REF] Cigolini | A new dynamic look-ahead scheduling procedure for batching machines[END_REF], Mason and Oey [START_REF] Mason | Scheduling complex job shops using disjunctive graphs: a cycle elimination procedure[END_REF] [START_REF] Oey | Scheduling batch processing machines in complex job shops[END_REF]. The work proposed by Mönch et al. [MFDP + 09] presents a scheduling problems in semiconductor manufacturing. Typical scheduling problems found in semiconductor manufacturing are identified in particular batching problems. The problem addressed in this paper deals with:

• Multiple machines at each stage. The authors in Mehta and Uzsoy [START_REF] Mehta | Minimizing total tardiness on a batch processing machine with incompatible job families[END_REF] presented a problem of minimizing total tardiness on a batch processing machine with incompatible job families. They proposed a dynamic programming to solve the problem. The paper of Balasubramanian et al. [START_REF] Balasubramanian | Genetic algorithm based scheduling of parallel batch machines with incompatible job families to minimize total weighted tardiness[END_REF] extends solutions of Mehta and Uzsoy [START_REF] Mehta | Minimizing total tardiness on a batch processing machine with incompatible job families[END_REF] to a batching problem with incompatible jobs on parallel machines which aims at minimizing weighted total tardiness. The authors in Kim et al. [START_REF] Kim | Scheduling wafer lots on diffusion machines in a semiconductor wafer fabrication facility[END_REF] focused on a problem of scheduling wafer lots on diffusion workstations in a semiconductor wafer fabrication facility. In their considered diffusion workstation, there are multiple identical machines, and each of them can process a limited number of wafer lots at a time. Their scheduling problem at the diffusion workstation is a problem of scheduling multiple job families on identical parallel batch-processing machines.

• Multiple stages. The described problem for the diffusion area takes into account all the different stages of this area which can reach 4. In the literature, the number of stages does not generally exceed 2. The author Su [START_REF] Su | A hybrid two-stage flow shop with limited waiting time constraints[END_REF] considered a hybrid two-stage flow shop with a batch processor in stage 1 and a single processor in stage 2. This is also the case in articles Sun and Min [START_REF] Sung | Scheduling in a two-machine flowshop with batch processing machines for earliness/tardiness measure under a common due date[END_REF] and Sung and Kim [START_REF] Sung | Minimizing due date related performance measures on two-batch processing machines[END_REF] which consider a batch problem on 2 stages. Oulamara et al. [START_REF] Oulamara | Flowshop scheduling problem with batching machine and task compatibilities[END_REF] studied a two-machine flow shop problem with conventional and batching machines in the first and second stage, respectively, and arbitrary job compatibilities.

• Multiple criteria. The goal of the paper is to simultaneously optimize three indicators. These indicators are the number of lots going through the line (to maximize), the number of wafer in a batch (to maximize) and the waiting time of lots (to minimize). Generally, articles in the literature tackled scheduling batch problems by considering one single indicator to optimize. Uzsoy in [Uzs95] tackled the problem with the objective of minimizing the completion time of jobs, Hung [START_REF] Hung | Scheduling of mask shop e-beam writers[END_REF] has an objective of maximizing the batch processing machine utilization. Perez et al. [START_REF] Perez | Minimizing total weighted tardiness on a single batch process machine with incompatible job families[END_REF] studied the problem of minimizing the total weighted tardiness on a single machine and Mathirajan and Sivakumar [START_REF] Mathirajan | Minimizing total weighted tardiness on heterogeneous batch processing machines with incompatible job families[END_REF] focused on the minimization of the total weighted tardiness on heterogeneous batch processing machines under condition of dynamic arrival of jobs, incompatible job families and nonidentical job sizes. Few articles deal with different criteria simultaneously. Furthermore, in our problem, setup times (loading and unloading) are not depending on the sequence and the complexity is added by the maximum times lags between batching operations.

Our problem is a complex variant of job shop problem. These types of problems can be found in Mason et al. [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF] and Ovacik and Uzsoy [START_REF] Ovacik | Decomposition methods for complex factory scheduling problems[END_REF]. These problems are frequently analyzed using a method known as the shifting bottleneck (SB) procedure see for example Adams et al. [START_REF] Adams | The shifting bottleneck procedure for job shop scheduling[END_REF]. Several aspects, like identifying appropriate sub problem solution procedures can be found in Demirkol and Uzsoy [START_REF] Demirkol | Decomposition methods for reentrant flow shops with sequence dependent setup-times[END_REF] and Uzsoy and Wang [START_REF] Uzsoy | Performance of decomposition procedures for job-shop scheduling problems with bottleneck machines[END_REF]. In Mason et al. [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF], a scheduling problem in semiconductor manufacturing close to the one tackled in this paper is solved by a modified SB heuristic. We propose to model solutions of the considered batching and scheduling problem through a variant of the disjunctive graph proposed in Mason et al. [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF].

Problem description

The diffusion area defines a batching and scheduling problem of wafer lots on two types of equipment: cleaning and furnaces machines. These resources are able to perform several lots simultaneously. Each lot requires one or more consecutive operations on the equipment and each operation has a recipe (specification on a process on how it should be executed on a tool; This pertains to requirements of maintaining proper temperature, pressure, and metal composition, among others) which determines its duration and the set of machines that are able to process it. On a 24-hour basis, each operation has to be assigned to an equipment, and included into a batch, i.e. a set of operations of the same recipe that are processed simultaneously by the equipment. The constraints in the diffusion area are divided into three types: equipment constraints, process constraints and line management constraints.

Constraints

Some of these constraints are common to the two types of resources while others are dedicated to furnaces.

Equipment constraints

• Common constraints

Dedicated equipment: Any equipment is able to process a limited set of recipes. Maximum batch size: Any equipment defines a maximum batch size corresponding to the capacity of that equipment.

Loading and unloading times: A time may be needed to load and unload a batch on the equipment.

Unavailability periods: The equipment may be unavailable during some periods (defined by time windows) due to qualification, repair, maintenance etc.

In process jobs: The equipment may be occupied at the beginning of the time horizon by in-process operations that have to be completed before the equipment becomes available.

• Specific furnace constraints

Minimum time between two batches on an furnace: Furnaces must be inspected after completing each batch.

Process constraints

• Common constraints Precedence: Operations must be performed following the manufacturing process of the lot. The operations of a lot are chained and no operation can start before the end of its predecessor, except for the first operation of each lot. Minimum time lag: There is a fixed handling and transport time between every two successive operations of a lot.

Release dates: They correspond to the arrival times of lots at the cleaning machines and furnaces. A lot cannot be scheduled before its release date. Because the diffusion area is a stage of the complex global production process, release dates are estimated by a simulation tool. Fixed recipe: Each operation of a lot is associated with a recipe, i.e. the lot should be processed on resources that are qualified for the corresponding recipe. This implies that all lots in the same batch must be processed with the same recipe. Process time: The process time of a batch on an equipment depends on the recipe. Maximum time lag: A time limit is given for two successive operations x and y of a lot. The difference between the starting time of y and the completion time of x cannot exceed this limit. Maximum time lags depend on the operations x and y. They can correspond to the maximum time limit between cleaning and furnace equipment.

Line management constraints

• Specific furnace constraints Minimum time between two batches of same recipe: There is a minimum time between the beginning of two batches of the same recipe on two different furnaces.

Objective

In semiconductor fabs, several indicators are used to measure the performance. The interested reader can refer to [START_REF] Montoya-Torres | Manufacturing performance evaluation in wafer semiconductor factories[END_REF] for more details. We describe thereafter the indicators that are relevant in our study.

• Number of moves (to maximize): It corresponds to the number of completed operations on the planning horizon, which can be compared to the target number fixed by the production managers.

• Batching coefficient (to maximize): Defined on the planning horizon, it is calculated as the number of moves divided by the sum of the number of batches performed on each machine, times the maximum capacity of that machine. Note that the denominator is the number of lots that could be performed if the equipment was loaded up to their maximum capacities.

• X-factor (to minimize): This indicator is used to evaluate the waiting times of lots in the diffusion area in order to reduce the cycle times. For a given lot, this factor is calculated as the total staying time of the lot in the diffusion area over its processing time.

These three indicators will be used to evaluate solutions in the fab. Other indicators could be relevant. Among them, the Work-In-Process (WIP), is defined as the number of wafers being in the fab at a given period, either in a production state or in a non-production state (e.g. transport and waiting). The law of Little [START_REF] Little | A proof for the queuing formula l = λw[END_REF] establishes that if a system is stable and stationary then the average WIP is proportional to the average cycle time. The throughput is defined as the outgoing number of wafers of the fab per unit of time. The evolution of the throughput in time makes it possible to know if the system is stable i.e. to know if there is no accumulation of lots in the fab or if the system evolves according to forecasts. Glassey and Resende [START_REF] Glassey | Closed-loop job release control for vlsi circuit manufacturing[END_REF], observed that there is a relation between the increasing of the throughput and the output of a fab. The cycle time drastically increases when throughput is close to the maximal capacity of the fab. Hence for both WIP and throughput indicators, the X-factor indicator is adequate.

The goal is to optimize the various performance measures, while taking into account the numerous complex constraints. The next subsection described the mathematical formulation of the problem.

Problem formulation

The scheduling problem can be formulated as follows. For sake of clarity, the above-described Unavailability periods, In-process jobs and Minimum distance between batches of the same recipe constraints are not included. In Section 4, we describe how we tackle these characteristics.

A set of jobs (lots) J = {J i |i = 1, . . . , n} has to be processed on a horizon T by a set of machines M = {M k |k = 1, . . . , m}. Each job J i is made of n i operations such that each operation O ij has a duration p ij > 0 and a set M ij ⊆ M of machines (the furnaces or the cleaning machines) able to process it. Let O k = {O ij ∈ O|M k ∈ M ij } denote the set of operations that can be assigned to machine M k . The value of p ij and the elements of the set M ij are determined by the recipe of operation O ij denoted ρ ij . In general we have M ij ⊂ M since each machine cannot be configured for all recipes. Each operation O ij has to be included in a batch on a resource k ∈ M ij . Each machine has a finite capacity R k which gives the maximal number of lots in the same batch. On each machine k, S k denote the setup time needed before starting a new batch, D k denote the removal time needed after the completion of a batch and s k denote the constant setup time needed between two different batches. s 0 k denotes the initial setup time on machine k, depending on the state of the resource at time 0. Two consecutive operations O ij and O i(j+1) of the same job are linked by minimal and maximal time lags. Once the batch of O ij is completed and removed from k, the setup for the batch of O i(j+1) cannot start before a minimal time lag τ min ij and has to start before a maximal time lag τ max ij . Let O = {O ij |i = 1, . . . , n; j = 1, . . . , n i } denote the set of all operations. Each job J i has a relative priority c i (c i < c j means that J i is more urgent than J j ). Each job corresponds to a number w i of wafers produced when the job is completed.

Table 1 summarizes the notations used.

Notations Description J = {Ji|i = 1, . . . , n} Set of jobs (lots) T Horizon length M = {M k |k = 1, . . . , m} Set of machines ni Number of operations of lot Ji O = {Oij |i = 1, . . . n; j = 1, . . . ni} Set of operations Oij Operation j of lot Ji pij Duration of operation Oij Mij Qualified machines to process Oij O k
Operations that can be processed on the machine

M k ρij Recipe of Oij R k Capacity of machine M k (batch max. size) S k
Setup before the operation on machine Decisions D1-D3 can all be represented by a family of batches B = {B kq } k∈{1,...,m},q∈{1,...,ν k } where B kq is the batch sequenced at position q on machine M k . ν k ∈ {0, . . . , |O k |} denote the number of batches assigned to machine M k . Decision D4 lead to a family of start times T = {t ij } Oij ∈O assigned to the operations. Once a solution {B, T } is determined, we have a machine assignment {m ij } Oij ∈O where m ij denote the machine O ij is assigned to, i.e. verifying that ∃q ∈ {1, . . . , ν k } such that O ij ∈ B mij q . To be feasible, a solution {B, T } and its corresponding assignment {m ij } Oij ∈O have to satisfy the following constraints. The operations of the same batch must have the same recipe, i.e.

ρ ij = ρ xy ∀B ∈ B, ∀O ij , O xy ∈ B (1)
Each operation must be assigned to a machine able to process its recipe.

m ij ∈ M ij ∀O ij ∈ O (2)
The batch capacity cannot be exceeded and each batch includes at least one operation.

1 ≤ |B kq | ≤ R k ∀B qk ∈ B (3)
An operation appears in only one batch, i.e.

B ∩ B = ∅ B, B ∈ B, B = B ( 4
)
All operations are included in a batch

∪ B∈B B = O (5)
The start time of the first operation of each lot cannot exceed the lot release date:

t i1 ≥ r i ∀J i ∈ J (6)
Each operation O ij , j > 1, cannot start before a minimal time lag after the end of its preceding operation O i(j-1) , which takes into account of the removal time of the batch of O i(j-1) , the minimal time lag τ min i(j-1)

and the setup time of batch of O ij :

t ij -t i(j-1) ≥ D m i(j-1) + p i(j-1) + τ min i(j-1) + S mij ∀i ∈ {1, . . . , n}, ∀j ∈ {2, . . . , n i } (7)
Each operation O ij , j > 1, has to start before a maximal time lag after the end of its preceding operation O i(j-1) , which takes account of the removal time of the batch of O i(j-1) , the maximal time lag τ max i(j-1) and the setup time of batch of O ij :

t ij -t i(j-1) ≤ D m i(j-1) + p i(j-1) + τ max i(j-1) + S mij ∀i ∈ {1, . . . , n}, ∀j ∈ {2, . . . , n i } (8)
The start times of two operations of the same batch must be equal:

t ij = t xy ∀B ∈ B, ∀O ij , O xy ∈ B (9)
An operation of a batch which is not at the first position on its machine cannot start before the end of the preceding batch on the machine, plus the necessary removal time of the preceding batch, plus the minimal setup time on the machine between two batches, plus the necessary setup time for the next batch.

t ij -t xy ≥ p xy + D k + s k + S k ∀B kq,q>1 ∈ B, ∀O ij ∈ B kq , ∀O xy ∈ B k(q-1) (10)
An operation of a batch in the first position on its machine cannot start before the initial setup time for this batch (we assume

S 0 k ≤ S k + D k + s k , ∀M k ∈ M): t ij ≥ s 0 mij ∀O ij ∈ O (11)
By definition, a feasible solution includes each operation inside a batch. In our problem, the scheduling horizon is limited to a scheduling horizon T . Hence only those batches released in the interval [0, T ] must be taken into account. Several criteria are used to measure the quality of a feasible solution. The number of moves is the number of wafers produced in [0, T ]:

f mov = θij ∈θ w i θ ij ( 12
)
where θ ij = 1 when the job is completed before T and θ i ∈ [0, 1[ denotes the completion ratio of job i before time T , otherwise.

θ ij = Oij ,tij +pij <=T p ij + Oij ,tij <T <tij +pij (T -t ij ) ni j=1 p ij
The batching coefficient is the average ratio of the actual size of each batch divided by its maximal size. Let B T = {B ∈ B|t ij < T, ∀O ij ∈ B} denote the set of batches started before time period T .

Batching coefficient = k B kq ∈B T |B kq |/R k |B T | (13) f batch = k B kq ∈B T |B kq |/(R k + Z k 100 ) |B T | ( 14
)
where Z k means the total number of qualified recipes on machine M k . The weighting of the batching coefficient by the number of qualified recipes on the concerned machine makes it possible to use the least general-purpose machines initially. The average X-factor is the average of the X-factor of each job weighted by the job priority. All the jobs do not have the same priority. Certain jobs are more important than others (important customers, tests to be carried out quickly, delay to catch up with, etc). Thus, a weight is assigned to the job to reflect its priority. Let

J T = {J i ∈ J |t ini + p ini ≤ T } denote the set of jobs completed before T . X-fac = Ji∈J T (t ini -r i )/p ini |J T | (15) f X-fac = Ji∈J T c i (t ini -r i ) |J T | (16)
The choice made together with the decision makers of the production unit is to combine these different objectives into a single one by maximizing the following weighted sum:

f = αf mov + βf batch -γf X-fac (17)
where α, β and γ are adjustable weights allocated to each objective function. The objective functions have been designed to integrate some requests of managers. In the calculation of f X-f ac , the delay of each lot is multiplied by its priority, in order to accelerate the urgent lots. In the calculation of the total batching coefficient, the batching coefficient of each equipment is multiplied according to the number of qualified recipes on this equipment. This leads to choose in priority the machines that are able to process less recipes, and thus to maintain the availability of the most flexible equipment. The objective of the problem is to search for a feasible selection {B, T } such that f is maximized. Note that, given a (feasible) solution {B, T }, f can be computed in O(N ) time where N = n i=1 n i is the total number of operations.

The disjunctive graph representation

The considered problem can be seen as an extension of the job-shop problem and, consequently, the disjunctive graph model can be used for batching and scheduling problems as proposed by Mason et al. [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF]. In their article, the authors consider a different objective function (total weighted tardiness). Sequence-dependent setup times and reentrant flows are also considered but there are no maximum time lags. Unfortunately, these constraints considerably increase the difficulty of the problem Gentner et al. [START_REF] Gentner | Batch production scheduling in the process industries[END_REF]. Let us explain how the problem is modeled using disjunctive graphs. We define the disjunctive graph G = (V, C, E) as follows.

• V is a set of nodes where there is one node per operation, denoted V ij , plus a dummy start node denoted 0.

• C is the set of conjunctive arcs representing the release dates and minimal and maximal time lags.

There is an arc from node 0 to node V i1 of each job J i . There is an arc from V ij to V i(j+1) and an arc from V i(j+1) to V ij , for each consecutive operations O ij and O i(j+1) of each job J i .

• Let B denote a partial or complete batching for the problem satisfying at least Constraints (1) through (4). B is a complete batching if Constraint (5) is also verified, otherwise it is a partial batching. Recall that B also defines the machine assignment m ij , for all

O ij ∈ ∪ B∈B B. We assume m ij = 0 if O ij is not batched in B, i.e. if O ij ∈ ∪ B∈B B.
B unambiguously defines a selection S as follows. For each distinct operations O ij and O xy such that • Each arc from 0 to the first operation O i1 is valued by max(r i , S 0 mij ), the maximal value between the release date of job i and the initial setup time for machine m i1 .

m ij = m xy = k = 0, select arc (V xy , V ij ) k if O xy is in a batch sequenced before the batch of O ij , select arc (V ij , V xy ) k if O ij is in
• Each arc from V ij and V i(j+1) is valued by

D mij + p ij + τ min ij + S m i(j+1)
, the value of the minimal time lag between O ij and O i(j+1) plus the setup and removal times linked to the assignment of the operations.

• Each arc from O i(j+1) and O ij is valued by -(D mij + p ij + τ max ij + S m i(j+1)
), the (negative) value of the maximal time lag between O ij and O i(j+1) plus the necessary setup and removal times.

• Each arc (V ij , V xv ) k is valued either by p ij + D k + s k + S k if the opposite arc is not selected or by 0 if the opposite arc is selected. Indeed, in the first case this arc represents the decision to sequence O xy after O ij on machine k whereas in the second case, both arcs (V ij , V xv ) k and (V xv , V ij ) k represent the synchronization of the operations included in the same batch.

We can state that the (partial) solution represented by the (partial) batching B and its selection S is feasible if and only if longest path problems from node 0 to each node V ij in G(S) have a solution. In the positive case and if B is complete, a feasible schedule T can be obtained by setting t ij to the length of the longest path from 0 to V ij . Furthermore T is the best schedule compatible with B one can obtain when the objective is to maximize f . The problem can be formulated as follows: Find the batching B verifying Constraints (1) through (5) such that the corresponding selection S is feasible and maximizes f .

As stated in the previous sections, the actual data issued from the fab have other characteristics that have been tackled, such as in-process jobs and machine down times. We can model easily these characteristics thanks to the disjunctive graph formulation. They can be both represented by operations with fixed start times on the involved machine. A start time t can be fixed by linking the node with node 0 by two opposite arcs valued by t and -t. As stated in Section 3, the actual problem also includes an important line management constraint: a minimum time between the start time of two batches of the same recipe scheduled on two different batches has to be respected. This can also be tackled through the disjunctive graph representation. A fictitious machine can be associated to each recipe and disjunctive arcs linking two operations of the same recipe can be defined. Then, whenever the two operations are batched on two different machines (furnaces), the disjunctive arc has to be oriented. Once oriented the arc is valued by the minimum required distance.

Solution Methods

We propose a two-phase heuristic method and a metaheuristic to solve the problem. The first phase of the heuristic is a constructive heuristic based on successive job insertions. The second phase is a local search method which aims at improving the initial solution. Both phases are based on the evaluation of a (complete or partial) selection through longest path calculations. For the metaheuristic, we propose a Simulated Annealing algorithm.

Evaluation of a partial or complete batching

Any partial or complete batching B and its selection S can be evaluated through the calculation of start time t ij , equal to the longest path from 0 to V ij in G(S) of each operation O ij . To compute such longest paths, since the graph includes arcs with negative weights, we use the Bellman-Ford algorithm which has a O(N |S ∪ C|) time complexity. If the algorithm finds a path of positive length, then the partial or complete solution is unfeasible. Otherwise the algorithm returns start times t ij and the objective function value f can be determined.

Computing an initial solution by a priority rule-based constructive heuristic

The initial solution (selection) is computed by a job insertion method. The jobs are first sorted in a list L according to the order: jobs involving maximal time lags first, then increasing release dates, then job priority.

The method starts with an empty batching. Then, the jobs are taken in the order given by the list and inserted one by one in the current batching. The insertion of J i is made as follows. Let B = {B kq } k∈M,q∈1,...,ν k denote the current batching including jobs located before J i in L. For each operation O ij of J i and for each resource k ∈ M ij , there are 2ν k + 1 insertion positions of O ij in the batch sequence of machine M k : indeed, for each batch B kq , we may insert O ij inside batch B kq or create a new batch at any position. Each of these insertion positions is evaluated with the algorithm described in the previous section and the one that maximizes f is kept to update B. If none of the insertion positions is feasible for an operation O ij , then all these partial solutions violate a maximal time lag and there exists insertion positions that violate only maximal time lag between O i(j-1) and O ij (the last position on each resource of M ij , for instance). Hence we have to remove O i(j-1) from B, the previous operation of job J i , and insert it at a later insertion position. If this is not feasible, they are not scheduled and delayed from the list of operations. Note that at most ni j=1 k∈Mij 2ν k + 1 insertion positions are tested per inserted lot. Thereafter, we will call by Priority-rule Based Insertion Algorithm (PBIA), the algorithm described above.

Improving the initial solution by iterative sampling

Iterative random sampling consists in iteratively applying the constructive algorithm presented in Section 5.2 and randomly changing each time the input operation list. Two different iterative sampling methods are proposed. In the first random sampling method the input list is initialized at each iteration as a random permutation. We call this algorithm Random Iterative Algorithm by Priority-rule Based Insertion (RIA-PBIA). In the second random sampling method, the list is changed by making few random modifications on the sorted list. It consists in dividing the list into intervals and, for each interval, in making a random modification following a uniform law. We call this algorithm Pseudo-Random Iterative Algorithm by Priority-rule Based Insertion (PRIA-PBIA).

Simulated Annealing

Simulated annealing (SA) belongs to the class of randomized local search algorithms and has been developed by [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] to handle hard combinatorial problems. The use of simulated annealing supposes the definition of neighborhood. The neighborhood we consider is defined by the following moves: 1. "Batch move". A batch is moved randomly (both equipment and position can change).

"Operation move". An operation is moved randomly in an existing batch or a new batch is created

(both equipment and position can change).

3. "Operation switch". Two batches with the same recipe are randomly selected. Then two randomly selected operations are switched.

In our implementation, 50% of the considered moves are batch moves, 25% are operation moves and 25% are operation switches. According to randomly generated moves, we proceed as follows. All random selections are made according to the uniform law. If a move is impossible due to a violation of a constraint, we restore the previous solution and try another move. Each neighbor solution is evaluated by means of the longest path computations presented in Section 5.1.

The algorithm starts with an initial solution s0 and then tries to find better solutions by searching in its neighborhood (obtained by the moves described above), and applying a stochastic acceptance criterion. When a neighbor (a solution s) of the current solution is selected we compare the difference ∆f = f (s) -f (s0). If the difference ∆f is negative, the neighbor replaces the current solution. Otherwise, the neighbor is accepted with a probability based upon the Boltzmann distribution P accept (∆f ) ≈ exp( -∆f kT ) where k is a constant and the temperature T is a control parameter. This temperature is gradually lowered following a geometrical law function g(T ) : g(T ) = αT with α < 1.

Computational experiments

We coded the proposed algorithms in Java and we tested them on a 1 GigaOctet RAM and 3.4 gigahertz processor computer. We conducted two types of experiment tests: the first type of tests is conducted on actual fab data and compared to the constructive initial solution obtained by PBIA solutions, with pseudo and random lists and the Simulated Annealing (SA) solutions. In the second type of tests, we compared the proposed heuristic PBIA to an algorithm for scheduling jobs on parallel batch processors with incompatible job families and unequal ready times problem proposed by Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF] as it is close to the one tackled in this paper.

Experimental tests on actual fab data

These tests have been performed on actual instances issued for two months of production. There are 700 jobs yielding a total number of 1400 operations with about 50 different recipes to schedule on 70 furnaces and 12 cleaning machines. Each furnace has a capacity of 4 or 6 lots and each cleaning machine of 2 or 4 lots. The time needed to load and unload a batch varies between 10 and 30 minutes. The minimum and maximum time lags varie from 10 minutes to 4 hours. The target time horizon is 24 hours. The weights of the components of the objective function have been computed both to normalize the different indicators and to take account of the user preferences. For the experiments, we selected α = 601, β = 1500001, γ = 41.

The number of replications is 5 for the simulated annealing and 20 for the pseudo and random lists. The computational time of the simulated annealing for one instance does not exceed 5 minutes and the constructive heuristic and the random and pseudo lists is limited to 1 minute. The percentage is calculated on the average relative improvement brought on the considered method over the reference solution obtained by the constructive heuristic (initial), i.e. M ethod solution-Initial list

Initial solution

From Table 2 and Table 3, we displayed the objective function solution values (note), the number of moves (f mov ), the batching coefficient (f batch ) and the X-factor (f X-f ac ) obtained from the three methods: Simulated Annealing (SA), Pseudo (PRIA-PBIA) and Random (RIA-PBIA) algorithms on the basis of solutions obtained by PBIA (Initial). It is important to note that the Simulated Annealing starts from the best solution of (PRIA-PBIA) and (RIA-PBIA).

We can notice that on average we improve solutions with the simulated annealing algorithm on the two considered months. An improvement of 47,05% on the number of moves f mov is noticed (resp. 6,6% for the second month), the batching coefficient f batch of 33,04% (resp. 4,15% for the second month), the f X-f ac increases by 41% (resp. is improved of 4,07% for the second month) and the objective value (note) is improved of 11,13% for the first month (resp. 12,61% for the second month). We improve on average initial solutions with PRIA-PBIA (Pseudo-Random) of 38,32% on the number of moves (resp. slightly deteriorated of 0,95% on second month), the batching coefficient of 28,23% (resp. 0,70% on the second month) and X-factor increases of 47,07% (resp. improved of 0,26% on the second month) and the objective function slightly decreases of 1,71% (resp. 1,42% on the second month). For the Random solutions (RIA-PBIA), the number of moves is improved of 37,26% (resp. decreased of 2,05% on the second month), the batching coefficient is improved of 29,96% (resp. 1,17% on the second month), the X-factor is deteriorated of 48,93% (resp. 1,16% on the second month) and the note slightly decreased of 3,86% (resp. 6,30% on the second month). In the Pseudo and random cases, we can notice that some indicators are sometimes slightly deteriorated (more in the case of Random lists) due to the fact we have a larger search space when using random lists than pseudo-random lists and then there is a capability of reaching better solutions but also worse see Table 4. (PBIA) has been implemented in this area and is actually running. We summarize some results obtained after the algorithm has been implemented. The comparison is done before and after using our algorithm. For confidential reasons, we are not allowed to describe the real values of the fab. This is why, with the agreement of the staff, we give results on equipment Type 1 and Type 2. They are respectively clean equipment and furnace equipment. Two indicators are presented: the X-factor and the batching coefficient. For the first month, the X-factor on the Type 1 equipment (resp. the Type 2 equipment) is improved of 26% (resp. 17% for the Type 2 equipment). The batching coefficient is also improved by 20% for the Oxidation and 29% for the Deposition equipment. For the second month, the X-factor is improved by 36% for the Type 1 and 23% for the Type 2 equipment. On the other hand, the batching coefficient is deteriorated by 6% for the Oxidation equipment and an improvement of 2% is noticed on the Deposition equipment. Substantial productivity improvements have been achieved in the diffusion area of the plant thanks to Priority-based Insertion algorithm (PBIA).

6.2 Experimental tests on data issued from Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF] We want to compare our algorithm (PBIA) to an algorithm proposed in the literature for a close problem. We decided to compare it to a scheduling of jobs on parallel machines with incompatible job families and unequal ready times Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF].

We show through experimental tests that the constructive heuristic PBIA we designed for the specific problem of batching and scheduling in a diffusion area can be easily adapted to solve the problem described by Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF]. Furthermore, our proposed heuristic outperforms the heuristic described in their article on some instances for the problem of minimizing Total Weighted Tardiness (TWT).

The aim of the problem tackled in Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF] consists in scheduling of parallel batch processing machines with incompatible job families and unequal ready times of jobs to minimize TWT include:

• Jobs of the same family have the same processing time.

• All the batch processing machines are identical in nature.

• Once a machine is started, it cannot be interrupted, i.e. no preemption is allowed.

To solve this NP-hard problem the authors proposed two different decomposition approaches. The first approach constructs fixed batches, then assigns these batches to the machines using a genetic algorithm (GA), and at last, sequences batches on each machine. The second approach first assigns jobs to machine using a GA, then constructs the batches on each machine for the jobs assigned to it and at last, sequences these batches. In conclusion, among the tested algorithms, the authors proved on their tests that the algorithm GA 2 BATC-II which is in the second type of approach provides in average better results. As the authors used a genetic algorithm, the After taking the 162-instances (randomly generated) from the authors, we have conducted experiment tests on the same computer as one described in the subsection 6.1. Let us recall that in this subsection the objective function is the TWT and the goal is to compare the TWT obtained from our proposed heuristic (PBIA) to algorithm GA 2 BATC-II proposed in Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF]. As it is known, the input of our proposed algorithm is a list. Then, to give a list to PBIA, we retain the best one list between two considered lists as follows. The first list is ordered by increasing order of due date and the second one is ordered by decreasing order of the priority of jobs. As the evaluation of a list does not require much time (less than 3 minutes), for each of the two lists we calculate the TWT and we retained the minimum one. With the solution obtained by the minimum one, we apply local search improvement (see Subsection 5.3). We have reported in Table 6 the assessment of this comparison.

From Table 6, we notice that on average as well as for the case of m = 3, m = 4 and m = 5, our heuristic (PBIA)is better in terms of objective function and computational times. In the considered article, the

M k D k Setup after the operation on machine M k s k Delay inter batch s 0 k

 0 Initial setup time τ min ij Min. delay between Oij and Oij+1 τ max ij Max delay between Oij and Oij+1 ci Priority of lot Ji wi Number of wafers of Ji B = {B kq } k∈{1,...,m},q∈{1,...,ν k } Batch at position q on machine M k ν k Number of batches on machine M k T = {tij } O ij ∈O Set of start times mij Machine which processes Oij θi Completion ratio of Oij B T Batches started before T Z k Number of qualified recipes on machine M k θ T Operations started before T N = n i=1 ni Total number of operations

Finding a feasible solution

  for the problem lies in making four types of decisions: D1 -Partition the operations into batches, D2 -Select a resource to process each batch, D3 -Order the batches on each resource, D4 -And assign a start time to each batch.

  E is the set of disjunctive arcs which represent the decisions of the problem. There are two opposite conjunctive arcs (V ij , V xv ) k and (V xy , V ij ) k for any machine k ∈ M and for any pair of operations O ij , O xy ∈ O k , O ij = O xy . These arcs represent the sequencing or the batching decision concerning O ij and O xy on machine k.

  a batch sequenced before the batch of O xy , and select both arcs (V ij , V xv ) k and (V xy , V ij ) k if O ij and O xy are assigned to the same batch. Once a selection is computed, we define a graph G(S) = (V, C ∪ S) where arcs C ∪ S are valued as follows (we assume s 0 0 = s 0 = S 0 = D 0 = 0):

  Table 2: Results on instances issued of the first month of the fab Table 3: Results on instances issued of the second month of the fab One objective of this study was to develop and propose an efficient scheduler algorithm to support fabrication operators and improve the principal indicators in the diffusion area. Our proposed algorithm

Table 1 :

 1 Summary of notations used to formalize the problem

Table 4 :

 4 Standart deviation values on the first and second month

Table 5 :

 5 Table5summarized the different parameters used in their algorithm (see Mönch et al.[START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF]). Parameter settings GA for approach 2 Mönch et al.[START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF]

	Parameter	Value
	Population	200
	Crossover probability	0.80
	Mutation probability	0.01
	Replacement probability	0.60
	Diversity	0.03
	Number of generations	1000

Acknowledgments

This work was part of the MEDEA+ European project HYMNE (High Yield driven MaNufacturing Excellence in sub 65 nm CMOS), partly funded by the "Ministère de l' Économie, de l'Industrie et de l'Emploi" (French Ministry of Economy, Industry and Employment).

We have used a heuristic which was not originally designed for the problem of minimization of TWT on parallel batch machines with incompatible job families and unequal ready times of jobs. The results of the tests showed that our heuristic outperformed in the majority of the cases the GA 2 BATC-II: we had better objective function values and less computation times.

Concluding remarks

We proposed a model and a method based on a disjunctive graph for batching and scheduling problem in a semiconductor manufacturing factory while taking into account complex constraints and optimizing multiple measures. A constructive algorithm has been proposed to solve the problem, local search improvements based on the disjunctive graph representation have been defined and a simulated annealing algorithm has been developed. The computational tests made on real instances of the factory plant showed that we obtain quality solutions in fast calculation times. For the industry, very substantial productivity improvements have been achieved in the diffusion area and also on the overall fab thanks to the use of the proposed algorithm named Priority rule-Based constructive algorithm (PBIA). With adjustments, the Priority-rule Based Insertion Algorithm (PBIA) has been embedded in a software used by the industry. The use of a disjunctive graph brings significant improvements for interactive scheduling at the fab level. A prototype software gives the graphical user interface of the software. It includes the off-line batching and scheduling phase and also an interactive module that allows the decision-makers to make modifications and to test options directly on the proposed plan.

Regarding, the simulated annealing algorithm, the procedures significantly improve the different criteria (number of moves, batching coefficient and X-factor) in comparison with a local search method. We also compared our algorithm to an algorithm designed for a problem of scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times Mönch et al. [START_REF] Mönch | Heuristic scheduling of jobs on parallel batch machines with incompatible job families and unequal ready times[END_REF] as it is close to the one tackled in this article. The computational tests showed that for more than half of the instances our heuristic called Priority-rule Based Insertion (PBIA) outperforms the algorithm proposed by Mönch et al. [MBFP05]. This research can be extended in many ways. Other types of moves such as the simultaneous move of two jobs linked by a maximum time lag should be tested. The maximal time lags are currently hard constraints. However in practice some of them can be relaxed and treated as soft constraints or objectives. Another important issue can be in a multicriteria analysis study. Computing several Pareto optimal solutions are useful particularly when the situation frequently changes as it is usually the case of semiconductor manufacturing.


Aldo Maalouf
email: maalouf@sic.sp2mi.univ-poitiers.fr

Mohamed-Chaker Larabi
email: larabi@sic.sp2mi.univ-poitiers.fr

A NO REFERENCE OBJECTIVE COLOR IMAGE SHARPNESS METRIC

. The multiscale structure tensor allows for accumulating multiscale gradient information of local regions of the color image. Thus, averaging properties are maintained while preserving edge structure. This structure tensor is capable of identifying edges in spite of the presence of noise. Once edges are identified, we define a sharpness metric based on the eigenvalues of the multiscale structure tensor. Particularly, we show that the difference of the eigenvalues of the multiscale structure tensor can be used to measure the sharpness of color edges. Based on this fact we formulate our no reference sharpness metric for color images. Experiments performed on LIVE database indicate that the objective scores obtained by the proposed metric agree well with the subjective assessment scores.

INTRODUCTION

Recently, there has been an increasing need to develop quality measurement techniques that can predict perceived image/video quality automatically. These methods are useful in various image/video processing applications, such as compression, communication, printing, display, analysis, registration, restoration, and enhancement. For example, a noise metric can be used to estimate the quantization error caused by compression without accessing the original pictures, while a sharpness metric can be used as a control parameter for sharpness enhancement algorithms applied to digital imagery; a sharpness metric can also be used to estimate the blur caused by image compression algorithms. Subjectively, to have a look is probably the best way to evaluate image quality, because human beings, e.g., end users, should make ultimate assessment on the performance of algorithms on digital images. However, subjective methods normally take/cost much time/resources as end users have to be highly involved, i.e., these methods cannot be routinely performed as there could even be a difference in assessment between different (groups of) end users. So, effective and efficient subjective IQA metrics are desirable but too hard to develop in real time systems. Therefore, objective IQA is more demanding. Depending on whether or not the original image is used or on how much information from the original image is used, objective IQA can be classified into three types: full-reference This work was supported by the French National Research Agency (ANR) under the QuIAVU project.

(FR), no-reference (NR), and reduced-reference (RR) metrics. FR metrics need full information of the original images and demand ideal images as references which can be hardly achieved in practice for some applications (such as broadcasting...). The traditional methods of FR (such as peak signal-to-noise-ratio PSNR) are based on pixel-wise error and have not always been in agreement with perceived quality measurement. On the other hand, RR metrics make use of a part of the information from the original images in order to evaluate the visual perception quality of the distorted ones. Particularly, the metric is not relative to a reference image, but rather an absolute value is computed based on some characteristics of the given image. NR metrics aim to evaluate distorted images without any cue from their original ones. Quality assessment without a reference is a challenging task; distinction between image features and impairments is often ambiguous. Of particular interest of this work is the noise resilient wavelet-based no reference (NR) objective sharpness metric for color and multispectral images. For that we used a wavelet-based multiscale multispectral tensor [START_REF] Maalouf | Regularization of Multivalued Images by Means of a Wavelet-Based Partial Differential Equation[END_REF] to detect the edges of multispectral images. This structure tensor is capable of detecting edges in spite of the presence of noise. Then, based on the eigenvalues of this structure tensor, we propose a NR sharpness metric. The rest of the paper is organized as follows: section 2 is a state-of-the-art of the previous work. In section 3 a review of the multiscale structure tensor is given. In section 4 we present our NR metric. The last two sections summarize the experimental results, the conclusions and directions for future work.

STATE OF THE ART

This section presents an overview of the NR sharpness metrics existing in the literature. The easiest way to assess the quality of an image from edge information is to use the variance [START_REF] Erasmus | An automatic focusing and astigmatism correction system for the sem and ctem[END_REF]. This method is based on the fact that, when an image is blurred, the transitions between the graylevels decrease. As a result the variance decreases. The problem of this method is that fine textures and edges in between small varying regions are regarded as blurred or smoothed edges. Another simple technique is to use the image histogram [START_REF] Firestone | Comparison of autofocus methods for automated microscopy[END_REF]. The histogram method consists of fixing a threshold which is usually the mean of the image, then, the metric is defined as a weighted sum of the histogram bins which are above the threshold. A sharp image contains a high number of bins above the threshold and thus, increases the metric. The problem of this metric is that the selection of a good threshold may not be obvious. Another histogram-based method has been proposed by Chern et al. in [START_REF] Chern | Practical issues in pixel-based autofocusing for machine vision[END_REF]. Chern et al. who defined their metric on the entropy computed from the histogram. The idea is that the entropy increases when the probability of occurrence of a gray level decreases and vice versa. Sharp images contain a larger number of gray levels, thus a lower probability and higher entropy. An autocorrelation metric was proposed by Batten in [START_REF] Batten | Autofocusing and astigmatism correction in the scanning electron microscope[END_REF]. Batten proposed to compute the difference between the autocorrelation values at two different distances along the horizontal and vertical directions. If the image is smoothed or blurred, the correlation between the neighboring pixels increases and the correlation increases. As a result, the sum of the difference metric decreases. Other NR metrics use differential methods such as the gradient and the laplacian [START_REF] Bovik | Handbook of image and video processing[END_REF] to assess the sharpness of an image. While these methods have a high accuracy, they are very sensitive to noise. Some other perceptual metrics such as the one proposed in [START_REF] Marziliano | Perceptual blur and ringing metrics: Applications to jpeg[END_REF] predict the sharpness of an image from the width of edges. In these metrics, an edge detection technique is first applied, then, the start and end positions of each edge are determined as well as the locations of the local maxima near the edges. Thereafter, sharpness measure is defined by the edge width. Another statistic-based sharpness measure is the well known kurtosis measure [START_REF] Zhang | A kurtosis-based statistitcal measure for two-dimensional processes and its application to image sharpness[END_REF] [START_REF] Caviedes | A new sharpness metric based on local kurtosis, edge and energy information[END_REF]. The kurtosis is a statistical measure of the peakness and flatness of a distribution and is inversely proportional to the sharpness. Recently, Ferzli et al. [START_REF] Ferzli | A no-reference objective image sharpness metric based on the notion of just noticeable blur (jnb)[END_REF] proposed a NR sharpness metric based on the notion of just noticeable blur (JNB). They showed that the human visual system will mask the blurriness around edges up to a certain threshold. They referred to this threshold as the JNB which is obtained by performing subjective tests. Then, they used a probability summation on the JNB to construct their metric. The above mentioned metrics do not deal with color images, i.e. they do not take into account the multispectral aspect of the geometry of a color image. Furthermore, they do not take into account the noise factor. In fact, the presence of noise may affect the performance of the sharpness metric. A solution to this problem could be to denoise the image before applying the metric. However, this may cause the degradation of edges and, therefore, the failure to predict the edge sharpness.

In this work, we address the problem of evaluating the sharpness of a color image with and without noise. For that, we dused a multiscale multisctructure [START_REF] Maalouf | Regularization of Multivalued Images by Means of a Wavelet-Based Partial Differential Equation[END_REF] tensor whose norm defines the edges of color images. Then, we define our sharpness metric based on the eigenvalues of this tensor.

REVIEW OF THE MULTISCALE WAVELET-BASED STRUCTURE TENSOR

The multiscale edge representation described in [START_REF] Mallat | Singularity detection and processing with wavelets[END_REF] and [START_REF] Mallat | Characterization of signals from multiscale edges[END_REF] is used to define the structure tensor. In this approach, x and y-directional wavelets are given by the partial derivatives of a separable, nonorthogonal scaling function θ (x, y) as follows: L 2 R 2 at scale j are defined by:

ψ 1 (x, y) , ψ 2 (x, y) = ∂ θ ∂ x (x, y) , ∂ θ ∂ y (x,
W 1 j (x, y)I W 2 j (x, y)I = I * ψ 1 j (x, y) I * ψ 2 j (x, y) = ∇ (I * θ j ) (1)
Where ψ l j (l = 1, 2) and θ j represent the wavelet and the scaling function at scale j, respectively, defined by: ψ l j (x, y) = ψ l x 2 j , y 2 j √ 2 j and θ j (x, y) = θ x 2 j , y 2 j √ 2 j . This stipulates that the wavelet transform of an image consists of the components of the gradient of the image, smoothed by the dilated smoothing function θ j . The direction of the gradient vector at a point (x 0 , y 0 ) indicates the direction along which the image I has the steepest slope. Therefore, a point (x 0 , y 0 ) is regarded as an edge point at scale j if the magnitude of the wavelet coefficient attains a local maximum along the gradient direction. Based on this theory of singularity detection, a multiscale multistructure diffusion tensor can be constructed for an mband image I (x): R 2 → R m with components for I i (x) : R 2 → R for i = 1, 2, 3, ....., m (m = 3 for color images) as follows:

G j w =    m ∑ i=1 W 1 n, j,i 2 m ∑ i=1 W 1 n, j,i W 2 n, j,i m ∑ i=1 W 1 n, j,i W 2 n, j,i m ∑ i=1 W 2 n, j,i 2    (2)
Where W k n, j,i is the undecimated wavelet coefficient computed at scale j and position n for the image channel i. The norm of G j w is defined in terms of its eigenvalues, G j w = λ +λ -, and it describes the total local derivative energy. Figure 1 shows the norms of the multistructure tensor of Di Zenzo defined in [START_REF] Zenzo | A note on the gradient of multi-images[END_REF] (Figure 1(b)) and the norm of the multiscale multistructure tensor defined in (2) (Figure 1(c)) of the noisy 'Lenna' image.

It is clear that the multiscale structure tensor provides a better characterization of the image edges. We formulate in the next section our NR sharpness metric for multispectral images based on the eigenvalues analysis of the multiscale structure tensor.

NO REFERENCE SHARPNESS METRIC

In this section, a noise-resilient sharpness metric for color images is presented. The proposed metric is based on the behavior of the eigenvalues of the wavelet-based structure tensor described in the previous section. Let λ j + and λ j -be the largest and the smallest eigenvalues of G j w computed at a wavelet scale j. λ j + is none other the derivative energy in the most prominent direction which is equal to the orientation in the image with maximum color change, while λ j -describes the amount of derivative energy perpendicular to the prominent local orientation. The difference λ j +λ j -describes the line energy, i.e., the derivative energy in the prominent orientation that is corrected for by the energy contributed by noise λ j -. As the scale j increases, λ j + increases at the edge pixels while λ j -decreases at noise pixels. Thus, the difference λ j +λ j can be used to predict image sharpness. Therefore, we define our sharpness metric as:

Score = J ∑ j=1 M ∑ m=1 N ∑ n=1 λ j + (m, n) -λ j -(m, n) (3)
where J is the largest wavelet scale (practically we set J = 3) and, M and N are the height and width of the image respectively.

EXPERIMENTAL RESULTS

A good sharpness metric is one which has a curve that declines slowly or does not decay abruptly when there is an increase in blurness. This ensures that the metric can perform robustly under a larger degree of blurness. To test the performance of our proposed metric we made use of the image set shown in figure 2.

The images of figure 2 are blurred with a Gaussian point spread function of different standard deviation σ values before separately introducing with two most common additive random noises, namely the zero-mean Gaussian white noise and salt-and-pepper noise. For each type of noise, two noise levels are introduced. For Gaussian white noise, its variance, ν is set equal to 0.010 and 0.020 while for salt-and-pepper noise, the noise density, d is set equal to 0.10 and 0.20, respectively. Examples of the noisy blurred images are shown in figure 3 for the Lenna image.

The scores obtained for each image of figure 2 in the presence of salt-and-pepper and Gaussian noises are shown in figure 4. It is clear that the proposed metric is quasi robust to different noise levels and types.

We also evaluate our wavelet-based sharpness metric using the LIVE database generated by [START_REF] Sabir | A statistical evaluation of recent full reference image quality assessment algorithms[END_REF], which presents the most recent and comprehensive survey of the performance of various image quality metrics available in the literature. Particularly, we used the dataset from the LIVE database selected by the authors of [START_REF] Ferzli | A no-reference objective image sharpness metric based on the notion of just noticeable blur (jnb)[END_REF] for evaluation. This dataset contains fifty images that were extracted from LIVE database. These images were subjectively evaluated for sharpness assessment. The sharpness mean opinion scores (MOS) as well as the source code for the objective metric proposed in [START_REF] Ferzli | A no-reference objective image sharpness metric based on the notion of just noticeable blur (jnb)[END_REF] are downloaded from [10]. Figure 5 shows the MOS vs the scores obtained by our metric as well as the scores obtained by the metric proposed in [START_REF] Ferzli | A no-reference objective image sharpness metric based on the notion of just noticeable blur (jnb)[END_REF]. Finally, we calculate the three measurements recommended by VQEG [START_REF] Vqeg | Final report from the video quality experts group on the validation of objective models of video quality assessment[END_REF] to test the consistency of our sharpness metric and the subjective perception, namely, the Pearson correlation coefficients (PCC), the Spearman coefficient correlation (SCC), the mean absolute prediction error (MAE) and the outlier ratio (OR). The results are shown in table 1 for our metric and the metric proposed in [START_REF] Ferzli | A no-reference objective image sharpness metric based on the notion of just noticeable blur (jnb)[END_REF]. From these results, we can see that the PCC and the SCC coefficients are high indicating that we can achieve a high accuracy with our metric.

CONCLUSION

This paper presented a wavelet-based noise-resilient sharpness metric for color images. First, a wavelet-based structure tensor is defined whose eigenvalues are used to define the metric. Experiment results showed that the proposed metric is resilient to noise and has a strong correlation with human judgment. The proposed metric can therefore eliminate the need to denoise the image before assessing the quality. For future work, we propose to conduct more extensive human experiment on high dynamic range HDR images and compare the consistency of our sharpness metric in evaluating the quality of this type of images.

Figure 1 :

 1 Figure 1: Norms of the Di Zenzo structure tensor and the multiscale structure tensor defined in (2)

Figure 2 :Figure 3 :Figure 5 :

 235 Figure 2: Overview of the test material. a-Lena, b-Iris, c-Lighthouse, d-Caster and e-Haifa.

Table 1 :

 1 Objective assessment scores

	Model	PCC	SCC	MAE	OR
	Ferzli et al. [7] 0.9477 0.9301 0.3428 0.3400
	Our Metric	0.960	0.951	0.3	0.201

18th European Signal Processing Conference (EUSIPCO-2010) Aalborg, Denmark, August 23-27, 2010 © EURASIP, 2010 ISSN 2076-1465


Y C Mombo Boussougou

T Le Gouguec

Y Quéré

F Huret

Study of Proximity-Coupling: Application to new RF-Interconnects

Keywords: radiofrequency, coupling, contactless interconnect, modulation, multiplexing, capacitive. I

The new concept of RF Interconnect presented in this paper is aimed at replacing conventional interconnections in digital circuits. The study was, thus, focused on the transmission channel composed of contactless transitions. After a brief recall of the theory, an RF interconnection was modeled through use of multiplexing digital CDMA and amplitude modulation and characterized by using a co-simulation tool. The signal forms observed along the transmission system and the frequency response of the RF channel are presented to highlight the capabilities of this RF interconnect.

INTRODUCTION

The performances of integrated circuits such as Systems on Chip (SoC), Systems in Package (SiP) or Systems on Package (SoP) are limited by the global interconnect layers. Indeed, the increase in the number of circuit features, despite the progress in terms of integration, leads to an increase of chip surface at the origin of a lengthening of global interconnects, which extends the delays in information transmission. To reduce these effects, new solutions for signal transmissions need to be developed. Among the solutions under consideration let us cite optical interconnects, 3D interconnects and radio frequency (RF) interconnects [START_REF]International Technology Roadmap for Semiconductors (ITRS)[END_REF][START_REF] Ho | The Future of Wires[END_REF][START_REF] Davis | Interconnect limits on gigascale integration (GSI) in the 21 st century[END_REF][START_REF] Havemann | High-Performance Interconnects: An Integration Overview[END_REF]. Optical interconnects allow important flows, but require a significant change in the chip-making technology, which implies a high implementation cost. By taking advantage of functional block superposition, 3D-interconnects allow reduction of line lengths. However, arranging such optimal system architectures may be very difficult or even impossible, and several long interconnects may still exist. There are two kinds of RF interconnects corresponding to the transport of information either by radio, or by a waveguide used as a transmission medium. But, the interest of using radio signals to distribute clock ones is decreased by the problems of electromagnetic compatibility (EMC) and the size of antennas.

Together these considerations led us to carry out an electromagnetic study of RF interconnects, i.e. microstrip or coplanar lines, compatible with the current microelectronic technologies.

Section II presents the concept of RF interconnects and its possible applications in the transmission of many binary signals. Then, Section III deals with the study and the modeling of different couplings aimed at realizing a contactless transmission to the microwave line. Finally, an example of digital/analog co-simulation for the characterization of RF-interconnections is given in Section IV. Our conclusions are drawn in Section V, and clues are given about the perspectives opened by this study.

II. RF-INTERCONNECT CONCEPT

RF interconnects can be used for the distribution of clock signals or to replace parallel data buses [START_REF] Ryujunwoo | RF Interconnect for Multi-Gbit/s Board-Level Clock Distribution[END_REF][START_REF] Mau-Chung | Advanced RF/Baseband Interconnect Schemes for Inter-and Intra-ULSI Communications[END_REF]. In this paper, Fig. 1 illustrates the principle of RF interconnect in the case where it replaces a classical 16-bit bus (Fig. 1). With such a parallel bus at 100 Mbits/s, information is transmitted at 16 x 100 Mbits/s =1.6 Gbits/s. But, the coupling between the signal lines makes difficult any increase of the speed. To be equivalent to the parallel 16-bit bus, the proposed serial interconnect (Fig. 2) must authorize a similar data rate (1.6 bits/s). By transmitting all bits simultaneously with a proper multiplexing technique, an RF-interconnect (Fig. 3) allows one to overcome the slowness of baseband transmission and to reduce the number of lines. To transmit signals to the RF-line many solutions are available, but here focus is on contactless transitions. Indeed, the low-frequency noise generated by digital circuits is filtered by using a contactless coupling between the digital circuits and the RF-line filters. Moreover, this is relevant in terms of achievement of reconfigurable circuits. The characteristics of a contactless coupling make that the digital signals have to be modulated for transmission via the microwave channel. In addition, because of the simultaneous transmission of the multiple signals, the use of multiplexing techniques such Frequency Division Multiple Access (FDMA) or Code Division Multiple Access (CDMA) is necessary. The next section will be devoted to the study of contactless coupling between two lines

RF Interconnect Interface Interface

III. CONTACTLESS TRANSITION BETWEEN TWO LINES

A. Capacitive Coupling

The replacement of metallic vias by a coupling contactless access to the RF interconnect would limit the baseband noise and enhance flexibility.

Let us study, at first, the case of a coupling through capacity as shown in Fig. 4 As a first approximation, the value of the coupling capacitance can be determined by the planar capacitance as follows:

e W W C RF R . 0 ε ε ≈ ( 1 )
where W RF is the width of the RF interconnect line, W is the width of the access line, e is the gap between the two lines, and ε R is the relative permittivity of the substrate between the two lines. More accurate expressions of the capacitance are available in [START_REF] Wong | An empirical three-dimensional crossover capacitance model for multilevel interconnect VLSI circuits[END_REF].

One of the simplest contactless couplings is a capacitive coupling between the digital and the RF line.

To determine the parameters that have a great influence on the transmission level, let us calculate the transfer function of a radiofrequency interconnect. The corresponding scheme is presented in Fig. 5; one should note that it is more complete than the one presented by Sun et al. [START_REF] Sun | Performance of inter-chip RF-Interconnect using CPW, capacitive coupler, and UWB Transceiver[END_REF] because of the addition of stubs with their loads to the sections of interconnects corresponding to ports. The stubs can be used to optimize both the transmission gain and the bandwidth, or to allow the simulation of another charge on the line.

Z c , d Z c Z c C t C r R 1 R 2 Z c1 ,L 1 Z c 2 ,L 2 Z c 4 ,L 4 Z c 3 ,L 3 R 4 R 3 V s V i V o A B I O Fig. 5

. Electric circuit for the transfer function calculation

For greater clarity the transfer function calculated here corresponds to the power gain of the structure depicted in Fig. 5 when the lengths, L 2 and L 4 , are both neglected. It represents the output-to-input power ratio expressed as follows:

I O P P H = ) (ω ( 2 ) 2 2 2 4 3 4 3 4 )) sinh( ) (cosh( ) // // ( ) // // ( ) ( d Z Z d Z Z Z Z Z R Z Z R Z R Z H in C C in C C in R C R P T R γ γ ω - × + + = (3) where
) tanh( ) tanh( d Z Z d Z Z Z Z L C C L C in γ γ + + = ( 4
) C C R L Z Z R Z Z R // ) // ( 4 3 + = (5) 1 // ) // ( R C C in P Z Z Z Z Z T + = (6)
Z R1 and Z R3 are, respectively, the equivalent impedances of stubs 1 and 3. This transfer function permits the assessment of the influence exerted by various parameters such as the coupling capacitance, C T or C R , the characteristic impedance of the RF line, Z C , and the access impedances, R 2 and R 4 . For example, Fig. 6 illustrates the evolution of the power gain for several values of coupling capacitance in the case where

C T = C R = C C .
As expected, it shows that the higher the capacitance value is, the higher the level of transmission is. So, the use of high coupling capacitances will necessitate a rather large surface area and a small gap between both levels.

At this stage, it is worth wondering about the impact of the characteristic impedance of the RF interconnect on transmission. It is illustrated by Fig. 7, where the energy transfer is maximized through use of high-impedance lines. However, in the case of a microstrip RF line of high impedance, its narrow width will lead to significant metal losses. It is, thus, paramount to find a compromise between the characteristic impedance values and metal losses. This study was conducted by considering the capacitive coupling between microstrip lines and/or coplanar waveguide. However, the availability of a large variety of microwave transmission lines can also allow the development of other forms of coupling.

The plots of variations of the transfer function versus the coupling capacitance show that high capacitance values are needed to enhance the transmission. But, according to technological considerations, this is uneasy. It is thus necessary to consider other solutions.

B. Microstrip line and slotline coupling

According to the investigations reported in [START_REF] Gupta | Microstrip Lines et Slotlines[END_REF] about the coupling between a microstrip line and a slotline, a broadband transmission is possible with high transmission levels. The principle of the transmission system is illustrated in Fig. 11: the signals are encoded by specific orthogonal codes, and then they are multiplexed and modulated on a highfrequency carrier because of the high pass channel. The information is retrieved through the reverse operation. Figure 12 illustrates various signal forms along the transmission system before and after the channel. With the co-simulation the digital and analog parts can coexist, and the information transfer is possible under certain constraints for the channel.

V. CONCLUSION

This study demonstrated the feasibility of the concept of RF interconnects through circuit and electromagnetic simulations. Then, we developed the expression of the transfer function so as to determine the parameters to be optimized. It highlighted not only the key-role of the capacitance value, but also the impact of the characteristic impedance value of the RF line.

At first, our investigations were focused on microstriplinemicrostripline transitions. Simulations of such a transition provided good values for the S-parameters. The use of this channel in a co-simulation tool enabled us to plot the signals before and after the channel.

There is no doubt that a more detailed study of the technological possibilities and constraints is worth being carried out. A more accurate assessment of such interconnects requires to gain more insight into the complexity of the electronic functions in emission and reception.

Fig. 1 .

 1 Fig. 1. Classical 16-bit bus interconnect

Fig. 2 Fig. 3 .

 23 Fig. 2. Serial transmission

Fig. 4 .

 4 Fig. 4. Coupling capacitance between two microstrip lines (a), equivalent circuit model (b)

Fig. 6 .Fig. 7 .

 67 Fig. 6. Frequency variations of the transfer function for different values of the coupling capacitance, C C

Fig. 8 .Fig. 9 .

 89 Fig. 8. 3D structure of a microstripline-slotline transition Figure 8 shows the 3D structure of a microstripline-slotline transition for a 2-port communication. Using the circuit model presented in Fig. 8, the impedances and stubs lengths structure can be optimized. Simulations run with optimized parameters, i.e. L MT = 552 um, L ST = 522 um, Z MT = 174 Ohms and Z ST = 411 Ohms, led to the optimized S-parameters curves presented in Fig. 9.

Fig. 10

 10 Fig. 10. S-parameters of a microstripline-slotline transitionOne should note that the device behaves like a high-pass function. Its frequency response is relatively flat and the losses induced by the double transition are low over a band from 30 to 70 GHz.

Fig. 11 .

 11 Fig. 11. Various signal forms along the transmission system

Fig. 11 .

 11 Fig. 11. Diagram of a CDMA transmission system with ASK modulation