André Rossi
email: andre.rossi@univ-ubs.fr

A Robustness Measure of the Configuration of Multi-Purpose Machines

Keywords: SCHEDULING, PARALLEL MACHINES Multi-purpose machines, Robustness, Scheduling, Configuration

This paper presents new results for assessing the robustness of a configuration for multi-purpose machines. The workshop manager is provided quantitative and meaningful information on how a configuration behaves when disturbances affect the demand. Thus, configurations can be compared on the basis of temporal performance, but also on robustness. The robustness measure of a configuration is returned by assessing the minimum magnitude of disturbances affecting the forecast demand that may lead to breaking the deadline provided by the decision-maker. First, it is shown how to minimize the completion time for a given demand. The demands such that the deadline can be met for a given configuration are characterized afterward, and two robustness measures for a configuration are also provided. The theoretical results are illustrated in detail through an application example.

Multi-purpose machines are often used to model industrial workshops composed of parallel machines subject to constraints which prevent the machines from processing all the product types [START_REF] Brucker | Complexity of scheduling problems with multi-purpose machines[END_REF]. The configuration of such workshops is a key choice as it has a significant impact on performance. However, determining a configuration in a real workshop is a complex issue because it also has an impact on its organization, flexibility and may involves costs in the whole production process.

This paper provides some new tools to assess the temporal performance of a given configuration, as well as two indicators of its robustness. While initially designed to address the problem of configuration in the field of semiconductor industry, the present work provides a theoretical framework that can help the decision-maker to quantitatively evaluate a configuration for any workshop modelled as a set of multi-purpose machines subject to deadline constraints and demand uncertainty.

In the semiconductor industry, the photolithography process is known to be a challenging step as it involves machines that are very expensive and subject to rapid obsolescence due to the constant advances in that field [START_REF] Akçali | Cycle-Time Improvements for Photolithography Process in Semiconductor Manufacturing[END_REF] [START_REF] Yoon | Deadlock-Free Scheduling of Photolithography Equipment in Semiconductor Fabrication[END_REF]. Moreover, the production flow is not easy to manage as products may have complex reentrant routes, and it is often observed that the photolithography workshop is a bottleneck. Finally, as the production processes are constantly updated to integrate new technologies, production volumes are hard to forecast because of unpredictable production losses. In such a context, robustness [START_REF] Kouvelis | Robust Discrete Optimisation and its Applications[END_REF] [4] is a highly desirable feature for a photolithography workshop configuration, this workshop being modelled as a set of multi-purpose machines. Previous works have focused on finding robust configurations using a branch-andbound approach [START_REF] Aubry | Optimisation pour la configuration robuste de systèmes de production de biens et de services[END_REF] and integer linear programming [START_REF] Aubry | Minimizing setup costs for parallel multi-purpose machines under load-balancing constraint[END_REF]. However, these approaches have appeared not to be satisfactory to practitioners. First because they return a single solution to a strongly multi-criteria problem [START_REF] Roy | On operational research and decision aid[END_REF] with some criteria being difficult to model, and second because in such complex contexts, decision-makers do not want a computer program with low situation awareness to make high-stakes decisions for them. That is the reason why the main purpose of this paper is to provide the decision-maker with tools allowing for the measurement of the robustness of a given configuration, by assessing the minimum magnitude of disturbances that may lead to break a given deadline. Then, the decision-maker can perform "What If" studies on several configurations in order to choose the most appropriate ones or integrate these tools into a complete decision making process. It should be stressed that the entire decision process is out of the scope of this article, since the proposed results can be applied to address the configuration to any workshop. Although the proposed approach is not always practicable (for very large of sparse configuration matrices), it has been successfully tested on real-word instances in the field of photolithography (see [START_REF] Rossi | Ordonnancement en milieu incertain, mise en oeuvre d'une démarche robuste[END_REF]).

This paper is organized as follows. The configuration and the underlying scheduling problem are presented in section 2. As solving the underlying scheduling problem may lead to unbalanced production plans, section 3 focuses on the characterization of the demands for which a balanced production plan can be found for a given configuration. It also shows the limits of the proposed approach through computational experiments. Then, these results are used for characterizing the demands that can be processed by a deadline in section 4. This characterization is the basis of the robustness measures presented in section 5, as it allows for the identification of the minimum magnitude of demand disturbances that may lead to breaking the deadline. The paper's contributions are illustrated through an application example that is dealt with in detail in section 6. Finally, the results of this work are discussed in the conclusion section.

2 Workshop configuration and the attached scheduling problem

Workshop configuration

The number of machines is denoted m, and n is the number of product types. The vector N models the demand, Ni being the quantity of products of type i that is to be processed by the workshop. The products have to be processed by one machine, provided the selected machine can process it. For some technological reasons, any machine cannot process any product type. The technological matrix is a boolean n × m matrix T , such that Tij = 1 if and only if the products of type i could possibly be processed by machine j. The processing actually becomes possible if some resources are added to the machine, and if some settings are performed. In that case, the machine j is said to be qualified for the products of type i. The configuration of the workshop is a boolean n × m matrix Q, such that Qij = 1 if and only if the machine j is qualified for the product type i. A valid configuration must be such that Qij ≤ Tij for all i in I = [1, n] and for all j in J = [1, m]. Furthermore, there should be at least one machine qualified for each product type, and each machine should be qualified for at least one product type. In a photolithography workshop, the decision-maker has to design a configuration matrix periodically in order to fit the forecast demand (that is also updated periodically), by ensuring that the underlying scheduling problem has a solution which meets a given deadline, even when the actual demand differs from the forecast demand.

Underlying scheduling problem

The underlying scheduling problem involves the speed matrix, that is an n × m matrix V such that Vij is the number of products of type i that machine j can process per unit of time (provided that Qij = 1). The machines are supposed to be uniform [START_REF] Brucker | Scheduling Algorithms[END_REF], i.e. V can be written as V = V P × V M , where V P is a strictly positive n × 1 column vector, and V M is a strictly positive 1 × m row vector.

Indeed, machine j has a processing speed V Mj , that is the amount of products that machine j can process per unit of time. However, processing one unit of a product of type i may take longer than processing one unit of a product of type i = i on the same machine. Thus, V Pi is defined as the speed factor for the products of type i. By convention, V P1 = 1: the product type 1 serves as a reference for expressing the speed factor of the other product types. More precisely, processing one unit of a product of type i can be achieved V Pi times faster than processing the same quantity of a product of type 1 on the same machine. Finally, machine j can process Vij = V Pi × V Mj units of products of type i per unit of time for all (i, j) in I × J. By construction, matrix V is a rank-one matrix: this characterizes uniform machines. Such a model is appropriate for most photolithography workshops.

Because the machines are very expensive and subject to rapid obsolescence, it makes sense to schedule the products to process so as to minimize the maximum completion time (this criterion is usually referred to as the makespan and denoted Cmax [START_REF] Pinedo | Planning and Scheduling in Manufacturing and Services[END_REF]). The solution for such a problem is a value for Rij for all (i, j) in I × J, where R is a production plan. The numerical value for Cmax is the problem objective, and can be deduced from R. More precisely, Rij is equal to the time spent by machine j processing products of type i. Preemption and splitting are allowed in that problem, that can be referred to as QM P M |split|Cmax using the three-field notation [START_REF] Graham | Optimization and approximation in deterministic sequencing and scheduling: a survey[END_REF]. This problem is polynomially solvable [START_REF] Aubry | Minimizing setup costs for parallel multi-purpose machines under load-balancing constraint[END_REF] and can be modelled as a linear program denoted LP .

(LP ) :

             Minimize Cmax (1) j∈J Qij VijRij = Ni (∀i ∈ I) (2) i∈I QijRij ≤ Cmax (∀j ∈ J) (3) Rij ≥ 0 (∀i ∈ I)(∀j ∈ J) (4)
(1) is the objective function. The set of constraints (2) enforces that the demand is met for each product type. The set of constraints (3) enforces that the time workload of each machine is less than Cmax. The set of constraints (4) ensures that the production plan is valid, as an amount of time cannot be negative. It should be stressed that Rij are the decision variables while the makespan Cmax is the problem objective. The configuration Q, the speed matrix V and the demand N are given.

LP makes possible the comparison between two configurations on the basis of temporal performance. Let Cmax(Q, N ) be the optimal value for Cmax after solving LP for configuration Q and demand N . Then, configuration Q1 performs better than configuration Q2 on demand N if and only if Cmax(Q1, N ) is smaller than Cmax(Q2, N ).

Balanced production plans

This section focuses on the characterization of the demands that are such that a balanced production plan can be found. Solving LP to optimality leads to one of the two following situations. The production plan R is either balanced, or not. More precisely, R is said to be balanced if all the machines share the same workload (i.e.

i∈I QijRij = Cmax for all j in J). In the rest of this paper, the workload of the machine j is denoted Lj , and C(N ), defined below, is a particular workload that will be proved to be equal to Lj for all j when the production plan is balanced (in [START_REF] Aubry | Minimizing setup costs for parallel multi-purpose machines under load-balancing constraint[END_REF] LP is shown to be a transportation problem [START_REF] Chvátal | Linear Programming[END_REF] when a balanced production plan exists).

C(N ) = i∈I Ni V Pi j∈J V Mj
Let us consider the following example with m = 2 machines and n = 3 product types. The machines are such that Vij = 1 for all i and j (identical machines), the configuration matrix Q and the demands N 1 and N 2 defined below are considered.

V = 1 1 1 1 1 1 , Q = 0 1 1 1 1 0 , N 1 = 3 7 2 , N 2 = 7 3 2
Solving LP to optimality for N 1 (respectively for N 2 ) leads to the production plans R 1 (respectively R 2 ). It should be stressed that the solutions of LP are generally not unique.

R 1 = 0 3 4 3 2 0 , C 1 max = C(N ) = 6, R 2 = 0 7 3 0 2 0 , C 2 max = 7
First it can be seen that R 1 is balanced: the two machines have the same workload, that is equal to C(N ). Second, solving LP to optimality for N 2 leads to a non-balanced production plan R 2 with a greater makespan. This small example shows that because of non-qualifications, it is sometimes impossible to get a balanced production plan even when LP is solved to optimality. Roughly speaking, it can be said that the configuration does not fit the demand when no balanced production plan can be found, because the completion time is strictly greater than C(N ) while it could be equal to C(N ) for another configuration (the one defined by Qij = 1 for all i and j for instance).

The aim of this section is to characterize the demands that are such that a balanced production plan can be found for a given configuration. This is useful for assessing the robustness of a configuration as shown in section 5. To do so, the following lemma is used. Lemma 1. Let R be an admissible solution to LP . Then,

• There exists jc in J such that Lj c ≥ C(N )

• There exists ju in J such that Lj u ≤ C(N )

Proof. The set of constraints (2) and (3) in LP can be rearranged as:

       j∈J Qij RijV Mj = Ni V Pi (∀i ∈ I) (2) i∈I Qij Rij = Lj ≤ Cmax (∀j ∈ J) (3)
Multiplying each inequality of (3) by V Mj and summing them up yields to:

j∈J i∈I Qij RijV Mj = j∈J V Mj Lj
Using (2), the last equality can be written as:

i∈I Ni V Pi = j∈J V MjLj
C(N ) = j∈J V Mj Lj j∈J V Mj (A)
Lemma 1 is now shown by contradiction using equation (A). First, it is assumed that Lj < C(N ) for all j in J. Then, Lj = C(N )αj where αj > 0 for all j. Thus, equation (A) is equivalent to:

j∈J αj = 0
This is a contradiction, as αj are strictly positive numbers for all j, and so, it is proved that there exists jc in J such that Lj c ≥ C(N ).

Second, it is assumed that Lj > C(N ) for all j in J. Then, Lj = C(N ) + αj where αj > 0 for all j. Equation (A) again leads to a contradiction, proving that there exists ju in J such that Lj u ≤ C(N ).

It can easily be shown that the configuration Q * defined by Q * ij = 1 for all (i, j) in I × J is such that any demand N can be processed in C(N ) units of time, and is associated with a balanced production plan. As a matter of fact, non qualifications (i.e. zero-entries in the configuration matrix Q) are responsible for non-balanced production plans. Given a configuration Q, the set of demands for which a balanced production plan can be found is characterized by a set of inequalities, as shown in the theorem 3.1 in the next section.

Characterization

If no balanced production plan exists for N , critical machines are defined as follows. A machine is said to be critical if its workload is equal to Cmax for any optimal production plan. For a given optimal production plan, if a machine happens to have a load equal to Cmax, it does not necessarily mean that it is critical, as it may exist a different production plan (which is also an optimal solution to LP ) for which this machine has a load that is strictly less than Cmax. However if a machine has a load which is strictly less than Cmax, then this machine is not critical. It should be mentioned that identifying the critical machines does not require the computation of all the optimal production plans. It can be achieved from any single optimal solution to LP , by using a polynomial algorithm that can be found in [START_REF] Rossi | Ordonnancement en milieu incertain, mise en oeuvre d'une démarche robuste[END_REF]. Furthermore, it can be noted that any product type processed by a critical machine cannot be processed by a non-critical machine. If it was the case, a critical machine would be able to transfer a non-zero quantity of product to a non-critical one, leading to a new production plan for which a critical machine would actually be non-critical which is a contradiction. For the same reason, any product type that can be processed either by critical or by non-critical machines has to be exclusively processed by non-critical machines for any production plan optimal for LP .

Example of critical machines

The following example aims at illustrating the definition of critical machines. The workshop configuration Q, the speed matrix V = V P × V M and the forecast demand N are given hereafter.

Q =    1 0 0 0 1 0 0 0 1 0 1 1    , V P =    1 1 1 1    , V M = 1 1 1 , N =    10 5 1 5

  

The following production plans R 1 and R 2 are two optimal solutions for LP , for which the makespan is Cmax = 10.

R 1 =    10 0 0 0 5 0 0 0 1 0 5 0    , R 2 =    10 0 0 0 5 0 0 0 1 0 0 5   
Production plan R 1 shows that machine 3 cannot be critical, since its load is equal to 1 (and so is strictly less than Cmax = 10). Besides, R 2 also shows that machine 2 is not critical, since there exists an optimal solution to LP (namely R 2 ) for which machine 2 has a load which is strictly less than Cmax = 10. It can be proved that machine 1 is critical because it is the only machine qualified for a product type (namely product type 1) that requires 10 units of time to be processed. As a consequence, there does not exist any other production plan which is an optimal solution to LP , and such that machine 1 has a load that is strictly less than 10. Theorem 3.1. For a given configuration Q, there exists a balanced production plan for the demand N if and only if the following set of constraints is satisfied.

(C k ) : i∈I k Ni V Pi j∈J k V Mj ≤ i∈ Īk Ni V Pi j∈ Jk V Mj (∀k ∈ K)
In the above set of constraints, I k is a non-empty strict subset of I and Īk is the complement of I k in I. Analogously, J k is a non-empty strict subset of J and Jk is the complement of J k in J. For all k in K, I k and J k are such that Qij = 0 for all (i, j) in I k × J k . The Cartesian product I k × J k is assumed to be a maximum rectangle of zeros in Q. A rectangle of zeros is said to be maximum if it is not included in another one, i.e. there does not exist i ∈ Īk such that (i ∪ I k ) × J k is a rectangle of zeros, and there does not exist j ∈ Jk such that

I k × (j ∪ J k ) is a rectangle of zeros.
All the maximum rectangles of zeros (denoted MRZ for short) in Q are computed using an algorithm shown in section 3.4, derived from Nourine and Raynaud algorithm [START_REF] Nourine | A fast algorithm for building lattices[END_REF]. They are numbered from 1 to |K|, where |K| can be up to 2 n -2 for m = n when Q is the identity matrix. Although such an enumerative algorithm is obviously non-polynomial [START_REF] Garey | Computers And Intractability, A Guide To The Theory Of NP-Completness[END_REF], it can be used for some real-life instances. Section 3.4 shows the limits the proposed approach.

The theorem 3.1 is shown below.

Proof. First, it is shown that if there exists a balanced production plan for N , then constraint (C k ) is satisfied for all k in K. This implication is proved by contraposition, i.e. it is shown that if there exists k in K such that constraint (C k ) is violated, then no balanced production plan can be found for N .

As I k × J k is a maximum rectangle of zeros, all the products of type i in I k are processed by some machines in Jk . These machines process at least all the products of type i in I k , as a consequence, lemma 1 ensures that there exists jc in Jk such that:

Lj c ≥ i∈I k Ni V Pi j∈ Jk V Mj I k × J k
is a maximum rectangle of zeros, which also means that the machines in J k only process products of type i in Īk . Then, lemma 1 ensures that there exists ju in J k such that:

Lj u ≤ i∈ Īk Ni V Pi j∈J k V Mj
Since constraint (C k ) is violated, it can be written that:

Lj u ≤ i∈ Īk Ni V Pi j∈J k V Mj < i∈I k Ni V Pi j∈ Jk V Mj ≤ Lj c
This clearly shows that no balanced production plan can be found for the demand N . Second, it is shown that if constraint (C k ) is satisfied for all k in K, then there exists a balanced production plan for N .

This implication is proved by contraposition, i.e. it is shown that if an optimal production plan for N is not balanced, then there exists k in K such that (C k ) is violated.

It is assumed that R is a non-balanced production plan for N that is also an optimal solution to LP . By definition, critical machines are processing products that non-critical machines cannot process (because they are not qualified for these product types). This means that I k × J k is a maximum rectangle of zeros in Q, where Jk is the set of the critical machines, I k is the set of product types that are processed by these machines, and J k is the set of the non-critical machines. Consequently, critical machines share the same load (this common load sets the value for Cmax)

= i∈I k Ni V Pi j∈ Jk V Mj = Cmax (∀j ∈ Jk )
Since there does not exist any balanced production plan for N , it can be written that Lj ≤ Cmax for any non-critical machine j. For any j in J k , multiplying Lj by V Mj and summing up all these inequalities yields to:

j∈J k V Mj Lj < j∈J k V Mj Cmax (B)
The last inequality is strict because, as no balanced production plan can be found for N , there exists (at least) one machine ju in J k , the load of which is strictly less than Cmax.

The machines j in J k are not qualified for the products of type I k . Consequently, Qij Rij = 0 for all i in I k . By definition of Lj , it can be written that:

Lj = i∈I QijRij = i∈ Īk QijRij (∀j ∈ J k )
The inequality (B) can then be written as follows:

i∈ Īk j∈J k QijRij V Mj < Cmax j∈J k V Mj
And so, the quantity

j∈J k QijRij V Mj can be replaced by Ni V Pi
for all i, as specified by the set of constraints (2) in LP .

Finally, replacing Cmax by its value leads to:

i∈ Īk Ni V Pi j∈ Jk V Mj < i∈I k Ni V Pi j∈J k V Mj
This clearly shows that (C k ) is violated.

Expressing C max for unbalanced production plans

Knowing that Cmax = C(N ) for any balanced production plan, this section shows that the makespan value can be computed analogously when the demand is such that any optimal solution to LP is not balanced. Theorem 3.1 shows that if an optimal solution to LP is not balanced, then there exists a maximum rectangle of zeros in Q denoted I0 × J0 such that the products of type in I0 are processed by critical machines only. This maximum rectangle of zeros can be found by using a polynomial algorithm (see [START_REF] Rossi | Ordonnancement en milieu incertain, mise en oeuvre d'une démarche robuste[END_REF]). Furthermore the machines in J0 only process the product types in I0 and these product types are only processed by the machines in J0.

The value of Cmax only depends on the critical machines that share the same workload on the product types they are processing. This situation can be modelled as a subproblem of the original problem, where the considered product types are in I0 and where the only considered machines are those in J0. This subproblem turns out to have a balanced production plan as the critical machines share the same workload. Thus, the makespan expression for a non-balanced production plan is:

Cmax = i∈I 0 Ni V Pi

Proposed algorithm

The computation of all the MRZ for a given configuration is a necessary step for characterizing the demands that can be processed with a balanced production plan. It can be achieved using Algorithm 1, that is derived from Nourine and Raynaud algorithm [START_REF] Nourine | A fast algorithm for building lattices[END_REF].

input : Q, a n × m configuration matrix output: I k and J k the h MRZ defined by I k × J k for all k ∈ [1, h] Initialization h = 1 is the number of MRZ currently found I0 = {1, 2, . . . , n} J0 = {∅} for b = 1 to n do A ← Set of column indexes of zeros entries located in the row b of matrix Q for k = 0 to h -1 do Building a candidate r = I × J I = A ∩ I k J = J k ∪ {b} Finding out if there exists q ∈ [0, h -1] such that Iq = I q = 0 stop = 0 while stop == 0 do if I == Iq then Jq = Jq ∪ J stop = 1 else q = q + 1 end if q ≥ h then h = h + 1 I h = I J h = J stop = 1 end end end end Delete I0 and J0

Algorithm 1: Enumerating the maximum rectangles of zeros in a configuration

This algorithm takes configuration Q as input, and aims at building all the MRZ in Q. The MRZ are found iteratively (it is an enumerative algorithm), and h is the number of MRZ currently found. Initially, h is set to one and a dummy MRZ denoted I0 × J0 is built. This MRZ does not exist in practice (J0 is the empty set) and will be deleted at the very end of the algorithm. Its purpose is to let actual MRZ being built (this step is called Building a candidate r = I × J ). The rows of Q serve as generators: element A (that is the set of indexes of zero elements in the row b of Q) is used to enlarge any preexisting MRZ r k by adding machine b to J k , and computing the intersection of I k and A. Thus, the dummy MRZ plays the role of a "neutral element" allowing to build candidates from the rows of Q. The resulting rectangle of zeros r = I × J (the so-called candidate) is looked for in the list of existing MRZ. If I already exists, the preexisting MRZ Iq × Jq is enlarged by including the machines of J in Jq. If I cannot be found in the preexisting MRZ, the list grows longer as the candidate becomes a new MRZ, leading h to be incremented by one. When the algorithm terminates, the dummy MRZ I0 × J0 is deleted and h = |K|.

It must be stressed that Algorithm 1 is a simplified version of the C program actually implemented, as many technical improvements can be added to decrease computation time. The purpose of this algorithm is to provide the reader with an overview of the necessary steps to compute all the MRZ. Thus, it can be observed that since the list of rectangles of zeros has to be searched for each new candidate, finding a new MRZ takes more and more time as the list grows longer. If each candidate is to be added to the list, the number of MRZ can be up to 2 n -2 [START_REF] Rossi | Ordonnancement en milieu incertain, mise en oeuvre d'une démarche robuste[END_REF]. Furthermore, as Nourine and Raynaud [START_REF] Nourine | A fast algorithm for building lattices[END_REF] say that finding h (i.e. the total number of MRZ for a given configuration) without enumerating them all is an open problem, it is not easy to forecast the algorithm running time. The next section reports computational experiments allowing to characterize the configurations that can be searched for MRZ in a reasonable amount of time.

Computational results

Since Algorithm 1 may return a non-polynomial number of maximum rectangles of zeros, some experiments have been conducted for determining the limits of the approach. The tested instances are randomly generated configurations characterized by the number of product types n, the number of machines m and the matrix density d, that is the probability for every single entry Qij to be set to one (i.e. a matrix with a low density is sparse). Algorithm 1 is run on twenty instances randomly generated for each triplet (n, m, d). Tables 3.4 n is greater than m without loss of generality because matrix Q and its transpose Q T have the same rectangles of zeros (I k × J k has only to be replaced by J k × I k ). Computational results for instances with n < 15 are not reported because their solving time is very close to zero whatever configurations density. As expected, it can be observed that computation time increases significantly with instance size. Thus, Table 3.4.2 displays results only up to n = 24 because solving larger instances requires more than one hour of computation time. However, this only highlights a single aspect of the limits of the proposed approach as Tables 3.4.2 and 3.4.2 show that instance size is far from being the most important parameter impacting computation time. Density turns out to have an even more significant impact as the average computation time for instances of size 24 × 48 is around half and hour for d = 0.25 (sparse configurations) but only 0.03 s for d = 0.75 (dense configurations). This is due to the fact that there are far more maximum rectangles of zeros to be found in sparse configurations than in dense ones (about 192,000 for d = 0.25 and around 400 for d = 0.75 in the case of 24 × 48 instances). Finally, density and instance size (in a lesser extend) appear to be the most relevant parameters when considering running Algorithm 1 on a given configuration. The proposed approach is practicable regardless of density for instances having less than 20 product types, and up to 30 product types for dense and mid-dense configurations. This approach turned out to be applicable in the context of photolithography workshops, as the number of product types does not exceed 15 as far as the author knows.

Characterization of deadline meeting

In this section, the demands that can be processed by a deadline for a given configuration are shown to be characterized by a set of inequalities based on the maximum rectangles of zeros. This result provides an efficient way of assessing the robustness of a configuration as the minimum magnitude of disturbances on the forecast demand. In the rest of this paper, τ denotes a deadline that is common to the whole demand. Such an assumption on the deadline is realistic because the demand is often a set of products (i.e. a lot) that is expected to leave the workshop at a given date. In other words, the purpose of this section is to find the set of demands N such that the optimal value for Cmax (returned by solving LP ) is less than or equal to τ .

       i∈I Ni V Pi ≤ τ j∈J V Mj (E0) i∈I k Ni V Pi ≤ τ j∈ Jk V Mj (E k ) (∀k ∈ K)
Proof. This theorem is proved by showing that the three following implications (1) ⇒ (2), ( 2) ⇒ (3) and ( 3) ⇒ (1) hold. First, it is shown that (1) ⇒ (2). It is assumed that R, a production plan for N meets the deadline τ . If R is balanced, then N is set to N and the property (2) is true. If R is not balanced, then every machine having a load that is strictly less than τ can be loaded up to τ by increasing the quantity of products of any type for which the machine is qualified. By proceeding this way, a new demand N satisfying N ≥ N and the associated balanced production plan meeting the deadline τ are built. Now, it is shown that (2) ⇒ (3). Demand N is supposed to be such that its corresponding production plan is balanced, and that it can be processed in τ units of time or less. According to the results presented in the previous section, it can be written that:

i∈I N i V Pi j∈ J V Mj ≤ τ
As N ≤ N the above inequality also holds for N . Thus, the constraint denoted (E0) in S is satisfied. Furthermore, as the production plan for N is balanced, constraint (C k ) holds for all k:

(C k ) : i∈I k N i V Pi j∈J k V Mj ≤ i∈ Īk N i V Pi j∈ Jk V Mj (∀k ∈ K)
Adding the quantity

i∈I k N i V Pi j∈ Jk
V Mj to both sides of the above inequalities leads to:

i∈I k N i V Pi j∈J V Mj ≤ i∈I N i V Pi j∈ Jk V Mj (∀k ∈ K)
The right-hand side can be upper-bounded by

τ j∈J V Mj j∈ Jk V Mj because it has been shown that i∈I N i V Pi
is less than or equal to τ j∈J V Mj , and then, by dividing the whole inequality by j∈J V Mj (which is a strictly positive quantity), it yields to the new inequality:

i∈I k N i V Pi ≤ τ j∈ Jk V Mj (∀k ∈ K)
Then, as N ≤ N , the above inequality also holds for N , which is the expression of constraint (E k ) for all k in K.

Finally, it is shown that (3) ⇒ (1). This implication is shown by contraposition, i.e. it is proved that N cannot be processed in less than τ units of time which implies that either constraint (E0) is violated, or there exists k in K such that constraint (E k ) is violated.

Thus, it is now assumed that any optimal production plan for N is such that Cmax > τ . If N is such that its production plan is balanced, then it can be deduced that:

Cmax = i∈I Ni V Pi j∈J V Mj > τ
This clearly shows that constraint (E0) is violated.

If N is such that its production plan is not balanced, then it has been shown in section 3.3 that there exists a maximum rectangle of zeros (i.e. an integer k in K such that I k × J k is a maximum rectangle of zeros) such that: Even in the field of scheduling, the term robustness may have different meanings as shown in [START_REF] Billaut | Flexibility and Robustness in Scheduling[END_REF] and in [START_REF] Kouvelis | Robust Discrete Optimisation and its Applications[END_REF]. Whereas its general meaning often refers to the ability to face disturbances or uncertainties, the application-dependent character of this term may not be surprising as the so-called disturbances and uncertainties may be very different from one problem to another, as well as the available means to face them, and the solution features that are relevant for protection from disturbances.

Cmax = i∈I k Ni V Pi j∈ Jk V Mj > τ Thus constraint (E k ) is violated.
In this paper, deadline satisfaction is assumed to be a key issue. This makes sense in the context of semiconductor manufacturing because the photolithography workshop is often a bottleneck. As a consequence, the robustness of a configuration is defined as its ability to maintain the deadline satisfaction despite the disturbances that may affect the forecast demand. The present section aims at assessing the minimum magnitude of disturbances that may lead to breaking the deadline. The higher this magnitude is, the more robust the configuration is.

The characterization of the demands that are such that the configuration can meet the deadline is a useful tool to assess the robustness because constraints (E0) and (E k ) for all k in K define the borders that should be watched to check if disturbances may compromise the ability to meet the deadline. To achieve this goal, it is necessary to evaluate the distance from any demand meeting the deadline to the borders of the characterized set. The distance is evaluated using the L1-norm because |N -N |1 returns the additional load between demands N and N when N ≥ N .

Distance from a demand to the constraints' frontiers

Distance to the frontier of (E 0 )

Demand N is on the frontier of constraint (E0) if and only if N satisfies:

i∈I Ni V Pi = τ j∈J V Mj
Where Q, V P , V M and τ are given.

The following lemma provides a handy formula for assessing the distance from any demand N to the frontier of constraint (E0). This is not a distance definition since the classical L1-norm distance is used. Lemma 2. The distance from any demand N to the frontier of constraint (E0) is denoted d(N, E0) and has the following expression:

d(N, E0) = τ j∈J V Mj - i∈I Ni V Pi V Pmin
Where i0 is defined by

V Pmin = min i∈I {V Pi}
The proof is not provided here for the sake of concision, but it is analogous to the proof used to show a very close lemma that can be found in [START_REF] Rossi | Ordonnancement en milieu incertain, mise en oeuvre d'une démarche robuste[END_REF].

Distance to the frontier of (E k )

For all k in K, demand N is on the frontier of constraint (E k ) if and only if N satisfies:

i∈I k Ni V Pi = τ j∈ Jk V Mj
It is assumed that Q, V P , V M , k and τ are given. The distance from any demand N to the frontier of the constraint (E k ) is denoted d(N, E k ) and has the following expression:

d(N, E k ) = τ j∈ Jk V Mj - i∈I k Ni V Pi V Pi k (∀k ∈ K) Where i k is defined by V Pi k = min i∈I k {V Pi}
The proof is not provided, but it is also derived from [START_REF] Rossi | Ordonnancement en milieu incertain, mise en oeuvre d'une démarche robuste[END_REF]. It is recalled that as the maximum rectangles of zeros cannot always be computed in a reasonable amount of time, assessing the distance from N to the frontier of all (E k ) may not always be practicable as shown in section 3.

Robustness measures

Because most workshops evolve in an uncertain environment, it is now considered that the actual demand Ñ that the workshop must process is not exactly the forecast demand N * . As this change may affect the workshop ability to meet the deadline τ , it is relevant to assess the robustness of the configuration, i.e. its ability to guaranty that the workshop can meet the deadline even when Ñ is different from N * .

More formally, the actual demand Ñ can be written as Ñ = N * + ∆N + -∆N -where ∆N + and ∆N -are nonnegative vectors satisfying ∆N + i × ∆N - i = 0 for all i in I. This means that ∆N + i and ∆N - i cannot both be strictly positive. ∆N + i is the additional quantity of product type i that is found in the actual demand Ñ and ∆N - i is the missing quantity of product type i in Ñ . As an example, N * 1 = 60 and Ñ1 = 72 lead to ∆N - 1 = 0 and ∆N + 1 = 12 because there are 12 more products of type 1 in Ñ than expected in N * . Thus the equality N * 1 = Ñ1 + ∆N + 1 -∆N - 1 holds. It would also hold with ∆N - 1 = 1 and ∆N + 1 = 13 but this set of values is not valid as ∆N + 1 × ∆N - 1 = 0. In the rest of this paper, it is assumed that the workshop can meet the deadline τ for the forecast demand N * . In other words, τ is supposed to be greater than or equal to the optimal Cmax value returned after solving LP for N * . It can easily be shown that if Ñ ≤ N * , the workshop ability to meet the deadline is not put in question. As a consequence, additional quantities of products are only taken into consideration when examining the differences between N * and Ñ . Thus, the robustness r(Q, N * , τ ) of a configuration Q is the minimum quantity of products added to the forecast demand N * (regardless of product types) that may lead to missing the deadline τ . More formally, r(Q, N * , τ ) is the minimum distance from N * to the frontier of (E0) and of (E k ) for all k:

r(Q, N * , τ ) = min d(N * , E0), min k∈K d(N * , E k )
Thus, the robustness is defined as a distance. In that sense, this definition is very close to the Sotskov stability radius [START_REF] Sotskov | On the Calculation of the Stability Radius of an Optimal or an Approximate Schedule[END_REF], however robustness is focused on additive perturbations only, as fewer products than expected cannot put in question the configuration ability to meet the deadline. Since the robustness measure is based on the maximum rectangles of zeros, its use may be unpracticable for some instances as shown in section 3.4.

Since Q is assumed to be such that a production plan meeting the deadline can be found for N * , r(Q, N * , τ ) is nonnegative. If a configuration does not satisfies this assumption, it makes no sense to assess its robustness because it is not even suitable for the forecast demand. The above expression clearly shows that the robustness value is either set by (E0), or by one or more constraints (E k ). These two cases have very different significations, as constraint (E0) does not depend on the configuration, while (E k ) is derived from the maximum rectangles of zeros in Q.

This distinction is meaningful as it provides the decision-maker with a key piece of information. If r(Q, N * , τ ) = d(N * , E0), then modifying the configuration cannot lead to increase robustness. The full workshop potential for facing demand uncertainty is already used. However, if there exists k in K such that d(N * , E k ) < d(N * , E0), then the configuration limits the robustness that the machines may offer. As a matter of fact, the rectangle of zeros I k × J k prevents the configuration from offering more robustness. Thus, the decision-maker can consider qualifying a machine in J k for a product type in I k for increasing the configuration robustness. It may also be relevant to let him or her know the current configuration robustness potential ρ(Q, N * , τ ) defined as follows:

ρ(Q, N * , τ ) = r(Q, N * , τ ) d(N * , E0)
It must be stressed that if d(N * , E0) is zero, then ρ(Q, N * , τ ) is set to one by convention. Indeed, d(N * , E0) = 0 means that N * is in the frontier of (E0). The robustness potential is zero, but no other configuration could improve that result.

The robustness potential ρ(Q, N * , τ ) is in [0, 1], and can be viewed as a fitness indicator. The higher it is, the more suitable the configuration is for facing uncertainties in the neighborhood of the forecast demand. ρ(Q, N * , τ ) = 1 means that the configuration takes advantage of the full workshop potential for facing disturbances.

Application example

All the theoretical framework presented in this paper is illustrated in this section through a small instance. The considered workshop is made of m = 4 machines for processing n = 5 product types. The technological matrix is supposed to be the n-by-m all-ones matrix (i.e. machines are not subject to technological constraints). The workshop configuration Q, the speed matrix V = V P × V M and the forecast demand N * are given hereafter. These results have been computed using Matlab R with the linprog command, but without the Simplex option. Thus, Matlab R uses an interior point based method [START_REF] Karmarkar | A polynomial-time algorithm for linear programming[END_REF], that is expected to run faster than the Simplex method on large instances. However, since solvers have a finite precision, numerical errors may affect the results. It is the case here, where it has been deliberately decided to display the results with a decimal precision limited to three digits. Such an imprecision is not a problem because determining critical machines allows to compute the Cmax value symbolically (with an infinite precision). It can be observed that the production plan is not balanced because machine 4 is underloaded whereas the other machines appear to be critical. Then, according to the results shown in section 3.3, the subproblem involving the product types that are not processed by machine 4 (i.e. I0 = {1, 2, 3, 5}) and all the machines but machine 4 (i.e. J0 = {1, 2, 3}) has a balanced production plan, and a completion time having a value identical to those of the original problem.

Q =      1 1 0 0 0 1 1 0 0 0 1 0 1 0 0 1 1 1 1 0      , V P =      1 3 3 5 2      , V M = 2 1 1 3 , N * =     
Cmax = i∈I 0 N * i V Pi j∈ J0 V Mj = 135 1 + 22 3 + 79 3 + 43 2 2 + 1 + 1 = 1141 24 ≈ 47.542
The Cmax value returned by the solver cannot be exact as 1141 is not divisible by 3. However, computational approximations do not prevent one from accessing the symbolical value for Cmax, provided that critical machines can be identified.

If the new configuration Q defined below is used instead of Q, solving LP leads to the following results for N * : It can be observed in R that the machines 2 and 3 are critical. This is due to the fact that in the configuration Q , the product types 1 and 3 are only processable on these machines. Consequently, it can be deduced that C max = 242 3 ≈ 80.667. This clearly shows that the configuration Q performs better than Q on the forecast demand N * as Cmax < C max . Since Q has 11 qualifications (i.e. one-entries) and Q has only 10, this result also shows that temporal performance is not directly related to the number of qualifications in the configuration matrix.

Q =      0 1 1 0 1 1 1 1 0 0 1 0 1 0 0 1 0 1 0 1      , R =      0

Robustness assessment

The set of demands that can be processed in less than τ units of time with configuration Q can now be determined. To do so, theorem 4.1 is applied. As it requires the computation of the maximum rectangles of zeros in Q, Algorithm 1 is run, returning the seven following maximum rectangles of zeros:

r1 = {2, 3} × {1, 4} r2 = {3, 4} × {2} r3 = {1, 4} × {3} r4 = {1, 2, 3, 5} × {4} r5 = {1} × {3, 4} r6 = {3} × {1, 2, 4} r7 = {4} × {2, 3}
V P2 + N3 V P3 ≤ τ V M2 + V M3 (E2) : N3 V P3 + N4 V P4 ≤ τ V M1 + V M3 + V M4 (E3) : N1 V P1 + N4 V P4 ≤ τ V M1 + V M2 + V M4 (E4) : N1 V P1 + N2 V P2 + N3 V P3 + N5 V P5 ≤ τ V M1 + V M2 + V M3 (E5) : N1 V P1 ≤ τ V M1 + V M2 (E6) : N3 V P3 ≤ τ V M3 (E7) : N4 V P4 ≤ τ V M1 + V M4
Besides constraints (E k ) for all k in K depending on the configuration Q, constraint (E0) must also be enforced:

(E0) : N1 V P1 + N2 V P2 + N3 V P3 + N4 V P4 + N5 V P5 ≤ τ V M1 + V M2 + V M3 + V M4
The set of demands that can be processed by τ units of time is denoted S in the statement of theorem 4.1 (it is recalled that S is made of constraints (E0) and (E k ) for all k in K). By construction, S is a bounded intersection of a finite set of half spaces. Consequently, S is a polyhedron. As inequality (E2) is equal to the sum of inequalities (E6) and (E7), it can be deduced that the proposed description for S is not minimal because some constraints may be redundant (where redundant should be understood as not facet-defining). Some techniques to get a minimal description of a polyhedron can be found in [START_REF] Wolsey | Integer Programming[END_REF]. Adapting these techniques to the present work may lead to saving memory space when storing the maximum rectangles of zeros in a configuration. However, it would be far more desirable to embed these techniques into the generation of maximum rectangles of zeros algorithm, as many constraints appear to be redundant in some particular cases. Such an improvement would be both time and space saving.

In the present example, it can be checked that the forecast demand N * does not satisfy all these constraints for τ = 45, as it cannot be processed in less than 47.542 units of time: (E4) is not satisfied for τ = 45.

However, constraints (E0) and (E k ) for all k in K are satisfied for τ = 50. It that case, the distances from N * to the frontier of these constraints are given below: It can be observed that constraint (E4) is the closest one to N * . For that reason, the robustness of configuration Q is r(Q, N * , 50) = d(N * , E4) = 9.8. This indicates that any demand N defined by N = N * + ∆N + + ∆N -can be processed in less than 50 units of time, provided that |∆N + |1 ≤ 9.8. The decision-maker is provided with a quantitative insight in the demand disturbances that the configuration can face while meeting the deadline. Besides, it should be stressed that this robustness measure relies on a worst-case scenario: some demands N defined by N = N * + ∆N + + ∆N -with |∆N + |1 > r(Q, N * , τ ) that can be processed in less than τ units of time may exist. In the present example, it can easily be checked that N 0 = N * + ∆N + + ∆N -with ∆N + and ∆N -defined below is processable in less than 50 units of time, while |∆N + |1 = 20. R 0 and C 0 max have been computed by solving LP with N 0 . In other words, |∆N + |1 ≤ 9.8 is a sufficient (but not necessary) condition for ensuring that the actual demand can be processed in less than 50 units of time.

It has been shown in section 5.2 that assessing the robustness of a configuration that does not allow to meet the deadline for the forecast demand is not relevant. This is clearly the case for Q since C max = 80.667 > τ . Consequently, this configuration is not suitable for N * (but Q is) and assessing its robustness is useless. However, this conclusion would be the opposite with the forecast demand N defined below. Where Cmax(Q, N ) is the optimal objective value for LP with configuration Q and demand N . Thus, configuration Q fits the forecast demand N whereas it is not the case for Q. This shows that robustness is not an intrinsic property of a configuration since it also depends on the forecast demand.

If the robustness offered by the current configuration is considered too low, the decision-maker may wish to add qualifications (i.e. replace some zero entries by one entries in Q). However, this can be a very costly operation (at least in the field of photolithography), and may not always be efficient for increasing robustness. In the present example, the configuration robustness potential value is ρ(Q, N * , τ ) = 6.9%. This means that Q is not really suitable for N * , as only 6.9% of its robustness potential is currently used. Then, the decision-maker can expect far better results by introducing new qualifications in the current configuration (this work does not address the effective choice of such additional qualifications). More precisely, by choosing the configuration Q * defined by Q * ij = 1 for all i and j, the robustness would be r(Q * , N * , τ ) = d(N * , E0) = 142.6. Thus, the robustness potential would be 100%. The comparison with Q * is rather extreme, as Q * is the most expensive configuration. This shows that the tools presented here provide useful information to the decision-maker that may wish to find the best possible tradeoff between configuration cost, performance and robustness.

Conclusion

The quantitative assessment of temporal performance along with robustness provides useful information to the decisionmaker who has to take into account many parameters and criteria to update a configuration for multi-purpose machines, or to design a new one. The proposed approach, while not polynomial, turns out to be practicable for real-life instances, provided they are dense enough and not too large. This work must be completed by a sensitivity analysis in order to assess the completion time deviation for any given value for the magnitude disturbance on the forecast demand. This would provide a valuable piece of information to the decision-maker, allowing him or her to estimate the trend for the completion time that the workshop may offer under unexpected overload. Future works should also address the problem of choosing the new qualifications to add to a configuration (if necessary) for increasing temporal performance and/or robustness. This could be achieved through a decision-making process mixing these tools to a configuration generator using a branch-and bound approach or an integer programming model. Then, the model could be extended to account for the qualification cost, that may depend on the machines as well as on the product types.

V

  Mj leads to the following equality, denoted (A):

4

 4 /mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.Computing the maximum rectangles of zeros

  .2, 3.4.2 and 3.4.2 show the average and maximum number of maximum rectangles of zeros (referred to as Avg. #MRZ and Max. #MRZ), as well as the average and maximum execution time in seconds for each couple (n, m), for d = 0.25, d = 0.5 and d = 0.75 respectively (referred to as Avg. CPU time and Max. CPU time in the Tables). The computer used is powered by a 3.2 GHz Intel R Pentium IV Processor with 1GB RAM under Microsoft R Windows XP. Algorithm 1 is implemented in C language and compiled with GCC.

Theorem 4 . 1 .

 41 The three following propositions are equivalent:1. The demand N can be processed in τ units of time or less 2. There exists a demand N such that: (a) N ≥ N (b) N can be processed in τ units of time or less (c) The production plan for N is balanced 3. The demand N satisfies (S), the following set of inequalities: (S) :

  The theorem is shown by transitivity of equivalence.

5

  Assessing the robustness of a configuration

  d(N * , E0) = 142.6 d(N * , E1) = 199.0 d(N * , E2) = 769.4 d(N * , E3) = 147.8 d(N * , E4) = 9.8 d(N * , E5) = 15.0 d(N * , E6) = 71.0 d(N * , E7) = 1164.0

  , N ) = 53.3 > τ, r(Q , N , 50) = 34

com/tprs Email: ijpr@lboro.ac.uk International Journal of Production Research

	Lj
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

http://mc.manuscriptcentral.

Table 1 :

 1 Maximum rectangles of zeros for low-density configurations (d = 0.25)

	Instances	Avg. #MRZ Max. #MRZ Avg. CPU time Max. CPU time
	n × m 15×30	4031	6683	0.274 s	0.484 s
	16×32	6659.40	10058	0.685 s	1.390 s
	17×34	9617	14496	1.730 s	4.749 s
	18×36	16200.35	29593	6.972 s	22.478 s
	19×38	25729	46072	27.676 s	99.691 s
	20×40	33424.25	45814	51.126 s	106.565 s
	21×42	55950.25	85555	140.927 s	267.348 s
	22×44	93573.55	140500	438.947 s	1036.250 s
	23×46	119337.35	208151	686.636 s	1492.820 s
	24×48	192475.50	295362	2001.847 s	4911.350 s
	25×50	-	-	> one hour	> one hour
	Instances	Avg. #MRZ Max. #MRZ Avg. CPU time Max. CPU time
	n × m 15×30	628.60	858	0.029 s	0.046 s
	16×32	824.85	1051	0.050 s	0.078 s
	17×34	1120.75	1723	0.064 s	0.109 s
	18×36	1487.90	2044	0.097 s	0.156 s
	19×38	2010.15	3271	0.158 s	0.297 s
	20×40	2419.55	3658	0.212 s	0.422 s
	21×42	3211.20	4325	0.337 s	0.517 s
	22×44	3804.55	5491	0.453 s	0.877 s
	23×46	4996.80	6692	0.726 s	1.284 s
	24×48	6191.60	8531	1.047 s	2.098 s
	25×50	8627.95	12335	2.096 s	3.801 s
	26×52	9939.65	15525	3.027 s	8.747 s
	27×54	12910.45	21516	5.630 s	18.255 s
	28×56	14583.40	21009	7.982 s	19.569 s
	29×58	18725.75	28465	15.454 s	40.599 s
	30×60	23950.70	34265	28.252 s	69.967 s
	31×62	26163.60	34319	36.464 s	64.189 s
	32×64	32463.50	39123	62.021 s	96.086 s

Table 2 :

 2 Maximum rectangles of zeros for mid-density configurations (d = 0.50)

Table 3 :

 3 Maximum rectangles of zeros for high-density configurations (d = 0.75)

	8

  Solving LP on the instance above leads to the following computational results for the production plan R and the completion time Cmax, but it must be stressed that this solution is not unique:

	International Journal of Production Research
	6.1 Solving the underlying scheduling problem
			47.542 39.917	0	0	
			0	0.271	7.063	0	
	R =	 	0	0	26.333	0	 	,	Cmax = 47.542
			0	0	0	5.733	
			0	7.354 14.146	0
	F o
	r
	P
		e e r
				R e v i e w
									O n
									l y
									135	
									22	
									79	 
									86	
									43
					12

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk

Acknowledgments

The author would like to thank the anonymous referees whose suggestions greatly improved the quality of this manuscript, and Joanna Ropers for correcting English grammar and spelling.


J C Nappé
email: jc.nappe@yahoo.fr.

M Benabdesselam

Ph Grosseau

B Guilhot

Effect of swift heavy ion irradiations in polycrystalline aluminum nitride

Keywords: AlN, irradiation, swift heavy ions, electronic interactions, thermoluminescence, optical absorption spectrophotometry PACS: 77.84.Bw, 61.82

Thanks to its high thermal conductivity, aluminum nitride may be a serious candidate as fuel coating for the Gas Fast Reactor. However, its behavior under irradiation is not entirely well understood. In order to catch a glimpse of this behavior, specimens were irradiated with swift heavy ions of different energies then characterised by both thermally stimulated luminescence and optical absorption spectrophotometry. With these techniques, the native defects, as well as those affected by irradiation, were identified: thus, by comparison to the virgin sample, no new defect detectable by these techniques is created by irradiations. Eventually, the fact that these techniques complement each other allowed to understand the effect of irradiation parameters on the defect concentration.

Introduction

Gas Fast Reactor (GFR) is one of the six new systems considered by Generation IV International Forum (GIF): it is designed to work under helium-pressure and high-temperature (1100-1300 K). These working conditions led to the selection of non-oxide refractory ceramics as fuel coating; the remarkable thermo-mechanical properties of carbides and nitrides make them great candidates for this application. However, their behavior under irradiation has to be studied in more details.

Among the studied materials, aluminum nitride is distinguished by its high thermal conductivity which should ensure the heat transfer between the coolant and the fuel. Actually, the values obtained for commercial aluminum nitride generally range between 160 and 200 W m -1 K -1 [START_REF] Couturier | Elaboration and characterization of a metal matrix composite: Al/AlN[END_REF][START_REF] Roulet | Corrosion behaviour of aluminium nitride in liquid aluminium: Influence of the microstructure[END_REF][START_REF] Jarrige | Effect of oxygen on the thermal conductivity of aluminium nitride ceramics[END_REF][START_REF] Virkar | Thermodynamic and Kinetic Effects of Oxygen Removal on the Thermal Conductivity of Aluminum Nitride[END_REF] with the highest value of 320 W m -1 K -1 measured on a single crystal [START_REF] Slack | Thermal expansion of some diamondlike crystals[END_REF].

The main drawback for the use of aluminum nitride in nuclear facilities is its ability to generate a radioactive long-lived isotope 14 C by transmutation of 14 N [6]. Nevertheless, this drawback can be reduced by the use of isotope tailored material, in which 15 N replaces 14 N [START_REF] Yano | Neutron irradiation effects on isotope tailored aluminum nitride ceramics by a fast reactor up to 2 × 1026 n/m2[END_REF].

The behavior of aluminum nitride ceramics is also studied in the framework of fusion reactors, in which this material, thanks to its high electrical resistivity (> 10 12  m at 300 K [START_REF] Hafidi | Contribution à l'étude des propriétés thermiques de céramiques à base de nitrure d'aluminium[END_REF]), could be used as electric insulator for toroidal current break, neutral beam injectors, and superconducting magnets. Thus, many works have been carried out to understand the behavior of aluminum nitride under neutron irradiation, especially by Yano et al. [START_REF] Yano | Neutron irradiation effects on isotope tailored aluminum nitride ceramics by a fast reactor up to 2 × 1026 n/m2[END_REF][START_REF] Yano | A HREM study of neutron-irradiation-induced dislocations in aluminum nitride[END_REF][START_REF] Yano | Effect of isochronal annealing on thermal diffusivity of neutron-irradiated AlN[END_REF][START_REF] Yano | High-resolution electron microscopy of neutron-irradiation-induced dislocations in SiC[END_REF][START_REF] Yano | Macroscopic length, lattice-parameter and microstructural changes in neutron-irradiated aluminum nitride due to annealing[END_REF][START_REF] Kim | Mechanical-properties of neutron-irradiated AlN[END_REF][START_REF] Yano | Neutron irradiation damage in aluminum oxide and nitride ceramics up to a fluence of 4.2 x 10(26) n/m(2[END_REF][START_REF] Yano | Swelling and microstructure of AlN irradiated in a fast-reactor[END_REF][START_REF] Yano | Thermal and mechanical-properties of neutron-irradiated aluminum nitride[END_REF][START_REF] Akiyoshi | Thermal conductivity of ceramics during irradiation[END_REF]. However, as far as we know, no study has been performed concerning the characterization of defects induced by irradiation.

Thus, the aim of this study is to characterize the evolution of native defects of aluminum nitride as a result of swift heavy ion irradiations, using thermally stimulated luminescence and optical absorption spectrophotometry.

Material and methods

Sample preparation and irradiation

The studied material was obtained by sintering of AlN powders with Y 2 O 3 (sintering aid: 3 mol%) at 2123 K in nitrogen atmosphere during 3 h. This polycrystalline aluminum nitride consists of 8 µm grains with 0.89 % of oxygen. Its thermal conductivity was measured between 198 and 214 W m -1 K -1 . Preparation of the samples before irradiation was already described elsewhere [START_REF] Nappé | Damages induced by heavy ions in titanium silicon carbide: Effects of nuclear and electronic interactions at room temperature[END_REF]. Briefly, they were polished with diamond suspensions down to 1 µm, and then subsequently irradiated with different swift heavy ions (SHI). This kind of irradiation aims to simulate the interactions of fission products without the drawback of neutrons activating the specimens. Thus, the specimens can be manipulated just after irradiations without any particular precautions. SHI were provided by GANIL facilities (Caen, France). In order to know the effect of electronic stopping power, samples were irradiated with different ion-energy couples. Table 1 shows the effect of three ion irradiations in AlN, to wit the average projected range (Rp), both the electronic (Se) and the nuclear (Sn) stopping powers. All these data were obtained with SRIM-2008 code [19]. In addition, this table contains the experimentally reached fluences and the irradiation temperatures. The fluences reached for 930 MeV Xe irradiations are lower than those reached with the two other irradiations. This is due to the highest energy, which may heat the samples during irradiation. Thus, in order not to heat the samples, the flux used for this irradiation was lower than for the two others. The required duration for equivalent fluences were hence much higher.

Thermally stimulated luminescence measurements

Principle of thermally stimulated luminescence (TSL) may be found elsewhere [START_REF] Porte | Thermoluminescence of aluminium nitride: Influence of the thermal treatment[END_REF]21]. For TSL measurements, specimens were excited at room temperature (RT) during five minutes by a deuterium UV-lamp (50 W) covering the 200-400 nm spectral range.

Two minutes after excitation, TSL readout was achieved by linear heating of the excited sample from RT to 673 K at a rate of 1 K s -1 . TSL intensity was recorded by means of a photomultiplier (Philips 2018 B) covering the 200-600 nm range. For spectral analyzis of TSL emission, spectra were caught by an optical multi-channel analyzer (Princeton Instruments) covering the 300-800 nm range.

Optical absorption measurements

Optical absorption spectrophotometry (OAS) was performed at room temperature using a Varian Carry 500 spectrophotometer equipped with a diffuse reflection sphere. The spectra were recorded on the 250-800 nm wavelength range with a scan rate of 600 nm min -1 . The absorption reference was made of polytetrafluoroethylene (PTFE). For all the samples, optical absorption was measured before and after ion irradiation in the same experimental conditions. The absorption bands affected by irradiation have been characterized by the difference curve, which is obtained by subtracting the absorption spectrum of the virgin sample from the spectrum of the same sample after ion irradiation.

Results

Thermally Stimulated Luminescence

Typical evolution of TSL glow curves as a function of fluence is shown in Figure 1. The TSL signal of the AlN sample before irradiation consists of one broad and structureless glow curve peaking at approximately 400 K. This broad peak is due to the presence of oxygen in the nitrogen sublattice. More precisely, it attests the presence of trapping centers (TC), certainly oxygen in substitution of nitrogen O N and/or bivalent 2O N centers [START_REF] Porte | Thermoluminescence of aluminium nitride: Influence of the thermal treatment[END_REF][START_REF] Rosa | Mechanism of thermoluminescence in Al:O, Czech[END_REF][START_REF] Trinkler | Thermoluminescence properties of AlN ceramics[END_REF]. Whatever the ion-energy couple or the irradiation temperature are, the TSL intensity decreases with increasing fluence. Moreover, neither shift of the main peak, nor presence of new TSL peaks induced by irradiation were observed. In order to compare the TSL signal for each irradiation, we calculated the area under each glow curve; Figure 2 shows the variation of TSL area as a function of fluence for each irradiation. In this Figure, we can note three different trends.

First, for 74 MeV Kr and 92 MeV Xe irradiations at RT, TSL area quickly decreases up to 10 13 cm -2 then tends to stabilize. Second, for irradiation at 773 K with 74 MeV Kr we can note a slight decrease up to 10 13 cm -2 followed by a sharp decline occuring between 10 13 and 10 14 cm -2 , then finally a stabilization of the signal. Third, whatever the temperature of 930 MeV Xe irradiation is, TSL area decreases more swiftly than for other irradiations up to stabilization. These trends suggest that the degradation process revealed by TSL is the same for different irradiations but of similar energies. noticeable for the highest fluences (see inset): the emissions at 3.17 ± 0.01 eV and 2.67 ± 0.17 eV have been respectively attributed to complex centers V Al -O N -3N and V Al -2O N -2N while the one at 3.60 ± 0.04 eV may be associated to transitions from V Al to O N [START_REF] Rosa | Mechanism of thermoluminescence in Al:O, Czech[END_REF]. Concerning the low bands appearing at lower energies, they might be linked to some impurities: the emission at 2.04 ± 0.05 eV is attributed to the presence of Mn 4+ in substitution of Al 3+

(Mn Al ) [21] whereas that centered at 1.79 ± 0.03 eV has been reported by Martin et al. as being due to the presence of Cr 3+ in substitution of Al 3+ (Cr Al ) [START_REF] Martin | Visible emission from amorphous AlN thin-film phosphors with Cu, Mn, or Cr[END_REF]. Inductively coupled plasma (ICP) analyzis confirmed the presence of both impurities in our samples, at concentrations of 1.9 ± 0.1 and 2.0 ± 0.1 ppm respectively. After irradiation, it is interesting to remark that no new RC is created: as shown for the glow curves, the intensity of the spectra decreases with the fluence, whatever the irradiation is. as explained in section 2.3., the induced absorption was obtained by subtracting the absorption spectrum of a virgin sample from that of the same sample after ion irradiation. After irradiation, the optical density of 7 bands increases more particularly; the energy of these bands (Table 2) was determined after Gaussian decomposition of the obtained spectra (not shown here). Each of these 7 bands was ever on the absorption spectrum of the virgin sample: irradiations do not induce any new absorbing defect in the 250-800 nm range. In the literature, few information are available concerning the identification of absorbing defects in AlN. The data presented in Table 2 has been performed using those from other materials, in particular from alumina. Indeed, according to the high chemical affinity of oxygen towards aluminum, it is likely that a thin layer of alumina is present on the surface of AlN samples. Finally, as for the impurities of both manganese and chromium observed in the TSL emission spectra, the presence of iron in AlN samples (Table 2) was clearly identified by ICP, at a ratio of 330 ± 20 ppm.

Optical absorption spectrophotometry

Discussion

The estimation of the concentration of each absorbing defect (AD) is complicated. However, for each irradiation parameter, we observed that when an absorption band increases, all the other bands do likewise, even if not with the same rate (Figure 4b); this observation was also right for decrease of the optical density of the irradiated sample compared with that of the virgin one. That is why, in order to understand the effects of SHI in AlN material, we decided to follow the variation of the overall absorption of each sample as a function of irradiation parameters by calculating, for each performed irradiation, the average of the optical density induced by irradiation; this global absorption was called average optical absorption induced by irradiation (Figure 5). Thus, whatever the irradiation carried out at RT is, the induced absorption remains constant and almost null up to 10 12 cm -2 , and then increases. For irradiations performed at 773 K with 930 MeV Xe ions, the variation is similar to that at RT except for highest fluence of 4.5x10 13 cm -2 for which the induced absorption decreases slightly at 773 K while it increases for RT. This behavior could be explained by an annealing of the defects created by irradiation from a critical fluence  c ranging between 10 13 and 4.5x10 13 cm -2 at 773 K. Thus, above  c the defects would be annealed, inducing a decrease of the overall concentration of defects. Finally, for the 74 MeV Kr irradiation at 773K, we observe a sharp drop in absorption compared to the virgin sample. This drop is followed by the same increase with the fluence as obtained for irradiation at RT: the average induced absorption hence reachs a comparable value to that before irradiation for 10 14 cm -2 fluence, and becomes higher for 10 15 cm -2 . Thus, if the electronic stopping power Se decreases (13.1 keV nm -1 for 74 MeV Kr and 21.8 keV nm -1 for 930 MeV Xe, see Table 1), the critical fluence of annealing  c decreases, being lower than 10 12 cm -2 for irradiation at 773 K with 74 MeV Kr. Moreover, the evolution of this average induced absorption reflects a balance between annealing by irradiation temperature and creation of defects by irradiation itself, which inclines toward the creation only for high fluences.

Thanks to OAS, we know that irradiations create defects that absorb the UV-visible radiations, and therefore some of the TSL emissions mainly located in this range. To better understand the irradiation effects due to SHI, we now compare the results obtained by both techniques.

First, by considering the annealing at 773 K of the absorbing defects induced by 74 MeV Kr ion irradiation (Figure 5a), one would expect an increase of the TSL signal for low 74 MeV Kr fluences, above the signal of the virgin sample.

However, according to Figure 2a, the TSL signal also decreases for this irradiation. Therefore, this ion-energy couple (74 MeV Kr) seems to induce the decrease of TSL center concentration. Moreover, regarding the fact that there is no difference between the TSL signals related to 930 MeV Xe irradiations performed at RT and 773 K, we can consider that an increase of the irradiation temperature does not induce any increase of TSL center concentration, but only the annealing of AD. Therefore, the differences between the TSL curves of 74 MeV Kr ion irradiations would be only due to both the creation of AD at RT and their annealing at 773 K.

Second, for the same fluences and for irradiations performed at RT, TSL signal is nearly the same for both 74 MeV Kr and 92 MeV Xe irradiations, and then decreases in the case 930 MeV Xe irradiation. As for these irradiations the AD concentration is the same whatever the electronic stopping power is, this latter has not any effect on TSL center concentration up to a threshold Se th above which this concentration strongly decreases; according to data reported in Table 1, this threshold would be between 17.7 and 21.8 keV nm -1 .

Eventually, when the fluence increases at RT, we noted that both the AD concentration increases in the same way whatever the Se is, and the TSL signal decreases as a function of Se. Therefore, the decrease of TSL signal as a function of the fluence cannot directly be attributed to the creation of AD. So we can conclude that when the fluence increases, the concentration of TSL centers decreases.

Thus, by comparing the TSL and OAS signals we now understand the effect of each irradiation parameter on the concentrations of both TSL centers and AD. But how to justify the fact that TSL center concentration decreases when either the electronic stopping power or the fluence increases? The first hypothesis is that, following the thermal spike model [START_REF] Toulemonde | Transient thermal-process after a high-energy heavy-ion irradiation of amorphous metals and semiconductors[END_REF][START_REF] Toulemonde | Thermal spike model in the electronic stopping power regime[END_REF], electronic stopping power provides heat in the material through electron-phonon interactions. This heat could induce some annealing effect as ever remarked by Benyagoub et al. in silicon carbide [START_REF] Benyagoub | Athermal crystallization induced by electronic excitations in ion-irradiated silicon carbide[END_REF]. Thus, the TSL center would be annealed, explaining the decrease of the TSL signal when increasing the fluence. Nevertheless, TSL centers observed in our specimens being due to the substitution of both Al and N atoms by impurities (oxygen, chromium, manganese, ...), their annealing appears to be difficult. Therefore, we rather think that electronic stopping power induces some modifications of these TSL centers. For instance, Mn 4+ which is in substitution of Al 3+ may be ejected from its site Mn Al and becomes an intersticial Mn i , or the complex centers V Al -O N -3N and V Al -2O N -2N may be dissociated by Se;

these modifications may also be met with trapping centers such as oxygen in substitution of nitrogen. This hypothesis would explain the decrease of the concentration of TSL centers all the while creating other (absorbing) defects in the material. Whatever the hypothesis is, further characterisations would be necessary; for example, photoluminescence measurements would allow us to know the variation of the recombination center concentration as function of the irradiation parameters.

Conclusion

In this study, we first identified both trapping and recombination centers by thermally stimulated luminescence, but not all the absorbing defects are affected by SHI irradiation. A major result is that whatever the analytical technique is, no new defect was detected after these irradiations.

Second, thanks to the fact that TSL and OAS complement each other, we discriminated the effect of fluence, electronic stopping power, and temperature on the overall concentration of both TSL centers and AD. Thus, the higher the fluence, the lower the TSL center concentration, and the higher the AD concentration. Moreover, we defined an electronic stopping power threshold Se th , which ranges between 17.7 and 21.8 keV nm -1 , above which the TSL center concentration strongly decreases. On the contrary, when Se < Se th , the electronic stopping power has no effect on both TSL centers and AD concentrations. Concerning the temperature, we concluded that this parameter has no effect on TSL center concentration, whereas it induces the decrease of the AD concentration if the fluence is higher than a critical fluence  c which depends on ion-energy couple. This result is quite encouraging for application of aluminum nitride as a component of fuel for the GFR, which should operate at temperatures well above 773 K.

Eventually, we tryed to know how the TSL signal decreases as a function of both electronic stopping power and fluence, and our hypothesis is that irradiations induce a change of the TSL centers, creating non-luminescent defects.

Figure 1 :

 1 Figure 1: TSL glow curves of polycrystalline AlN sample as a function of 92 MeV Xe ion fluence.

Figure 2 :

 2 Figure 2: TSL area under glow curve as a function of ion fluence for (a) low energy irradiations, and (b) high energy irradiations; lines are only to guide the eye.

Figure 3

 3 Figure3shows the typical evolution of TSL emission spectra of the main peak (400 K) as a function of fluence; position of the bands was determined by Gaussian decomposition (not shown here). The spectrum of the virgin sample shows the main recombination centers (RC) in AlN already identified byBenabdesselam et al. [21]. First, the main broad band centered at around 3.2 eV is the convolution of the emission of probably three RC, which are more easily

Figure 3 :

 3 Figure 3: Spectral distribution of the main TSL peak (400 K) as a function of 92 MeV Xe ion fluence in polycrystalline AlN samples; the inset is a zoom of the main emission band for the highest fluences.

Figure 4a shows the

  Figure 4a shows the optical absorption spectrum of AlN before irradiation. It presents several absorption bands in the 250-800 nm range. In this Figure, one can see the spectrum of the sample irradiated at RT with 92 MeV Xe to 10 15 cm -2 ;

Figure 4 :

 4 Figure 4: Optical absorption of polycristalline AlN in the 250-800 nm wavelength range: (a) comparison between the virgin sample and the same sample irradiated at room temperature with 92 MeV Xe to 10 15 cm -2 , and (b) optical absorption bands induced by 92 MeV Xe irradiation.

Figure 5 :

 5 Figure 5: Average optical absorption induced by (a) low energy irradiations, and (b) high energy irradiations within the 250-800 nm wavelength range; lines are only to guide the eye.

Table 1 :

 1 Summary of irradiation conditions; RT stands for room temperature.

	Ion irradiation Rp (µm) Se (keV nm -1 ) Sn (keV nm -1 ) Fluences (cm -2 )	Temperatures
	74 MeV Kr	9.7	13.1	0.054	10 12 , 10 13 , 10 14 , 10 15	RT & 773 K
	92 MeV Xe	9.4	17.7	0.132	10 12 , 10 13 , 10 14 , 10 15	RT
	930 MeV Xe 46.2	21.8	0.018	10 11 , 10 12 , 10 13 , 4.5x10 13 RT & 773 K

Table 2 :

 2 Identification of absorption bands affected by 92 MeV Xe irradiation in polycristalline AlN sample; F centers are defects based on anionic vacancies.

	 (nm) E (eV)	Related defect	References
	678 ± 7 1.83 ± 0.02 Cr 3+ (zero-phonon)	[25-29]
	630 ± 6 1.97 ± 0.02 Fe Al	[30]
	511 ± 6 2.43 ± 0.03 ?
	455 ± 9 2.73 ± 0.05 V N or Al i or F 2 2+	[31-36]
	406 ± 8 3.06 ± 0.06 V Al or Cr Al	[37,38]
	370 ± 6 3.35 ± 0.05 F 2 +	[39]
	290 ± 8 4.28 ± 0.12 O related center (F 2 ?) or Al i [36,40-42]

Acknowledgment

Authors would like to gratefully thank Isabelle Monnet (GANIL) for her great help during the irradiation of the samples. They also would like to thank Yves Jorand (INSA, Lyon, France) for providing the virgin materials. This study was partly funded by the French research group MATINEX.