Robust scheduling optimization for flexible manufacturing systems with replenishment under uncertain machine failure disruptions

Control Engineering Practice 92 (2019) 104094

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Robust scheduling optimization for flexible manufacturing systems with replenishment under uncertain machine failure disruptions✩ Zhiguo Wang a , Chee Khiang Pang b ,∗, Tsan Sheng Ng a a b

Department of Industrial Systems Engineering and Management, National University of Singapore, Singapore 117576, Singapore Engineering Cluster, Singapore Institute of Technology, 10 Dover Drive, Singapore 138683, Singapore

ARTICLE

INFO

Keywords: Flexible manufacturing systems (FMSs) Machine failure disruption Distributionally robust optimization Mixed-integer linear program (MILP)

ABSTRACT This paper studies the scheduling problem for the flexible manufacturing systems (FMSs) under uncertain machine failure disruptions, where machine allocations and job schedules need to be determined to achieve a set of production due-date requirements as well as possible. A robust scheduling optimization model is proposed based on the concept of threshold scenario, bounded by which the due-dates are guaranteed to be achieved. It is shown that the associated stochastic scheduling problem can be equivalently solved by computing the solution of a mixed-integer linear program (MILP). Computational results show that our proposed model performs well in achieving the planned due-dates under uncertainty when compared to various standard approaches. The practical applicability of our approach is verified using a real stamping industry application, in which the stamping parts are various types of voice coil motor yokes used in commercial hard disk drive actuators. Apart from FMSs, the proposed approach can also be applied to various other industries including project scheduling, airline scheduling, transportation scheduling.

1. Introduction Flexible manufacturing systems (Pach, Berger, Sallez, & Trentesaux, 2015; Wang, Pang, & Ng, 2018; Yue, Xing, Hu, Wu, & Su, 2016) are commonly seen in modern automated production systems (Yang, Gao, Simon, Zhu, & Su, 2018) which possess great flexibility in machine allocation and part routing. Each machine may be capable of processing multiple jobs, and each job may be processed by alternative machines (Hu, Zhou, & Li, 2011; Sindicic, Bogdan, & Petrovic, 2011). A good scheduling method is shown to be able to further improve the energy efficiency, production effectiveness, and environmental sustainability (Chen, Zhang, Arinez & Biller, 2013; Li et al., 2013). Uncertainties commonly exist in the manufacturing system in practice, for instance, the uncertain processing times (Wang, Xia, Meng & Li, 2017), the uncertain processing cost (Wang et al., 2018), and the uncertain machine failure disruptions. In many existing scheduling methods (Blazewicz, Ecker, Pesch, Schmidt, & Weglarz, 2013; Jin, Zhao, Han, & Wang, 2018; Wang, Liao, Tang & Ning, 2017), the machines are assumed to be reliable and not subject to failure throughout the entire scheduling process. The reasons why machine failure was not previously considered are mainly twofold. First is that many traditional manufacturing systems like flow shop do not possess the flexibility of FMSs. The machine failure will cause the termination of the manufacturing system, and the production process is not able to

resume until the malfunctioned machine is replenished. Second is that the task allocations in traditional manufacturing systems are rather static. The machine tool wear rate can be forecasted to implement preventive maintenance in advance. As a contrast, due to the flexible and dynamic characteristics of modern FMSs, machine failure has become a considerably frequent disturbance and mot importantly the disruption can be uncertain in practice (Chen, Wei & Hu, 2013; Li & Ierapetritou, 2008). As a result, FMSs scheduling model which takes the machine failure disruption into consideration is in great demand (Pinedo, 2016). To the best of our knowledge, the existing methods which consider about machine failure disruptions can be classified into several types. One type of methods focus on developing deadlock avoidance supervisory control strategies such that the failure of any machine does not propagate through blocking to stall other portions of the system (Chew & Lawley, 2006; Wang, Chew, & Lawley, 2008). Luo, Xing, Zhou, Li, and Wang (2017) presents a new scheduling method by combining a robust supervisory control policy and hybrid heuristic search for the automated manufacturing system. The other type of methods are known as preventive scheduling (Janak, Lin, & Floudas, 2007; Li & Ierapetritou, 2008). These methods make use of historical data and forecasting techniques to predict the behavior of uncertainty in the future. This type of method requires accurate and sufficient information prior to realization of the uncertainties. They usually serve as the

✩ This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. ∗ Corresponding author. E-mail addresses: [email protected] (Z. Wang), [email protected] (C.K. Pang), [email protected] (T.S. Ng).

https://doi.org/10.1016/j.conengprac.2019.07.012 Received 1 February 2019; Received in revised form 6 June 2019; Accepted 19 July 2019 Available online 26 August 2019 0967-0661/© 2019 Published by Elsevier Ltd.

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

2. Problem description

planning basis and a periodic update is required as time passes to obtain an accurate forecast. Another type of method is reactive scheduling (Janak, Floudas, Kallrath, & Vormbrock, 2006; Pach et al., 2015), which is a re-scheduling technique to modify the baseline schedule during the process to adapt to the change in disruption. The common techniques include partial re-scheduling, and full re-scheduling.

In this section, the scheduling of the flexible manufacturing systems (FMSs), under uncertain machine failure disruptions will be described and modeled. The notations in this work follow the conventions as follows: a random variable is denoted by the tilde sign, for instance, 𝑥. ̃ Bold letters such as 𝒙 represent vectors or matrices. A vector of ̃ The set of non-negative real numbers random variables is denoted as 𝒙. is denoted as R+ . The expression 𝒙̃ ≤ 𝒚̃ refers to the relationship that 𝒚̃ stochastically dominates 𝒙̃ in the first-order component-wise.

These approaches either require an excessive amount of information regarding the disruptions or only passively react when the random disruption is realized. This paper introduces a proactive scheduling optimization model based on the concept of threshold scenario, over which the planned production due-dates are guaranteed to be achieved. In addition, most existing works assume that only a single or a set of machines are failure-prone. In this work, the most general case of uncertainty, namely every machine is subject to disruption while processing any job, will be considered. In addition, machine replenishment which is a common practice in the industry has rarely been considered in the existing literature of FMSs scheduling. Motivated by this, our formulation allows the disrupted machine to be replenished to process jobs after repair. This gives the model more flexibility and it therefore has wider range of applications. The developed max–min scheduling model belongs to the class of distributionally robust optimization models (Goh & Sim, 2010; Mak, Rong, & Zhang, 2014) in the literature. These models assume partial distributional information of the uncertain parameters, namely the probabilities of machine failure scenarios in this context. The decisions are chosen to optimize the worst-case objective value among all possible distributions with the specified probabilities of each individual machine. Compared to estimating the joint probabilities of a large number of possible machine failure scenarios, it is less burdensome to estimate the failure probability of each individual machine. The failure probability of a machine can be evaluated by the cumulative operation time and the operating conditions such as speed, temperature, and mode (Liu, 2000). Several studies (Alidaee & Li, 2014; Jia, Nadarajah, & Guo, 2018) have studied about the machine failure probability in scheduling and modeled it using Weibull distribution.

Nomenclature 𝛱 |⋅| 𝜋𝑎 𝑣𝑎,𝑏 𝑖 𝑛 [𝑛] R R𝑖 𝑅𝑟 V𝑐 (V𝑛𝑐 ) 𝑥𝑖𝑟 𝑠𝑖 𝒔 𝑑𝑖𝑟 𝛥𝑑𝑖𝑟 𝑐𝑖𝑟 P Q 𝑤𝑖𝑗 𝑧𝑖𝑗

The computational complexity is another motivation of this work. In practice, the occurrence of the machine disruption is uncertain (Chew, Wang, & Lawley, 2011; Hsieh, 2004). Given a number of machines in the system, the possible combination of machine disruption and replenishment scenarios can be enormous. The decisions on the machine allocation and the task starting time are dependent on the disruption scenarios, which results in a multi-stage stochastic scheduling problem. This can be challenging to solve. Built upon the concept of threshold scenario, the optimal solution to the proposed model can be equivalently computed by solving a mixed integer linear program (MILP) with a reasonable size. This can be solved using the existing commercial solver IBM ILOG CPLEX. The model performance and practical applicability are evaluated through simulated computational experiments as well as real stamping industry case study. The results show that the proposed model has a more comparable performance in terms of both success rate in achieving the planned due-dates and computational time. The rest of the paper is presented as follows: Section 2 describes the discussed FMSs scheduling problem as well as its mathematical formulation. Section 3 presents the proposed robust scheduling model under uncertain machine disruptions uncertainties. Section 4 illustrates the use of the proposed model with a simulated case study and evaluates the model performance across a wide range of full size FMSs. Section 5 tests the practical applicability of the proposed model using one real stamping industry application. The model performance is evaluated against various benchmark approaches. Section 6 concludes the important findings of this work and some future work directions.

𝑀 𝝐̃ 𝝐ℎ 𝝐 𝑖𝑟 𝑡𝑖 𝜶 𝜶 𝑖𝑟 𝑃̂𝑖𝑟 ℎ(𝒔,𝒅)

ℎ𝑖(𝒔,𝒅)

Set of part types to be manufactured. The cardinality of a set. Job sequence of 𝑎th part type. The 𝑏th job in part type 𝜋𝑎 . The index of job 𝑖. The total number of jobs to be processed. The index set of all jobs to be processed. Set of available manufacturing resources. Set of resources capable of processing job with index 𝑖. The 𝑟th resource in R. Set of choice (non-choice) jobs. Binary variable indicating if machine 𝑅𝑟 is assigned to job 𝑖. Starting time of job 𝑖. Vector of the start time for all jobs. Processing time for job 𝑖 using machine 𝑅𝑟 . Replenishment time for machine 𝑅𝑟 before the start of job 𝑖. Unit time manufacturing cost for job 𝑖 using Set of precedence pairs among all jobs. Set of job pairs without precedence. Binary variable indicating if the same machine is assigned to job 𝑖 and job 𝑗. Binary variable indicating if the same machine is assigned to job 𝑗 before job 𝑖. An arbitrarily large number. The random matrix which indicates the machine failure status. The binary matrix which indicates the status of the ℎth machine failure scenario. The (𝑖, 𝑟)th component of 𝝐. The planned due-date for job 𝑖. The threshold scenario of machine disruptions. The (𝑖, 𝑟)th element of 𝜶. Probability of machine failure in 𝑅𝑟 while processing job 𝑖. The set of available historical information regarding job start time 𝒔 and durations 𝒅. The set of available historical information regarding (𝒔, 𝒅) up to but does not include job 𝑖.

2.1. Descriptions of FMSs A general FMS with shared machines and flexible parts routing will be presented in this section. The typical assumptions that are commonly used in manufacturing systems modeling are as follows: • No pre-emption is allowed: once allocated, a machine cannot be assigned to another job until it completes the current job. • Pre-requisite: a job will hold the machines already allocated to it until all the required pre-requisite machines are ready. 2

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

• No multitasking: each machine can process at most one job at a time. • All jobs are released at the beginning of the manufacturing process. • All machines are available when the jobs are released to FMSs. • The travel time between machines is negligible. • There is no set-up time for processing a job. • A machine is subject to stochastic failure only when processing jobs. There is no other type of uncertainties which may disturb job processing except for machine failure. • The processing of jobs on a machine terminates when machine failure occurs and resumes immediately after replenishment.

Otherwise if |R𝑖 | = 1, job 𝑖 is called a non-choice job. As such, the set of jobs V can be partitioned into two mutually exclusive sets V = V𝑐 ∪ V𝑛𝑐 , where V𝑐 and V𝑛𝑐 represent the set of choice and non-choice jobs respectively. This kind of flexibility results in a vast number of scheduling combinations. As a matter of fact, a small number of jobs and machines may lead to a complex scheduling problem. Example 1. Fig. 1 is the block diagram of an example FMS. This example consists of six machines 𝑅1 − 𝑅6 and produces two part types 𝜋1 , 𝜋2 . The manufacturing process starts with the input of raw materials and ends with the delivery of final products. It has the operational flexibility in the sense that two types of products are produced in the system using two different processes and some jobs can be processed by alternative machines. Part type 1 consists of three jobs 𝑣1,1 , 𝑣1,2 , 𝑣1,3 . The first job 𝑣1,1 is a choice job which can be processed by either 𝑅2 or 𝑅3 . The rest of the jobs are non-choice jobs which are processed by 𝑅5 , 𝑅6 . Part type 2 consists of four jobs 𝑣2,1 , 𝑣2,2 , 𝑣2,3 , 𝑣2,4 . Job 𝑣2,2 is a choice job which can be processed by either 𝑅2 or 𝑅4 . Job 𝑣2,1 , 𝑣2,3 , 𝑣2,4 are non-choice jobs, each of which is processed by 𝑅1 , 𝑅5 , and 𝑅6 .

The key characteristics of the class of typical FMSs are summarized as (i) There are multiple part types to be produced. Each part type contains a pre-defined sequence of jobs. (ii) There exists a common pool of available production machines. Each machine is distinct and may be capable of processing a few different jobs. (iii) A single machine is assigned to process a job each time. (iv) The job processing time and machine production cost in general can vary with respect to different jobs and machines. (v) The machine allocation and parts routing are flexible, which is subject to the decision maker’s strategy. The set of part types is defined as 𝛱 = {𝜋𝑎 ∶ 𝑎 = 1, … , |𝛱|},

In network modeling, the FMS system networks are commonly represented by one of the most widely used network representations, namely activity-on-node (AON) (Klastorin, 2004). The nodes (circles) denote the jobs and the arcs (arrows) denote the precedence relationships. As mentioned in (4), the complete set of jobs can be indexed as [𝑛]. Fig. 2 illustrates the AON representation of the FMS network in Example 1. Without loss of generality, let 𝑛 be a final (possibly ‘dummy’) job with zero duration. To denote the start of the manufacturing process, let 0 be a first (possibly ‘dummy’) job with zero duration. The execution of a job takes a certain amount of time, which is called the processing time or job duration. The duration of job 𝑖 using machine 𝑅𝑟 is denoted as 𝑑𝑖𝑟 . This implies that the machine occupied by a job becomes available again only after at least a period of 𝑑𝑖𝑟 .

(1)

where | ⋅ | represents the set cardinality. As mentioned, each part type 𝜋𝑎 is produced by following a production sequence of jobs as 𝜋𝑎 ≜ 𝑣𝑎,1 … 𝑣𝑎,𝑏 … 𝑣𝑎,|𝜋𝑎 | ,

(2)

where 𝑣𝑎,𝑏 is the 𝑏th job in manufacturing part type 𝜋𝑎 . For instance, the job of processing raw metal ores into useful metal parts requires a sequence of jobs like crushing, smelting, and refining. In practice, these required jobs as well as their production procedures are specified by the specialist with subject expertise in the relevant industrial application. This is beyond the scope of this work. For the matter of schedule planning, the job sequences are assumed to be known to the decision maker. Each part type is considered as completed only when all the intermediate jobs are completed. The complete set of all jobs to be processed across all part types in the FMS can be summarized as the set {𝑣𝑎,𝑏 ∶ 𝑎 = 1, … , |𝛱|, 𝑏 = 1, … , |𝜋𝑎 |}.

2.2. Machine failure disruption Machine failure disruptions like breakdown are very common in practice due to reasons like machine degradation, tool wear, heavy workload. Define 𝝐 𝑛×|R| to indicate the machine status, where each entry of the matrix, denoted as 𝝐 𝑖𝑟 , 𝑖 ∈ [𝑛], 𝑟 ∈ R𝑖 , takes value either 0 or 1. Under uncertainty of machine failure disruptions, 𝝐 is random and thus is denoted as 𝝐̃ . Similarly, 𝝐 𝑖𝑟 is denoted as 𝝐̃ 𝑖𝑟 under uncertainty and is defined as { 1, if machine 𝑅𝑟 fails at the start of job 𝑖, 𝝐̃ 𝑖𝑟 = (6) 0, otherwise.

(3)

Alternatively, these jobs can be re-indexed as job 𝑖 = 1, … , 𝑛, 𝑛 = ∑|𝛱| |𝜋 | and let 𝑎=1 𝑎 [𝑛] = {1, … , 𝑖, … , 𝑛}

The domain set 𝜩 𝑛×|R| is defined by {0, 1}𝑛×|R| and 𝜩 𝑛×|R| contains all the possible values that the machine status 𝝐̃ can take, namely 𝝐̃ ∈ 𝜩. Alternatively, it can be denoted as 𝟎𝑛×|R| ≤ 𝝐̃ ≤ 𝟏𝑛×|R| . In the next section, this set will be effectively formulated as 𝟎 ≤ 𝝐̃ ≤ 𝜶 when upper bound 𝜶 is introduced. It is assumed that the disrupted machine can be replenished to continue the job processing. The replenishment can refer to the repair of the existing malfunctioning machine or the purchase of a new machine. The machine replenishment is modeled by extending the normal processing duration of the corresponding job by the replenishment period 𝛥𝑑𝑖𝑟 . As a consequence, the general job duration with replenishment is denoted as 𝑑̃𝑖𝑟 = 𝑑𝑖𝑟 + 𝝐̃ 𝑖𝑟 𝛥𝑑𝑖𝑟 , 𝑟 ∈ R𝑖 .

(4)

denote the index set of jobs. The purpose of doing this is mainly to simplify the notations in later model formulation. As can be seen later, this does not affect the results while it actually helps make the model formulation neater and readable. The set of available manufacturing machines is defined as R = {𝑅1 , … , 𝑅𝑟 , … , 𝑅|R| }.

(5)

In FMSs, the machine commonly refers to, but is not limited to, the manufacturing machines. Generally speaking, it can refer to any machine involved in the manufacturing process, like equipment, vehicles, or even technicians in a more general setting. Each of the machines may be capable of performing one or a few different jobs depending on its mechanical properties. The subset of machines that are capable of performing job 𝑖 is defined as R𝑖 , where R𝑖 ⊆ R and R𝑖 ≠ ∅. This leads to the concept of the choice job. A choice job is defined as the job that can be processed by at least two different machines. A non-choice job is defined as the job that can be only processed by a single machine. Due to the flexible routing, choice jobs are ubiquitous in FMSs. Mathematically, if |R𝑖 | ≥ 2, job 𝑖 is called a choice job.

2.3. Scheduling formulation for FMSs The schedule planning in FMSs is mainly about determining the machine-job allocation as well as the scheduled start time of each job. As such, to represent the machine-job allocation plan, the binary decision variables 𝑥𝑖𝑟 , 𝑖 ∈ [𝑛], 𝑟 ∈ R𝑖 , are introduced as { 1, if machine 𝑅𝑟 is allocated to job 𝑖, 𝑥𝑖𝑟 = (7) 0, otherwise. 3

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

Fig. 1. Layout of a flexible manufacturing system example.

Fig. 2. AON representation of Example 1.

In the class of FMSs discussed, the job is processed by one machine each time. Mathematically, it is formulated as the constraint, ∑ 𝑥𝑖𝑟 (𝝐) = 1, ∀𝑖 ∈ [𝑛], ∀𝝐 ∈ 𝜩, (8)

another set of scheduling decisions to be made would be the start times of each job, which are denoted as,

𝑟∈R𝑖

where 𝑠𝑖 represents the start time of job 𝑖. By default, it is assumed that the start of each part type is not delayed. In addition, it is assumed that no setup or lap time is required before and after each job. As such, the start time of a succeeding job may also represent the completion time of its preceding job. As a convention, two dummy jobs with zero duration are artificially created to denote the start and the end of the whole scheduling process. Without loss of generality, let 𝑛 be the final (possibly ‘dummy’) job with zero duration, and hence 𝑠𝑛 represents the completion time of the entire manufacturing process, also known as the production makespan. The job sequence mentioned earlier actually specifies one important relationship among jobs in a manufacturing system, that is the precedence relationships. Job 𝑖 precedes job 𝑗 if job 𝑗 can only start when job 𝑖 has completed. For instance, in automobile manufacturing, the wheel and engine can be installed only when the chassis is molded. Since each part type refers to the production of a different product, it

𝒔 = (𝑠1 , … , 𝑠𝑖 , … , 𝑠𝑛 ) ∈ R𝑛+ ,

where the machine allocation decision 𝒙(𝝐) is formulated as a policy affected by the machine failure disruption 𝝐 ∈ 𝜩. This is because the choice of machine allocation can be affected by the machine failure disruptions. Since 𝝐̃ is uncertain, this constraint (8) has to hold for all 𝝐 ∈ 𝜩 such that all possible disruption scenarios are covered. When dealing with machines, the decision maker has to consider about the usage cost of machines in practice. It is assumed that there is a cost budget constraint on the machine usage, where 𝑐𝑖𝑟 is the cost of using machine 𝑅𝑟 to process job 𝑖, and 𝐵 is the total available machine budget. The constraint is formulated as ∑ ∑ 𝑐𝑖𝑟 𝑥𝑖𝑟 (𝝐) ≤ 𝐵, ∀𝝐 ∈ 𝜩. (9) 𝑖∈[𝑛] 𝑟∈R𝑖

A scheduling plan consists of both the machine allocation and the job start time. The allocation has been represented by 𝒙. As a result, 4

(10)

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

is assumed that the part types are independent of each other. That is to say, the jobs across different part types do not possess any precedence relationship. A job pair is denoted as (𝑖, 𝑗) ∶= 𝑖 → 𝑗 if job 𝑖 precedes job 𝑗. The set of precedence pairs among all jobs is denoted as the set P. The implication of the precedence relationship is reflected by the constraint in the start time. The start time of successor job should be no earlier than the completion time of the predecessor job. Mathematically, it is expressed as, ∑ 𝑥𝑖𝑟 (𝝐)𝑑̃𝑖𝑟 , ∀(𝑖, 𝑗) ∈ P, ∀𝝐 ∈ 𝜩, 𝑠𝑗 (𝝐) ≥ 𝑠𝑖 (𝝐) + (11)

constraint (16) will thus become active by controlling the start time of job 𝑖 to be later than job 𝑗 by at least the period of 𝑑̃𝑗𝑟 . In the case that jobs 𝑖, 𝑗 are processed by different machines, i.e. 𝑤𝑖𝑗 = 0, there is no restriction regarding the start time. In this case, constraints (15) and (16) will just become void. The ultimate goal in schedule planning is to decide on the optimal schedule that can meet certain manufacturing production due-dates, which is formulated as ∑ 𝑠𝑖 (𝝐) + 𝑥𝑖𝑟 (𝝐)𝑑̃𝑖𝑟 ≤ 𝑡𝑖 , 𝑖 ∈ [𝑛], ∀𝝐 ∈ 𝜩, (17)

where the job completion time is computed as the sum of the start time 𝑠𝑖 and the job duration 𝑑̃𝑖𝑟 with respect to the allocated machine 𝑅𝑟 . In the job scheduling, it is assumed that pre-emptions are not allowed. Hence, once a machine is allocated to a job, it is considered as not available for other jobs, until it is released from its current job. Controlled by (11), no machine clash between a pair of jobs will happen if they are required to comply with the precedence relationships. For activity pairs (𝑖, 𝑗) with no precedence relationship, machine clash can be avoided by imposing the constraint that either, job 𝑖 and job 𝑗 do not share the same machine, or that they share the same machine with no overlapping period. To facilitate this, set Q is defined to contain all the pairs of jobs without precedence relationship, and the following set of binary variables 𝑤𝑖𝑗 are defined as { 1, if the same machine is assigned to job 𝑖 and 𝑗, 𝑤𝑖𝑗 = (12) 0, otherwise.

where 𝑡𝑖 stands for the planned due-date of job 𝑖. A common example of due-date 𝑡𝑛 is the production makespan, namely completion time of the entire manufacturing process. To present a more general setting, there can be due-date requirement 𝒕 = (𝑡1 , … , 𝑡𝑖 , … , 𝑡𝑛 ) for each job or only a subset of the jobs.

𝑟∈R𝑖

𝑟∈R𝑖

3. Robust scheduling optimization A robust scheduling optimization model will be formulated in this section to deal with the uncertain machine failure disruptions. 3.1. Uncertain machine failure disruptions Given the set of jobs [𝑛] and the set of machines R, each of the machines may be subject to random failure disruption at the start of any job. The number of possible scenarios of disruption is thus exponential. It is not plausible to literally test each of the possible scenarios and determine whether it returns a feasible schedule. As a result, the strategy we propose is to set up an imaginary threshold scenario of machine failure disruption, bounded by which the due-dates can be achieved as well as possible. Define 𝜶 𝑛×|R| as the threshold scenario that is selected to define the range of machine failure disruptions, 𝟎 ≤ 𝝐̃ ≤ 𝜶, over which the duedates are guaranteed to be achieved. All the components in 𝜶 take values either 0 or 1. Define the set

For all the job pairs (𝑖, 𝑗) in Q, the machine allocation schedule has to be controlled by the following constraint, 𝑥𝑖𝑟 (𝝐) + 𝑥𝑗𝑟 (𝝐) − 1 ≤ 𝑤𝑖𝑗 (𝝐), ∀(𝑖, 𝑗) ∈ Q, 𝑅𝑟 ∈(R𝑖 ∩ R𝑗 ), ∀𝝐 ∈ 𝜩.

(13)

First, the resource set is defined by (R𝑖 ∩ R𝑗 ), since the machines that are subject to potential allocation clash between jobs 𝑖 and 𝑗 are only those which are capable of processing both jobs. Next, in the case that 𝑤𝑖𝑗 = 1, jobs 𝑖 and 𝑗 are allowed to share the same machine by definition, and both 𝑥𝑖𝑟 and 𝑥𝑗𝑟 can take value 1 in (13). On the contrary, if 𝑤𝑖𝑗 = 0, (13) is constrained to guarantee that at most one of 𝑥𝑖𝑟 and 𝑥𝑗𝑟 can take value 1 for the same machine 𝑅𝑟 . In this way, job 𝑖 and job 𝑗 are controlled not to share the same machine. As for the job start time, the only possibility of assigning the same machine to two different jobs without clash is to use the machine with no overlapping time period. In other words, one job has to wait for the machine until the other job finishes using it. In this case, there exists an order that indicates which job is able to be processed by the machine first. The binary variables 𝑧𝑖𝑗 are therein defined to indicate the order, { 1, if the machine is assigned to job 𝑗 before job 𝑖, 𝑧𝑖𝑗 = (14) 0, if the machine is assigned to job 𝑖 before job 𝑗.

L(𝜶) ∶= {𝝐 ∶ 𝝐 ≤ 𝜶, ∀𝝐 ∈ 𝜩}

to represent the set of disruption scenarios which are guaranteed to be feasible. The larger the set L(𝜶) is, the wider the range of uncertain 𝝐̃ that can be protected over. As a consequence. the objective is to find the optimal 𝜶 which maximizes the probability of random 𝝐̃ falling into the set L(𝜶), namely 𝝐 ≤ 𝜶). Since each element in 𝝐̃ is binary, there are ∑𝑛 P(̃ in total 𝐻 = 2 𝑖=1 |R𝑖 | possible scenarios of 𝝐̃ . Denote the ℎth scenario as 𝝐 ℎ , where ℎ = 1, … , 𝐻. The (𝑖, 𝑟)th element value of scenario 𝝐 ℎ is represented as 𝝐 ℎ𝑖𝑟 . Further, the probability function can be evaluated ∑ ℎ ℎ as 𝐻 ℎ=1 P(𝝐 )1(𝝐 ≤ 𝜶). Many existing scenario optimization models (Dembo, 1991) require specifying the occurrence probabilities of each individual scenario as a prior. This can be very difficult since each scenario involves the status of multiple machines and multiple jobs. It is not easy to compute the joint probability and there are exponential number of them to be estimated. Motivated by this, we propose a distributionally robust optimization model (Goh & Sim, 2010; Mak et al., 2014) which makes use of only partial distributional information as specified by a probability measure defined as

The order of machine allocation will affect the start time which is formulated as: ∑ 𝑠𝑖 (𝝐) − 𝑠𝑗 (𝝐) + 𝑥𝑖𝑟 (𝝐)𝑑̃𝑖𝑟 ≤ 𝑀(1 − 𝑤𝑖𝑗 (𝝐)) 𝑟∈(R𝑖 ∩R𝑗 ) (15) + 𝑀𝑧𝑖𝑗 (𝝐), ∀(𝑖, 𝑗) ∈ Q, ∀𝝐 ∈ 𝜩, 𝑠𝑗 (𝝐) − 𝑠𝑖 (𝝐) +

∑

F(̃𝝐 ) ∶= {

𝑥𝑗𝑟 (𝝐)𝑑̃𝑗𝑟 ≤ 𝑀(1 − 𝑤𝑖𝑗 (𝝐))+

𝐻 ∑ ℎ=1

(16)

𝑟∈(R𝑖 ∩R𝑗 )

(18)

P(𝝐 ℎ ) = 1,

𝐻 ∑

𝝐 ℎ𝑖𝑟 P(𝝐 ℎ ) ≤ 𝑃̂𝑖𝑟 , P(𝝐 ℎ ) ≥ 0}.

(19)

ℎ=1

In this formulation, only the marginal machine failure probability, denoted as 𝑃̂𝑖𝑟 , ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 , is required. Each 𝑃̂𝑖𝑟 represents the probability of machine failure while processing a corresponding job. The efforts for estimating this partial distributional information, although significant, are less burdensome. The marginal machine failure probability 𝑃̂𝑖𝑟 can be computed through reliability analysis. Various studies (Fernandez & Vazquez, 2012; Jia et al., 2018) have discussed

𝑀(1 − 𝑧𝑖𝑗 (𝝐)), ∀(𝑖, 𝑗) ∈ Q, ∀𝝐 ∈ 𝜩, where 𝑀 is an arbitrarily large number. When the same machine is allocated to job 𝑖 first and then allocated to job 𝑗, constraint (15) takes into effect by controlling the start time of job 𝑗 to be later than job 𝑖 by at least a period of 𝑑̃𝑖𝑟 . Similarly, when the same machine is allocated to job 𝑗 first and then allocated to job 𝑖, 5

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

Thus, it can be obtained that 𝛿𝝐 (1) + (1 − 𝛿)𝝐 (2) ≤ 𝜶 and A is proved to be convex. As B: suppose there exist any two vectors 𝜷 (1) , 𝜷 (2) ∈ B. As 𝜷 (1) ∈ B, it means that 𝜷 (1) is not dominated by 𝜶. In other words, at least one element of the vector 𝜷 (1) is strictly greater than the corresponding element of 𝜶. Without loss of generality, we assume that there exists some 𝑖 ∈ [𝑛], 𝑟 ∈ R𝑖 such that 𝜷 (1) 𝑖𝑟 > 𝜶 𝑖𝑟 . Given that 𝜶 𝑖𝑟 ∈ {0, 1} and (1) 0 ≤ 𝜷 (1) 𝑖𝑟 ≤ 1, it is therefore easy to find that 𝜶 𝑖𝑟 = 0, 𝜷 𝑖𝑟 > 0 in this case. As a consequence, it can be obtained that 𝛿𝜷 (1) > 𝛿𝜶 𝑖𝑟 as well as 𝑖𝑟 ≥ (1 − 𝛿)𝜶 for any 𝛿 ∈ (0, 1]. (1 − 𝛿)𝜷 (2) 𝑖𝑟 𝑖𝑟 Similarly, there will exist some 𝑖∗ ∈ [𝑛], 𝑟∗ ∈ R𝑖 such that 𝜷 (2) > 𝑖∗ 𝑟∗ 𝜶 𝑖∗ 𝑟∗ . By analogy, it can be obtained that 𝜂𝜷 (2) > 𝜂𝜶 ∗ 𝑟∗ as well as ∗ ∗ 𝑖 𝑖 𝑟 (1 − 𝜂)𝜷 (1) ≥ (1 − 𝜂)𝜶 𝑖∗ 𝑟∗ for any 𝜂 ∈ (0, 1]. 𝑖∗ 𝑟∗ (1) (1) If (𝑖, 𝑟) = (𝑖∗ , 𝑟∗ ), 𝜷 𝑖𝑟 and 𝜷 𝑖∗ 𝑟∗ refer to the same entry, which can be just denoted as (𝑖, 𝑟) without loss of generality. Based on the (2) relationships derived, it can be obtained that 𝛿𝜷 (1) 𝑖𝑟 + (1 − 𝛿)𝜷 𝑖𝑟 > (2) (1) 𝜶 𝑖𝑟 , ∀𝛿 ∈ (0, 1] and 𝛿𝜷 𝑖𝑟 + (1 − 𝛿)𝜷 𝑖𝑟 ≥ 𝜶 𝑖𝑟 with 𝛿 = 0. That is to say, 𝛿𝜷 (1) + (1 − 𝛿)𝜷 (2) ∈ B. (1) (2) If (𝑖, 𝑟) ≠ (𝑖∗ , 𝑟∗ ), at least one of 𝛿𝜷 (1) 𝑖𝑟 + (1 − 𝛿)𝜷 𝑖𝑟 > 𝜶 𝑖𝑟 or 𝜂𝜷 𝑖∗ 𝑟∗ + (1 − (2) 𝜂)𝜷 𝑖∗ 𝑟∗ > 𝜶[𝑘] will hold with 𝜂 ∈ (0, 1], 𝜶[𝑘] denotes the 𝑘th element (2) of the vector 𝜶. Similarly, at least one of 𝛿𝜷 (1) 𝑖𝑟 + (1 − 𝛿)𝜷 𝑖𝑟 ≥ 𝜶 𝑖𝑟 or (2) 𝜂𝜷 (1) + (1 − 𝜂)𝜷 ≥ 𝜶[𝑘] will hold with 𝜂 = 0. It can be obtained that 𝑖∗ 𝑟∗ 𝑖∗ 𝑟∗ 𝛿𝜷 (1) + (1 − 𝛿)𝜷 (2) ∈ B. In both cases, B is proved to be convex. As such, it can be concluded that both A and B are convex. This completes the proof. □

about ways in conducting the reliability analysis. Since this is not the main focus of this work, it is assumed that 𝑃̂𝑖𝑟 is known in the model formulation by now. In practice, the manager may adopt the existing methods developed in the area of reliability to estimate the machine failure probability first and then apply our proposed model. As a proper probability measure, the probability of each scenario P(𝝐 ℎ ) should be non-negative and the sum of the probabilities of all scenarios should be one. Under uncertainty, the decisions are chosen to optimize the worstcase objective value among all possible distributions within the specified distribution family (19). The proposed scheduling model can therefore be formulated as a max–min problem as follows, max inf 𝜶,𝒙,𝒔 𝝐̃

𝐻 ∑

s.t.

ℎ=1 𝐻 ∑

𝐻 ∑

P(𝝐 ℎ )1(𝝐 ℎ ≤ 𝜶)

ℎ=1

P(𝝐 ℎ ) = 1,

(20)

𝝐 ℎ𝑖𝑟 P(𝝐 ℎ ) ≤ 𝑃̂𝑖𝑟 , ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 ,

ℎ=1 ℎ

P(𝝐 ) ≥ 0, ∀ℎ, or it can be written compactly as max

𝐻 ∑

inf

𝜶,𝒙,𝒔 P(̃𝝐 )∈F(̃𝝐 )

P(𝝐 ℎ )1(𝝐 ℎ ≤ 𝜶).

(21)

ℎ=1

The dual model of this problem can be formulated as max sup 𝜙 + 𝜶,𝒙,𝒔 𝜙,𝜃

𝑖𝑟

s.t.

𝜙+

∑∑

𝑛 ∑ ∑

𝑃̂𝑖𝑟 𝜃𝑖𝑟

𝝐 ℎ𝑖𝑟 𝜃𝑖𝑟 ≤ 1(𝝐 ℎ ≤ 𝜶), ∀ℎ,

(22)

𝑟

𝑖

Proposition 2. The optimal solution to problem (24) can be equivalently computed by solving the following problem, ∑∑ max sup 𝜙 + 𝑃̂𝑖𝑟 𝜃𝑖𝑟

𝑖=1 𝑟∈R𝑖

𝜶,𝒙,𝒔 𝜙,𝜃

𝝐 ℎ𝑖𝑟 ∈ {0, 1}, ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 , ∀ℎ, 𝜃𝑖𝑟 ≤ 0, ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 ,

− 𝑀 ⋅ 𝜶 𝑖𝑟 ≤ 𝑦𝑖𝑟 − 𝜃𝑖𝑟 ≤ 𝑀 ⋅ 𝜶 𝑖𝑟 , ∀𝑖, 𝑟, 𝜃𝑖𝑟 ≤ 0, ∀𝑖, 𝑟. Proof. In order to guarantee that the condition holds for all 𝝐 ℎ , it is sufficient just to make sure the worst case is satisfied. As a consequence, the first constraint in (24) can be equivalently formulated as ∑∑ sup {𝜙 + 𝝐 ℎ𝑖𝑟 𝜃𝑖𝑟 } ≤ 1, 𝝐 ℎ𝑖𝑟 ∈ [0, 1], ∀ℎ, (28)

The formulation can thus be re-formulated as: 𝑛 ∑ ∑ max sup 𝜙 + 𝑃̂𝑖𝑟 𝜃𝑖𝑟 𝜶,𝒙,𝒔 𝜙,𝜃

𝑖𝑟

𝜙+

𝑛 ∑ ∑ 𝑖=1 𝑟 𝑛 ∑ ∑

𝑖=1 𝑟∈R𝑖

𝝐 ℎ𝑖𝑟 𝜃𝑖𝑟

𝝐 ℎ ∶𝝐 ℎ ≤𝜶

ℎ

≤ 1, 𝝐 ≤ 𝜶, ∀ℎ,

𝝐 ℎ𝑖𝑟

∈ [0, 1],

𝝐 ℎ ∶𝝐 ℎ ≰𝜶

𝑖=1 𝑟

Proposition 1. It can be proved that both A = {𝝐 ∶ 𝝐 ≤ 𝜶, 𝝐 𝑖𝑟 ∈ [0, 1], ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 } and B = {𝝐 ∶ 𝝐 ≰ 𝜶, 𝝐 𝑖𝑟 ∈ [0, 1], ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 } are convex. Proof. As for the convexity of A, suppose there exist any two vectors 𝝐 (1) , 𝝐 (2) ∈ A. That is to say, the following relationships hold: 𝝐

≤ 𝜶.

(25)

As a consequence, for any 𝛿 ∈ [0, 1], it is true that 𝛿𝝐 (1) ≤ 𝛿𝜶, (1 − 𝛿)𝝐 (2) ≤ (1 − 𝛿)𝜶.

𝑟

𝑖

(29)

𝑟

Since 𝜃𝑖𝑟 ≤ 0, the worst case happens when 𝝐 ℎ takes the least number of non-zero values in its elements. It is easy to understand this from a concrete example. For instance, in R3 , if 𝜶 = [0 0 1]𝑇 , then A and B can be defined as

𝜃𝑖𝑟 ≤ 0, 𝜙 ∈ R.

(2)

𝑖

and the second constraint can be equivalently formulated as ∑∑ sup {𝜙 + 𝝐 ℎ𝑖𝑟 𝜃𝑖𝑟 } ≤ 0, 𝝐 ℎ𝑖𝑟 ∈ [0, 1], ∀ℎ.

(24)

𝝐 ℎ𝑖𝑟 𝜃𝑖𝑟 ≤ 0, 𝝐 ℎ ≰ 𝜶, ∀ℎ, 𝝐 ℎ𝑖𝑟 ∈ [0, 1],

𝝐 (1) ≤ 𝜶,

(27)

− 𝑀 ⋅ (1 − 𝜶 𝑖𝑟 ) ≤ 𝑦𝑖𝑟 ≤ 0, ∀𝑖, 𝑟,

Since 𝝐 ℎ takes value only 0 or 1, the constraint of 𝝐 ℎ𝑖𝑟 ∈ {0, 1} can be relaxed to [0, 1]. This will not affect the optimal solution. In fact, the optimal solution will lie on the corner point or the boundary of the convex hull. The indicator function is defined as { 1, if 𝝐 ℎ ≤ 𝜶, 1(𝝐 ℎ ≤ 𝜶) = (23) 0, if 𝝐 ℎ ≰ 𝜶.

𝜙+

𝑟

𝜙 + 𝑦𝑖𝑟 ≤ 𝜶 𝑖𝑟 , ∀𝑖, 𝑟,

𝜙 ∈ R.

s.t.

𝑖

𝑖𝑟

s.t. 𝜙 ≤ 1,

⎛⎡0⎤ ⎡0⎤⎞ A = ⎜⎢0⎥ , ⎢0⎥⎟ , ⎜⎢ ⎥ ⎢ ⎥⎟ ⎝⎣0⎦ ⎣1⎦⎠

(30)

⎛⎡0⎤ ⎡1⎤ ⎡0⎤ ⎡1⎤ ⎡1⎤ ⎡1⎤⎞ B = ⎜⎢1⎥ , ⎢0⎥ , ⎢1⎥ , ⎢1⎥ , ⎢0⎥ , ⎢1⎥⎟ . ⎜⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎟ ⎝⎣0⎦ ⎣0⎦ ⎣1⎦ ⎣0⎦ ⎣1⎦ ⎣1⎦⎠

(31)

As can be seen, the worst 𝝐 ℎ ∈ A will be [0 0 0]𝑇 and the worst ∈ B will be one of [0 1 0]𝑇 or [1 0 0]𝑇 . Now we explore the condition in the general situation. 𝝐ℎ

(26) 6

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

∑ For a general vector 𝜶, suppose that out of the total 𝑖 |R𝑖 | ele∑ ′ ′ ments, 𝑚 of them are 0 𝑠 and 𝑘 of them are 1 𝑠 (𝑚 + 𝑘 = 𝑖 |R𝑖 |). The index of the 𝑚 zero elements in 𝜶 is collated as set M. Mathematically, it is represented as M = {(𝑖, 𝑟) ∶ 𝜶 𝑖𝑟 = 0, 𝑖 ∈ [𝑛], 𝑟 ∈ R𝑖 }. Similarly, the index of entries with value 1 is denoted in the set K = {(𝑖, 𝑟) ∶ 𝜶 𝑖𝑟 = 1, 𝑖 ∈ [𝑛], 𝑟 ∈ R𝑖 }. In order for a vector 𝝐 ℎ to be in the set A, it should satisfy the condition that 𝝐 ℎ𝑖𝑟 = 0 for all (𝑖, 𝑟) ∈ M. For other element values 𝝐 ℎ𝑖𝑟 , (𝑖, 𝑟) ∈ K, they can have the value of either 0 or 1. As a result, the size of set A is 2𝑘 . Among all 𝝐 ℎ ∈ A in (28), the scenario 𝟎𝑇 returns ∑ ∑ the maximum value of 𝜙+ 𝑖 𝑟 𝝐 ℎ𝑖𝑟 𝜃𝑖𝑟 and it leads to the condition that 𝜙 ≤ 1. On the other hand, the number of elements in the set B will then be 2𝑚+𝑘 − 2𝑘 . This is the case since 𝝐 ℎ will either lie in A or B but not both. In order to satisfy the condition that 𝝐 ℎ ≰ 𝜶 for 𝝐 ℎ ∈ B, at least one of the elements 𝝐 ℎ𝑖𝑟 , (𝑖, 𝑟) ∈ M has to take value 1. Under this circumstance, the scenarios with the least number of non-zero values in its vector elements will be the ones with only non-zero value 1 in ( ) one of the indexes (𝑖, 𝑟) ∈ M. There are 𝑚1 possible choices. As a result, the constraint (29) can be re-formulated as 𝜙 + 1 ⋅ 𝜃𝑖𝑟 ≤ 0, if 𝜶 𝑖𝑟 = 0, (𝑖, 𝑟) ∈ M, 𝜙 + 0 ⋅ 𝜃𝑖𝑟 ≤ 1, if 𝜶 𝑖𝑟 = 1, (𝑖, 𝑟) ∈ K.

of ℎ(𝒔,𝒅) , which refers to all start times and processing durations up to but does not include the 𝑖th job, namely, ℎ𝑖(𝒔,𝒅) = {(𝑠1 , 𝑑1 ), … , (𝑠𝑖−1 , 𝑑𝑖−1 )}.

Each admissible sub-sequence ℎ𝑖(𝒔,𝒅) is assumed to be consistent with the set of precedence relationships P. That is to say, (𝑠𝑗 , 𝑑𝑗 ) does not appear before (𝑠𝑖 , 𝑑𝑖 ) in the sequence if 𝑖 → 𝑗. Let H be the set of all admissible sub-sequences observed at the start of schedule planning. The non-anticipativity constraints of start times are 𝑖 𝑠𝑖 (ℎ𝑖(𝒔,𝒅) ) ≥ 𝑠𝑗 (ℎ𝑗(𝒔,𝒅) ), ℎ𝑗(𝒔,𝒅) ∈ 𝐻(𝒔,𝒅) , ∀𝑖 ∈ [𝑛],

{ } | 𝑖 𝐻(𝒔,𝒅) = ℎ𝑗(𝒔,𝒅) ∈ H | ℎ𝑗(𝒔,𝒅) ⊂ ℎ𝑖(𝒔,𝒅) ∈ H . |

− 𝑀 ⋅ 𝜶 𝑖𝑟 ≤ 𝑦𝑖𝑟 − 𝜃𝑖𝑟 ≤ 𝑀 ⋅ 𝜶 𝑖𝑟 , ∀𝑖, 𝑟,

(32)

(33)

(34)

where 𝑀 is an arbitrarily big number. As can be seen from (34), when 𝜶 𝑖𝑟 = 0, the first condition is void and the second condition becomes 0 ≤ 𝑦𝑖𝑟 − 𝜃𝑖𝑟 ≤ 0 which can be further translated as 𝑦𝑖𝑟 − 𝜃𝑖𝑟 = 0. On the contrary, when 𝜶 𝑖𝑟 = 1, the second condition is void and the first condition will lead to 𝑦𝑖𝑟 = 0. This completes the proof. □

(39)

The conventional modeling approach of multi-stage problems is to consider the time horizon. Nevertheless, as the time horizon increases, the decision tree will grow quite bushy with a large number of combinations, given the flexibility of systems. As such, in our formulation, we consider the full horizon scenario 𝝐̃ directly without specifying the history of the process. A potential difficulty is that, when we have split up the scenarios, the non-anticipativity of the decisions may have been lost since they now would include information of the outcomes up to the end of the horizon. In the next section, we enforce non-anticipativity by adding constraints explicitly in the formulation.

All other decision variables are modeled as policies which depend on the past information ℎ(𝒔,𝒅) . The current formulation (8), (9), (11), (13), (15), (16), (17), with non-anticipativity is presented in (39). In the current formulation, the scheduling model is in the format of mixed-integer, two-stage stochastic optimization program, where the set of decision variables themselves are adjustable based on the historical information ℎ(𝒔,𝒅) . The next proposition will show that the optimal solution to the current formulation in fact can be obtained by solving a mixed-integer linear program.

3.2. Non-anticipativity

Proposition 3. The optimal solution of problem (39) can be equivalently computed by solving (42).

The non-anticipativity can be interpreted as ‘‘the decision maker cannot anticipate every possible outcome so our decisions are nonanticipative of future outcomes’’. Simply, the decision maker cannot completely observe or anticipate the knowledge of the future outcomes while making the decision now. Let ℎ be the set of available historical information regarding job start times and processing durations ℎ(𝒔,𝒅) = {(𝑠1 , 𝑑1 ), … , (𝑠𝑛 , 𝑑𝑛 )},

(38)

| ∑ ⎧ ⎫ | 𝑥𝑖𝑟 (ℎ𝑖(𝒔,𝒅) ) = 1, ∀𝑖 ∈ [𝑛], ⎪ ⎪ | | 𝑟∈R𝑖 ⎪ ⎪ | | ∑ ∑ ⎪ ⎪ | 𝑖 ⎪ ⎪ 𝑐𝑖𝑟 𝑥𝑖𝑟 (ℎ(𝒔,𝒅) ) ≤ 𝐵, | | ⎪ ⎪ | 𝑖∈[𝑛] 𝑟∈R𝑖 | ⎪ ⎪ ∑ | 𝑗 𝑖 𝑖 ̃ ⎪ ⎪ | 𝑠𝑗 (ℎ(𝒔,𝒅) ) ≥ 𝑠𝑖 (ℎ(𝒔,𝒅) ) + 𝑥𝑖𝑟 (ℎ(𝒔,𝒅) )𝑑𝑖𝑟 , | ⎪ ⎪ | 𝑟∈R𝑖 | ⎪ ⎪ | ⎪ ⎪ | ∀(𝑖, 𝑗) ∈ P, | ⎪ ⎪ | 𝑗 𝑗 𝑖 𝑖 | 𝑥 (ℎ ) + 𝑥 (ℎ ) − 1 ≤ 𝑤 (ℎ , ℎ ), ⎪ ⎪ 𝑗𝑟 (𝒔,𝒅) 𝑖𝑗 (𝒔,𝒅) | 𝑖𝑟 (𝒔,𝒅) (𝒔,𝒅) ⎪ ⎪ | | ⎪ ⎪ | ∀(𝑖, 𝑗) ∈ Q, 𝑅𝑟 ∈ (R𝑖 ∩ R𝑗 ), | ⎪ ⎪ 𝐻 | ∑ ∑ ⎪ ⎪ | inf P(𝝐 ℎ )1(𝝐 ℎ ≤ 𝜶) | 𝑠𝑖 (ℎ𝑖(𝒔,𝒅) ) − 𝑠𝑗 (ℎ𝑗(𝒔,𝒅) ) + 𝑥𝑖𝑟 (ℎ𝑖(𝒔,𝒅) )𝑑̃𝑖𝑟 ⎬ ⎨max | 𝜶,𝒙,𝒔 P(̃𝝐 )∈F(̃𝝐 ) ⎪ ⎪ ℎ=1 | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) | ⎪ ⎪ | 𝑗 𝑖 | ⎪ ⎪ | ≤ 𝑀(1 − 𝑤𝑖𝑗 (ℎ(𝒔,𝒅) , ℎ(𝒔,𝒅) )) ⎪ ⎪ | | ⎪ ⎪ + 𝑀𝑧𝑖𝑗 (ℎ𝑖(𝒔,𝒅) , ℎ𝑗(𝒔,𝒅) ), ∀(𝑖, 𝑗) ∈ Q, | | ⎪ ⎪ | ∑ ⎪ ⎪ | 𝑠 (ℎ𝑗 ) − 𝑠 (ℎ𝑖 ) + 𝑗 𝑥𝑗𝑟 (ℎ(𝒔,𝒅) )𝑑̃𝑗𝑟 ⎪ | 𝑗 (𝒔,𝒅) 𝑖 (𝒔,𝒅) ⎪ | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) | ⎪ ⎪ | ⎪ ⎪ | 𝑗 𝑖 | , ℎ ≤ 𝑀(1 − 𝑤 (ℎ )) ⎪ ⎪ 𝑖𝑗 (𝒔,𝒅) (𝒔,𝒅) | | ⎪ ⎪ | + 𝑀(1 − 𝑧𝑖𝑗 (ℎ𝑖(𝒔,𝒅) , ℎ𝑗(𝒔,𝒅) )), ∀(𝑖, 𝑗) ∈ Q, ⎪ ⎪ | | ⎪ ⎪ | | 𝑠 (ℎ𝑖 ) + ∑ 𝑥 (ℎ𝑖 )𝑑̃ ≤ 𝑡 , 𝑖 ∈ [𝑛]. ⎪ ⎪ | 𝑖 (𝒔,𝒅) 𝑖𝑟 (𝒔,𝒅) 𝑖𝑟 𝑖 ⎪ ⎪ | 𝑟∈R𝑖 | ⎩ ⎭

Define 𝑦𝑖𝑟 = 𝜃𝑖𝑟 (1 − 𝜶 𝑖𝑟 ). Clearly, 𝑦𝑖𝑟 = 𝜃𝑖𝑟 if 𝜶 𝑖𝑟 = 0 and 𝑦𝑖𝑟 = 0 if 𝜶 𝑖𝑟 = 1. As such, the above constraint can be further re-formulated as − 𝑀 ⋅ (1 − 𝜶 𝑖𝑟 ) ≤ 𝑦𝑖𝑟 ≤ 0, ∀𝑖, 𝑟,

(37)

where

Based on the definition of M, (32) can be equivalently written as 𝜙 + 𝜃𝑖𝑟 (1 − 𝜶 𝑖𝑟 ) ≤ 𝜶 𝑖𝑟 , ∀𝑖, 𝑟.

(36)

Proof. As discussed earlier, the support of 𝝐̃ 𝑖𝑟 ∈ {0, 1} can be relaxed to [0, 1]. In order to guarantee that the schedule is feasible under random failure scenario 𝝐̃ , all the conditions affected by 𝝐̃ are constrained to hold for all possible values taken in the entire support [0, 1]. This is known as a semi-infinite linear programming problem, which is not tractable. As a result, instead of optimizing all the scenarios, it is sufficient to guarantee the worst case is satisfied. If a feasible solution can be obtained for this worst case, it should also work for all the other cases. In this way, a dynamic constraint is transformed to a deterministic one. Based on our definition, 𝜶 is the worst-case scenario in the formulation.

(35)

where (𝑠𝑖 , 𝑑𝑖 ) denotes respectively the start time and processing duration of the 𝑖th job in the sequence. The information of start times and processing durations may not be available for all jobs 𝑖 ∈ [𝑛] in practice. Let ℎ𝑖(𝒔,𝒅) denote a sub-sequence 7

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094 Table 1 System specifications of the studied FMS.

In particular, sup 𝑥𝑖𝑟 𝑑̃𝑖𝑟 = sup 𝑥𝑖𝑟 (𝑑𝑖𝑟 + 𝝐̃ 𝑖𝑟 𝛥𝑑𝑖𝑟 ) 𝝐̃ 𝑖𝑟

𝝐̃ 𝑖𝑟

(40)

= 𝑥𝑖𝑟 𝑑𝑖𝑟 + 𝑥𝑖𝑟 𝛥𝑑𝑖𝑟 sup 𝝐̃ 𝑖𝑟 For a maximization problem, the worst while still feasible result provides a lower bound to the optimal solution. On the other hand, instead of controlling the results to be feasible for dynamic constraints under infinitely many possible realizations, constraining them to be feasible for one particular set of realizations is considered as a relaxation of the original problem. Once the realization of the random 𝝐̃ is given, the dynamic scheduling problem will become deterministic since all the second-stage decision variables can be decided. In particular, 𝜶 is one possible value of 𝝐̃ , which results in a relaxation of the original problem, and the returned results should serve as an upper bound to the optimal solution. To obtain a tractable formulation, we replace 𝑥𝑖𝑟 𝜶 𝑖𝑟 by 𝑞𝑖𝑟 and introduce the following constraints,

Machines

Jobs

R𝑖

𝑑𝑖𝑟

𝛥𝑑𝑖𝑟

8

𝑅1 𝑅2 𝑅3

1 2 3 4 5

{1, 2, 3} {1, 2, 3} {1, 2, 3} {1, 2, 3} {1, 2, 3}

{2, 3, 3} {3, 2, 3} {2, 4, 3} {2, 3, 3} {1, 1, 2}

{10, 8, 8} {8, 9, 7} {10, 6, 9} {11, 8, 9} {10, 12, 7}

Intel Core i5 processor and 8-gigabyte random-access memory (RAM). This configuration is widely available for most desktop computers nowadays. ILOG CPLEX is used to solve the scheduling model. The random network generator 𝑅𝑎𝑛𝐺𝑒𝑛1 (Demeulemeester, Vanhoucke, & Herroelen, 2003) is used to generate diverse network topologies with different number of jobs and machines. An 𝑀 = 10,000 is used in all simulation studies in this section.

𝑞𝑖𝑟 ≤ 𝜶 𝑖𝑟 , ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 , 𝑞𝑖𝑟 ≤ 𝑀𝑥𝑖𝑟 , ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 ,

𝑡6

4.1. Simulation setup (41)

The designed case studies are tested on the full system. A full system refers to the case that all jobs are choice jobs and each machine is capable of performing all jobs. The reason for using the full system is to exhaust the maximum number of possible combinations of sub-problems and solution space. The FMS system plotted in Fig. 3 will be discussed in this case study. Job 1 to job 5 are the actual manufacturing jobs. Job 0 and job 6 are dummy jobs artificially created to denote the start and end of the manufacturing process as a common modeling practice. Both of them have duration zero and no machine is allocated to them. The set of precedence pairs is P = {(0, 1), (1, 2), (2, 3), (2, 4), (3, 5), (4, 5), (5, 6)} and the set of non-precedence pairs is Q = {(3, 4)}. The system specifications are presented in Table 1. There are five machining jobs 1, … , 5 and three manufacturing machines 𝑅1 , 𝑅2 , 𝑅3 . Each of the machines is capable of performing all jobs. As such, R𝑖 = {1, 2, 3}, ∀𝑖 = 1, … , 5. The durations of processing different jobs using different machines are listed as 𝑑𝑖𝑟 and machine replenishment time in the case of failure is listed under 𝛥𝑑𝑖𝑟 .

𝑞𝑖𝑟 ≥ 𝜶 𝑖𝑟 − 𝑀(1 − 𝑥𝑖𝑟 ), ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 . All in all, the solution to (42) is both a lower bound and an upper bound to the optimal solution of the original problem. This completes the proof. □ As can be seen, apart from the motivation of solving an important industrial problem, introducing the setup of threshold failure scenario 𝜶 also has another important practical implication. That is, it manages to largely release the computational complexity existing in the dynamic multi-stage scheduling problem. The original dynamic scheduling problem under uncertainty is known as the semi-infinite linear programming problem, which has been shown to be NP-hard by various studies. This above proposition shows that its optimal solution can be equivalently computed from solving a deterministic scheduling problem, which can be solved using a commercial solver like IBM ILOG CPLEX. | ∑ ⎧ ⎫ | 𝑥𝑖𝑟 = 1, ∀𝑖 ∈ [𝑛], | ⎪ ⎪ | 𝑟∈R𝑖 | ⎪ ⎪ | ∑ ∑ ⎪ ⎪ | 𝑐 𝑥 ≤ 𝐵, | 𝑖𝑟 𝑖𝑟 ⎪ ⎪ | 𝑖∈[𝑛] 𝑟∈R | ⎪ ⎪ 𝑖 | ∑ ⎪ ⎪ | 𝑠 ≥ 𝑠 + (𝑥 𝑑 + 𝑞 𝛥𝑑 ), ∀(𝑖, 𝑗) ∈ P , | 𝑖 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖 ⎪ ⎪ | 𝑗 𝑟∈R𝑖 | ⎪ ⎪ | ⎪ ⎪ | | 𝑥 + 𝑥 − 1 ≤ 𝑤 , ∀(𝑖, 𝑗) ∈ Q, ⎪ ⎪ 𝑗𝑟 𝑖𝑗 | 𝑖𝑟 | ⎪ ⎪ | 𝑅𝑟 ∈ (R𝑖 ∩ R𝑗 ), ⎪ ⎪ | | ⎪ ⎪ ∑ | | 𝑠𝑖 − 𝑠𝑗 + ⎪ ⎪ (𝑥𝑖𝑟 𝑑𝑖𝑟 + 𝑞𝑖𝑟 𝛥𝑑𝑖𝑟 ) 𝐻 | ∑ ⎪ ⎪ | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) inf P(𝝐 ℎ )1(𝝐 ℎ ≤ 𝜶) | ⎨max ⎬ | 𝜶,𝒙,𝒔 P(̃𝝐 )∈F(̃𝝐 ) ⎪ ⎪ ℎ=1 | ≤ 𝑀(1 − 𝑤𝑖𝑗 ) + 𝑀𝑧𝑖𝑗 , ∀(𝑖, 𝑗) ∈ Q, | ⎪ ⎪ | ∑ | ⎪ ⎪ | 𝑠𝑗 − 𝑠𝑖 + (𝑥 𝑑 + 𝑞 𝛥𝑑 ) 𝑗𝑟 𝑗𝑟 𝑗𝑟 𝑗𝑟 ⎪ ⎪ | | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) ⎪ ⎪ | | ⎪ ⎪ | ≤ 𝑀(1 − 𝑤 ) + 𝑀(1 − 𝑧 ), ∀(𝑖, 𝑗) ∈ Q, ⎪ ⎪ 𝑖𝑗 𝑖𝑗 | | ⎪ ⎪ ∑ | | 𝑠𝑖 + ⎪ ⎪ (𝑥 𝑑 + 𝑞 𝛥𝑑 ) ≤ 𝑡 , 𝑖 ∈ [𝑛], 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖 | ⎪ ⎪ | 𝑟∈R𝑖 | ⎪ ⎪ | | 𝑞𝑖𝑟 ≤ 𝜶 𝑖𝑟 , ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 , ⎪ ⎪ | ⎪ ⎪ | | ⎪ ⎪ | 𝑞𝑖𝑟 ≤ 𝑀𝑥𝑖𝑟 , ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R𝑖 , | ⎪ ⎪ | 𝑞 ≥ 𝜶 − 𝑀(1 − 𝑥 ), ∀𝑖 ∈ [𝑛], ∀𝑟 ∈ R . ⎭ ⎩ | 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖

4.2. Results discussion The resulting schedule from solving the proposed scheduling model is 𝒔 = (0, 0, 2, 5, 4, 7, 8). The chosen machine for jobs 1 to 5 are 𝑅1 , 𝑅2 , 𝑅1 , 𝑅2 , 𝑅2 respectively, which are highlighted in green in the plot with the indicated job duration. The optimal threshold machine failure scenario 𝜶 = (0, 0, 1). The feasible set of machine failure scenarios is L = {(0, 0, 0), (0, 0, 1)}. The threshold scenario (0, 0, 1) refers to the failure of machine 𝑅3 . As can be seen, 𝑅3 is not chosen for any job. In fact, 𝑅3 takes the longest processing durations for most of the jobs. As such, the failure of this machine is ‘feasible’ in the sense that the planned due-date can still be achieved. 4.3. Computational complexity The CPU times for solving the proposed model in various FMS systems with a wide range of sizes are reported in Table 2. The activity durations 𝒅̃ 𝒊𝒓 are randomly generated using RanGen and the machine replenishment time in the case of failure 𝜟𝒅 𝒊𝒓 are set as 5𝒅̃ 𝒊𝒓 . As can be seen, for the FMS system with 100 jobs and 10 machines, it takes only 0.47 s to solve. For a big FMS system with 1000 jobs and 50 machines, it takes about 60.46 s (≈1 min) to solve. For completeness, the extremely large cases which are usually rare in practice are also tested. For systems with 2000 jobs, it takes about 3 to 6 min to compute which is considered as quite comparable in solving a FMS scheduling problem. For systems up to 4000 jobs, it takes about 15 min to solve. All in all, the results from the demonstrative case study show that the proposed model indeed has the desired performance as predicted.

(42)

4. Simulation studies In this section, one demonstrative case study will be conducted to test the performance of the proposed scheduling model. All computational tests are carried out on a digital computer equipped with 8

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

Fig. 3. AON representation of the simulated 5-job FMS. Table 3 Machine profile parameters.

Table 2 Computational complexity. Number of jobs

Number of machines

CPU time (seconds)

Machine

Capability

Duration (seconds)

100 200 400 600 800 1000 1000 2000 2000 4000 5000

10 10 20 20 20 30 50 30 50 30 40

0.47 1.21 5.24 9.97 16.07 33.82 60.46 212.42 370.85 897.62 >1000

𝑅1

𝑣11 , 𝑣12

4.76, 3.73

𝑅2

𝑣21

3.89

𝑅3

𝑣21 , 𝑣31

3.48, 4.26

𝑅4

𝑣31 , 𝑣41

5.68, 4.46

𝑅5

𝑣31

4.37

𝑅6

𝑣12

4.30

𝑅7

𝑣22

5.27

𝑅8

𝑣22

5.52

5.2. Model comparison Two benchmark approaches in the literature are formulated to compare the performance with our proposed model. They are known as the deterministic model and the sample average model. The deterministic model assumes that the machine disruption is known beforehand as a static scenario 𝝐 ∗ . The model is presented as (43). The deterministic model provides a quick schedule to the decision maker in practice whereas it does not take the uncertainty factor into consideration.

5. Industrial application A real industrial case study is carried out to evaluate the practical applicability of the proposed FMS scheduling model. An application in the stamping industry is selected for this experiment. Stamping includes a number of metal forming processes like coining, punching, piercing, bending. In this case study, input raw materials are metal sheets and output products are various types of voice coil motor (VCM) yokes used in commercial hard-disk drive actuators. An example of VCM yokes is shown in Fig. 4.

⎫ ⎧ | ∑ | 𝑥𝑖𝑟 = 1, ∀𝑖 ∈ [𝑛], ⎪ ⎪ | | 𝑟∈R𝑖 ⎪ ⎪ | ⎪ ⎪ | ∑ ∑ | ⎪ ⎪ | 𝑐 𝑥 ≤ 𝐵, 𝑖𝑟 𝑖𝑟 | ⎪ ⎪ | 𝑖∈[𝑛] 𝑟∈R 𝑖 ⎪ ⎪ | ∑ | ⎪ ⎪ ∗ | 𝑠𝑗 ≥ 𝑠𝑖 + 𝑥 (𝑑 + 𝜖 𝛥𝑑 ), ∀(𝑖, 𝑗) ∈ P , 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖 𝑖𝑟 | ⎪ ⎪ | 𝑟∈R ⎪ ⎪ 𝑖 | | ⎪ ⎪ | 𝑥 + 𝑥 − 1 ≤ 𝑤 , ∀(𝑖, 𝑗) ∈ Q, 𝑅 ∈ (R ∩ R ), 𝑗𝑟 𝑖𝑗 𝑟 𝑖 𝑗 | 𝑖𝑟 ⎪ ⎪ | ⎪ ⎪ ∑ | ∗ min max 𝑡 | ⎨ 𝑥𝑖𝑟 (𝑑𝑖𝑟 + 𝜖𝑖𝑟 𝛥𝑑𝑖𝑟 ) ≤ 𝑀(1 − 𝑤𝑖𝑗 ) + 𝑀𝑧𝑖𝑗 , ⎬ 𝑖 𝑠 − 𝑠𝑗 + 𝑖∈[𝑛] | 𝑖 ⎪ ⎪ | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) | ⎪ ⎪ | | ∀(𝑖, 𝑗) ∈ Q, ⎪ ⎪ | ⎪ ⎪ | ∑ | ⎪ ⎪ ∗ | 𝑠𝑗 − 𝑠𝑖 + 𝑥𝑗𝑟 (𝑑𝑗𝑟 + 𝜖𝑗𝑟 𝛥𝑑𝑗𝑟 ) ≤ 𝑀(1 − 𝑤𝑖𝑗 ) | ⎪ ⎪ | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) ⎪ ⎪ | | ⎪ ⎪ | + 𝑀(1 − 𝑧 ), ∀(𝑖, 𝑗) ∈ Q, | 𝑖𝑗 ⎪ ⎪ | ⎪ ⎪ ∑ | ∗ |𝑠 + ⎪ ⎪ 𝑥𝑖𝑟 (𝑑𝑖𝑟 + 𝜖𝑖𝑟 𝛥𝑑𝑖𝑟 ) ≤ 𝑡𝑖 , 𝑖 ∈ [𝑛]. | 𝑖 | ⎪ ⎪ 𝑟∈R𝑖 | ⎭ ⎩

5.1. System specifications This stamping system can be characterized by the class of FMSs scheduling problem in this paper. Each stamping part type has a predetermined sequence of jobs, where some jobs can be processed by more than one machine and some machines can process more than one type of jobs. Currently, the scheduling decision is made by the subjective judgment in the factory. This studied stamping system consists of two part types 𝜋1 , 𝜋2 and eight stamping machines 𝑅1 − 𝑅8 . Part type 𝜋1 has a job sequence 𝜋1 ≜ 𝑣11 𝑣21 𝑣31 𝑣41 and part type 𝜋2 has a job sequence 𝜋2 ≜ 𝑣12 𝑣22 . Each of the machines has a different efficiency profile. Their performance profile parameters are summarized in Table 3. Jobs 𝑣21 , 𝑣31 , 𝑣12 , 𝑣22 are choice-jobs which can be processed by more than one machine. Represented by AON, the jobs 𝑣11 − 𝑣41 in part type 1 are indexed as jobs 1–4 and the jobs 𝑣12 − 𝑣22 in part type 2 are indexed as jobs 5–6. Job 0 denotes the start of the manufacturing process and job 7 marks the completion of the process. The makespan due-date is set as 𝑡7 = 20.

(43) 9

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

Fig. 4. Example of VCM yokes.

The planned due-dates 𝑡𝑖 are now defined as decision variables, and the objective function max𝑖∈[𝑛] 𝑡𝑖 then represents the makespan of the manufacturing process. (43) yields a simple mixed-integer linear programming formulation.

as 𝑠∗𝑛 . Taking each outcome in the testing set as a deterministic scenario, the resulting makespan can be solved as 𝑠𝑙𝑛 , 𝑙 = 1, … , 𝑁. Then, ∑𝑁 1(𝑠∗𝑛 ≤ 𝑠𝑙𝑛 ) × 100%. (45) 𝜓 = 𝑙=1 𝑁

⎧ ⎫ | ∑ | ⎪ ⎪ 𝑥𝑖𝑟 = 1, ∀𝑖 ∈ [𝑛], | ⎪ ⎪ | 𝑟∈R𝑖 | ⎪ ⎪ | ∑ ∑ | ⎪ ⎪ | 𝑐𝑖𝑟 𝑥𝑖𝑟 ≤ 𝐵, ⎪ ⎪ | | 𝑖∈[𝑛] 𝑟∈R𝑖 ⎪ ⎪ | ∑ | ⎪ ⎪ | 𝑠 ≥ 𝑠 + 𝑥 (𝑑 + 𝝐 𝛥𝑑 ), ∀(𝑖, 𝑗) ∈ P , ∀𝝐 ∈ 𝜩, 𝑖 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖𝑟 𝑖 ⎪ ⎪ | 𝑗 | 𝑟∈R𝑖 ⎪ ⎪ | | ⎪ ⎪ | ⎪ ⎪ | 𝑥𝑖𝑟 + 𝑥𝑗𝑟 − 1 ≤ 𝑤𝑖𝑗 , ∀(𝑖, 𝑗) ∈ Q, 𝑅𝑟 ∈ (R𝑖 ∩ R𝑗 ), | ⎪ ⎪ ∑ | | 𝑠𝑖 − 𝑠𝑗 + ⎪ ⎪ 𝑥𝑖𝑟 (𝑑𝑖𝑟 + 𝝐 𝑖𝑟 𝛥𝑑𝑖𝑟 ) | ⎪ ⎪ | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) 𝑡𝑖 (̃𝝐 )) | ⎨min EP (max ⎬ | 𝑖∈[𝑛] ⎪ ⎪ | ≤ 𝑀(1 − 𝑤 ) + 𝑀𝑧 , 𝑖𝑗 𝑖𝑗 | ⎪ ⎪ | | ⎪ ⎪ ∀(𝑖, 𝑗) ∈ Q, ∀𝝐 ∈ 𝜩, | ⎪ ⎪ | ∑ | ⎪ ⎪ | 𝑠𝑗 − 𝑠𝑖 + 𝑥𝑗𝑟 (𝑑𝑗𝑟 + 𝝐 𝑗𝑟 𝛥𝑑𝑗𝑟 ) | ⎪ ⎪ | 𝑅𝑟 ∈(R𝑖 ∩R𝑗 ) ⎪ ⎪ | | ⎪ ⎪ | ≤ 𝑀(1 − 𝑤 ) + 𝑀(1 − 𝑧 ), | 𝑖𝑗 𝑖𝑗 ⎪ ⎪ | ⎪ ⎪ | ∀(𝑖, 𝑗) ∈ Q, ∀𝝐 ∈ 𝜩, | ⎪ ⎪ | ∑ | ⎪ ⎪ | 𝑥𝑖𝑟 (𝑑𝑖𝑟 + 𝝐 𝑖𝑟 𝛥𝑑𝑖𝑟 ) ≤ 𝑡𝑖 , 𝑖 ∈ [𝑛], ∀𝝐 ∈ 𝜩. ⎪ ⎪ | 𝑠𝑖 + | ⎪ ⎪ 𝑟∈R𝑖 | ⎩ ⎭

Clearly, the higher the success rate is, the better the model performance is in achieving the due-dates. 5.3. Results discussion The comparison of model performance is plotted in Fig. 5. Five cases named as test cases 1, … , 5 are studied. The difference across these cases is the level of uncertainties in the data set of random machine disruptions. From test case 1 to test case 5, the level of variation in the uncertain machine disruption scenarios gradually increases. In particular, test case 1 is chosen as the least uncertain case, which is a deterministic case. In other words, the machine disruption is fixed in test case 1. From test case 2 onwards, the simulated disruption scenarios become more and more random and diverse. The success rate is computed across all three models under the five test cases. As can be seen from the results, when there is no randomness in test case 1, all three models perform effectively as like a deterministic model and all achieve the same success rate. When the randomness in the machine disruptions increases, the success rate of all three models will gradually drop. This is expected as the more uncertain the disruption is, the more difficult it is to achieve the due-dates successfully. Besides, the deterministic model has the worst performance since it does not take the uncertainty into consideration at all. The sample average model outperforms the deterministic model by evaluating the sample average over the uncertain disruptions. The proposed robust optimization model has the best performance as it effectively considers the threshold scenario, in which the due-dates are guaranteed to be achieved. Notably, although the success rates across all models decrease from test case 1 to test case 5, the magnitude of the drop for the robust optimization model is much smaller compared to the other two models. This shows some kind of robustness in the proposed model. The robust optimization model can still maintain a reasonably good success rate in test case 5 whereas the deterministic model and sample average model are comparable only when there is little randomness.

(44) The sample average model considers about the uncertainty of the machine failure scenarios 𝝐̃ . In practice, empirical data of machine disruption scenarios can be obtained from the historical observations. The sample average model then tries to solve the scheduling optimization problem which minimizes the expected makespan of the manufacturing process in the data sample. The model is presented as (44). As for the test, 2𝑁 random failure scenarios 𝝐̃ are generated. The data set is divided into half, where the first half is used for training the models while the second half is used for testing. The performance metric success rate 𝜓 is defined as the probability of achieving the targeted makespan of the manufacturing process. The resulting objective value, namely the makespan, obtained from the respective models is denoted 10

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094

Fig. 5. Performance comparison.

6. Conclusion

Chen, G., Zhang, L., Arinez, J., & Biller, S. (2013). Energy-efficient production systems through schedule-based operations. IEEE Transactions on Automation Science and Engineering, 10, 27–37. Chew, S. F., & Lawley, M. A. (2006). Robust supervisory control for production systems with multiple resource failures. IEEE Transactions on Automation Science and Engineering, 3, 309–323. Chew, S. F., Wang, S., & Lawley, M. A. (2011). Resource failure and blockage control for production systems. International Journal of Computer Integrated Manufacturing, 24, 229–241. Dembo, R. S. (1991). Scenario optimization. Annals of Operations Research, 30, 63–80. Demeulemeester, E., Vanhoucke, M., & Herroelen, W. (2003). Rangen: A random network generator for activity-on-the-node networks. Journal of Scheduling, 6, 17–38. Fernandez, A., & Vazquez, M. (2012). Improved estimation of Weibull parameters considering unreliability uncertainties. IEEE Transactions on Reliability, 61, 32–40. Goh, J., & Sim, M. (2010). Distributionally robust optimization and its tractable approximations. Operations Research, 58, 902–917. Hsieh, F.-S. (2004). Fault-tolerant deadlock avoidance algorithm for assembly processes. IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, 34, 65–79. Hu, H., Zhou, M., & Li, Z. (2011). Supervisor optimization for deadlock resolution in automated manufacturing systems with petri nets. IEEE Transactions on Automation Science and Engineering, 8, 794–804. Janak, S. L., Floudas, C. A., Kallrath, J., & Vormbrock, N. (2006). Production scheduling of a large-scale industrial batch plant. II. Reactive scheduling. Industrial & Engineering Chemistry Research, 45, 8253–8269. Janak, S. L., Lin, X., & Floudas, C. A. (2007). A new robust optimization approach for scheduling under uncertainty: II. Uncertainty with known probability distribution. Computers & Chemical Engineering, 31, 171–195. Jia, X., Nadarajah, S., & Guo, B. (2018). Exact inference on weibull parameters with multiply type-i censored data. IEEE Transactions on Reliability. Jin, F., Zhao, J., Han, Z., & Wang, W. (2018). A joint scheduling method for multiple byproduct gases in steel industry. Control Engineering Practice, 80, 174–184. Klastorin, T. (2004). Project management: Tools and trade-offs. NJ: Wiley Hoboken. Li, Z., & Ierapetritou, M. (2008). Process scheduling under uncertainty: Review and challenges. Computers & Chemical Engineering, 32, 715–727. Li, J., Morrison, J. R., Zhang, M. T., Nakano, M., Biller, S., & Lennartson, B. (2013). Automation in green manufacturing. IEEE Transactions on Automation Science and Engineering, 10, 1–4. Liu, P. H. (2000). A comparative study of three tool replacement/operation sequencing strategies in a flexible manufacturing system. Naval Research Logistics, 47, 479–499. Luo, J. C., Xing, K. Y., Zhou, M. C., Li, X. L., & Wang, X. N. (2017). Scheduling of deadlock and failure-prone automated manufacturing systems via hybrid heuristic search. International Journal of Production Research, 55, 3283–3293. Mak, H.-Y., Rong, Y., & Zhang, J. (2014). Appointment scheduling with limited distributional information. Management Science, 61, 316–334. Pach, C., Berger, T., Sallez, Y., & Trentesaux, D. (2015). Reactive control of overall power consumption in flexible manufacturing systems scheduling: A potential fields model. Control Engineering Practice, 44, 193–208. Pinedo, M. L. (2016). Scheduling: Theory, algorithms, and systems. Springer. Sindicic, I., Bogdan, S., & Petrovic, T. (2011). Resource allocation in free-choice multiple reentrant manufacturing systems based on machine-job incidence matrix. IEEE Transactions on Industrial Informatics, 7, 105–114.

This paper studies the scheduling problem for the flexible manufacturing systems (FMSs) with replenishment, where machine allocations and manufacturing job schedules need to be decided to achieve some planned production due-dates under uncertain machine failure disruptions. A robust scheduling model is proposed. Based on the concept of threshold scenario, the optimal solution to the multi-stage dynamic scheduling problem can be equivalently computed as solving a mixedinteger linear program (MILP) in a reasonable amount of time. The results from the computational experiments and real stamping industrial case study show that the proposed model indeed has a comparable performance in terms of achieving higher and more stable success rate of meeting the due-dates. One possible future work direction is to introduce dynamic machine disruption rate. Currently, the machine disruption rate is assumed to be static and independent of the machine usage in the formulation. This is true when the effect of degradation is not significant. Nevertheless, the more jobs allocated to a machine and the longer the machine is used, the higher the chance that the machine will get faulty. One possibility is to add a feedback control between the usage rate of the machines and their corresponding disruption rate, which as a consequence, links to the machine allocation decisions. In addition, due to reasons like machine degradation, tool wear, the power consumption of machines in processing the job can be uncertain. Further, the duration of processing a job can be affected and become uncertain as well. Another future work direction is to consider about more uncertain factors in the scheduling like the machine cost uncertainty, job processing duration uncertainty. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References Alidaee, B., & Li, H. (2014). Parallel machine selection and job scheduling to minimize sum of machine holding cost, total machine time costs, and total tardiness costs. IEEE Transactions on Automation Science and Engineering, 11, 294–301. Blazewicz, J., Ecker, K. H., Pesch, E., Schmidt, G., & Weglarz, J. (2013). Scheduling computer and manufacturing processes. Springer Science & Business Media. Chen, X., Wei, T., & Hu, S. (2013). Uncertainty-aware household appliance scheduling considering dynamic electricity pricing in smart home. IEEE Transactions on Smart Grid, 4, 932–941. 11

Z. Wang, C.K. Pang and T.S. Ng

Control Engineering Practice 92 (2019) 104094 Wang, B., Xia, X., Meng, H., & Li, T. (2017). Bad-scenario-set robust optimization framework with two objectives for uncertain scheduling systems. IEEE/CAA Journal of Automatica Sinica, 4, 143–153. Yang, F., Gao, K., Simon, I. W., Zhu, Y., & Su, R. (2018). Decomposition methods for manufacturing system scheduling: a survey. IEEE/CAA Journal of Automatica Sinica, 5, 389–400. Yue, H., Xing, K., Hu, H., Wu, W., & Su, H. (2016). Petri-net-based robust supervisory control of automated manufacturing systems. Control Engineering Practice, 54, 176–189.

Wang, S., Chew, S. F., & Lawley, M. A. (2008). Using shared-resource capacity for robust control of failure-prone manufacturing systems. IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, 38, 605–627. Wang, Y., Liao, Z., Tang, T., & Ning, B. (2017). Train scheduling and circulation planning in urban rail transit lines. Control Engineering Practice, 61, 112–123. Wang, Z., Pang, C. K., & Ng, T. S. (2018). Data-driven scheduling optimization under uncertainty using renyi entropy and skewness criterion. Computers & Industrial Engineering, 126, 410–420.

12