Emna Mhiri
email: emna.mhiri@grenoble-inp.fr

Fabien Mangione

Mireille Jacomino

Philippe Vialletelle

Guillaume Lepelletier

Heuristic algorithm for a WIP projection problem at finite capacity in semiconductor manufacturing

Keywords: WIP projection, finite capacity planning, semiconductor industry, mixed integer programming, iterative algorithm

In this paper, we propose a heuristic approach for fixing work-in-progress (WIP) projection issues in the semiconductor industry especially for High Mix Low Volume (HMLV) facilities. The considered problem consists of estimating the start and end dates for each remaining process step of the production lots in the WIP and anticipating the fab loading taking into account the constraints of the maximum throughput of machines considered as capacity constraints and customer delivery commitments. The objective being to guarantee on-time delivery, we focus on minimizing the total weighted tardiness (TWT). We have formulated the problem into a mixed-integer programming (MIP) and we have empirically shown its computational intractability. Due to the computational intractability using actual production data, a heuristic algorithm is proposed. It is an iterative finite capacity planning system that considers as inputs lots due dates and equipment capabilities and capacities. The performance of the heuristic is assessed using industrial instances. It turns out that it achieves targeted objectives with satisfactory results in terms of quality of the solution and computation time.

I. INTRODUCTION

W ORK-in-progress (WIP) projection is a mid-term ca- pacity planning activity. The objective is to compute a mid-term target schedule in order to drive factory execution, to anticipate production issues and to calculate net demand and net resource capacities. In our study, the outcome is a weekly-released schedule that depicts the start and end dates of each remaining processing step as well as the expected workload accumulated on each equipment per time bucket over the planning horizon.

In this study, the WIP projection problem is considered in one of the most dynamic industries in the world, the semiconductor industry. The semiconductor manufacturing process is extremely complex and constantly innovating. The considered wafer production plant is a High Mix Low Volume (HMLV) production line: there are several hundreds of products, different technologies and heterogeneous toolsets i.e. collections of nonidentical multi-purpose parallel machines (or tools). Moreover, typical semiconductor fabrication processes require several hundreds of different steps. As, for obvious reasons, HMLV fabs cannot multiply machines, their production flows are re-entrant: the same machine can process products at different stages of their fabrication. This also means that according to the decisions taken on the production line, products may experience various cycle times depending on the priority given either to a given product, to satisfy customer demand, or to a certain technological level for the purpose of line balancing. Hence, the capacity planning issue is difficult to solve and it is particularly more complex than in other industries [START_REF] Mönch | Production planning and control for semiconductor wafer fabrication facilities[END_REF].

Semiconductor manufacturing is composed of four major phases: wafer fabrication (fab), wafer probe, assembly, and final test. Wafer fabrication, often referred to as "front end", represents the most complicated, expensive and time-consuming phase of all four stages [START_REF] Uzsoy | A review of production planning and scheduling models in the semiconductor industry part I: system characteristics, performance evaluation and production planning[END_REF]. In this phase, hundreds of circuits are layered through successive operations on a silicon wafer. The manufacturing process in wafer fabs involves a highly complex sequence of processing operations which can be classified into various types, as for example: oxidation and thermal treatment, film deposition, planarization, photolithography, etching and ion implantation. These operations are repeated for each layer of circuitry on the wafer. Figure 1 presents a simplified view of the wafer fabrication process. Fig. 1: Wafer fabrication Process. [START_REF] Mönch | Production planning and control for semiconductor wafer fabrication facilities[END_REF] Each operation shown in Figure 1 can include multiple elementary steps (cleaning, process, measurement). The total number of steps per flow typically ranges between 400 and 800 for current production technologies and up to 1200 steps for latest generations. Some of the processing steps in a flow are performed on individual wafers, others on groups of wafers (lots), and still others on groups of lots (batches).

Steps performed on individual wafers or lots of wafers are referred to as serial steps, while those performed on groups of lots are called batch steps. A lot is generally composed of 25 wafers, while a typical batch contains up to six lots. In the considered case study, the lot requires 8 to 10 weeks to be processed. Steps are executed on more than one hundred workstations called "toolset" [START_REF] Mönch | Production planning and control for semiconductor wafer fabrication facilities[END_REF]. Due to flow re-entrance, lots visit the same toolset more than once during the manufacturing process.

For each step, the wafer has to be processed on various types of tools using a well-defined recipe. The recipe contains the detailed instructions to be used at the machine level in order to proceed the intended physical transformations or measurements. The identification of the candidate tools to be used is made through qualification of recipes on the tools. However, in HMLV fabs, because of multiple differences in hardware and software configurations, hence variety of recipes to be used, it is not possible to qualify all recipes on every machine. Qualification is one of the characteristics of the HMLV semiconductor manufacturing. It determines the processing authorization of a product on a machine. It acts like an eligibility constraint that allows production volume allocation of a product to a machine. It is known also as the process capability constraint [START_REF] Rowshannahad | Qualification management with batch size constraint[END_REF].

Besides, each toolset has an identified throughput considered as its capacity which refers to its upper loading threshold under a given product mix condition. To establish a feasible production schedule over a planning horizon of several weeks, thus requires to consider capability and capacity constraints. Moreover, as for other industries, semiconductor manufacturing facilities must respect customers delivery commitments to survive in competitive business environments. For HMLV fabs, actual cycle time is widely spread and skewed due to large variability of numerous sources: equipment heterogeneity, product priorities, low redundancy, steps qualifications, etc. It is then crucial to consider also variable cycle times while defining a production plan. In practice, fab's historical data and various applications of the queuing theory are often used.

In this paper, a mixed integer program (MIP) and a heuristic algorithm are proposed to project current inventory and new wafer starts throughout the remaining processing sequence, taking into account all the cited constraints. The objective is to establish a feasible midterm schedule, in a fast execution time (less than 5 minutes, the required computation time of capacity planners of the industrial partner), while minimizing lots delivery delays and optimizing workload balance among all toolsets. This study is applied to the Crolles 300 mm wafer fab of STMicroelectronics. Thus, data from actual production process are collected and used to evaluate the performance of the developed approaches. This paper is organized as follows. This section introduces the main characteristics of the considered industrial environment. In section 2, some background on existing related capacity planning problems is provided. In section 3, the problem is stated and the MIP formulation is presented. The proposed iterative heuristic algorithm is explained in section 4. Next, in section 5, experiments conducted and analyses carried out are discussed. Finally, section 6 draws conclusions and provides suggestions for future work.

II. PREVIOUS RELATED WORK

As the semiconductor industry is considered as one of the most complex manufacturing processes, many researchers have paid attention to the capacity planning problems encountered in this environment.

The various problems investigated have considered different phases of the manufacturing process of integrated circuits, different constraints, different methods and techniques used for capacity planning and different performance measures. Mönch et al. [START_REF] Mönch | Production planning and control for semiconductor wafer fabrication facilities[END_REF], Uzsoy et al. [START_REF] Uzsoy | A review of production planning and scheduling models in the semiconductor industry part I: system characteristics, performance evaluation and production planning[END_REF], [START_REF] Uzsoy | A review of production planning and scheduling models in the semiconductor industry part II: shop-floor control[END_REF] and Gupta et al. [START_REF] Gupta | Operational planning and control of semiconductor wafer fabrication[END_REF] have mentioned in their reviews different capacity planning techniques used in the semiconductor environment which can be classified in infinite and finite capacity planning techniques. They can also be divided, according to the length of the planning horizon, into long-term (strategic), mid-term (tactical) and short-term (operational) planning tools.

Among the methods used for capacity planning, classical techniques were successfully used in many industries especially for tactical and operational production planning, such as Material Requirement Planning (MRP) developed by Orlicky [START_REF] Orlicky | Material requirements planning[END_REF], Manufacturing Resource Planning (MRPII) [START_REF] Rondeau | Evolution of manufacturing planning and control systems: from reorder point to enterprise resource planning[END_REF], Just In Time (JIT) [START_REF] Golhar | The just-in-time philosophy: a literature review[END_REF] and Theory Of Constraints (TOC) [START_REF] Goldratt | Theory of constraints: What is this thing called Theory of Constraints and how should it be implemented[END_REF]. The application of these traditional techniques for capacity planning in semiconductor industry presents some shortcomings. Indeed, it is proven that MRP method can be inefficient and may produce unrealistic production schedules when used in field applications. It ignores capacity constraints and assumes fixed cycle times ( [START_REF] Rossi | A simulation-based finite capacity mrp procedure not depending on lead time estimation[END_REF], [START_REF] Jodlbauer | Material and capacity requirements planning with dynamic lead times[END_REF], [START_REF] Sun | Comparing dynamic risk-based scheduling methods with mrp via simulation[END_REF], [START_REF] Aouam | Zero-order production planning models with stochastic demand and workload-dependent lead times[END_REF]). However, in semiconductor facilities, cycle times depend on many factors, such as machine utilization rate, lot size, inventory and dispatching rules, and are thus variable. Either shortcoming above leads to infeasible production schedules, fluctuating workloads over time and significant users effort to adjust the plans. The JIT technique proves its strengths [START_REF] Levitt | Just-In-Time methods for semiconductor manufacturing[END_REF]; however, it presents some limitations in the high-mix low-volume production systems. It seems to be more suitable for a repetitive production environment with stable demand and low product mix [START_REF] Carlson | Mixed model assembly simulation[END_REF]. The TOC seems an efficient capacity planning technique in semiconductor industry [START_REF] Rippenhagen | Implementing the theory of constraints philosophy in highly reentrant systems[END_REF] but it considers only bottleneck resources and it can not deal with changes in the bottlenecks.

In addition to these classical industrial methods, authors use discrete event simulation models, queueing theory, linear programming and heuristics for capacity planning applied to semiconductor industry. Discrete event simulation is often used for capacity planning decisions in wafer fabs [START_REF] Resende | A program for simulation of semiconductor wafer fabrication[END_REF] in order to evaluate the performance of production planning strategies ( [START_REF] Tullis | Successful modeling of a semiconductor R & D facility[END_REF], [START_REF] Thompson | Using simulation-based finite capacity planning and scheduling software to improve cycle time in front end operations[END_REF], [START_REF] Fowler | Measurable improvements in cycle-time-constrained capacity[END_REF], [START_REF] Weintraub | A simulation-based finite capacity scheduling system[END_REF], [START_REF] Grewal | Integrating targeted cycle-time reduction into the capital planning process[END_REF]). Indeed, discrete-event simulation is considered as the only practical method that explicitly calculates the cycle time as a function of resource availability and production rate. The simulation model can be used also to determine bottlenecks under a given product mix and to make strategic decisions concerning equipment purchase [START_REF] Potti | Using simulation to improve semiconductor manufacturing[END_REF]. However, simulation models used for capacity planning in the semiconductor industry present some severe limitations. Their set-up is very time-consuming due to the volume and often complexity of the data required for the models involved. Moreover, these models do not provide a means for optimization of the plan ( [START_REF] Pritsker | Production scheduling using FAC-TOR[END_REF], [START_REF] Ignizio | Fab simulation and variability[END_REF]).

Concerning queueing network models, Shanthikumar et al. [START_REF] Shanthikumar | Queueing theory for semiconductor manufacturing systems: A survey and open problems[END_REF] presented a survey of the different applications of queueing theory for semiconductor manufacturing systems. They recognized that in spite of fast computing time compared with simulation models, the accuracy of classical queueing models is not satisfactory due to the complexity of the semiconductor manufacturing process.

The linear programming (LP) approach is widely applied to specific issues encountered in capacity planning for the semiconductor industry. Mixed-integer programming (MIP) models are developed for strategic planning in order to maximize the profit ( [START_REF] Bermon | Capacity optimization planning system (CAPS)[END_REF], [START_REF] Swaminathan | Tool capacity planning for semiconductor fabrication facilities under demand uncertainty[END_REF], [START_REF] Barahona | Robust capacity planning in semiconductor manufacturing[END_REF]) or to minimize the machine tool operating costs, new tool acquisition costs, and inventory holding costs taking into account capacity constraints [START_REF] Erengüc | Tool capacity planning in semiconductor manufacturing[END_REF]. A good source of previous work related to more strategic capacity planning is provided by Geng and Jiang [START_REF] Geng | A review on strategic capacity planning for the semiconductor manufacturing industry[END_REF].

LP (sometimes in combination with discrete-event simulation) is also used to solve medium-term finite capacity planning problems. The work of Hung and Leachman [START_REF] Hung | A production planning methodology for semiconductor manufacturing based on iterative simulation and linear programming calculations[END_REF] is an example of such an approach. Leachman [START_REF] Leachman | Modeling techniques for automated production planning in the semiconductor industry[END_REF] used LP for production planning and presents a corporate capacity planning model, which includes multiple facilities integrated with the production process. Habla et al. [START_REF] Habla | A finite capacity production planning approach for semiconductor manufacturing[END_REF] suggested a MIP formulation to determine completion time targets for the lots on bottleneck steps. Bermon et al. [START_REF] Bermon | Capacity analysis of complex manufacturing facilities[END_REF] introduced a linear programming model to analyze the capacity of large and complex manufacturing production lines.

Due to the intractability of LP models, they are generally combined with heuristics such as genetic algorithms [START_REF] Hsiung | Tool planning in multiple product-mix under cycle time constraints for wafer foundries using genetic algorithm[END_REF] or decomposition techniques as Benders [START_REF] Bard | Midterm planning to minimize deviations from daily target outputs in semiconductor manufacturing[END_REF] or Lagrangian relaxation ( [START_REF] Erengüc | Tool capacity planning in semiconductor manufacturing[END_REF], [START_REF] Habla | A finite capacity production planning approach for semiconductor manufacturing[END_REF]) to reduce execution time. Besides, approximate methods have also been widely used to develop either infinite or finite capacity planning systems for the semiconductor industry.

Infinite capacity planning systems are developed to estimate the future loading of equipment in order to identify bottleneck resources and to balance the loading of each production resource over the planning horizon ( [START_REF] Chen | Capacity planning with capability for multiple semiconductor manufacturing fabs[END_REF], [START_REF] Chen | Capacity requirements planning for twin fabs of wafer fabrication[END_REF], [START_REF] Chen | Infinite capacity planning for IC packaging plants[END_REF]).

Bearing in mind the importance of capacity constraints, many authors developed finite capacity planning systems using algorithmic approaches. Fargher et al. [START_REF] Fargher | A planner and scheduler for semiconductor manufacturing[END_REF] used a beam-search algorithm in combination with backtracking steps for lot release and for the determination of schedules in an aggregated sense.

Horiguchi et al. [START_REF] Horiguchi | Finitecapacity production planning algorithms for a semiconductor wafer fabrication facility[END_REF] proposed an algorithm that estimates the start and finish date of each job scheduled on each critical resource: their algorithm considers the available time for all the feasible combinations of time bucket and critical resource, and it reduces the available time whenever a new production order is added to the schedule. This approach, due to the high aggregation level in modeling resources and relationships, might lead to orders overlapping on the same resource in the same time bucket (i.e. infeasible plans).

Habenicht and Mönch [START_REF] Habenicht | A finite-capacity beam-search-algorithm for production scheduling in semiconductor manufacturing[END_REF] used also a beam-search algorithm to determine planned start and completion dates for the macro operations (sets of consecutive process steps) of a lot.

Chua et al. [START_REF] Chua | An intelligent multi-constraint finite capacity-based lot release system for semiconductor backend assembly environment[END_REF] developed an intelligent multi-constraint finite capacity-based lot release system. This system has been designed, developed and implemented to solve the lot release problems in a discrete manufacturing environment with a huge product mix and multiple capacity constraints.

In this study, we are interested in operations research related (or mathematical) optimization approaches and we consider a medium-term finite capacity planning problem, applied to a semiconductor production line. So far, the literature review has pointed out that the debate about this problem is still open, and the proposed approaches by several authors still have some limits. Table II presents a taxonomy of studies considering the same problem and using operations research solving tools.

Steps cycle time can be either defined as a fixed input by the proposed approach or a variable output of the procedure. Capability constraints are relevant, since they can be embedded or not in the proposed procedure. Lots due dates, relevant as well, can be considered as input parameters or not. Finally, the model can be tested via data generated by authors (random instances) or through data from real-life production systems (real case).

Whilst capacity and cycle time are tightly linked one another through the Little's law [START_REF] Little | A proof for the queuing formula: L= λ w[END_REF], cycle time is considered by most approaches as a fixed input parameter. Moreover, some methods ignore capability constraints thus leading to infeasible production plans. Finally, the applicability in field to real-life companies has not been reported for all the anterior studies.

Furthermore, Table II presents, for each study, its algorithmic and operational objectives. As one can notice, the main issues treated in the existing studies are generally limited to dispatching rules and release control policies which are outside the scope of this paper.

In the literature, there are few studies considering the WIP projection problem in the semiconductor industry ( [START_REF] Kim | Decomposition method application to a large scale linear programming wip projection model[END_REF], [START_REF] Lee | Production planning in semiconductor wafer fab considering variable cycle times[END_REF], [START_REF] Govind | Setting performance targets in a 300mm wafer fabrication facility[END_REF]). In these works, authors consider different objectives and do not take into account all the cited constraints.

Kim and Leachman [START_REF] Kim | Decomposition method application to a large scale linear programming wip projection model[END_REF] proposed a LP formulation and a decomposition heuristic method to determine net demand and net resource capacities taking into account capacity constraints. They tested their approaches using random data. Lee et al. [START_REF] Lee | Production planning in semiconductor wafer fab considering variable cycle times[END_REF] employed deterministic linear programming techniques for the WIP projection problem in the wafer fab, that explicitly considers the variable cycle time. Govind and Fronckowiak [START_REF] Govind | Setting performance targets in a 300mm wafer fabrication facility[END_REF] consider WIP projection problem to measure production performance at IBM's 300 mm wafer fab by computing productivity and WIP targets at infinite capacity. As one can see, even if some papers tackle similar planning problems, none of the already proposed models explicitly address our specific problem.

The research work outlined here tried to overcome some of the limits above: the proposed finite capacity planning algorithm does not consider fixed steps cycle time, it takes into account lots due dates and it has been tested in a real-life industrial context. Furthermore, it meets the key requirement of semiconductor industrials, consisting on fast computing of feasible production plans (in five minutes at most on a personal computer) to facilitate "what-if" analysis.

III. MATHEMATICAL FORMULATION OF THE PROBLEM

A. Problem description

The WIP projection problem may be defined as follows. A set of lots l ∈ {1, . . . , L}, composed of Q l wafers each, is considered. For each lot l, it remains an identified number of steps S l to be processed on a time horizon discretized in T periods t ∈ {1, . . . , T } of equal length P t . Each lot of a weight w l indicating its priority, has a release date r l and a due date d l .

The performance measurement to be minimized in this problem is total weighted tardiness TWT. TWT is a measure that incurs a penalty for each lot that finishes processing after its promised delivery date. This penalty increases with the magnitude of the tardiness, and therefore schedules that minimize the weighted (by lot priority) sum of penalties provide good on-time delivery performance, whereas higher values of total weighted tardiness indicate that many important lots are not being delivered on time. Indeed, a processing schedule will provide a completion time, C l , for each lot. The tardiness, T l , of lot l is then defined as

T l = max(0, C l -d l ).
The weighted tardiness of the lot l (W T l ) is defined as W T l = (w l × T l ). Total weighted tardiness computes the weighted sum of tardiness values: T W T = l W T l . Each remaining step s l ∈ {1, . . . , S l } of the lot l is processed on one or several qualified toolsets i ∈ {1, . . . , I}. The quantity of wafers of a lot l assigned to the toolset i, processing the step s l during the period t, is denoted a s l ,l,i,t . It has a waiting time wt s l ,l and it consumes a unit processing time p s l ,l,i on each of its qualified processing toolsets. It also has a start date s s l ,l and an end date e s l ,l (Figure 2). Each toolset i has a finite capacity C i,t , which gives the maximal loading L i,t over a period t.

B. Mixed-Integer Program

In this subsection, an appropriate MIP formulation is presented for the multi-product, multi-period and multi-resource capacity planning problem. The proposed MIP is similar to LP capacity planning models that can be found in standard textbooks, with some variations and extensions. Using the Loading of toolset i in period t y s l ,l,t = s s l ,l if the step s l of lot l is released in period t, 0 otherwise x s l ,l,t =1 if step s l of lot l is processed in period t, 0 otherwise notation presented above, the MIP is as follows:

min l w l T l (1)
s.c. s 1,l ≥ r l l = 1, . . . , L (2)
s s l ,l + i p s l ,l,i × a s l ,l,i,t × x s l ,l,t = e s l ,l s l = 1, . . . , S l , l = 1, . . . , L (3) s s l ,l ≥ e s l -1,l s l = 2, . . . , S l , l = 1, . . . , L

t y s l ,l,t = s s l ,l s l = 1, . . . , S l , l = 1, . . . , L (5)

t x s l ,l,t = 1 s l = 1, . . . , S l , l = 1, . . . , L (6)
C l = e S l ,l l = 1, . . . , L (7)
T l ≥ C l -d l l = 1, . . . , L (8)
T l ≥ 0 l = 1, . . . , L (9) t × P t × x s l ,l,t ≤ y s l ,l,t s l = 1, . . . , S l , l = 1, . . . , L, t = 1, . . . , T (10)
(t + 1) × P t × x s l ,l,t > y s l ,l,t s l = 1, . . . , S l , l = 1, . . . , L, t = 1, . . . , T -1 (11) L i,t = l s l p s l ,l,i × x s l ,l,t × a s l ,l,i,t i = 1, . . . , I, t = 1, . . . , T (12)
L i,t ≤ C i,t i = 1, . . . , I, t = 1, . . . , T (13)
x s l ,l,t = {0, 1} s l = 1, . . . , S l , l = 1, . . . , L, t = 1, . . . , T The objective function (1) minimizes the total weighted tardiness (TWT). The MIP constraints can be classified in two kinds: temporal constraints ((2). . . [START_REF] Jodlbauer | Material and capacity requirements planning with dynamic lead times[END_REF]) and cumulative constraints (constraints ( 12)-( 13)). Constraints (2) define the start date of the first remaining step for each lot. The end date of each remaining step of each lot is computed using constraints [START_REF] Rowshannahad | Qualification management with batch size constraint[END_REF]. Constraints (4) present precedence constraints of processing steps. Constraints [START_REF] Gupta | Operational planning and control of semiconductor wafer fabrication[END_REF] guarantee that each remaining step of each lot is released once. Constraints [START_REF] Orlicky | Material requirements planning[END_REF] verify that each remaining step of each lot is processed once during the planning horizon. Constraints [START_REF] Rondeau | Evolution of manufacturing planning and control systems: from reorder point to enterprise resource planning[END_REF] define the lots completion date. Constraints ( 8) and ( 9) compute the tardiness for each lot. Constraints [START_REF] Rossi | A simulation-based finite capacity mrp procedure not depending on lead time estimation[END_REF] and [START_REF] Jodlbauer | Material and capacity requirements planning with dynamic lead times[END_REF] indicate that each remaining step of each lot is processed in one period. Constraints [START_REF] Sun | Comparing dynamic risk-based scheduling methods with mrp via simulation[END_REF] calculate the workload accumulated by each toolset over each period taking into account the qualification of the toolset to the processed step and the quantity of wafers assigned to the considered toolset. Constraints [START_REF] Aouam | Zero-order production planning models with stochastic demand and workload-dependent lead times[END_REF] are the capacity constraints. Constraints ( 14) are the binary constraints for the decision variable.

The mathematical model presented above has been solved by ILOG CPLEX solver. Experiments were run on an Intel Core TM i5 PC running a 2.7 GHz processor and 4 GB of RAM. Tests have been performed on 30 randomly generated instances of the problem in order to highlight the main characteristics of the industrial data and to maintain a certain degree of generality in order to preserve all the difficulty of the problem. Indeed, based on the observation made in the literature, we identified seven important problem parameters which could affect the performance of the proposed approach: number of lots (L), maximum number of remaining steps for each lot (max S l ), number of toolsets (I), length of the planning horizon (T ), lots steps unit processing times (p s l ,l,i ), lots due dates (d l ) and machines capacities (C i,t ).

We consider the cases of 3, 5, 10, 20, 100 and 300 parallel toolsets with a fixed capacity corresponding to the maximum equipment utilization rate which is equal to 100%. The lots weights w l are chosen from a uniform distribution over (0,1). The lots release dates r l and lots quantity of wafers Q l are supposed equal to 0 and 25 for all lots, respectively. The range of lots due dates d l and steps unit processing times p s l ,l,i is extracted from real data. d l are ranging from 1 to 210 days relative to the release date and p s l ,l,i range between 0.0005 and 0.5 hours. The planning horizon is set to 24 periods (weeks). Table III presents the different tests parameters generating 30 instances.

Optimal results were obtained in reasonable execution time while testing the MIP on instances of reduced size. Further increasing the size of the tested instances (up to about 4000 steps plan), the resolution of MIP was halted as it required a very large amount of time and computer memory (Figure 3). Indeed, the whole real problem presents 70 742 400 constraints and 69 371 200 variables. It corresponds to a WIP composed of 2000 lots, each lot having a maximum of 680 remaining steps to process on 300 toolsets over a planning horizon composed of 24 periods (weeks). Thus, the size of real instances is obviously too large to be solved using the proposed MIP (Figure 3). From the empirical evidence on the computational difficulties in getting optimal schedule considering lots due dates and capacity constraints, it is obvious that the problem of WIP projection applied to the real case study will be computationally intractable. Furthermore, Garey and Johnson [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NPCompleteness[END_REF] highlighted in their study that production planning, capacity planning and scheduling problems in complex job shops like semiconductor manufacturing are known as strongly NP-hard problems. This has motivated us to develop a heuristic algorithm for the research problem considered in this study to provide near optimal solutions and/or efficient solution in a reasonable time.

The proposed heuristic algorithm is presented in the next section.

IV. HEURISTIC ALGORITHM An alternative methodology for the above problem should be accurate and, at the same time, fast and small enough to be stored and implemented in a mainframe or work station computer system. Bearing this in mind, a heuristic approach for WIP projection problem in HMLV semiconductor manufacturing line has been developed. It is an iterative algorithm composed of three main modules: (i) WIP projection at infinite capacity, (ii) workload accumulation and capacity analysis and (iii) workload and capacity balancing. The algorithm is executed by iterations on periods of the planning horizon. The principle of iterative running of the algorithm is inspired from the literature [START_REF] Winkler | Modeling of complex decision making using forward simulation[END_REF] and the detailed scheduling in the commercial ERP/APS. For each defined period, WIP projection module estimates the evolution of the WIP, lot by lot, based on lots due dates. Then, workload accumulation module calculates the expected equipment loading. In case of toolsets oversaturation i.e. the loading of toolsets exceeds their maximal capacity, workload and capacity balancing module is employed to reduce toolsets loading by shifting their affected steps to subsequent periods. The following sections will detail the three major modules. Figure 4 depicts the flow of the developed system. Fig. 4: Finite Capacity Planning Algorithm Flow.

A. WIP projection module

The objective of this module is to push lots, one by one, forward along their routes, from their current positions up to their due dates. It also aims to compute over a period the activity required by step to ensure the delivery plan.

For each selected time bucket t of the planning horizon, this module requires the following data input:

1) WIP status and wafer starts at the beginning of the considered projection period (position r l , quantity Q l ), 2) Routing information, including a partition of each route into consecutive steps, 3) Steps unit processing times p s l ,l,i , 4) Lots due dates d l and weights w l , 5) Flow factor Xf actor s l ,l that reflects possible waiting times between consecutive process steps to achieve the target cycle time, extracted from historical data. It is defined as the step mean cycle time divided by the step raw processing time RPT [START_REF] Martin | Key factors in designing a manufacturing line to maximize tool utilization and minimize turnaround time[END_REF]. WIP projection module includes three steps. Step 1 computes, for each lot, from its position in the route, four parameters which are remaining process time RemP T l , remaining reference cycle time RemRef CT l , remaining expected cycle time RemExpCT l and cycle time coefficient CT Coef f l . RemP T l is equal to the sum of lot remaining steps unit process time multiplied by lot quantity of wafers Q l .

RemP T l = S l s l =1 I i=1 p s l ,l,i × Q l ( 15
)
RemRef CT l corresponds to the sum of the reference cycle times of lot remaining steps Ref CT s l ,l .

RemRef CT l = S l s l =1 Ref CT s l ,l (16)
In the industrial context considered, each step has a reference cycle time, extracted from historical data, named Ref CT s l ,l .

Ref CT s l ,l corresponds to the product of the unit step process time with the quantity of lot wafers Q l and the flow factor Xf actor s l ,l . It is the maximum amount of time that a lot would spend at that step, including waiting and processing times.

Ref CT s l ,l = I i=1 p s l ,l,i × Q l × Xf actor s l ,l (17)
RemExpCT l is equal to the maximum between the difference between the due date and the current time t and RemP T l .

RemExpCT l = max(d l -t, RemP T l )

The lot cycle time coefficient CT Coef f l identifies the necessary and sufficient speed for lots to achieve their due date according to the reference cycle time. It is equal to the ratio between lot remaining expected cycle time RemExpCT l and lot remaining reference cycle time RemRef CT l :

CT Coef f l = RemExpCT l RemRef CT l ( 19
)
In step 2, the RemExpCT l is split on the elementary steps of each lot l to compute an expected cycle time per step ExpCT s l ,l which is equal to the product of ObjCT s l ,l and CT Coef f l .

ExpCT s l ,l = Ref CT s l ,l × CT Coef f l (20)
Equation [START_REF] Fowler | Measurable improvements in cycle-time-constrained capacity[END_REF] gives a rough estimation of queuing time at each step. Hence, waiting time by step wt s l ,l can be computed:

wt s l ,l = ExpCT s l ,l - I i=1 p s l ,l,i × Q l (21)
In step 3, having the waiting time by step, start dates and end dates for all lots remaining steps, decision variables x l,s l ,t , lots completion date and tardiness are computed. Figure 5 illustrates an example of 2 lots with different due dates, having 3 remaining steps each. The first one has an earlier due date i.e. a higher priority and less RemExpCT l than the second.

Using the classical projection based on historical data, the two lots have the same distribution of remaining steps waiting times over the planning horizon, independently of their due dates, because they have the same remaining process time RemP T l . However, the proposed projection module allocates steps expected cycle times taking into account lots priorities i.e. due dates. Indeed, there are multiple priority levels of production lots. Production priorities can be divided into two levels according to the urgency of delivery: hot and standard. So, to respect these priorities, the projection module shrinks steps waiting times in order to satisfy the hot lot's due date. However, for a standard lot, it extends steps waiting times respecting the lot due date.

To further explain the concept of WIP projection, a simple random instance is tested. The considered WIP consists of 10 lots of 25 wafers each, following different routes, and having different due dates. Table IV presents, for each lot, the number of remaining steps, RemP T l , RemRef CT l , RemExpCT l and CT Coef f l .

Figure 6 illustrates projection results of the 10 lots during the first period of the planning horizon. For some lots, a sequence of steps is repeated twice (lots 1,4,5,7,8,9) i.e. lots visit the same toolset twice which illustrates the re-entrant flows. Figure 6 shows start and end dates, waiting time and processing time for each remaining process step during the considered period. Some steps (step 4.5 and step 10.4) start in the first period and finish in the subsequent periods of the planning horizon. This figure demonstrates that the projection engine allows the extension of steps waiting times for lots having a far due date which is the case of lots 1 and 6 and it shrinks steps cycle times in case of close due date for lots 2, 4 and 8. Lots 2, 4 and 8 are not delivered on time. Their due dates are not reachable so their fab-out dates are equal to the sum of the current date (t = 0) and the remaining process time RemP T l .

B. Workload accumulation and capacity analysis module

After WIP projection, the loading of toolsets, over each considered period L i,t , is computed based on the assumption of infinite capacities.

To optimize the computation time, toolsets are distributed in balancing groups. A balancing group is a set of toolsets x s l ,l,t Decision variables values, results of WIP projection module in period t v s l ,l,r =1 if recipe r corresponds to step s l of lot l, 0 otherwise a s l ,l,i Quantity of wafers of lot l in step s l processed by the toolset i p r,i

Processing time of recipe r on toolset i Q l

Quantity of wafers of lot l Decision variables: Loading, for a given recipe r, of the least loaded toolset among those on which r is qualified Using the above parameters, and decision variables, the linear program formulation can be represented as follows:

L i,
       Minimize α • Lmax -β • Lmin + γ • r Lmax r -δ • r Lmin r + δ • ( I b i L i,t -Lmin) with α = I 2 b , β = I b , γ = 1, δ = 1/I b s.t. L i,t = r p r,i × W r,i i = 1, . . . , I b (22)
Ir i=1 W r,i = l s l x s l ,l,t × v s l ,l,r × a s l ,l,i r = 1, . . . , R b (23)
L i,t ≥ Lmin r r = 1, . . . , R b , i = 1, . . . , I r (24) L i,t ≤ Lmax r r = 1, . . . , R b , i = 1, . . . , I r (25) L i,t ≥ Lmin i = 1, . . . , I b (26) L i,t ≤ Lmax i = 1, . . . , I b (27)
The linear program seeks to :

• Minimize the workload of the most loaded toolset in the balancing group Lmax. • Maximize the workload of the least loaded toolset in the balancing group Lmin. • Minimize the total workload of toolsets

I b
i L i and maximize the total workload of the least loaded toolset per recipe r Lmin r , with the same degree of priority.

• Minimize the total workload of the most loaded toolsets per recipe r Lmax r .

For the example cited above, the remaining steps of 10 lots are considered to be processed by 6 toolsets {M 1, M 2, M 3, M 4, M 5, M 6}. These toolsets are classified into 4 balancing-groups {M 1, M 6}, {M 2, M 4}, {M 3} and {M 5}. Figure 7 illustrates the saturation percentage of toolsets i.e. the ratio of the loading to the available capacity during the first period of the planning horizon ( Li,1

Ci,1 , i = 1..6) while processing the remaining steps ordered in increasing order of the start date. In this example, the capacity of all the considered toolsets (C i,1 , i = 1..6) is equal to 24 hours/day. Figure 7 shows that there are two over-saturated toolsets (M 2 and M 6) which workloads exceed saturation threshold. Fig. 7: Workload accumulation at infinite capacity for the first period of the planning horizon.

C. Workload/Capacity Balancing Module

As a result of the workload accumulation module, loading of some toolsets may exceed their maximal capacities i.e. constraints [START_REF] Aouam | Zero-order production planning models with stochastic demand and workload-dependent lead times[END_REF] are not satisfied. In this case, the toolset is unable to process all its affected steps during the considered period so its loading should be balanced over subsequent periods. The principle of this module is to postpone additional lots in order to bring back workload of over-saturated toolsets below their maximal saturation and to smooth the activity over the planning horizon. The algorithm for workload/capacity balancing module is as follows:

1) Sort toolsets in decreasing order of saturation.

2) Select lots executed on over-saturated toolsets.

3) Sort selected lots in increasing order of a computed ranking coefficient (rankingCoef f l ). The rankingCoef f l illustrates the priority of the lot in terms of its position in the process sequence of the considered toolset and the urgency of delivery. The position of a lot in the process sequence of a toolset is determined by the processing date of its last remaining step treated by the considered toolset denoted s S l ,l,t . The urgency of delivery is defined by the lot cycle time coefficient (CT Coef f l ). To compute the rankingCoef f l , the lot position in the process se-quence is normalized by the period length P t . Hence, the rankingCoef f l is equal to :

rankingCoef f l = 1 CT Coef f l + s S l ,l,t P t (28)
4) For the first selected lot in the sorted list, shift the last step executed in the considered over-saturated toolset and its successors to the next period. 5) Remove the processing time of shifted steps s l from the loading of its qualified processing toolsets while considering the quantity of wafers a s l ,l ,i processed by each toolset. The new value of the loading of the considered toolsets L i,t is, then, equal to :

L i,t = L i,t - l s l (a s l ,l ,i × p s l ,l ,i × x l ,s l ,t ) (29)
6) Repeat steps 2, 3, 4 and 5 for all toolsets until the saturation criterion is satisfied for all toolsets over the period t. Hence, this module modifies steps projection at period t as well as the WIP for the beginning of period t + 1. For instance, to balance the capacity and the workload of the over-saturated toolsets M 2 and M 6 in the considered example, the balancing module selects M 2 as the most over-saturated toolset ( L2,1 C2,1 = 109.3%). Then, it selects lots 2, 4, 5, 7, 8 and 10 processed by this resource (Figure 7). These lots are sorted in increasing order of rankingCoef f l as it is mentioned in Table V. In order to decrease the loading of the toolset M 2, steps 5.3 and 7.2 and its successors are shifted to the next period of the planning horizon. Hence, the loading of M 2 becomes less than its maximum capacity: L6,1 C6,1 = 96.4%. M 4 is also qualified for step 5.3, so its loading decreases by 5.06%. Step 5.4 projected in the first period is also postponed as it is the successor of the shifted step 5.3. Thus, the loading of M 1 processing step 5.4 becomes equal to 58.83%. Shifting the successors of step 7.2 (steps 7.3, 7.4 and 7.5) leads to decreasing the loading of toolsets M 1, M 4 and M 5. The same algorithm is applied to the toolset M 6 by shifting step 9.1 and its successor step 9.2. So, its loading decreases to 72%.

The toolsets workload obtained after steps shifting is illustrated in Figure 8. Table VI presents the WIP and the computed parameters (RemP T l , RemObjCT l , RemExpCT l and CT Coef f l ) in the beginning of the next period.

The proposed approach is tested over a five-day planning horizon. Indeed, as mix variations were present in industrial dataset used for this study, it was decided to focus on a very short planning horizon to evaluate the proposed approach. Clearly, the shorter the length of the period, the more accurate the results of the approach. The final obtained schedule for this instance is illustrated in Figure 9. For this instance, the TWT is equal to 1.46 days and we have five delayed lots.

Fig. 9: The obtained schedule using heuristic approach.

V. RESULTS AND DISCUSSION

The proposed algorithm is coded in Java and it is tested on a 4 GigaOctet RAM and 2.7 GigaHertz processor computer. We conducted two types of experiments to evaluate the performance of the proposed approach. The first type corresponds to a comparison between the exact method and the heuristic using a set of randomly generated instances. In the second type of experiment, we compare the projected schedule obtained by the proposed approach using real data with what is really going on in the wafer fab following this schedule.

A. Evaluation of the proposed heuristic algorithm in comparison with optimal solution

For this evaluation, random instances were generated and solved using the MIP and the proposed heuristic algorithm. The parameter, number of lots (L), assumes only five levels (L=10, 15, 20, 50 and L=100). The parameters generated for the proposed instances are presented in Table VII. So, three random problem instances for each fixed parameter combination are obtained, giving a total of 270 test problems. Each of the 270 problem instances generated has been solved using the ILOG CPLEX solver and the proposed heuristic algorithm. The results on TWT obtained for each instance using MIP model and using the proposed iterative algorithm are recorded. Based on these results, the heuristic solution matched exactly with the optimal solution 53 times. Furthermore, for each instance with a size L × max S l × I, we compute:

TABLE VII: Test data parameters

• The absolute deviation= |TWT value from a heuristic algorithm -optimal TWT value| • The relative deviation= absolute deviation value optimal TWT value Figure 10 shows the relative deviation over 270 instances plotted against the absolute deviation.

In this figure, we can define four zones or classes according to the size of the instance:

B. Experimental tests on real fab data

The aim of this section is to evaluate the ability of the proposed approach to tackle real world problems. The test of the real instance (L=2000, max S l =680, I=300, T =24), unsolved in reasonable execution time using the MIP approach in Section 3, is treated. The execution time of this instance with the proposed algorithm is around 30 seconds. In the calculated production schedule, 80 % of projected lots are delivered on time. Furthermore, the saturation of toolsets is kept below the pre-defined saturation threshold while minimizing lots lateness. Figure 11 illustrates the obtained weekly saturation at infinite and finite capacity of a photo-lithography toolset considered as a bottleneck. In semiconductor fabs, several indicators are used to measure performance [START_REF] Montoya-Torres | Manufacturing performance evaluation in wafer semiconductor factories[END_REF]. Jointly with evaluate the waiting times of lots for processing. In this section, we compare the cited indicators of performance of the heuristic solution with the indicators determined in the real production line. To ensure this experiment, six tests have been performed on actual instances issued over four months of production: September, October, November and December 2015. We have made projections in six different periods (week1, week2, week3, week4, week5 and week6) and we have determined the three indicators for each projection. For confidential reasons, we are not allowed to provide the real values of the fab. This is why, we compute the relative deviation between the predicted value and the real one for each period of the planning horizon:

Relative deviation = |Estimated value-real value| real value 1) Analysis based on the performance measure: number of moves: Figure 12 shows relative deviations of number of moves over 15 time buckets (weeks) of the planning horizon. It illustrates that in the first 6 periods for the different instances, the relative deviation between the real number of processing steps and the calculated value is low. The average of the average relative deviations over six periods for the different tests is equal to 12.7%, reflecting a small difference between the estimated number of moves and the achieved one. Further being away from the beginning of the projection, the relative deviation between the obtained solution and the real number of moves increases which is explained by the variability of the process. Hence, there is a convergence between what is estimated and what is achieved in terms of periodic activity for a short-term planning horizon.

2) Analysis based on the performance measure: number of moves by usage: To evaluate how the heuristic solution anticipates the fab loading, we compute the absolute deviation of the number of moves by set of toolsets sharing the same qualifications named "usage" over the six instances for each period of the planning horizon. Figure 13 shows the difference between the total number of completed steps processed by two types of bottleneck usages (photolithography and etching). For this indicator also, we observe a convergence between the planning and the real process for the first 6 periods with an average of the average relative deviations over these periods equal to 6.5% for the usage of photolithography and 12.3% for the usage of etching. Therefore, the heuristic provides good estimations of the tools loading close to the real workload while respecting capacity constraints. 3) Analysis based on the performance measure: TWT: To compare between the real total weighted tardiness and the obtained value of this indicator using the iterative algorithm for the six tests, absolute and relative deviations are computed and reported in Table VIII. From Table VIII, we note that the estimated value of TWT is close to the real tardiness while respecting lots due dates. Indeed, the average of relative deviations over six instances is equal to 4%.

VI. CONCLUSIONS AND PERSPECTIVES

This paper has examined the problem of WIP projection at finite capacity to minimize the TWT, and has proven empirically the computational complexity in obtaining optimal solution and suggested a simple, fast and efficient heuristic. The motivation for this research is to compute a feasible production plan to drive the execution of wafer fabs. This problem is of considerable practical value because the heuristic, proposed in this paper, can be used in planning of a large number of production lots while respecting lots due dates and capacity and capability constraints. The computational tests, made on real production instances, showed that acceptable solutions are obtained in reasonable execution time. Indeed, the TWT could be minimized and the average tool utilization rate could be balanced significantly by using the developed system. Besides, the computation for real instances is achieved in around 30 seconds which is efficient for planning problems with a horizon of weeks up to months in real situations. Hence, this decision support tool outperforms simulation and analytic models for establishing a feasible production schedule rapidly. Finally, it is observed (as well as statistically verified) from the results of the comparison of the different criteria (total number of moves, number of moves by usage and TWT) an obvious convergence between what is predicted using the developed approach and what is achieved in the real process over a short-term planning horizon. These results show that the implementation of the finite capacity planning system in real fabs seems very interesting to minimize lots lateness and to establish a feasible production schedule. There are a number of interesting extensions of the problems that can be pursued. The first important issue would be to perform a more thorough multi-criteria analysis while shifting lots to balance toolsets loadings. Besides, it is necessary to implement the developed finite capacity planning system in the production plant to guarantee the performance of the solution. Considering other specificities of semiconductor industry such as batching or sequence dependent setup times may be interesting to enhance the accuracy of the developed system.

Fig. 2 :

 2 Fig. 2: Problem description.

Fig. 3 :

 3 Fig. 3: Limits of MIP resolution.

Fig. 5 :

 5 Fig. 5: Principle of cycle time computation for each processing step of a lot.

Fig. 6 :B

 6 Fig. 6: Simple instance: Production schedule at infinite capacity.

Fig. 8 :

 8 Fig. 8: Workload accumulation at finite capacity.

  Number of lots (L) 10, 15, 20, 50, 100 5 Maximum number of remaining steps of lot l (max S l ) times p s l ,l,i 0.0001× uniform(5

Fig. 10 :

 10 Fig. 10: Comparaison between the optimal and the approximate solution.

Fig. 11 :

 11 Fig. 11: Weekly saturation of a photo-lithography toolset at infinite and finite capacity.

Fig. 12 :

 12 Fig. 12: Number of moves comparison actual versus forecast.

Fig. 13 :

 13 Fig. 13: Total number of moves processed by photolithography and etching usages comparison actual versus forecast.

Table I

 I

summarizes the notation.

TABLE I :

 I Summary of problem notation Unit processing time of step s l of lot l on qualified toolset i, 0 on non-qualified toolset i C i,t Capacity of toolset i in period t a s l ,l,i Quantity of wafers of lot l in step s l processed by the toolset i Decision variables Description s s l ,l Start date of step s l of lot l e s l ,l End date of step s l of lot l C l

	Indices	Description
	L	Number of lots
	l = 1..L	Lot index
	S l	Number of remaining steps of lot l
	s l = 1..S l	Lot's step index
	I	Number of toolsets
	i = 1..I	Toolset index
	T	Number of time buckets
	t = 1..T	Period index
	Parameters	Description
	Pt	Length of period t
	Q l	Quantity of wafers of lot l
	r l	Release date of lot l
	w l	Weight of lot l
	d l	Due date of lot l
	p s l ,l,i
		Completion date of lot l
	T l	Tardiness of lot l
	L i,t

TABLE II :

 II Taxonomy of the literature of finite capacity planning approaches applied to semiconductor industry.

	Test				Real case study			Real case study		Example		Real case study		Example	Example	Real case study	Real case study
		Variable	cycle times
	Constraints and assumptions	Capability Due dates constraints
		Capacity	constraints
	Objectives	Algorithmic Operational	Maximize the profit Determine wafer release quantities	Generate capacity-	Maximize the profit feasible start and out	schedules	Minimize total weighted Determine completion time targets for the tardiness bottleneck steps	Reduce cycle time and Determine the work to release into the factory at the variance of cycle time any time	Estimate the start and Improve delivery finish date of each job performance and system scheduled on each critical predictability resource	Determine the start date Establish a feasible and the end date of each production plan operation of the lot	Compute orders release	dates for semiconductor Solve lot release problem	back end assembly	Minimize total weighted Establish a feasible target	tardiness production plan
	Reference		Hung and	Leachman	[32], Bermon	et al. [35]	Leachman [33]		Habla et al.	[34]	Fargher et al.	[41]	Horiguchi et	al. [42]	Habenicht and	Mönch [43]	Chua et al. [44]	Our study
	Approach		Linear	programming	based					Algorithms	/heuristics	based

TABLE III :

 III Summary of tests parameters

TABLE IV :

 IV Data for simple instance

	Lot l	Weight	Number	RemP T l	Rem-	Rem-	CT-
		w l	of re-	in	Ref CT l	ExpCT l	Coef f l
			maining	days	in days	in days
			steps S l
	Lot 1	0.33	6	1.1	1.6	5	3.125
	Lot 2	1	4	0.8	1.1	0.5	0.45
	Lot 3	0.5	2	0.25	0.41	1.5	3.65
	Lot 4	0.5	8	1.7	2.3	1.5	0.65
	Lot 5	0.5	6	1	1.4	1.5	1.07
	Lot 6	0.33	4	0.75	1.02	5	4.9
	Lot 7	0.5	8	0.86	1.05	1.5	1.43
	Lot 8	1	4	0.8	1.05	0.5	0.48
	Lot 9	0.5	4	0.8	1.05	1.5	1.43
	Lot 10	0.5	6	1.4	1.9	1.5	0.79

TABLE V :

 V Order of lots processed on M 2 according to rankingCoef f l

	Lot l	CTCoeff l	Step s l	s S l ,l,t	RankingCoeff l
	Lot 5	1.07	Step 5.3	0.68	1.25
	Lot 7	1.43	Step 7.2	0.265	1.43
	Lot 10	0.79	Step 10.2	0.23	2.04
	Lot 4	0.65	Step 4.2	0.2	2.34
	Lot 8	0.48	Step 8.1	0	3.08
	Lot 2	0.45	Step 2.1	0	3.22

TABLE VI :

 VI WIP parameters in the beginning of the second period

	Lot l	Weight	Number	RemP T l	Rem-	Rem-	CT-
		w l	of re-	in	Ref CT l	ExpCT l	Coef f l
			maining	days		in days	in days
			steps S l
	Lot 1	0.33	6	1.1		1.6	4	2.5
	Lot 3	0.5	2	0.25		0.41	0.5	1.22
	Lot 4	0.5	3	0.58		0.76	0.7	0.92
	Lot 5	0.5	4	0.6		0.83	0.5	0.6
	Lot 6	0.33	4	0.75		1.02	4	3.92
	Lot 7	0.5	7	0.76		0.91	0.5	0.55
	Lot 9	0.5	4	0.8		1.05	0.5	0.47
	Lot 10	0.5	2	0.3		0.41	0.5	1.22

TABLE VIII :

 VIII TWT comparison actual versus forecast

	Instance	Absolute deviation	Relative deviation
		(days)	(%)
	Instance 1	228.255	6
	Instance 2	98.305	2.62
	Instance 3	108.86	3.39
	Instance 4	47.77	1.83
	Instance 5	50.13	2.08
	Instance 6	146.23	7.92

ACKNOWLEDGMENT

This work is supported by the ENIAC European Project INTEGRATE. The authors also gratefully acknowledge STMicroelectronics for their support on the knowledge of the semiconductor industry.


Maxime Reynal
email: maxime.reynal@isae.fr

Yvanne Colineaux
email: yvanne.colineaux@aviation-civile.gouv.fr

André Vernay
email: andre.vernay@aviation-civile.gouv.fr

Frédéric Dehais
email: frederic.dehais@isae.fr

Pilot Flying vs. Pilot Monitoring during the approach phase: an eye-tracking study

Keywords: Flight safety, Eye-tracking, Human error, Pilot activity

The adequate monitoring of the flight parameters in the cockpit is a critical issue for flight safety. However, little is known about how the crew supervises the flight deck. In this paper, the preliminary results of a project dedicated to analyze pilot flying and pilot monitoring eyes movements collected in full flight simulator during approach phases are presented. First analyses were conducted over 32 approach phases (8 different crews performing 4 approaches each). The results revealed that the pilot flying and the pilot monitoring exhibited similar ocular behavior during the approach. Moreover, the findings suggested that the pilot monitoring's attentional allocation may not be optimal especially during the short final with low percentage of dwell time on the speed indicator and high percentage of dwell time out of the window.

INTRODUCTION

Operating an aircraft is a complex activity that requires efficient teamwork between the Captain and the First Officer (F/O). During flight operations, the two crewmembers can be alternatively either "Pilot Flying" (PF) or "Pilot Monitoring" (PM). Their roles as PF and PM are defined by the Standard Operating Procedures (SOPs) [START_REF]Roles and Responsibility for PF and PM[END_REF] :

-The PF is responsible for managing the aircraft flight path and gives orders to the PM.

-The PM is responsible for monitoring the current and projected trajectory, the status (e.g. flight path, energy) of the aircraft and all external hazards (e.g. cumulonimbus).

-The PM must inform the PF (or intervenes if necessary) of any deviation of the flight parameters, and executes the PF's orders.

This distribution of roles emphasizes on the adequate monitoring of the cockpit, thus challenging pilots' attentional abilities. However, it is well admitted that attentional abilities are bounded. Recent studies have shown that the occurrence of unexpected events is likely to capture human operator's attention on a single task to the detriment of the supervision of the flight [START_REF] Régis | Formal Detection of Attentional Tunneling in Human Operator-Automation Interactions[END_REF] despite auditory alarms [START_REF] Dehais | Failure to detect critical auditory alerts in the cockpit: Evidence for inattentional deafness[END_REF]. Operational fatigue and drowsiness may lead to the state of mind wandering, leaving the pilot's ill equipped to face complex and challenging situations when a problem arises [START_REF] Durantin | Using near infrared spectroscopy and heart rate variability to detect mental overload[END_REF]. As stated by [START_REF] Casner | Vigilance impossible: Diligence, distraction, and daydreaming all lead to failures in a practical monitoring task[END_REF]: "Real-world monitors may be caught between a continuous vigilance approach that is doomed to fail, a dynamic environment that cannot be fully controlled, and what may be an irresistible urge to let one's thoughts drift." The National Transportation Safety Board (NTSB) and the International Civil Aviation Organization (ICAO) stated that poor monitoring issues were involved in most of major civilian accidents, such as Colgan Air Flight 3407 [START_REF]Loss of Control on Approach Colgan Air, Inc. Operating as Continental Connection Flight[END_REF], Asiana Air Flight 214 [START_REF]Descent Below Visual Glidepath and Impact With Seawall, Asiana Airlines Flight 214 Boeing[END_REF], Turkish Airlines Flight 1951 [START_REF]Dutch Safety Board, Crashed during approach, Boeing 737-800[END_REF], Air Inter Flight 148 or more recently the UPS Airlines Flight 1354 [START_REF]Crash During a Nighttime Nonprecision Instruments Approach Landing UPS Flight 1354 Airbus A300-600[END_REF] to name a few [START_REF]Loss of Control Action Group, Monitoring Matters -Guidance on the Development of Pilot Monitoring Skills[END_REF]. The growing awareness of the need to better understand these events motivated the creation of the Active Pilot Monitoring Group. This group identified several contributive factors such as time pressure, human limitation, poor mental models, automated flight deck issues and a "corporate climate that does not emphasis on monitoring". This work resulted in the "Practical Guide for Improving Flight Path Monitoring" that proposes several countermeasures to crew's poor monitoring [START_REF]The Active Pilot Monitoring Working Group, A Practical Guide for Improving Flight Path Monitoring[END_REF].

However, little is known on how the crewmember actually monitors the flight deck and cross-check the flight parameters especially during critical phases such as the approach [START_REF] Causse | Affective decision making under uncertainty during a plausible aviation task: An fMRI study[END_REF]. To that end, measuring eyes movements with eye-tracking technique offers promising perspective to undercover pilots' voluntary or not attentional strategies. Several studies revealed the suitability of the eye tracking technique for understanding attentional vulnerabilities of pilots interacting with highly automated flight deck [7][17].

In order to better understand why trained pilots fail sometimes to adequately monitor flight parameters, the DGAC/DSAC initiated the Pilot Vision project. This project aimed at analyzing eye tracking data collected by ISAE in different full flight simulators (Airbus-A330 and Boeing-777 full flight simulators) during approach phase preceding a go-around procedure.

In this paper, a first analysis conducted over 32 approaches is detailed. The aim of this study was to address the following questions: 1) In which proportion do the PF and the PM glances at the different areas of interests (AOIs) of the cockpit flight deck during the approach? 2) Does the monitoring in line with the SOPs? 3) Does the pilots' status (Captain versus First Officer) interfere with the pilots' role in the cockpit (PM versus PF)?

Figure 1 -Heatmaps examples of a Captain when he is PM (on top) and a First Officer when he is PF (at the bottom), between 1500 and 500 ft (the two pictures come from an Airbus A330 cockpit). The 'hot' zones where the users focused their gaze with a higher frequency are indicated by yellow/red colors whereas less attended zones are depicted in dark/blue colors.

MATERIAL AND METHOD Participants

Eye-tracker and Areas of Interest

Eye tracking data were collected with two synchronized Pertech eye-trackers (0.25° -0.5° of accuracy). Head movements were corrected by an alignment of three infrared emitters to map participants' fixations on an image of reference (see Figure 1 for a graphical view). The 15 following AOIs were created (see Figure 2 Average percentage of dwell time on each AOI was calculated for each participant, for each approach.

Procedure

Data were collected during the ASAGA project (Aeroplane State Awareness during Go-Around) leaded by the Bureau d'Enquêtes et d'Analyses (BEA), the French safety board. This project was dedicated to the study of PM and PF's eye movements during go-around procedures (GA) in Boeing 777 and Airbus A330 full flight simulators, equitably distributed over the crews (please note that these two aircraft simulators were not equipped with a head-up display) [START_REF]Bureau d'Enquêtes et d'Analyses pour la sécurité de l'aviation civile (BEA), Study of Aeroplane State Awareness during Go-Around[END_REF]. During the experiment, the crews executed four different approach procedures: the first three leading to a GA during the short final and the last one leading to a landing. We unfortunately did not counterbalance the order of the landings across the participants due to implementation constraints in the full flight simulators. In the present paper, we focused our analyses on the four approaches. Each approach was segregated into three phases: 2500 -1500 ft, 1500 -500 ft and 500 ft -TOGA or touch down (this last phase was temporally delimited when the PF put his hands and the thrust levers in order to execute the GA). In the two first scenarios, the Captain was PF and the First Officer was PM. In the last two scenarios, the Captain was PM and the First Officer was PF. approach, with tailwind will increase from 15 to 20 kt. ATC effectively announced this wind conditions change only during short final approach. The crew was expected to perform a GA.

Approach 3) GA due to IMC (Captain = PM, First Officer = PF): Still above Marseille airport, the crew had to make a breakthrough LOC DME 13 left under radar guidance. However, poor visibility led the crew to execute a GA.

Approach 4) and landing (Captain = PM, First Officer = PF): After the third GA, the crew managed to land on Marseille Airport.

RESULTS

Descriptive analysis

The next figures illustrate our descriptive results. We first merged AOIs related to the primary and secondary instruments for both PF and PM (see Figure 3). Then, the average percentage of dwell time on the different AOIs in two different ways was plotted to highlight the possible interactions between pilot's status and role (see Figure 4). Eventually, the last figure illustrates PF and PM average percentage of dwell time during the three segments of the approach.

Statistical analyses Primary vs. secondary instruments

Approaches 1 and 2 were averaged into a first single situation in which Captain was PF and First Officer was PM, and approaches 3 and 4 were averaged into a second single situation in which Captain was PM and First Officer was PF.

AOIs were first merged into two main groups, which are Primary instruments and Secondary instruments (see Figure 3). These two groups were taken as grouped AOIs then compared with External view, NZ and OZ.

A first inferential analysis was performed using Statistica Right MCDU. The statistical analysis revealed a main effect of role [F(1, 14) = 5.04, p < 0.5, η² p = .27], with PM (M = 1.60, SD = 4.01) spending more time on it than PF (M = .53, SD = 1.36, p < .05). A significant main effect of status was found [F(1, 14) = 5.80, p < 0.5, η² p = .29], with Captains (M = .34, SD = 1.13) glancing less at it than First Officers (M = 1.79, SD = 4.02, p < .05). A significant Role × Status interaction was also found [F(1, 14) = 8.99, p < 0.1, η² p = .39]. Post hoc analysis revealed that First Officers glanced more at it when they were PM (M = 3.04, SD = 5.30) than when they were PF (M = .53, SD = 1.16, p < .05). Moreover, First Officers glanced more at right MCDU (M = 3.04, SD = 5.30) than Captains (M = .16, SD = .34, p < .01) where they were both PM.

No statistical result was found for the following instruments: FMA, ND, FCU/MCP, Left MCDU, Flaps panel and External view.

To summarize, these statistical analyses revealed that:

(i) PF fixed more the Primary instruments (Attitude and Altitude indicators, Speed and Heading) than PM (see Table 1);

(ii) PF and PM exhibited almost similar average percentage of dwell time on the Airspeed, Altitude indicator, Heading (i.e. three out of four primary instruments), FMA, ND, FCU/MCP, left MCDU, Flaps panel (the most part of secondary instruments), and moreover External view AOIs;

(iii) The PF glanced more at the Attitude Indicator than the PM;

(iv) The status of the pilot led them to glance differently at primary instruments Airspeed (more fixed by First Officers) and Heading (more fixed by Captains);

(v) There is also a Role × Status interaction on Airspeed: with PM-First Officers glancing more at the Airspeed than PM-Captains.

DISCUSSION

The objective of the present paper was to present the preliminary findings conducted over eight crews performing four stabilized approaches (i.e. 32 analyses).

The main motivation was to investigate PM and PF's ocular behavior (in terms of dwell times) with regard to SOPs. To the authors' knowledge, this study was the first to measure both PM and PF's eye movements in the context of realistic operational scenario, with multi-crew environment in full flight simulator and during approach phases.

Firstly, the results confirmed that both pilots glanced more at the primary flight parameters than the secondary ones. This is consistent with the SOPs as it is mandatory for pilots to particularly monitor the parameters such as the aircraft attitude, the speed, the altitude and the heading. We also note that the rank ordering of the different AOIs is closely correlated with the scan data of PF observed by [START_REF] Sarter | Pilots' monitoring strategies and performance on automated flight decks: An empirical study combining behavioral and eye-tracking data[END_REF], for commercial PF and for pilots [START_REF] Dehais | Embedded Eye-Tracker in a Real Aircraft: New Perspectives on Pilot/Aircraft Interaction Monitoring[END_REF].

Secondly, the results showed that the PF and the PM exhibited the same ocular behavior at the exception of the Attitude indicator. Indeed, the PF glanced more at the Attitude indicator than the PM. Though this latter result is consistent with the SOPs (i.e. the PF is in charge of monitoring the trajectory), one could have expected more pronounced differences between the PF and the PM's ocular behavior. Indeed, these two pilots are supposed to behave differently as defined by the SOPs. We found no effect of the status on the pilots' behavior to the exception of the speed indicator. That is, whatever their status and experience (i.e. Captains are more experienced than First Officers), all PMs and PFs behaved the same way. Only the Captains, when they were PM, glanced more at the airspeed indicator than when they were PF. However, Captains fixated less this critical AOI than First Officers.

Thirdly, results suggested that PMs' visual behavior was not optimal regarding the prioritization of the flight parameters. For instance, PMs spent few time monitoring the speed during the three sequences of the approach. PMs progressively redistributed their attention across the three segments to particularly focus on the external view. This was particularly true during the last sequence (500 ft to GA/touch down) when PMs spend approximately 35% of the time on the external view. This finding is surprising as long as PMs are supposed to keep their head down to monitor critical flight parameters related to the aircraft state. This possible misallocation of attention could find two explanations. The first one is that the pilots can monitor the speed using peripheral vision/covert attention thanks to the "speed trend vector" [START_REF] Dehais | Embedded Eye-Tracker in a Real Aircraft: New Perspectives on Pilot/Aircraft Interaction Monitoring[END_REF] (i.e. an arrow indicating the speed value in the next 10 seconds). This indicator can lead to less eye movements to integrate "trends" of the speed rather than the exact value of the speed itself. A second explanation may rely on automation bias issues [START_REF] Parasuraman | Performance consequences of automation-induced "complacency[END_REF] leading the pilots to over rely on the autothrust to manage the speed, and thus to pay less attention to this parameter. These issues on the speed indicator could provide explanation to recent accidents such as Asiana Air Flight 214 for example, when the pilots failed to identify a critical drop of the speed [START_REF]Descent Below Visual Glidepath and Impact With Seawall, Asiana Airlines Flight 214 Boeing[END_REF].

Eventually a last result concerned the right MCDU: the PMs glanced more at it compared to PFs, and the First Officers fixated more this AOI when they were PM than when they were PF, and they also glanced more than Captains in general. This has to do with the fact that the PM used this user interface to check wind conditions during approach 3.

CONCLUSION, PERSPECTIVES AND LIMITATION OF THE STUDY

We believe that these first results demonstrate the need to conduct eye tracking studies to undercover both PM and PF eye movements during critical phases such as landing. These first findings show that there is a need to establish standards on visual pattern especially for PM and PF to be consistent with SOPs. These eye tracking results are consistent with pilots' training purposes recently recommended by the Federal Aviation Administration (FAA), stated that "by March 2019, air carriers must include specific training pertaining to improve monitoring".

However, this study has several limitations that need to be considered. First, the eye tracking results have to be taken with care. Indeed, our sample was composed of only 16 pilots (eight crews) and the accuracy of eye-tracking techniques still remains a challenge, especially in ecological conditions. Secondly we did not counterbalance properly the order and the design (i.e. PF, PM) of the four approaches. This design was imposed by the ASAGA project and the complexity of using full flight simulators. Thirdly, this study was conducted with French pilots and thus it does not take into account cultural effect. Actually, it's known that several differences exist between countries concerning the conduct of checklists and procedures [START_REF] Merritt | Culture in the cockpit do Hofstede's dimensions replicate?[END_REF].

Eventually, our study will continue with the analyses of the remaining data collected during the ASAGA project (i.e. four more crews -16 approaches), and other data collected during unstabilized approach (4 landings × 28 PFs) in Boeing-737 full flight simulator.

Focus should be placed on the analysis of pilots' ocular dynamics via their scan path, and other eye metrics such as saccades, number of fixations, or time duration between two fixations, for example [START_REF] Dehais | Automation Surprise in Aviation: Real-Time Solutions[END_REF].

  Eight voluntary French crews (i.e., one Captain and one First Officer) coming from different airlines, 16 pilots in total, took part to the experiment. The roles of the Captains and the First Officers (i.e., PM and PF) varied during the flight (see next section). PFs' mean age was 46.5 years (SD = 7.7) with a mean flight experience of 14450 hours (SD = 2192). PMs were on average 41.3 years (SD = 4.3), with a mean flight experience of 7991 hours (SD = 2795).

  ): 1) Airspeed (Speed), 2) Attitude indicator (AI), 3) Altitude indicator (Alt.), 4) Heading (HDG), 5) Flight Mode Annunciator (FMA), 6) Navigation Display (ND), 7) Electronic Centralized Aircraft Monitor (ECAM) for Airbus-A330 and Engine Indicating and Crew Alerting System (EICAS) for Boeing-777, 8) left Multipurpose Control Display Unit (MCDU), 9) right MCDU, 10) External view (Ext.), 11) Auto Flight Control Panel: Flight Control Unit (FCU) for Airbus-330 and Multi Control Panel (MCP) for Boeing-777, 12) Flaps control panel (Flaps), 13) landing gears control panel (Gears), 14 ) No Zone (NZ; i.e.all what is being viewed and which does not correspond to an AOI), and 15) Out of Zone (OZ; i.e. including all the data that was not captured by the device; this is not an AOI but a non-captured quantity of data).

Figure 2 -

 2 Figure 2 -Example of Captain's AOIs delimiting flight deck instruments (Airbus A330).

Approach 1 )

 1 GA requested from ATC (Captain = PF, First Officer = PM): While the crew began an approach on runway 18 at Lyon airport, ATC announced a change in the selected runway then requested an Instrument Landing System (ILS) break in 36 left under radar guidance to occupy the crew during approach. Under 200 ft, runway was occupied and ATC requested for an unexpected GA. This first GA was ordered by ATC and disrupted by a change in the aimed altitude. Approach 2) GA due to tailwind conditions (Captain = PF, First Officer = PM): The crew had to change aircraft course to Marseille airport in accordance with the flight record. They had to execute a standard ILS 31 right Z

Figure 3 -

 3 Figure 3 -Average dwell time for PF (on the left) and PM (on the right), regarding Primary instruments (i.e. Airspeed, Heading, Attitude and Altitude indicators), Secondary instruments (i.e. FMA, ND, ECAM/ECAIS, FCU/MCP, MCDUs, Flaps panel, and landing Gears), External view, Out of Zone (OZ) and No Zone (NZ).

Figure 4 -

 4 Figure 4 -Average dwell time per AOI, per role then per status (on the left), and per status then per role (on the right), for all the four approaches.

Figure 5 -

 5 Figure 5 -Average dwell time per AOI, per phases and roles, over all approaches.

  Tukey's HSD post hoc analysis for Role × Grouped AOIs interaction showed that PF (M = 51.59, SD = 12.11) fixed more at primary instruments than PM (M = 35.74, SD = 7.09, p < .001). No other significant result was found.Considering all AOIs separatelyA second general mixed 14 × 2 × 2 (AOIs [Speed, Attitude indicator, Altitude indicator, Heading, FMA, ND, ECAM/ECAIS, FCU/MCP, left MCDU, right MCDU, Flaps, External view, NZ, OZ] × Role [PF, PM] × Status [Captain, First Officer]) ANOVA with both AOI and Role implemented as within factors and Status implemented as between factor, was performed (see

Table 1 -

 1 This first ANOVA revealed a main effect of the AOIs, [F(4, 56) = 115.00, p < .001, η² p = .89], and a Role × Grouped AOIs interaction [F(4, 56) = 17.00, p < .001, η² p = .55] (see next Table1for a summary of the results).

	General ANOVA (with grouped AOIs)
	Variable	ddl	F	p	η² p
	Grouped AOIs	4, 56	115.00	< .001* .89
	Role	1, 14	1.00	.29	.08
	Status	1, 14	1.00	.44	.04
	Role × Status	1, 14	1.00	.30	.08
	Role × gr. AOIs	4, 56	17.00	< .001*	.55
	Status × gr. AOIs	4, 56	1.00	.43	.07
	Role × gr. AOIs × Status 4, 56	1.00	.60	.05

10. A general mixed 5 × 2 × 2 (Grouped AOIs [Primary instruments, Secondary instruments, Ext., OZ, NZ] × Role [PF, PM] × Status [Captain, First Officer]) analysis of variance (ANOVA) with both Grouped AOIs and Role implemented as within factors and Status implemented as between factor. Results from 5 AOIs (Primary instruments, Secondary instruments, Ext., OZ and NZ) × 2 Roles × 2 Status ANOVA.

Table 2

 2

	General ANOVA
	Variable	ddl	F	p	η² p
	AOIs	13, 169	66.27	< .001*	.84
	Role	1, 13	46.24	< .001*	.78
	Status	1, 13	.14	.72	.01
	Role × Status	1, 13	.36	.56	.03
	Role × AOIs	13, 169	15.23	< .001*	.54
	Status × AOIs	13, 169	1.75	.06	.12
	Role × AOIs × Status	13, 169	1.52	.12	.11

for the results). Data collected for Gears AOI being not balanced because they were missing for one session, the study does not integrate them into the statistical analyses. This ANOVA showed a main effect of the role [F(1, 13) = 46.24, p < .001, η² p = .780]. The analysis revealed a second main effect of the AOIs, [F(13, 169) = 66.27, p < .001, η² p = .84] and a Role × AOIs interaction [F(13, 169) = 15.23, p < .001, η² p = .54].

Table 2 -

 2 Results from 14 AOIs × 2 Roles × 2 Status general ANOVA. < .05, η² p = .31]. Post hoc analyses revealed that Captains gazed more at the airspeed when they were PF (M = 7.35, SD = 6.24) than when they were PM (M = 4.57, SD = 3.07, p < .05). When endorsing the role of PM, First Officer glanced significantly more at the Airspeed indicator (M = 8.53, SD = 5.68) than Captains (M = 4.57, SD = 3.07, p < .05). Finally, no difference on this AOI where found when First Officers acted as PF or PM.

	Attitude and Altitude indicators, Heading, FMA, ND,
	ECAM/ECAIS, FCU, MCDUs, Flaps panel and External
	view).
	For each AIO, a mixed 2 × 2 (Role [PF, PM] × Status
	[Captain, First Officer]) ANOVA was run with Role
	implemented as within factor and Status implemented as
	between factor. Tukey honest significant difference
	(Tukey's HSD) was used for post hoc testing.
	Airspeed. The analysis revealed a main effect of status
	[F(1, 14) = 5.60, p < .05, η² p = .29], with Captains (M =
	5.96, SD = 1.42) glancing significantly less at the Airspeed
	indicator than First Officers (M = 8.38, SD = 2.53). A
	trend was found for the role, [F(1, 14) = 4.03, p = .06, η² p
	= .22], with PF (M = 7.79, SD = 6.17) having high values
	for this AOI than PM (M = 6.55, SD = 4.96). The analysis
	showed a significant Role × Status interaction [F(1, 14) =
	6.14, p

To simplify statistical analyses, twelve different mixed ANOVAs were launched, one per AOI (i.e. Airspeed, Attitude indicator. The statistical analysis revealed a main effect of role [F(1, 14) = 31.42, p < .001, η² p = .69], with PFs spending more time staring at the attitude indicator (M = 30.96, SD = 18.57) than PMs (M = 16.19, SD = 8.14). No effect of the status, nor Role × Status interaction were found.

Altitude indicator. The analysis revealed a significant Role × Status interaction [F(1, 14) = 11.61, p < .01, η² p = .45], but the post hoc analysis revealed no significant differences. No PM/PF difference and no effect of the status were found.

Heading. The analysis revealed a significant main effect of status [F(1, 14) = 10.33, p < .01, η² p = .43], with Captains glancing more (M = 5.19, SD = 2.92) at the Heading indicator than First Officers (M = 1.53, SD = 1.37). The analysis revealed no significant effect between PM and PF, or Role × Status interaction.

ACKNOWLEDGMENTS

We would like to thank Guillaume Adam (BEA), leader of the ASAGA project, and all the pilots who took part in this experiment.