Scheduling Jobs with Stochastic Processing Time on Parallel Identical Machines

1 Introduction

This paper addresses the scheduling of jobs with stochastic processing times on parallel identical machines. The aim of this problem is to assign jobs to machines such that the assignment maximizes the probability that all the jobs are completed before a given common due date. Due to the structure of the objective, the problem can be classified as a β-robust scheduling problem [Daniels and Carrillo, 1997], i.e., approaches that maximize the probability that the objective is below some threshold.

The scheduling problem is given by a set of n independent jobs J and a set of m identical machines M. In order to exclude trivial instances, we assume that n > m. Each job j ∈ J is characterized by a stochastic processing time πj . The processing time is defined by a normal distribution function N (µj , σ2 j ) with mean µj ∈ N and variance σ2 j ∈ N. The last parameter of the problem is a common due date δ ∈ N0. In the classical α|β|γ notation of [Graham et al., 1979], the problem can be denoted by P |πj ∼ N (µj , σ2 j )|Pr [Cmax ≤ δ]. Therefore, the solution of the problem is related to the partitioning of jobs between the machines such that machines are similarly loaded.

To solve the problem, one can utilize the property that the sum of independent normally distributed variables has a normal distribution. Then by introducing a binary decision variable zi,j, which is 1 if and only if job j ∈ J is assigned to machine i ∈ M, the problem can now be formulated [Ranjbar et al., 2012] as the following non-linear model

∗Corresponding Author

subject to

where Φ is normalized cumulative normal distribution function with µMi and σ2Mi being the mean and the variance of the total processing time on machine i respectively.

1.1 RelatedWork

The problem stated above is a fundamental problem for many real-world applications where deterministic processing time cannot be assumed. A good example is the duration of surgeries in the operating room scheduling problem [Sagnol et al., 2018]. Many problems with stochastic processing time are modeled as β-robust. For example, P|πj ∼ N (µj , σ2 j )|Pr [Cmax ≤ δ] introduced by Ranjbar et al. [Ranjbar et al., 2012], is based on β-robustness. In the same paper, the authors devised two Branch-and-Bound algorithms that differ in their branching schemes. The main limitation of this approach are symmetries induced by the branching scheme. The authors mitigate this issue by a suitable solution encoding and by a dominance rule; however, the results indicate that the CPU time sharply grows with a increasing number of machines. The same problem with a different objective is addressed in [Pishevar and Tavakkoi-Moghaddam, 2014]. The authors maximize the probability that the sum of job completion times will be lower than a given bound. The problem is slightly simpler since the objective function does not require the product operator. The paper proposes a mixedinteger nonlinear programming model which is approximated by a mixed-integer programming model. The same β-robust problem where the sum of job completion times in the objective is substituted by the sum of job completion flow times is studied by Alimoradi et al. [Alimoradi et al., 2016]. The authors propose a Branch-and-Bound algorithm that outperforms an earlier work [Wu et al., 2009] addressing the singlemachine problem.

A robust scheduling problem minimizing the number of identical machines required while completing all the given jobs before a deadline is investigated in [Song et al., 2018]. Unlike the above-mentioned works, the processing time is not defined via a distribution function but using intervals, similarly as in [Hamaz et al., 2018]. The "price" of robustness, i.e., the number of jobs that can deviate from their nominal processing times, is defined as a parameter. The problem is solved by a Branch-and-Price procedure. Another class of robust scheduling problems, where job processing times are stochastic without any assumed form of the distribution, is known as distributionally robust. Chang et al. [Chang et al., 2019] studies parallel machine scheduling problem which minimizes the worst-case expected total flow time out of all given probability distributions. This min-max problem is solved by a reduction to an integer second-order cone program. Unrelated machine scheduling problem with stochastic processing times was addressed by [Skutella et al., 2016].

1.2 Contribution and Paper Outline

In this paper, we study an exact method for solving the problem P|πj ∼ N(µj, σ2 j )|Pr [Cmax ≤ δ]. The devised algorithm is based on the Branch-and-Price approach. In comparison with [Ranjbar et al., 2012] our approach is simpler to implement and has better scalability since there are no symmetries in the exploration of the solution space. The experimental results show better results compared to Branch-andBound algorithms from [Ranjbar et al., 2012]. Furthermore, we prove that the Pricing Problem is at least weakly NP-hard and we show how it can be solved by mixed-integer linear programming. Last but not least, we substantially improved the initial heuristic described in [Ranjbar et al., 2012]. The proposed heuristic is used in the Branch-and-Price algorithm to initialize the Master Problem.

The paper is organized as follows. Description of the Branch-and-Price algorithm is divided into three sections: decomposition, Pricing Problem and branching scheme. Section 5 compares our results with the algorithm from [Ranjbar et al., 2012]. The last section concludes the work.

2 Branch-and-Price Decomposition

Branch-and-Price [Barnhart et al., 1998] is one of the prominent exact methods for solving combinatorial problems. It is a branch-and-bound algorithm where each node of the search tree consists of a Linear Program (LP), which is called Master Problem. The Master Problem typically contains an exponential number of variables that serves as indicators, to denote whether the particular configuration is a part of the solution. A configuration represents a fragment of the complete solution. All variables are not explicitly enumerated in the LP, but rather lazily generated during the solution as needed. Generation of new configurations is an iterative process, where at each step, the dual solution of the LP with the current set of configuration is considered. With this, we ask for a cut in dual LP which is violated by this dual solution. The question which configuration (i.e., a cut in dual) might improve the current primary objective value is solved by so-called Pricing Problem. This method is usually known as column generation [Desaulniers et al., 2006]. Most importantly, it can be shown [L¨ubbecke and Desrosiers, 2005], that it is not always necessary to generate all the variables (configurations) in order to prove the optimality of the full model. When no new configuration improving the solution can be found, a branching procedure is used to ensure an integer solution.

In the next section, we show how to reformulate the formulate problem (1)-(4) to fit the Branch-and-Price method.

2.1 Master Problem

The Master Problem utilizes the fact that maximizing the sum of logprobs maximizes the product of probabilities, hence, linearizing the objective. The configurations P in the Master Problem correspond to all possible machine configurations, i.e., a subset of jobs that are scheduled on a machine. Each configuration k ∈ P consists of a vector of column coefficients ak for the constraint matrix and a single objective coefficient log pk. The ak is a characteristic vector, i.e., aj,k = 1 if and only if job j ∈ J is scheduled on some machine with configuration k ∈ P . The objective coefficient log pk is given as

where µ = (µ1, . . . , µn)T and σ = (σ2 1 , . . . , σ2n )T . Hence, pk is the probability that the makespan on a machine with configuration k ∈ P is less or equal to δ. Then, the problem (1)-(4) can be reformulated as the Master Problem as follows

subject to

where yk is an indicator, whether configuration k ∈ P is selected for some machine. The constraint (6) states that each job is scheduled on some machine while (7) ensures that we use at most |M| machines.

In practice, the set of all configurations P is replaced with a smaller set P ′ ⊆ P that initially starts with |M| configurations, providing an initial feasible solution. The set P ′ is iteratively enlarged by new configurations generated by the Pricing Problem.

2.2 Initial Heuristics

The Branch-and-Price starts with the set of configurations P ′ that corresponds to a feasible solution found heuristically. Our initial heuristics, called Large Job Allocated First

Algorithm 1 LJAF heuristics

(LJAF), is shown in Algorithm 1. It first sorts jobs in a nonincreasing order of products of means and variances µj · σ2 j . The ties are broken by larger µj. In each step, a job is taken and assigned to the machine i ∈ M with the currently largest probability pi.

We note that according to our experiments in Section 5, Algorithm 1 produces better initial solutions than the lower bound heuristics proposed by [Ranjbar et al., 2012].

3 Pricing Problem

3.1 Problem Statement and Complexity

From the dual problem of (5)-(8) it can be derived, that the Pricing Problem takes the form of

The formulation is based on binary decision variable xj determining whether job j will be selected into the new configuration or not. Furthermore, ψ ≥ 0 are negative dual prices associated with constraints (6). Since ψj ≥ 0 and µj ≥ 0, it can be seen, that the Pricing Problem is balancing the gain ψj from taking job j ∈ J into the new configuration and the loss, proportional to µj . We note that without loss of generality we may assume that σT x > 0 as in any optimal solution at least one job is allocated on each machine.

In the rest of the subsection we show, that the Pricing Problem is at least weakly NP-hard by the reduction from KNAPSACK PROBLEM [Kellerer et al., 2004]:

Definition (KNAPSACK PROBLEM). The instance of the problem and the solution is given as follows:

INPUT: S = {1, 2, . . . , n}, v = (v1, v2, . . . , vn) ∈ Nn , w = (w1, w2, . . . , wn) ∈ Nn and C, k ∈ N.

OUTPUT: Is there S′ ⊆ S such that ∑i∈S′ vi ≥ k and ∑i∈S′ wi ≤ C?

Proposition 1. The Pricing Problem is at least weakly NPhard.

Proof. The general idea for the reduction is the following. We set ψ = v, µ = w and δ = C + ∆ where ∆ > 0 is some small non-negative number. Let Φσ(·) denote a cumulative normal distribution with zero mean and variance σ2 . Then it can be seen, as σ2 → 0, the logΦσ(δ - µT x) approaches to zero, when δ > µT x; see Figure 1. Hence, if we would set σ → 0, logΦ(·) serves as a penalty function reflecting whether µT x ≤ δ is violated.

However, since σ > 0 has to hold, we will set it close to zero such that logΦ(·) is almost zero within the controllable error ϵ for µT x ≤ C + ∆ for some small ∆ > 0 and a large negative number for C + ∆ < µT x such that every optimal solution of the Pricing Problem represents a feasible solution for KNAPSACK PROBLEM.

Assuming σi being a constant equal to o, we look for 1 > o > 0, such that

where δ = C + ∆ with ∆ = 10-1T w and ϵ = 10-1T v for any x such that wT x ≤ C. Hence, solving for o we get

where Φ-1 is quantile function of N(0, 1). Therefore, we set σ = o · 1 for the reduction. Note that the assumption of setting o < 1 is not crucial for the reduction; instead of it, we could multiply δ and µ by a factor ⌈o-1 ⌉ and let σ = 1 to achieve the same controllable error ϵ.

Next, we ensure, that σ is so small, such that every optimal solution of the Pricing Problem represents a feasible solution for KNAPSACK PROBLEM. Hence, we need to ensure, that the smallest violation of wT x ≤ C dominates the largest possible gain of vT x. That can be formulated as

from which we get

Therefore, we can see that the choice

suffices.

It is easy to see that the decision version of the Pricing Problem in NP. Hence, we are left to show that the Pricing Problem has a solution greater than some k′ > 0 if and only if KNAPSACK PROBLEM has solution greater or equal than k > 0, such that k′ = k - ϵ.

• KNAPSACK PROBLEM→ PRICING PROBLEM Assume knapsack solution S′ given as characteristic vector s′ . We have that vT s′ = ψT s′ ≥ k. Further-

more, since µT s′ = wT s′ ≤ C, we have that

where the first inequality follows from that S′ is a feasible solution, the second from the fact that o < 1. The last inequality follows from the choice of σ. Therefore, the the Pricing Problem has solution with objective at least k′ = k - ϵ.

• PRICING PROBLEM→KNAPSACK PROBLEM

Let us assume a solution x⋆ of the Pricing Problem with objective greater than some k′ > 0. Let us define a solution of KNAPSACK PROBLEM as S′ = {i | ∀x⋆i = 1} with the corresponding characteristic vector s′ = x⋆ . We need to show that (i) wT s′ ≤ C and (ii) vT s′ ≥ k.

• (i) By contradiction. Assume that wT s′ ≥ C + 1. Then by (11) we have that

which contradicts the assumption k′ > 0.

• (ii) By (10) we have that the log term is bounded by

Therefore

However, since ψT x⋆ = vT s′ ∈ N, k ∈ N and 0 < ϵ < 1, then it follows that vT s′ ≥ k.

3.2 Exact Algorithm

To find an exact solution of the Pricing Problem, we adapt a trick used for solving stochastic knapsack problems [Morton and Wood, 1998]. We observe that solving (9) is equivalent to the following

where V = {v, v + 1, . . . , v - 1, v} such that v and v is a lower and upper bound on σT x respectively. The advantage is that for any fixed v, the expressions δ/√v and µ/√v become constants. Hence, for every v ∈ V , we optimize

which can be done by formulating it as a Mixed-Integer Programming (MIP) where logΦ(·) is approximated by a piecewise linear function. The formulation can be solved by a special simplex method for piecewise linear functions implemented in Gurobi solver [Gurobi, 2019]. Furthermore, we impose the constraint on the minimum and the maximum number of jobs that can be scheduled on a machine derived by [Ranjbar et al., 2012] to speed up the solution.

The objective value (9) of an optimal solution x⋆ is compared with γ, which is the dual price for constraint (7). When the objective (9) is greater than γ, the configuration is added to the Master Problem and resolved. When it is less or equal to γ, the Master Problem is solved optimally.

However, we note that this approach for solving the Pricing Problem is efficient only when bounds on σT x are tight, i.e., |V | ≈ 20. For other cases, we fall back to a linearized model of (1)-(4) for a single machine, except two differences:

• the termψT z is included in the objective
• the constraint (2) relaxed.

In the next section, we describe a branching procedure and the rest of the Branch-and-Price algorithm.

4 Branching Scheme

When the optimal solution of the Master Problem is not integer, a branching procedure is used. Unlike the original problem formulation (1)-(4) the branching scheme in this Branchand-Price is chosen such that it does not consider the assignment to machines. The reason is not to introduce symmetries caused by the fact that the machines are identical. Therefore the branching directly uses decision variables of the Pricing Problem. The scheme selects a couple of jobs (a, b) for which one branch enforces that both jobs must always be scheduled on the same machine (i.e., xa = xb), while the other branch defines the opposite (i.e., xa + xb ≤ 1).

The couple of jobs for branching is chosen with respect to y⋆ , i.e., the last result of theMaster Problem. Indices of jobs j are sorted in ascending order with respect to max yk : aj,k = 1. Then the algorithm first tries the first job as a with the second one as b. If this couple is not selected, then it tries the first one with the third one, etc. When the second index gets to the end the second one is taken as a and the third one as b, etc. Each couple (a, b) is tested whether it was not selected before in the active branch and whether some of its branches does not cause a conflict with the previous branching decisions stored in set B. An element of the set is defined as tuple (a, b, e), where e is true if jobs a and b must be scheduled on the same machine and false otherwise. The branching scheme takes the first non-conflicting couple (a, b).

Finally, it remains to explain the way the algorithm tests whether a set of branching decisions is non-conflicting. If B is a set of non-conflicting branching decisions, it means that there must be at most |M| configurations that covers all jobs in J. If there are no such configurations the Master Problem cannot be feasible. Then the test can be formulated with the following constraints

where xj,i is a binary decision variable equal to one if job j is in configuration i and zero otherwise. Constraints (14) guarantee that all jobs are selected at least once and the rest of the constraints represent the branching decisions formB. In fact, the structure of the above formulation shows that it is a decision version GRAPH COLORING problem. The vertices of the graph are individual jobs. Two vertices a and b are merged into a single one if there is a branching decision (a, b, true) in B. Branching decisions (a, b, false) correspond to edges of the graph. The problem tests whether the graph can be colored with |M| colors. Even though the decision problem is NP-complete for |M| ≥ 3, the experiments show that its computational difficulty is negligible compared to the difficulty of the Pricing Problem.

The test whether jobs (a, b) are admissible for branching is carried out as follows. The algorithm creates two sets B+ = B ∪ {(a, b, true)} and B- = B ∪ {(a, b, false)}. Both are tested using formulation (14)-(16). Then, (a, b) is admissible for braching if both B+ and B- are feasible. Moreover, the configurations found by formulation (14)-(16) are added to the Master Problem in order to preserve its feasibility after removing configurations conflicting with an appropriate branching decision.

After branching, constraints defined in the new B need to be propagated to the Pricing Problem such that generated configurations are consistent with respect to the given branching decisions. Since we solve the Pricing Problem as a MIP, these constraints are naturally incorporated by putting equivalent constraints into the model (12)-(13) without affecting the structure of the problem.

5 Experimental Results

5.1 Instances and Testing Environment

All the instances were generated according to the methodology proposed by [Ranjbar et al., 2012]. We chose 2, 4 and 6 parallel machines with 14,16,18 and 20 jobs to schedule. Parameter c, which influences the range in which the variances were generated, had two possible values, specifically 0.25 and 0.75. For each combination of the parameters, 10 instances were created, resulting in a total of 240 instances.

Both the Branch-and-Price and Branch-and-Bound [Ranjbar et al., 2012] were implemented in C++. The Master and Pricing problems were solved using Gurobi 8.0 solver. The program was run on a server with two Intel Xeon E5 2.60 GHz processors, 252 GB RAM memory and a 64-bit operating system. The Branch-and-Price method was allowed to use up to 16 CPU cores. Since the Branch-and-Bound method proposed in [Ranjbar et al., 2012] is not designed to use more than single core, it was run only on a single CPU. Time limit for each instance was set to 3600 seconds.

5.2 Branch-and-Price and LJAF Heuristics

The experimental evaluation is focused on the comparison of our algorithm with the best of the Branch-and-Bound algorithms introduced by [Ranjbar et al., 2012]. The results for different benchmark sets are summarized in Table 1. Each benchmark set is characterized by the number of jobs n and machines |M|. Each set has 10 instances generated with c = 0.25 and 10 instances with c = 0.75, i.e., 20 all together.

For the Branch-and-Price algorithm, we provide the mean running time in seconds, mean number of expanded nodes in the branching scheme and percentage of instances that did not finish before the specified time limit. Note that if the value of expanded nodes for some instance is 1, it means that no branching was necessary and the instance was solved in the initial node. For both the runtimes and expanded nodes we also provide standard deviation. The Branch-and-Bound algorithm is also described by the mean runtime in seconds, then by the the number of expanded nodes in total, and again by percentage of instances that were not solved before the given time limit. Again, runtimes and the number of expanded nodes are given with their corresponding standard deviations.

Comparing the results for instances with two parallel machines (i.e., |M| = 2) we observe that although our algorithm (i.e., Branch-and-Price) almost did not branch, it performed significantly worse than the reference Branch-and- Bound. The reason is discussed later in this section. However, we still manage to solve most of the instances below the given time limit.

For 4 and 6 machines, our algorithm outperforms the Branch-and-Bound in most cases. We observe that the Branch-and-Bound is faster for instances with 14 jobs; however, the difference is already minimal. For instances with a higher number of jobs, our approach becomes dominant. Compared to the Branch-and-Bound, we were able to drastically reduce the running times even in cases with 18 and 20 jobs, where the other algorithm had 90% timeout rate or higher. Our algorithm suffered from several timeouts on these instances as well, but only for cases with 20 jobs and even in those we did not have more than 3 unsolved instances.

Regarding the instances with two parallel machines, it was revealed that the Pricing Problem generates new configurations for the Master Problem that assign approximately n/|M| jobs to each configuration and rarely adds a few more. This implies that instances with a lower number of available machines and a higher number of jobs (i.e., where the number of jobs to schedule to a single machine is also higher) have more possible combinations of configurations. Further inspection had indeed shown that for the instances with two machines, the Branch-and-Price generated in general more configurations per instance than for instances with more parallel machines with the same number of jobs. Moreover, the mean time spent for generation of one such configuration in the Pricing Problem was also higher. Therefore, the combination of both these factors is responsible for longer runtime. However, a possible improvement would be not to solve the Pricing Problem to optimality in every iteration of the column generation and thus save some computational time. Hence, a heuristic method for Pricing Problem can be employed to improve the performance.

Finally, the comparison of the initial heuristic solutions by [Ranjbar et al., 2012] and our LJAF approach is shown in Figure 2. The histogram depicts the distribution of initial lower bounds produced by these two heuristics, the higher, the better. It can be seen that our approach is more likely to provide better initial bound than the former one. This is most likely due to the similarity with the LPT (largest processing time first) rule used for deterministic problem P||Cmax, where it acts as a 4/3-approximation algorithm.

5.3 Non-linear Mixed-Integer Model

For further comparison, we have also developed a non-linear model, i.e., the equivalent of model (1)-(4). However, since cumulative normal distribution function Φ does not have an analytical form, for x ≥ 0 we use the approximation

To further improve the performance of the model, we have added a symmetry breaking constraint in form of

to break the symmetry between identical machines. The results show that for instances with |M| = 4 machines and n = 10 jobs, the solution with SCIP 6.0.0 [Gleixner et al., 2018] non-linear mixed-integer solver takes about 3000 s. Therefore, we have excluded this method from further experiments. During the experiments, we observed that results were similar for different approximations of Φ.

6 Conclusions

This research work reflects a growing interest in stochastic scheduling problems. In this paper, we concentrated on a stochastic variant of P||Cmax where the processing time of each job is given by a normal distribution function. This variant is formulated as a β-robust scheduling problem where the objective is to maximize the probability the schedule is completed before a given common due date.

The main contribution of our work is a new algorithm significantly outperforming a Branch-and-Bound algorithm proposed in [Ranjbar et al., 2012]. The main reason for the better performance is a better elimination of symmetrical solution in the solution search space. Nevertheless, the Branch-andBound algorithm of Ranjbar et al. performs better on the data set with two machines. This case would probably need a different approach to the pricing problem since the experimental results indicate that this is the bottleneck of the algorithm. Other data sets show that our algorithm scales significantly better and can solve larger instances in a shorter time.

For future work, we suggest generalizing our approach to other families of probability distributions. For distributions that have known distributions of their sum, such as the ones that are closed under convolution (e.g., Poisson, Cauchy, Gamma) the current form of the Pricing Problem can be modified to handle them. For others, such as lognormal distributions, approximations of their sum exist [Mehta et al., 2007]. To handle general cases, sampling methods from stochastic knapsack problems [Morton and Wood, 1998] might be applicable to modify the structure and solution of the Pricing Problem.

Acknowledgements

This work was supported by the EU and the Ministry of Industry and Trade of the Czech Republic under the Project OP PIK CZ.01.1.02/0.0/0.0/15 019/0004688 and by the Grant Agency of the Czech Technical University in Prague, grant No. SGS19/175/OHK3/3T/13.