A NEW AND GENERAL STOCHASTIC PARALLEL MACHINE SCHELOC PROBLEM WITH LIMITED LOCATION CAPACITY AND CUSTOMER CREDIT RISK

. Scheduling-Location (ScheLoc) problem considering machine location and job scheduling simultaneously is a relatively new and hot topic. The existing works assume that only one machine can be placed at a location, which may not be suitable for some practical applications. Besides, the customer credit risk which largely impacts the manufacturer’s profit has not been addressed in the Sch-eLoc problem. Therefore, in this work, we study a new and general stochastic parallel machine ScheLoc problem with limited location capacity and customer credit risk. The problem consists of determining the machine-to-location assignment, job acceptance, job-to-machine assignment, and scheduling of accepted jobs on each machine. The objective is to maximize the worst-case probability of manufacturer’s profit being greater than or equal to a given profit (referred to as the profit likelihood). For the problem, a distributionally robust chance-constrained (DRCC) programming model is proposed. Then, we develop two model-based approaches: (i) a sample average approximation (SAA) method; (ii) a model-based constructive heuristic. Numerical results of 300 instances adapted from the literature show the average profit likelihood proposed by the constructive heuristic is 9.43% higher than that provided by the SAA, while the average computation time of the constructive heuristic is only 4.24% of that needed by the SAA.


Introduction
Recent information and automation technologies provide foundations for Industry 4.0 [2,15].The core idea of Industry 4.0 is to use the emerging technologies in a way that business and engineering processes are deeply integrated [24].Scheduling problem and location problem play important roles in manufacturing, and both of them have received much attentions from academia [1,14,26,27].In practice, these two problems closely intertwine and impact the production performance together.For example, in container port operations, decision makers must determine the positions of ships in berth (i.e., location problem) and the schedule for loading and unloading containers (i.e., scheduling problem), simultaneously [11].Another example comes from the mining industry.The planners must jointly determine which candidate locations should be selected to deploy crushers (i.e., space problem) and how to sequence minerals on the crushers at those locations (i.e., scheduling problem) [10].These facts indicate the necessity of investigating the integrated scheduling-location problem.
The scheduling-location problem is referred to as ScheLoc problem which is first introduced by Hennes and Hamacher [9].The deterministic single machine ScheLoc problem is well addressed by Hennes and Hamacher [9] and Kalsch and Drezner [11].Since parallel machine manufacturing environment is quite common in practice, recent works focus on the deterministic parallel machine ScheLoc problem [8,13,28].In reality, unexpected events often occur.Especially, job processing time can be uncertain as it may be affected by various factors such as machine and equipment conditions, manufacturing environment and operator skills [3].Therefore, Liu et al. [19] study a stochastic parallel machine ScheLoc problem with uncertain job processing time.Moreover, Liu and Liu [16] address the parallel machine ScheLoc problem with only mean vector and covariance matrix of job processing time known.These works promote the studies of ScheLoc problems.However, all of them assume that only one machine can be placed at a location, which may not be suitable for some practical applications.In reality, some locations are larger to accommodate more machines than the others.Therefore, this study extends the existing works and investigate a new and general parallel machine ScheLoc problem.For this novel problem, we assume a location can place several machines and its holding capacity is limited.
Usually, manufacturers complete required jobs and get the corresponding payment timely.However, according to Liu et al. [20], many manufacturers have experienced late payment from their customers.Manufacturers prefer jobs from customers with low credit risk to curb the adverse consequences of high customer late payment probability and maintain competitiveness.Tsosie and Nicastro [22] point out that customer credit has been applied to evaluate the customer payment probability.Liu et al. [20] indicate further that around 34% manufacturers have started to assess the customer payment probability based on their credit history.They are also the first to consider customer credit risk for flow shop scheduling problem.However, there is no research investigating the ScheLoc problem under customer credit risk.
Motivated by the above facts, this work addresses a novel and general stochastic parallel machine ScheLoc problem with limited location capacity and customer credit risk.The customer credit risk and job processing time are both assumed to be ambiguous, and only their partial distributional information is given.The studied problem consists of determining the machine-to-location assignment, job acceptance, job-to-machine assignment, and scheduling of accepted jobs on each machine.The objective is to maximize the profit likelihood, i.e., the worst-case probability of manufacturer's profit being greater than or equal to a given level.Besides, the worst-case probability of job completion time respecting the due date (referred to as the service likelihood) is also guaranteed.For the problem, a distributionally robust chance-constrained (DRCC) programming model is constructed.Then, we develop two model-based approaches: (i) a sample average approximation (SAA) method; (ii) a constructive heuristic.The contributions of this paper mainly include: (1) A new and general stochastic parallel machine ScheLoc problem with limited location capacity and customer credit risk is investigated.Especially, the assumption in the literature that only one machine can be placed at a location is relaxed.(2) For the problem, a DRCC programming model is proposed.
(3) Two model-based approaches, i.e., a SAA method and a constructive heuristic, are developed to efficiently solve the problem.
The remainder of this paper is laid out as follows.A literature review is presented in Section 2. In Section

Literature review
This paper addresses a general stochastic parallel machine ScheLoc problem with limited location capacity and customer credit risk.For the problem, a DRCC programming model is constructed.In the following, we review the related works about ScheLoc problems and distributionally robust (DR) approaches.
The ScheLoc problem, first introduced by Hennes and Hamacher [9], has received much attention from academia [8,12,13,16,19,25,28].Most existing works concern the deterministic ScheLoc problems [8,9,12,13,25,28].Hennes and Hamacher [9] first study a deterministic single machine ScheLoc problem, to minimize the makespan.A polynomial-time algorithm based on earliest release date (ERD) first rule is developed to solve the problem.Since parallel machine environment is quite common in practice, Heßler and Deghdak [8] investigate a deterministic parallel machine ScheLoc problem.The objective is to minimize the makespan.For the problem, an mixed integer programming (MIP) model is first proposed.Then, several constructive heuristics are developed for it.1450 instances are tested to demonstrate the effectiveness of the developed algorithms.The same problem is then studied by Wang et al. [25], Kramer and Kramer [12] and Li et al. [13].Wang et al. [25] first propose a new MIP model.Then, three constructive heuristics are designed to solve the problem.Numerical experiment results on benchmark instances proposed by Heßler and Deghdak [8] show that both the new MIP model and the proposed heuristics perform better than that of Heßler and Deghdak [8].Kramer and Kramer [12] propose a new MIP model from the perspective of arc-flow formulation.Then, a column generation method, two constructive heuristics and an iterated local search metaheuristic are designed.Numerical experiments show that the optimal solutions of the same benchmark instances can be obtained.Li et al. [13] propose three new mixed integer linear programming (MILP) models and a novel logic-based Benders decomposition (LBBD) method.Numerical results demonstrate the efficiency of the new formulations and the algorithm.Zhang et al. [28] investigate a new variant of the parallel machine ScheLoc problem considering delivery times and due dates.The work aims to minimizing the total cost.They first develop MILP models.A tailored LBBD method and a matheuristic are then developed to solve the problem.
In practice, unexpected events such as machine breakdown often occur.Therefore, some recent works focus on the ScheLoc problem with uncertain job processing time.Liu et al. [19] first study a stochastic parallel machine ScheLoc problem.The objective is to minimize the weighted sum of the expected total completion time and the location cost.For the problem, a two-stage stochastic programming model is established.Then, a SAA method, a genetic algorithm and a scenario-based heuristic are developed to solve the model.Liu and Liu [16] consider the stochastic parallel machine ScheLoc problem, in which only the mean and covariance matrix of job processing time are known.The objective is to minimize the machine location cost.The worst-case probability of job completion time respecting the due date is also guaranteed.For the problem, a distributionally robust (DR) formulation is first proposed.Then, an approximated mixed integer second order cone programming model is developed for it.
The above works provide beneficial support for decision-making in the ScheLoc problem.However, all of them restrict that only one machine can be placed at a location, and ignore the customer credit risk which has a non-negligible impact on the manufacturers' profit [20].Note that Liu et al. [20] are the first one to study a stochastic flow shop scheduling problem with customer credit risk.Nevertheless, there is no research work integrating it into the ScheLoc problem, to the best of our knowledge.
Besides, due to data scarcity, it is usually impractical to accurately estimate the probability distribution of random parameters [5].DR optimization, first introduced by Scarf [21], is an appropriate tool to handle uncertainty in the case where its distributional knowledge cannot be fully obtained and optimize the worst-case performance measure over an uncertainty set, namely ambiguity set.Wagner [23] and Delage and Ye [5] use moment-based ambiguity sets for deriving tractable reformulations of DR models.Referring to their works, we also consider a moment-based ambiguity set which is constructed based on the given mean vector and covariance matrix of uncertain parameters for this study.
Based on above facts, in the work, we study a new and general stochastic parallel machine ScheLoc problem with limited location capacity and customer credit risk.Especially, we assume that customer credit risk and job processing time are ambiguous and only their partial distributional information is given due to data scarcity.
For the problem, a DRCC programming model is constructed.Table 1 summarizes main differences between our work and the closely related literature.

Problem description and formulation
In the section, the detail of the studied problem is first described.Then, a new mathematical programming model is constructed.Besides, a complexity analysis is conducted for the studied problem.

Problem description
For the studied problem, consider a set of candidate locations , a set of machines  and a set of jobs .A location  ∈  can be selected to place at most   machines, where   is called the capacity of location  ∈ .When location  ∈  is selected, there is a location selection cost    .A cost    is incurred when a machine  ∈  is assigned to a location  ∈ .A job  ∈  with a customer payment probability   ∈ [0, 1] corresponding to a credit risk can be rejected, but results in penalty   .If job  ∈  is accepted to process, the manufacturer can obtain an expected revenue   •   , where   denotes the revenue without credit risk [20].In addition, job  ∈  is located initially at a storage area and has a release time   that corresponds to its moving time from the storage area to a location  ∈ .Each job  ∈  has a due date   .The uncertain processing time of job  ∈  is denoted as   , for which only the mean and covariance are known.Besides, this work follows basic assumptions: (1) A location can place several machines, and its capacity is limited.
(3) Job  ∈  can only be assigned to exactly one machine if it is accepted to process.(4) A machine can process at most one job at a time.
(5) Any two uncertain parameters are independent of each other.
To illustrate the studied problem, an example with 3 candidate locations, 2 machines and 8 jobs is presented.The job storage areas and candidate locations are dispersed on a network in Figure 1.The release time of job  to location  is set to be their Euclidean distance (e.g.,  21 = 2) where we assume that the common moving speed is 1 per time unit.The means of job processing times and customer payment probabilities are set to be [5,4,6,3,3,2,9,7] and [0.6, 0.8, 0.3, 0.7, 0.5, 0.8, 0.1, 0.3], respectively.The job rejection penalties, the revenues for processing job without credit risk, and the due dates are set to be [4, 2, 4, 2, 1, 5, 1, 2], [8,10,7,11,12,14,6,10], and [12,7,19,6,18,8,23,19], respectively.The location selection costs are set to be [18,  30, 10], and the machine-to-location assignment costs and the capacity of each location are uniformly set to be 1.The problem consists of determining the machine-to-location assignment, job acceptance, job-to-machine assignment and scheduling of accepted jobs on each machine.A feasible solution is shown in Figure 2. It can be observed that machines 1 and 2 are assigned to locations 1 and 3, respectively.Job 7 is rejected.Jobs 2, 1 and 3 are assigned to the machine 1 at location 1, while jobs 4, 6, 8 and 5 are assigned to the machine 2 at location 3. The expected profit is calculated as 13.8.

A new DRCC programming model
In the section, problem parameters and decision variables are first defined.Then, a DRCC programming model is constructed.The completion time of job  ∈ .

𝐾
In the following, we present the DRCC programming model, named as P1.
[P1]: max inf The objective function ( 1) is to maximize the profit likelihood, i.e., the worst-case probability of manufacturer's profit being greater than or equal to the given level .Notably, the manufacturer's profit consists of four parts, i.e., the revenue by processing accepted jobs, the job rejection penalty, the cost of selecting locations, and the cost of assigning machines to these selected locations.
Constraint (2) denotes the number of machines placed at location  cannot exceed its capacity.Constraint (3) indicates a machine cannot be placed at multiple locations.Constraint (4) implies that if a job has been accepted for processing, it must be assigned to a machine.Constraint (5) denotes a machine can process jobs only when it has been placed at a location.Constraint (6) indicates there is no processing precedence relationship between job  and job  if they are not assigned to the same machine  at location .Constraint (7) expresses that either job  is processed before job , or job  is processed before job  provided that they are both assigned to the same machine  at location .Constraint (8) ensures that the service likelihood is greater than or equal to a certain level 1 − .Constraint (9) denotes the domain of decision variables.

The complexity analysis
In this subsection, we conduct a complexity analysis for the parallel machine ScheLoc problem with limited location capacity and customer credit risk.
Theorem 1.The concerned parallel machine ScheLoc problem with limited location capacity and customer credit risk is NP-complete.
Proof.The proof is based on the following fact that PARTITION reduces to the parallel machine ScheLoc problem with limited location capacity and customer credit risk.Given a set of integers  1 ,  2 , . . .,  || and an auxiliary integer , such that  = 1 2 ∑︀ || =1   .The PARTITION problem is to decide whether the set can be partitioned into two subsets such that the sum of the integers of the two subsets are both equal to .It is known that the PARTITION problem is NP-complete, which is a decision problem in terms of complexity theory.
In the following, we construct a special case of the studied problem (i.e., a decision problem).The number of candidate locations, machines and jobs are set to be 2, 2, ||, respectively.The capacities of both locations are set to be 1.The location selection cost, machine-to-location assignment cost, and job rejection penalty are uniformly set to be 0. The customer payment probability of each job and the revenue for processing each job are uniformly set to be 1.The release time of each job from its storage area to any location is set to be 0. The processing time of job  ∈  is set to be   , and all jobs' due dates are uniformly set to be .To sum up,   = 1,    = 0,    = 0,   = 0,   = 1,   = 1,   = 0,   =   ,   = ,  ∈ {1, 2},  ∈ .In the case, a schedule with profit of || exists if and only if all jobs have been accepted and processed no later than the common due date .This can be done if and only if the jobs can be partitioned into two sets with equal total processing time, which is the PARTITION problem.

Solution methods
Due to the objective function (1) in probabilistic form and DR chance Constraint (8), model P1 cannot be directly solved via commercial solvers.Thus, we first propose a widely applied SAA method to solve the problem approximately.However, from the computation results in Section 5, we can see that no feasible solution can be found within the time limit for some large-scale instances.Therefore, we further propose a model-based constructive heuristic to solve the problem efficiently.

The SAA method
The principle of the SAA method is to replace the unknown probability distribution by an empirical distribution, which corresponds to a set of scenarios.In recent years, it has been widely applied to transform the DRCC programming model into a deterministic one [7,18,20].Referring to these works, in the following, we apply this method to convert the model P1 into a tractable model P2.  2) − ( 7), ( 9)

Ω
The objective function (10) is to maximize the proportion of scenarios in which the manufacturer's profit is greater than or equal to the given level .Constraints (11) express whether the manufacturer's profit being greater than or equal to a given level  in scenario  ∈ Ω. Constraints ( 12) and ( 13) calculate the completion time of job  ∈  in scenario  ∈ Ω. Constraints (14) indicate whether the completion time of each job is no later than its due date.Constraints (15) represent the number of scenarios in which at least one job's completion time exceeds its due date is no greater than the give level (i.e., |Ω| • ).Constraints ( 16) denote the domain of decision variables.The SAA method transforms the original DRCC programming model P1 into a scenario-based mixed-integer linear programming (MILP) model P2, which can be directly solved by commercial solvers such as Gurobi.However, the computation time of the method will grow exponentially with the size of the problem.Therefore, to solve the large-scale problem efficiently, we propose a model-based constructive heuristic.

The model-based constructive heuristic
In this subsection, we propose a new model-based constructive heuristic to solve the large-scale problem efficiently.Before presenting the heuristic, we first employ some common methods to transform the DRCC programming model P1 into a deterministic mathematical programming model.

An equivalent transformation
In the first step, we introduce an auxiliary continuous decision variable  ∈ [0, 1] to equivalently transform the original objective function (1) into formulas ( 17)-( 19) as follows. min

Bonferroni inequality
In the second step, we apply Bonferroni inequality to break the joint chance constraint (8) into a set of individual chance constraints.Specifically, joint chance constraint ( 8) can be conservatively approximated as individual chance constraints ( 20) and ( 21), where  ′  and  ′′  must satisfy . For convenience, we divide the risk tolerance  evenly in line with Do Chung et al. [6], i.e.,

A deterministic mathematical programming model
In the third step, we refers to Theorem 2.2 in Wagner [23] to equivalently transform the DR individual chance constraints into deterministic constraints (i.e., without probability).Before applying this theorem, we first introduce some new parameters.

𝑗
The variance of   , where  ∈ ; The expectation of   , where  ∈ ;  2

𝑗
The variance of   , where  ∈ .
After that, a new deterministic mathematical programming model P3 is presented as follows.
Model P3 still cannot be directly solved via commercial solvers due to the nonconvex constraint (22).However, note that this nonconvex constraint become a convex (i.e., second-order conic) one if the decision variable  is fixed.Motivated by this fact, we design a constructive heuristic based on model P3.

The constructive heuristic
The basic idea of the heuristic is trial-and-error.That is to say, given a certain , we can check whether model P3 is feasible by calling commercial solvers such as Gurobi.If the condition is true, the value of  is reduced via the golden-section search method, and vice versa.Besides, to offset the impact of conservative approximation caused by Bonferroni inequality, we enlarge the solution search space via modifying the preset profit level and job due date, i.e.,  ′ = / 1 and  ′  =  2 •   , with tuning parameters  1 ,  2 > 1.The time complexity of this constructive heuristic is O(M N T ), where M denotes the maximum number of updates for  and   , N represents the number of iterations for , and we use O(T ) to denote the time complexity of checking the feasibility of P3 with fixed  (Note that the feasibility version of model P3 with fixed  is a decision problem in terms of complexity theory.Such a decision problem is NP-complete, because Constraints (24) are Knapsack-like constraints.Therefore, for simplicity, we use O(T ) to denote the time complexity of solving this decision problem).Algorithm 1 presents the details of the constructive heuristic.)︂ ; Output: The minimum  in  * and its corresponding solution.

Numerical experiments
In this section, the illustrated example in Section 3 is first solved to demonstrate the validity of the SAA method and the constructive heuristic.Then, the performance of the two developed methods, and one adapted version of the best heuristic in stochastic ScheLoc area from the literature (i.e., the scenario-based heuristic proposed by Liu et al. [19]), are further compared by 300 instances adapted from Heßler and Deghdak [8].All numerical experiments are performed on a personal computer with Core I5, 2.11 GHz processor, 16 GB RAM, and Windows 10 Operating System, and the two solution methods are both coded in Python 3.6 combined with Gurobi 9.0.The computation time of each method is limited to 3600 s.

Evaluation standard
Solutions obtained by different methods are compared by the out-of-sample test in line with Liu et al. [20].Specifically, 10 000 scenarios are first randomly generated, which represents various realizations of uncertain parameters.Then, the profit likelihood (denoted as Υ 1 ) and service likelihood (denoted as Υ 2 ) are measured as follows: (1) The profit likelihood Υ 1 is measured by ( 1 /10, 000), where  1 denotes the number of scenarios in which manufacturer's profit no less than the given level; (2) The service likelihood Υ 2 is measured by ( 2 /10, 000), where  2 denotes the number of scenarios in which no job's completion time exceeds its due date.
Accordingly, the performance of different methods is compared by the values of Υ 1 and Υ 2 .

The solutions obtained by the two methods for the example
In this part, the example in Section 3 is solved to demonstrate the effectiveness of the two methods in Section 4. The number of scenarios (i.e., |Ω|) for the SAA method is set to be 30, and the values of parameters M and N for the constructive heuristic are set to be 5 and 10, respectively.The given profit level is set to be 6, and the risk tolerance  is set to be 0.2.
The schedules obtained by the two methods are presented in Figure 3. From Figure 3, we can see that the selected locations of two methods are identical, i.e., both locations 1 and 3 are selected.Recall that the capacity of each location is set to be 1.Therefore, two machines are assigned to two locations, respectively.For the SAA method, it can be observed that jobs 2 and 4 are rejected.Jobs 1, 5 and 3 are assigned to machine 1 at location 1, while jobs 6, 8 and 7 are distributed to machine 2 at location 3.In contrast, for the constructive heuristic, we can see that jobs 3 and 4 are rejected.Jobs 2, 1 and 5 are distributed to machine 1 at location 1, while jobs 6, 8 and 7 are assigned to machine 2 at location 3.
The expected profits of solutions obtained by the two methods are calculated as 6.7 and 10.6, respectively.Moreover, out-of-sample test shows the profit likelihood proposed by the constructive heuristic is (0.9259 − 0.5905)/0.5905× 100% = 56.8%higher than that provided by the SAA, while the service likelihood of the constructive heuristic is only (0.8537 − 0.8014)/0.8537× 100% = 6.1% lower than that of the SAA.Note that the service likelihood of the two methods are both higher than the given level (i.e., 1 −  = 0.8).

Numerical experiments on 300 instances
In this part, the performance of the two developed methods and the scenario-based heuristic from the literature are compared by 300 adapted instances.In the following, the instance generation rule is first illustrated.Then, computation results of different methods for these instances are analyzed.

Instance generation rule
The instances are adapted from the large-scale data set for the ScheLoc problem provided by Heßler and Deghdak [8].This data set includes the following information (i) the number of candidate locations; (ii) the number of machines; (iii) the number of jobs; (iv) job processing time; (v) job release time.]︁ , in which both  1 and  2 are adjustable parameters [17].In this paper,  1 is set to be 0.2, and  2 is set to be 0. ]︁ , in which the job processing revenue without credit risk and the job rejection penalty are randomly generated following  [10,15] and  [1,5], respectively [20], and the location selection cost is randomly generated following  [10,50].The capacity of each location is uniformly set to be 1, and the machine-to-location assignment cost is uniformly set to be 0. The standard deviation of each uncertain parameter is set to be 0.3 times its expectation [20].
Note that we adopt the first sixty combinations in Heßler and Deghdak [8]'s data set, and for each problem setting, we randomly generate five instances.Thus a total of 300 instances are tested.

Computation results
The computation results obtained by different methods for the 300 instances are shown in Table 2. Specifically, the first column denotes the sequence number of the tested problem setting, and the second column reports its size which includes the number of candidate locations, machines and jobs.The third to fifth columns give the profit likelihood, the service likelihood and the computation time of the SAA method, respectively.The sixth to eighth and the last three columns show the corresponding results of the scenario-based heuristic and the constructive heuristic, respectively.The average value of each indicator is reported in the last row of the table, which is calculated based on the instances for which all the methods can find feasible solutions within 3600 s.
From Table 2, we can see that the scenario-based heuristic and our proposed constructive heuristic can obtain solutions for all the combinations, while the SAA method can only obtain solutions for 53 of them.Moreover, the average computation time of the constructive heuristic is only 51.4/1212.7 = 4.24% and 51.4/231.9= 22.2% of those needed by the SAA and the scenario-based heuristic, respectively.These results demonstrate the efficiency of the constructive heuristic for solving the studied problem.

Figure 1 .
Figure 1.An illustrated example for the problem.

Figure 2 .
Figure 2. A feasible solution of the example.

Figure 3 .
Figure 3.The schedules obtained by the two methods for the example in Section 3.1.

Table 1 .
Comparison of ScheLoc related works.
Release time of job  ∈  to location  ∈  that corresponds to its moving time from 's storage area to location ;   Due date of job  ∈ ;   Uncertain processing time of job  ∈ ; ∈   + max ∈ {  }, which is taken as a large enough positive number; P [( − )( − )  ] = Γ  Equals 1 if location  ∈  is selected, and 0 otherwise; ℎ  Equals 1 if job  ∈  is accepted to process, and 0 otherwise;   Equals 1 if machine  ∈  is placed at location  ∈ , and 0 otherwise;   Equals 1 if job  ∈  is assigned to machine  ∈  at location  ∈ , and 0 otherwise;   Equals 1 if job  ∈  is processed after job  ∈  on machine  ∈  at location  ∈ , and 0 otherwise; Set of scenarios, which is indexed by ;   Customer payment probability of job  ∈  in scenario  ∈ Ω;   Processing time of job  ∈  in scenario  ∈ Ω;   The completion time of job  ∈  in scenario  ∈ Ω;   Equals 1 if at least one job's completion time exceeds its due date in scenario  ∈ Ω, and 0 otherwise.
∈  +max ∈ {  }, which is taken as a large enough positive number in scenario  ∈ Ω.