Robust Scheduling of Two-Agent Customer Orders with Scenario-Dependent Component Processing Times and Release Dates

Abstract

Although some uncertainty factors can occur in many practical environments, customer order scheduling problems involving two agents in such uncertain environments have not received attention in the current literature. Motivated by this observation, we address a two-agent customer order scheduling problem where various customer orders have scenario-dependent component processing times and release dates in order to find an appropriate schedule to minimize the maximum of the total completion time of the customer orders that belong to one agent and are subject to a constraint with the other agent. In order to solve this problem, a lower bound and six dominant properties are derived and used to propose a branch-and-bound algorithm to find an exact optimal solution. Afterward, three local search heuristics and two variants of a simulated annealing hyper-heuristic are proposed and empirically evaluated in order to find approximate solutions. Finally, we conclude the paper with a summary of our findings and some directions for future research.

1. Introduction

Recent calls to enhance the practical relevance of scheduling research has resulted in the consideration of the following three problem types: Customer Order Scheduling Problems (COSPs), first introduced in the studies by Ahmadi and Bagchi [1,2] and Gupta et al. [3]; multi-agent scheduling problems (MASPs), first defined by Gupta et al. [3], Baker and Smith [4], Wu et al. [5], and Agnetis et al. [6]; and scheduling problems with scenario-dependent processing times (SDPTSPs), first formalized by Daniels and Kouvelis [7,8]. More recently, Framinan et al. [9] provided a comprehensive and unified picture of COSPs. Agnetis et al. [10] provided a comprehensive survey of the multi-agent scheduling problems (MASPs); meanwhile, Perez-Gonzalez and Framinan [11] reviewed a survey paper on multi-criteria order scheduling problems. More recently, Braga-Santos et al. [12] proposed an efficient size-reduction algorithm to solve a COSP in order to minimize the total tardiness. Taking the total completion time as the measurement criterion, de Athayde Prata et al. [13] applied an innovative discrete differential evolution algorithm with differential mutations to solve a COSP with sequence-dependent setups. Pinto Antonioli et al. [14] addressed a mixed-integer linear programming as well as two adaptation heuristics and two metaheuristics to minimize the total tardiness COSP with a sequence-dependent setup time.

However, scheduling problems with scenario-dependent processing times (SDPTSPs) have not received much attention in the literature. SDPTSPs deal with a situation where the uncertain job processing times depend on the specific scenario that may arise in practice, thus giving rise to discrete job processing times called scenario-dependent processing times. In such situations, because of the uncertainty associated with the realization of any specific scenario, optimizing the given objective function among all possible scenarios is desired (Kouvelis and Yu [15]). Applications of SDPTSPs in several manufacturing environments, flexible manufacturing situations, and other practical circumstances are discussed by Daniels and Kouvelis [8]). For more examples, we refer readers to a few existing research efforts by Yang and Yu [16], Monch et al. [17], Aloulou and Della Croce [18], Aissi et al. [19], Kasperski and Zielinski [20], Wang et al. [21], Liu et al. [22], Wang et al. [23,24], Wu et al. [25], Xing et al. [26], Ren et al. [27], and Wu et al. [28], among others.

A review of the existing scheduling literature reveals that the research on COSPs involving multiple facilities and MASPs considered the job processing times or release dates as fixed integers. This is at odds with real-life environments. In fact, there are many significant uncertainty factors in practical production environments. For example, employees might become unable to work, machines might break, the working environments might change, and some external complex factors might provoke job cancellations or a changed tool quality. In such situations, the worst-case criterion for a system might be more important than the average-case criterion (Kouvelis and Yu [15]). Thus, inspired by these observations and in order to extend the application of scheduling research to the solving of practical problems, we address a two-agent COSP with this paper, for the first time in the scheduling literature. More specifically, we address scenario-dependent component processing times and release dates that minimize the maximum of the total completion time of the customer orders that belong to one agent and which are subject to the following constraint—that the maximum total completion time of the customer orders of the other agent should be no greater than an upper bound for all scenarios and all possible schedules.

Table 1. Summary of most recent contributions published concerning the COSP.

Problem Setting	Algorithm	Objective Form	References
COSP without setup times	A size-reduction algorithm	${{\displaystyle \sum }}_{i=1}^n{T}_i\left(\sigma \right)$	Braga-Santos et al. [12]
COSP with sequence-dependent setups	A differential evolution	${{\displaystyle \sum }}_{i=1}^n{C}_i\left(\sigma \right)$	de Athayde Prata et al. [13]
COSP with sequence-dependent setups	Hybrid matheuristics	${{\displaystyle \sum }}_{i=1}^n{T}_i\left(\sigma \right)$	Pinto Antonioli et al. [14]
COSP without setup times	Iterated greedy; BB	$ma{x}_s\left\{{{\displaystyle \sum }}_{i=1}^n{T}_i^s\left(\sigma \right)\right\}$	Wu et al. [28]
COSP without setup timeswith two agents	Simulated Annealing Hyper-heuristic; BB	$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}\mathrm{s}.\mathrm{t}.$$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}\le U$	This paper (2022)

summarizes the most recent contributions published concerning the COSP.

In view of the fact that the proposed problem is an NP-hard one, we propose a branch-and-bound method to find an exact solution. For more details on branch-and-bound, the reader could refer to Zegordi et al. [29] and Wu et al. [25]. On the other hand, we also propose some simulated annealing algorithms (SA) and two variants of simulated annealing hyper heuristics to find near-optimal solutions. For details on SA or SAHH, readers can refer to Kirkpatrick et al. [30], Bıçakcı et al. [31], Azimi [32], Anagnostopoulos and Koulinas [33], Bai et al. [34], and Liang et al. [35].

The several contributions of this article include: (1) The presentation of a new COSP with two agents as well as scenario-dependent component processing times and release dates. (2) The proposal of a branch-and-bound algorithm along with six properties and a lower bound to find the optimal solution. (3) The proposal of three heuristics and two variants of simulated annealing hyper heuristics in order to find near-optimal solutions. (4) The execution of simulation tests to evaluate the performances of all proposed methods.

The rest of this study is organized as follows. Section 2 describes the notations, defines the problem, proposes lower bound and six dominance properties, and outlines a branch-and-bound algorithm to optimally solve the problem. Section 3 discusses three local search heuristics and two variants of simulated annealing hyper heuristics to obtain optimal or near-optimal solutions. In Section 4, we empirically investigate the performance of the proposed heuristics in generating optimal or near-optimal solutions for the considered problem. Finally, we conclude the paper in Section 5.

2. Problem Definition and Properties

Throughout this paper, we use the following notations.

Notations

$N=\left\{1,2,\dots ,n\right\}$: represents a set of n jobs (orders);

$\Phi =\left\{x, y\right\}$: denotes the set of two agents, x and y;

$O_i^x(O_i^y)$: denotes a customer order i from agent x (y), i$\in N$;

$\Omega$_N = $\left\{O_1^{\Phi }, O_2^{\Phi }, \cdots ,O_n^{\Phi }\right\}$: presents a set of n customer orders, where $O_i^{\Phi }$ (or in brief, order i) presents a customer order i from agent x or agent y;

$O^x={\bigcup }_{i\in N}{O}_i^x, O^y={\bigcup }_{i\in N}{O}_i^y$;

$M=\left\{1, 2, \dots , m\right\}$: represents a set of indexes of m parallel machines,M₁, M₂, .., M_m;

s: scenario of customer order parameters, s = 1, 2;

$\sigma ,{\sigma }^{\prime }$: denote two full schedules of n customer orders;

$\delta ,{\delta }^{\prime }$: denote two partial schedules of n customer orders;

$t_{i\upsilon }^{\left(s\right)x}$(or $t_{i\upsilon }^{\left(s\right)y}$): is the component processing time of a customer order $O_i^{\Phi }$ of agent x (or agent y) on a machine $M_v,\upsilon \in M$, s = 1, 2;

$r_i^{\left(s\right)x} \left(or r_i^{\left(s\right)y}\right)$: denotes the ready (or release) time of a customer order i of agent x (or agent y), i$\in N$;

[ ]: presents the position of a customer order in a schedule;

$C_{\left[k\right]}^{\left(s\right)x}\left(\sigma \right)$ (or $C_{\left[k\right]}^{\left(s\right)\mathrm{y}}\left(\sigma \right)$): denotes the completion time of a customer order, say i, of agent x (or agent y) scheduled in the k-th position in $\sigma$, i$\in$N, where

$C_{\left[k\right]}^{\left(s\right)x}\left(\sigma \right)={\max }_{\upsilon \in M}\left\{\max \left\{C_{\left[k-1\right]}^{\left(s\right)x}\left(\sigma \right),r_i^{\left(s\right)x}\right\}+t_{i\upsilon }^{\left(s\right)x}\right\}$ or

$C_{\left[k\right]}^{\left(s\right)y}\left(\sigma \right)={\max }_{\upsilon \in M}\left\{\max \left\{C_{\left[k-1\right]}^{\left(s\right)y}\left(\sigma \right),r_i^{\left(s\right)y}\right\}+t_{i\upsilon }^{\left(s\right)y}\right\}$, s = 1, 2.

2.1. Problem Definition

Let DPm denote the m dedicated machines used to produce m components of each job, and sdpt represent the scenario-dependent processing times of the components. Following the standard three-field classification of Perez-Gonzalez and Framinan [11] and Framinan et al. [9], to define the DPm → 0|r_i, 2-agent, sdpt|($ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}$ /$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}\le U$) problem considered in this paper, consider a set of n customer orders belonging to n different customers, in which a customer order has m components to be processed on m parallel machines, and where a given component of each customer order is processed on the pre-specified dedicated machine. These orders are from two agents, agent x and agent y. Due to the fact that some factors of significant uncertainties are present, it is assumed that a customer order has scenario-dependent component processing times $t_{iv}^{\left(s\right)x}$($t_{iv}^{\left(s\right)y}$) on machine $M_v,v\in M$, and a scenario-dependent ready time $r_i^{\left(s\right)x}$(or $r_i^{\left(s\right)y}$), s = 1, 2. Furthermore, the orders of agent y have a common upper limit U on their completion time.

With the above description and notations, the problem considered in this paper, the DPm → 0|r_i, 2-agent, sdpt| ($ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}$/$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}\le U$) problem, is one of finding an appropriate schedule, among all possible schedules, to minimize the maximum of the total completion time of all customer orders that belong to agent x and are subject to the following constraint—that the maximum of the total completion time of the customer orders belonging to agent y should be no greater than an upper bound (U) for scenario s = 1, 2. In other words, we wish to determine a robust optimal schedule to minimize $ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}$, and which is subject to $ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}\le U$. Given that Wu et al. [5] and Lenstra et al. [36] respectively showed that the problems $1\left|r_j\right|\sum C_j^x:\sum C_j^y$ and $1\left|r_j\right|\sum C_j\le U$ are NP-hard, it follows that the problem considered in this paper is also an NP-hard problem. This can be summarized as follows.

Theorem 1.

The DPm → 0|r_i, 2-agent, sdpt|($ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}$/$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}\le U$) problem is NP-hard.

Below is an illustrative example with a feasible solution to the problem under study. Consider a two-agent order schedule $\sigma$, with two orders and two scenarios to be processed on two machines; then, let $(t_{11}^{\left(1\right)x},t_{12}^{\left(1\right)x},r_1^{\left(1\right)x})=\left(5, 6,0\right)$ $(t_{11}^{\left(2\right)x},t_{12}^{2\left(1\right)x},r_1^{\left(2\right)x})=(7,3,$1) for order x, while $(t_{11}^{\left(1\right)y},t_{12}^{\left(1\right)y},r_1^{\left(1\right)y})=\left(2,1,2\right)$ $(t_{11}^{\left(2\right)y},t_{12}^{\left(2\right)y},r_1^{\left(2\right)y})=(3,2,2$) for order y. Let $\sigma$ = ($O_1^x,O_2^y$), and let U = 12.

The completion times for each order scheduled in $\sigma$ for scenario 1 are calculated as follows:

$$C_{\left[1\right]}^{\left(1\right)x}\left(\sigma \right)=max\left\{r_1^{\left(1\right)x}+t_{11}^{\left(1\right)x}, r_1^{\left(1\right)x}+t_{12}^{\left(1\right)x} \right\}=max\left\{0+5, 0+6\right\}=6$$

$$C_{\left[2\right]}^{\left(1\right)y}\left(\sigma \right)=max\left\{\max \left\{5,2\right\}+2,\max \left\{6,2\right\}+1\right\}=max\left\{7, 7\right\}=7$$

The completion times for each order scheduled in $\sigma$ for scenario 2 are calculated as follows:

$$C_{\left[1\right]}^{\left(2\right)x}\left(\sigma \right)=max\left\{r_1^{\left(2\right)x}+t_{11}^{\left(2\right)x}, r_1^{\left(2\right)x}+t_{12}^{\left(2\right)x} \right\}=max\left\{1+7, 1+3\right\}=8$$

$$C_{\left[2\right]}^{\left(2\right)y}\left(\sigma \right)=max\left\{\max \left\{8,2\right\}+3,\max \left\{4,2\right\}+2\right\}=max\left\{11, 6\right\}=11$$

Since $ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}=\max \left\{C_{\left[2\right]}^{\left(1\right)y}\left(\sigma \right),C_{\left[2\right]}^{\left(2\right)y}\left(\sigma \right)\right\}=\max \left\{7,11\right\}\le 12=U$, the objective function is given by the following equation:

$$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}=\max \left\{C_{\left[1\right]}^{\left(1\right)x}\left(\sigma \right),C_{\left[1\right]}^{\left(2\right)x}\left(\sigma \right)\right\}=\max \left\{6,8\right\}=8$$

2.2. Proposed Lower Bound

To solve this intractable problem, we will derive several properties and a lower bound to be used in a branch-and-bound algorithm (Ogan and Azizoglu [37]) to accelerate the process finding an exact solution. Assume that $\delta$ is a specified partial schedule with q customer orders, while ${\delta }^{\prime }$ denotes the undetermined set containing (n − q) customer orders. For the simplification of the notation, we set (n − q) = $n_x+n_y$, where $n_x$($n_y)$ denotes the number of unscheduled customer orders for agent x (agent y). Furthermore, we set ${\theta }_v^{\left(s\right)}$ to be the completion time of the last customer order on $M_v,v\in M$, for $s=1,2,$ in $\delta$, and ${\theta }^{\left(s\right)}={\max }_{\upsilon \in \mathrm{M}}\left\{{\theta }_v^{\left(s\right)}\right\}$, i.e., $C_{\left[q\right]}^{\left(s\right)}\left(\delta \right)={\theta }^{\left(s\right)}$. To obtain the proposed lower bound, we only choose $n_x$ x-agent orders to schedule on the (q + 1)-th, (q + 2)-th, …, and (q + $n_x$)-th positions appended immediately after the q-th position. From the previous definitions, therefore, it follows that:

$$C_{\left[q+1\right]}^{\left(s\right)x}\left(\sigma \right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_{\left[q\right]}^{\left(s\right)},r_{\left[q+1\right]}^{\left(s\right)x}\right\}+t_{\left[q+1\right]v}^{\left(s\right)x}\right\}>{\theta }_{}^{\left(s\right)}+{{\displaystyle \sum }}_{\upsilon \in \mathrm{M}}^{}{t}_{\left[q+1\right]v}^{\left(s\right)x}/m$$

$$C_{\left[q+2\right]}^{\left(s\right)x}\left(\sigma \right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_{\left[q+1\right]}^{\left(s\right)},r_{\left[q+2\right]}^{\left(s\right)x}\right\}+t_{\left[q+2\right]v}^{\left(s\right)x}\right\}>{\theta }_{}^{\left(s\right)}+\left\{{{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{\left[q+1\right]v}^{\left(s\right)x}+{{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{\left[q+1\right]v}^{\left(s\right)x}\right\}/m.$$

Similarly, we have the following:

$$C_{\left[q+n_x\right]}^{\left(s\right)x}\left(\sigma \right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_{\left[q+n_x-1\right]}^{\left(s\right)},r_{\left[q+n_x\right]}^{\left(s\right)x}\right\}+t_{\left[q+n_x\right]v}^{\left(s\right)x}\right\}>{\theta }_{}^{\left(s\right)}+{{\displaystyle \sum }}_{j=1}^{n_x}{{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{\left[q+j\right]v}^{\left(s\right)x}/m.$$

Therefore, the total completion times for $n_x$x-agent orders can be summarized as follows:

$${{\displaystyle \sum }}_{j=1}^{n_x}{C}_{\left[q+j\right]}^{\left(s\right)x}\left(\sigma \right)>n_x{\theta }_{}^{\left(s\right)}+{{\displaystyle \sum }}_{j=1}^{n_x}\left(n_x-j+1\right){{\displaystyle \sum }}_{\upsilon \in \mathrm{M}}^{}{t}_{\left[q+j\right]v}^{\left(s\right)x}/m,\mathrm{for}s=1,2.$$

For the schedule $\sigma$, we have the equation below:

$${{\displaystyle \sum }}_{j=1}^n{C}_{\left[j\right]}^{\left(s\right)x}\left(\sigma \right)>{{\displaystyle \sum }}_{j=1}^q{C}_{\left[j\right]}^{\left(s\right)x}\left(\sigma \right)+n_x{\theta }_{}^{\left(s\right)}+{\displaystyle \sum }_{j=1}^{n_x}\left(n_x-j+1\right){{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{\left[q+j\right]v}^{\left(s\right)x}/m,\mathrm{for}s=1,2.$$

Therefore, we have the following equation:

$$\begin{gathered} ma{x}_{s=1,2}\left\{{{\displaystyle \sum }}_{j=1}^q{C}_{\left[j\right]}^{\left(s\right)x}\left(\sigma \right)+{{\displaystyle \sum }}_{j=1}^{n_x}{C}_{\left[q+j\right]}^{\left(s\right)x}\left(\sigma \right)\right\}>\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {{\displaystyle \sum }}_{s=1}^2{{\displaystyle \sum }}_{j=1}^q{C}_{\left[j\right]}^{\left(s\right)x}\left(\sigma \right)+n_x{\mathsf{\theta }}_{}^{\left(s\right)}+{{\displaystyle \sum }}_{j=1}^{n_x}\left(n_x-j+1\right){{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{\left[q+j\right]v}^{\left(s\right)x}/2m. \end{gathered}$$

Using the fact that the maximum value is larger than the mean value and a lemma, as defined by Hardy et al. [38], the proposed lower bound (lbdd) can be summarized as follows.

Proposition 1.

$lbdd={\displaystyle \sum }_{s=1}^2{{\displaystyle \sum }}_{j=1}^q{C}_{\left[j\right]}^{\left(s\right)x}\left(\sigma \right)+n_x{\theta }_{min}^{\left(s\right)}+{\displaystyle \sum }_{j=1}^{n_x}\left(n_x-j+1\right)t_{\left(q+j\right)*}^{\left(s\right)x}/2m$, where $t_{q+j*}^{\left(s\right)x}={{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{\left[q+j\right]v}^{\left(s\right)x}$, and $t_{\left(q+1\right)*}^{\left(s\right)x}\le \dots \le t_{\left(q+n_x\right)*}^{\left(s\right)x}$ is the non-decreasing order of $\{t_{q+1*}^{\left(s\right)x}$, …, $t_{q+n_x*}^{\left(s\right)x}\}$.

2.3. Dominance Properties

Consider two order sequences $,\sigma =(\delta ,O_i^{\Phi }$, $O_j^{\Phi },{\delta }^{\prime })$ and ${\sigma }^{\prime }=(\delta ,O_j^{\Phi }$, $O_i^{\Phi },{\delta }^{\prime })$, in which $\delta ,{\delta }^{\prime }$ are two subsequences. To prove that “$\sigma$ is no worse than sequence ${\sigma }^{\prime }$”, we need to confirm that for s = 1, 2, $C_i^{\left(s\right)x}\left(\sigma \right)+C_j^{\left(s\right)x}\left(\sigma \right)<C_j^{\left(s\right)x}\left({\sigma }^{\prime }\right)+C_i^{\left(s\right)x}\left({\sigma }^{\prime }\right)$ and $C_j^{\left(s\right)x}\left(\sigma \right)\le C_i^{\left(s\right)x}\left({\sigma }^{\prime }\right)$. Furthermore, recall that ${\mathsf{\theta }}_v^{\left(s\right)}$ is the completion time of the last job in the δ subsequence of schedule $\sigma \mathrm{or} {\sigma }^{\prime }$ on machine $M_v,v\in M$, s = 1, 2. With these defined conditions, we describe the following six properties to curtail the search for an optimal solution. The details of proofs of Property 1 are provided, while other properties are omitted because they are similar to those of Property 1.

Property 1.

For any two customer orders $O_i^{\Phi }$ and $O_j^{\Phi }$ from agent x to be scheduled adjacent to each other, if $\forall$ s = 1, 2, $ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\}\le ma{x}_{\upsilon \in M}\left\{t_{jv}^{\left(s\right)x}\right\}$, $r_j^{\left(s\right)x}>r_i^{\left(s\right)x}\ge ma{x}_{\upsilon \in M}\left\{{\theta }_v^{\left(s\right)}\right\}$, and $r_j^{\left(s\right)x}\le r_i^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\}$, then sequence $\sigma$ is not worse than sequence ${\sigma }^{\prime }$.

Proof. 

Suppose that there are q customer orders scheduled in the subsequence δ of $\sigma$(or ${\sigma }^{\prime }$), and that the complete time of the last job is $C_{\left[q\right]}^{\left(s\right)}$. According to the definition:

$$C_i^{\left(s\right)x}\left(\sigma \right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_{\left[q\right]}^{\left(s\right)x},r_i^{\left(s\right)x}\right\}+t_{i\upsilon }^{\left(s\right)x}\right\},$$

(1)

$$C_j^{\left(s\right)x}\left(\sigma \right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_i^{\left(s\right)x}\left(\sigma \right),r_j^{\left(s\right)x}\right\}+t_{j\upsilon }^{\left(s\right)x}\right\},$$

(2)

$$C_j^{\left(s\right)x}\left({\sigma }^{\prime }\right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_{\left[q\right]}^{\left(s\right)x},r_j^{\left(s\right)x}\right\}+t_{j\upsilon }^{\left(s\right)x}\right\},$$

(3)

$$C_i^{\left(s\right)x}\left({\sigma }^{\prime }\right)=ma{x}_{\upsilon \in M}\left\{\mathit{max}\left\{C_j^{\left(s\right)x}\left({\sigma }^{\prime }\right),r_i^{\left(s\right)x}\right\}+t_{i\upsilon }^{\left(s\right)x}\right\}.$$

(4)

Applying the given condition $r_j^{\left(s\right)x}>r_i^{\left(s\right)x}\ge ma{x}_{\upsilon \in M}\left\{{\theta }_v^{\left(s\right)}\right\}$ to (1) and (3), we have the following:

$$C_i^{\left(s\right)x}\left(\sigma \right)=r_i^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\},C_j^{\left(s\right)x}\left({\sigma }^{\prime }\right)=r_j^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{jv}^{\left(s\right)x}\right\},\mathrm{for}s=1,2.$$

Further applying the given condition $r_j^{\left(s\right)x}\le r_i^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\}$ to $C_i^{\left(s\right)x}\left(\sigma \right)$:

$$C_j^{\left(s\right)x}\left(\sigma \right)=r_i^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\}+ma{x}_{\upsilon \in M}\left\{t_{jv}^{\left(s\right)x}\right\},$$

(5)

and further applying $r_j^{\left(s\right)x}>r_i^{\left(s\right)x}$ to $C_j^{\left(s\right)x}\left({\sigma }^{\prime }\right)$, we have the following equation:

$$C_i^{\left(s\right)x}\left({\sigma }^{\prime }\right)=r_j^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{jv}^{\left(s\right)x}\right\}+ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\}.$$

(6)

Hence, from Equations (5) and (6), we have the following:

$$C_i^{\left(s\right)x}\left({\sigma }^{\prime }\right)-C_j^{\left(s\right)x}\left(\sigma \right)=r_j^{\left(s\right)x}-r_i^{\left(s\right)x}>0.$$

Combine Equations (1), (2), and (5) with inequality (6), we have the following:

$[C_j^{\left(s\right)x}\left({\sigma }^{\prime }\right)+C_i^{\left(s\right)x}\left({\sigma }^{\prime }\right)\left]-\right[C_i^{\left(s\right)x}\left(\sigma \right)+C_j^{\left(s\right)x}\left(\sigma \right)]=2(r_j^{\left(s\right)x}-r_i^{\left(s\right)x}$) + $\left({\max }_{\upsilon \in \mathrm{M}}\left\{t_{jv}^{\left(s\right)x}\right\}-{\max }_{\upsilon \in \mathrm{M}}\left\{t_{iv}^{\left(s\right)x}\right\}\right)$ > 0, indicating that sequence $\sigma$ is not worse than sequence ${\sigma }^{\prime }$. □

Property 2.

For any two customer orders $O_i^{\Phi }$ and $O_j^{\Phi }$ from agent x to be scheduled adjacent to each other, if $\forall$ s = 1, 2, $r_i^{\left(s\right)x}\ge ma{x}_{\upsilon \in M}\left\{{\theta }_v^s\right\}$, and $r_j^{\left(s\right)x}>r_i^{\left(s\right)x}+ma{x}_{\upsilon \in M}\left\{t_{iv}^{\left(s\right)x}\right\}$, then sequence $\sigma$ is not worse than sequence ${\sigma }^{\prime }$.

Property 3.

For any two customer orders $O_i^{\Phi }$ and $O_j^{\Phi }$ from agent x to be scheduled adjacent to each other, if $\forall$ s = 1,2, $max\left\{r_i^{\left(s\right)x}, r_j^{\left(s\right)x}\right\}\le mi{n}_{\upsilon \in M}\left\{{\theta }_v^s\right\}$, and $ma{x}_{\upsilon \in M}\left\{{\theta }_v^s+t_{iv}^{\left(s\right)x}\right\}<ma{x}_{\upsilon \in M}\left\{{\theta }_v^s+t_{jv}^{\left(s\right)x}\right\}$, then sequence $\sigma$ is not worse than sequence ${\sigma }^{\prime }$.

Property 4.

For any two customer orders $O_i^{\Phi }$ and $O_j^{\Phi }$ from an agent, x to be scheduled adjacent to each other, if $\forall s=1,2, max\left\{r_i^{\left(s\right)x}, r_j^{\left(s\right)x}\right\}\le mi{n}_{\upsilon \in M}\left\{{\theta }_v^s\right\}$, and $ma{x}_{\upsilon \in M}\left\{{\theta }_v^s+t_{iv}^{\left(s\right)x}\right\}<ma{x}_{\upsilon \in M}\left\{{\theta }_v^s+t_{jv}^{\left(s\right)x}\right\}$, then sequence $\sigma$ is not worse than sequence ${\sigma }^{\prime }$.

Property 5.

If $\exists O_i^y\in {\delta }^{\prime }$ and $\forall s=1, 2, {{\displaystyle \sum }}_{k\in \delta }{C}_{\left[k\right]}^{\left(s\right)y}>U,$ then $\sigma =\left(\delta ,{\delta }^{\prime }\right)$ can be eliminated from the search for an optimal (robust) schedule. (Note that $\delta$ denotes the scheduled part, while ${\delta }^{\prime }$ denotes the unscheduled part.)

Property 6.

If $\exists O_i^y\in {\delta }^{\prime }$ and $\forall s=1, 2,$${{\displaystyle \sum }}_{k\in \delta }{C}_{\left[k\right]}^{\left(s\right)y}+ma{x}_{\upsilon \in M}\left\{max\{{\theta }_v^{\left(s\right)},r_i^{\left(s\right)y}\right\}+t_{iv}^{\left(s\right)y}\}>U,$ then $\sigma =\left(\delta ,{\delta }^{\prime }\right)$ can be eliminated from the search for an optimal (robust) schedule.

2.4. Proposed Branch-and-Bound Algorithm

The lower bound described above, the six dominance properties, and the best near-optimal solution obtained from the proposed heuristics in the next section are used to describe the proposed branch-and-bound algorithm (B&B) by using the depth-first method for the small-sized job case in this problem. The details of the B&B algorithm are provided as follows.

Algorithm 1 A branch-and-bound method

00:  Input a set $M=\left\{1, 2, \dots , m\right\}$ of m parallel machines and a set $N=\left\{1,2,\dots ,n\right\}$ of n orders with scenario-dependent component processing times $t_{iv}^{\left(s\right)x}$($t_{iv}^{\left(s\right)y}$) and ready times $r_i^{\left(s\right)x}$ (or $r_i^{\left(s\right)y})$, i = 1,2,…,n}; criterion function: Minimize$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}$, subject to$ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)y}\left(\sigma \right)\right\}\le U$. 01:  Input the best heuristic solution with an upper bound. 02:  Start to branch from the root node by appending each order to create a new node. 03:  For each active node, (i) find the lower bound based on proposition 1 in ${\pi }^c$; then, (ii) evaluate if the lower bounds $\ge$ the incumbent upper bound and cut those nodes and all nodes below them in the branching tree. 04:  Delete the unwanted nodes from the branching tree by property 1 to property 6. 05:  To determine if the node is complete or not, we do the following: (i) find its criterion function equal to $ma{x}_{s\in \left\{1,2\right\}}\left\{\sum C_i^{\left(s\right)x}\left(\sigma \right)\right\}$; (ii) evaluate if this criterion function is smaller than the upper bound, then update it with the new one and keep the corresponding schedule; (iii) For the remaining nodes, branch from the node with the minimum lower bound to create a set of new active nodes. 06:  Repeat step 02 through step 05 until all nodes have been visited, and set the final complete schedule as σ. 07:  Output the final complete schedule σ as an optimal solution.

3. Simulated Annealing Hyper-Heuristic

In this section, we propose three heuristics and two variants of a simulated annealing hyper-heuristic with different cooling temperatures for finding optimal or near-optimal solutions to the problem. To investigate the effect of the customer orders’ component scenario-dependent processing times, we consider three types of combinations of the processing times $t_{i\upsilon }^{\Phi \mathrm{s}}$ for each customer order, s = 1, 2; i.e.:

$\left(\mathrm{i}\right)0.5ma{x}_{\upsilon \in M}\{r_i^{\left(1\right)\Phi }+t_{i\upsilon }^{\left(1\right)\Phi }\}+0.5ma{x}_{\upsilon \in M}\{r_i^{\left(2\right)\Phi }+t_{i\upsilon }^{\left(2\right)\Phi }\}$;

$\left(\mathrm{ii}\right)0.5mi{n}_{\upsilon \in M}\{r_i^{\left(1\right)\Phi }+t_{i\upsilon }^{\left(1\right)\Phi }\}+0.5mi{n}_{\upsilon \in M}\{r_i^{\left(2\right)\Phi }+t_{i\upsilon }^{\left(2\right)\Phi }\}$;

and

$\left(\mathrm{iii}\right)0.5\left(r_i^{\left(1\right)\Phi }+{{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{iv}^{\left(1\right)\Phi }\right)+0.5(r_i^{\left(2\right)\Phi }+{{\displaystyle \sum }}_{\upsilon \in M}^{}{t}_{iv}^{\left(2\right)\Phi })$.

Based on these three formulas, we schedule customer orders using the smallest processing time (SPT) first rule to obtain three initial schedules, S_MM, S_mm, and S_mean. To ensure good quality for the approximate solutions, we consider the two following policies: (1) We improve each of the S_MM, S_mm, and S_mean by using a pairwise interchange method and record these three heuristics as Mpi, mpi, and meanpi, respectively. (2) We use these three improved heuristics as three initial seeds in a simulated annealing super-heuristic algorithm. They are termed as SAHM, SAHm, and SAHmean. We propose a simulated annealing hyper-heuristic algorithm (SAH), in which the low-level heuristics may be based on the proposed eight candidate mutation operators. They are termed as HL₁, HL₂, ..., and HL₈.

The frequency (coded as $f_i$) of the improvement of a low-level heuristic HL_i is recorded when the accumulated performance of HL_i is obtained in each cycle. The selection probability P(HL_i) for HL_i is defined as $f_i/{{\displaystyle \sum }}_{i=1}^8{f}_i$. These low-level heuristics are selected randomly based on the values of the selection probabilities (P(HL₁), P(HL₂), …, P(HL₈)). For the first cycle (or I_cmax=1), the selection probabilities of these low-level heuristics are set to be equal, i.e., $f_i=1$, i = 1,…, 8. For I_cmax > 1, the probabilities are determined according to the corresponding low-level heuristic additive performance. To ensure the diversity of all low-level heuristics in the list, we set the minimum value $f_i=\max \left\{f_i,1\right\}$. With these definitions, the details of the eight proposed low-level heuristics are described below:

HL₁: 

Two-order swap heuristic.

HL₂: 

One step to the right heuristic.

HL₃: 

Two steps to the right heuristic.

HL₄: 

Two randomly selected orders, O_i, O_j, where 1 ≤ I < j < n, and swap both orders with their immediately succeeding orders.

HL₅: 

Two randomly selected orders, O_i, O_j, where 1 < i < j ≤ n, and swap both orders with their closest preceding orders.

HL₆: 

Two randomly selected orders, O_i, O_j, where 1 ≤ i < i + 1 < j ≤ n, swap the order in front O_j with its immediately succeeding order, and swap another O_i with its closest preceding order.

HL₇: 

Two randomly selected orders O_i, O_j, where 1 < i < j < n, swap the order in front O_j with its closest preceding order, and swap another O_i with its immediately succeeding order.

HL₈: 

The current order of jobs is reversed.

For any two schedules σ and ${\sigma }_t$, we set $RTC\left(\sigma \right)=\underset{s=1,2}{\max }{\displaystyle \sum }_{i\in O^x}{C}_{\left[i\right]}^{\left(s\right)x}\left(\sigma \right)$ and $RTC\left({\sigma }_t\right)=\underset{s=1,2}{\max }{\displaystyle \sum }_{i\in O^x}{C}_{\left[i\right]}^{\left(s\right)x}\left({\sigma }_t\right)$ as the values of their objective function. Let I_cmax denote the maximum number of cycles for SAHs, and ($T_i$,$T_{\mathrm{f}}$, Cf, Nr) denote the initial temperature, the final temperature, the cooling factor, and the maximum number of iterations per cycle. Then, we describe below the steps of the SAHs algorithms.

Algorithm 2 The steps of the SAHs procedure

01: Input Icmax, Nr, Ti, Tf, and Cf 02: Input each initial sequence σ and its objective function RTC(σ) 03: Set T = Ti; fk = 1, k = 1,…, 8; i_no = 0 04: Do while {i_no ≤ Icmax and T > Tf} 05:   num = 0 05:   Do while {num ≤ Nr} 06:    Randomly select an HLk by the roulette wheel according to their probabilities 07:    Utilize it to form a new sequence σt and to find the RTC(σt) of σt 08:    If Lt = f(st), RTC(σt) < RTC(σ), replace σ with σt; otherwise, replace σ with σt using the probability $\exp \left(-\left(RTC\left({\sigma }_t\right)-RTC\left(\sigma \right)\right)/RTC\left(\sigma \right)T\right)$ 09:     num = num + 1 and $f_k=f_k+1$ 10:   End do 11:   Update $f_k/{{\displaystyle \sum }}_{k=1}^8{f}_k$ for each HLk; $T=T\times \mathrm{Cf}$; i_no = i_no +1 12: End do 13: Output final sequence σ and its objective function RTC(σ).

In addition to the three heuristics, Mpi, mpi, and meanpi, developed to solve the studied problem, we employ two variants of the SAH algorithms labeled as (SAHM, SAHm, SAHmean) and (SAHMb, SAHmb, SAHmeanb), respectively. Both variants use the three initial schedules S_MM, S_mm, S_mean, respectively, similar to the seeds in the above SAH algorithm procedure. However, SAHMb, SAHmb, and SAHmeanb use the modified temperature formula from Bai et al. [34] as $T=T/\left(1+\beta T\right)$, where $\beta =\left(T_i-T_f\right)*Nr$/(I_cmax *Nr*$T_i*T_f)$, to fit the studied problem.

4. Computational Results

This section presents the results of extensive simulation experiments, performed to determine the effectiveness of the branch-and-bound algorithm, three local heuristics (Mpi, mpi, meanpi), and six SAHs algorithms (SAHM, SAHm, SAHmean, SAHMb, SAHmb, SAHmeanb) in solving the considered problem. The instances generated are described as follows. The processing times of customer orders were generated from the uniform distribution U(1, 100) for scenario 1 and from U(1, 200) for scenario 2. Following the design of Reeves [39], the order release (ready) times were generated from the uniform distribution U(0, 100nλ) for scenario 1 and from the uniform distribution U(0, 200nλ) for scenario 2, where n is the number of orders and λ is a control variable. Three different types of problem instances were created according to the values of λ at 0.25, 0.5, and 0.75. To generate a feasible problem instance, we first generated $n_y$ orders of agent y to put before $n_x$ orders, computed the total completion times for scenario 1 and scenario 2, and then set the value of the upper bound at $U=2.5\times \underset{s=1,2}{\max }\left\{{{\displaystyle \sum }}_{i\in O^y}^{}{C}_{\left[i\right]}^{\left(s\right)y}\right\}$, where 2.5 was a tested number that ensured a feasible instance.

4.1. Tuning the Parameters in SAH Algorithms

To obtain appropriate values of related parameters in the SAH algorithms, the number of customer orders was set at n = 10 and the number of machines was set at m = 3, where $\left(n_x,n_y\right)=\left(5,5\right)$. The component processing time ($t_{i\upsilon }^{\left(1\right)\Phi })$ of an order was generated from a uniform distribution U(1, 50), while the processing time $t_{i\upsilon }^{\left(2\right)\Phi }$ was generated from another uniform distribution U(1, 100). Meanwhile, the ready time ($r_i^{\left(1\right)\Phi })$ of order I for scenario one (s = 1) was generated from a uniform distribution U(1, 50·λ·n), while the ready time ($r_i^{\left(2\right)\Phi })$ of order I for scenario two (s = 2) was generated from another uniform distribution U(1, 100·λ·n), where $\lambda$ was taken as 0.25. The value of the upper bound was set at $2.5\times \underset{s=1,2}{\max }\left\{{{\displaystyle \sum }}_{i\in O^y}^{}{C}_{\left[i\right]}^{\left(s\right)y}\right\}$, where 2.5 was a tested number that ensured a feasible instance. One hundred problem instances were tested for each combination of parameters. The final temperature was set at T_f = 10⁻⁸ to tune all the values of the parameters in the SAH algorithms.

The average error percentage (AEP) was used as the index of performance; i.e., AEP = [(HA − O*)/O*] × 100[%], where HA is the total completion time of the orders belonging to agent x obtained by using three heuristics (Mpi, mpi, meanpi) or six SAHs (SAHM, SAHm, SAHmean, SAHMb, SAHmb, SAHmeanb), and O* is the total completion time of the orders belonging to agent x obtained by using the branch-and-bound algorithm. For simplicity, all the determined parameters of SAHM will be used in the other SAHs.

The proposed SAHs have four parameters, namely initial temperature (Ti), cooling factor (Cf), the runs of low-level loops for local improvement (Nr), and the runs of top-level loops (I_cmax). With Cf = 0.95, Nr = 1000, and I_cmax = 15, a simulation was conducted under a variation in Ti, from 10⁻⁵ to 10⁶, with multiples of 10 times each. As shown in Figure 1a, the AEP reached the lowest point when Ti was equal to 0.1; hence, the optimal setting for Ti was 0.1.

With Ti = 0.1, Nr = 1000, and I_cmax = 15, a simulation was conducted under a variation in Cf, from 0.90 to 0.99, with increments of 0.1 each. As shown in Figure 1b, the AEP approximated zero (below 0.01%) when Cf was equal to 0.91, indicating an effective reduction of AEP; hence, the optimal setting for Cf was 0.91.

With Ti = 0.1, Cf = 0.91, and I_cmax = 15, a simulation was conducted under a variation in Nr, from 200 to 1000, with increments of 50 each. As shown in Figure 1c, the AEP reached its minimum when Nr was equal to 850, also indicating an effective reduction of AEP with an increase in Nr; hence, the optimal setting for Nr was 850.

Finally, with Ti = 0.1, Cf = 0.91, and Nr = 850, a simulation was conducted under a variation in I_cmax, from 10 to 30, with increments of two each. As shown in Figure 1d, the AEP was the lowest when I_cmax was equal to 18, 22, 24, 26, 28, and 30; hence, the optimal setting for I_cmax was 22.

After previous tested results, we adopted the parameters Ti = 0.1, Cf = 0.91, Nr = 850, and I_cmax = 22 used in SAH algorithms for the following tested problems.

4.2. Small-n Experimental Results Analysis

For a small number of orders, the number of jobs is set at n = 8 and 10, and the machine number is set at m = 2, 3, and 4. To evaluate their impacts on the algorithms, the numbers of x-agent order and y-agent order $(n_x,n_y)$ are set at (2, 6), (4, 4), and (6, 2) for n = 8, and at (3, 7), (5, 5), and (7, 3) for n = 10. A set of 100 instances were randomly generated for each combination of n, m,$n_x$, and λ. Consequently, a total of 5400 problem instances were tested. The algorithms were set to skip to the next set of data if the number of nodes exceeded 10⁸.

The impacts of the number of orders (n), machine number (m), the number of x-agent orders (n_x), and the control parameter of the range of release dates (λ) on the performance of the branch-and-bound algorithms, three heuristics, and the six SAH algorithms are shown in

Table 2. A summary of the performance of the branch-and-bound algorithm.

n	m	Node	CPU_Time
8	2	3540	0.0299
	3	4649	0.0442
	4	4707	0.0489
10	2	141,867	0.7474
	3	140,400	0.8877
	4	187,833	1.3859
	nx
8	2	3913	0.0409
	4	2891	0.0331
	6	6092	0.0491
10	3	127,915	0.8616
	5	117,367	0.7827
	7	224,819	1.3766
	$\lambda$
8	0.25	4047	0.0406
	0.50	4001	0.0388
	0.75	4848	0.0438
10	0.25	112,313	0.7410
	0.50	144,859	0.9204
	0.75	212,928	1.3595

.

To measure the performance of the branch-and-bound algorithm, we recorded the average number of nodes and the average execution times (in seconds) for n = 8 and 10. Table 2 presents its capability, and the mean nodes and CPU times increase as n dramatically increases from 8 to 10 (columns 3 and 4 of Table 1). As demonstrated in Table 2, the mean nodes and CPU times increase as m increases regardless of the number of jobs n. When the number of x-agent orders (n_x in Table 2) becomes equivalent to half of the number of total orders (n), the mean nodes and CPU times are lower than other values of n_x independent of the number of jobs n. As the range of release dates becomes larger (λ increases from 0.25 to 0.75), the search space for the studied problem becomes wider; thus, the nodes and CPU times increase, especially for n = 10.

For the three heuristics, Mpi, mpi, and meanpi, and the six SAH algorithms, we record the AEP of the objective function; the AEP = [(HA − O*)/O*] × 100[%], where HA and O* are the values of the objective function obtained by running heuristics/SAHs and the branch-and-bound algorithm, respectively. The results are summarized in

Table 3. Performance of heuristics and SAH algorithms for a small n.

		AEP(%)
n	m	Mpi	mpi	meanpi	SAHM	SAHm	SAHmean	SAHMb	SAHmb	SAHmeanb
8	2	10.33	10.96	11.41	0.0000	0.0064	0.0000	0.0297	0.0315	0.0238
	3	10.89	11.64	11.76	0.0001	0.0009	0.0004	0.0332	0.0435	0.0312
	4	11.21	11.71	11.89	0.0000	0.0001	0.0000	0.0311	0.0261	0.0298
10	2	15.84	16.29	16.59	0.0029	0.0046	0.0026	0.0482	0.0603	0.0538
	3	16.35	16.38	16.74	0.0007	0.0055	0.0062	0.0435	0.0320	0.0498
	4	16.83	17.37	17.31	0.0048	0.0014	0.0048	0.0783	0.0515	0.0545
	nx
8	2	8.23	8.77	8.88	0.0001	0.0001	0.0004	0.0698	0.0717	0.0694
	4	14.03	15.21	15.56	0.0000	0.0073	0.0000	0.0242	0.0294	0.0154
	6	10.17	10.33	10.62	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
10	3	12.94	13.54	13.48	0.0077	0.0115	0.0136	0.1567	0.1314	0.1442
	5	18.17	18.93	19.09	0.0007	0.0000	0.0000	0.0133	0.0124	0.0139
	7	17.92	17.57	18.07	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
	$\lambda$
8	0.25	17.62	18.76	19.25	0.0001	0.0073	0.0004	0.0698	0.0700	0.0641
	0.50	9.78	10.27	10.51	0.0000	0.0001	0.0000	0.0175	0.0306	0.0172
	0.75	5.04	5.28	5.30	0.0000	0.0000	0.0000	0.0068	0.0006	0.0035
10	0.25	26.45	27.02	27.60	0.0069	0.0097	0.0108	0.1340	0.0967	0.1247
	0.50	14.07	14.49	14.37	0.0002	0.0018	0.0014	0.0244	0.0349	0.0219
	0.75	8.50	8.53	8.68	0.0012	0.0000	0.0014	0.0116	0.0123	0.0115
mean		13.58	14.06	14.28	0.0014	0.0032	0.0023	0.0440	0.0408	0.0405

.

Concerning the impact of the machine number (m) on the three heuristics and six SAH algorithms, the results in Table 3 show that there is not much difference in the performance of the AEP for heuristics or SAHs. As for the impacts of the number of x-orders (n_x) on the three heuristics and six SAH algorithms, it was observed that the three heuristics performed the worst when n_x = n_y, and the AEPs decreased (except for the SAHm at n = 8) as n_x increased for the six SAHs algorithms. It can also be observed in Table 2 that the AEP declines as the value of λ increases from 0.25 to 0.75 for almost all nine heuristics/algorithms (except for the SAHM at n = 10). This means that the wider the release dates of job orders, the better the solution quality resulting from these nine heuristics/algorithms. When the three heuristics and six SAH algorithms are divided into three groups, the best group is (SAHM, SAHm, SAHmean) with the lowest AEPs at (0.0014, 0.0032, 0.0023), the second is (SAHMb, SAHmb, SAHmeanb) with AEPs at (0.0440, 0.0408, 0.0405), and the worst is (Mpi, mpi, meanpi) with the AEPs at (13.58, 14.06, 14.28). Figure 2 and Figure 3 depict this peculiarity.

To further analyze the statistically significant differences of the solution quality among the three heuristics and six SAH algorithms, we made use of the SAS software in conducting an ANOVA (analysis of variance) on the AEPs for the variables, algorithms, n, m, n_x, and λ. However, we found that the normality assumption of the linear model did not hold. Therefore, the nonparametric Friedman test was utilized to differentiate the performance of the nine heuristics/SAH algorithms. Based on the sum of ranks of the AEP calculated for each of the 54 (n·m·n_x·λ = 2·3·3·3) blocks of the 100 test problem instances, the Friedman test showed that the p-value was less than 0.0001 (with a chi-square value of 374.4 and eight degrees of freedom). These test results verify that the AEP samples do not follow the same distribution at the level of significance α = 0.05. To further compare the pairwise differences among the three heuristics and six SAH algorithms, we applied the WNMT test (Wilcoxon–Nemenyi–McDonald–Thompson procedure; Holland et al. [40]).

Table 4. The rank-sum of heuristics and SAHs algorithms.

Heuristic/ Algorithm	No. of Obs.	Rank Sum (Ri, i = 1, 2, …, 9)
		Small n	Large n
Mpi	54	410.0	453.0
mpi	54	432.0	461.0
meanpi	54	454.0	382.0
SAHM	54	148.5	157.0
SAHm	54	152.5	138.0
SAHmean	54	154.5	140.0
SAHMb	54	228.5	249.0
SAHmb	54	232.5	223.0
SAHmeanb	54	217.5	227.0

Note that the critical value of the WNMT test is approximated at 88.4. (Holland et al. [40]); i.e., two algorithms are significantly different if |Ri − Rj| > 88.4.

(column 3) shows the sums of the rank of the AEPs across the 54 blocks for the nine heuristics/SHA algorithms. The WNMT test shows that (SAHM, SAHm, SAHmean, SAHMb, SAHmb, SAHmeanb) is a better performing group, with lower rank sums at (148.5, 152.5, 154.5, 228.5, 232.5, 217.5), and the group (Mpi, mpi, meanpi) is a worse performing group, with rank sums at (410.0, 432.0, 454.0). The test confirms that the performances of the two groups among the heuristics and SAH algorithms are statistically different at α = 0.05. SAHM is the best overall performer for the problem instances containing a small number of customer orders.

Concerning the usage of the eight low-level heuristics and the variation of probabilities of calling them for the SAHM, as displayed in Figure 4, HL₇ was called the most often, followed by HL₆. However, HL₈ was rarely called.

4.3. Large-n Experimental Results Analysis

In the set of experiments with a large number of orders, we set n = 80 and 100, and the machine number at m = 10, 20, and 30. We set $(n_x,n_y)$ at (20, 60), (40, 40), and (60, 20) for n = 80, and (30, 70), (50, 50), and (70, 30) for n = 100 to evaluate their impacts on the performance of the heuristics and the SAH algorithms. A set of 100 instances were randomly generated for each case. Consequently, 5400 problem instances were tested. We recorded the mean relative percentage deviance (RPD) for each of the three heuristics and six SAH algorithms. The RPD was calculated as $\mathrm{RPD}=100\left[\left(H_k-A_{*}\right)/A_{*}\right]\left[\%\right]$, where $H_k$ was obtained from each heuristic/algorithm, and $A_{*}$ was the smallest value among the H_ks from the three heuristics and six SAH algorithms.

Table 5. Summary of performance of heuristics and SAH algorithms for large n.

					RPD (%)
n	m	Mpi	mpi	meanpi	SAHM	SAHm	SAHmean	SAHMb	SAHmb	SAHmeanb
80	10	51.8	51.8	51.6	0.4	0.7	0.4	0.6	0.6	0.6
	20	52.5	52.5	52.1	0.4	0.4	0.4	0.6	0.6	0.6
	30	52.4	52.4	52.2	0.4	0.4	0.4	0.7	0.6	0.6
100	10	54.0	54.0	53.9	0.7	0.7	0.7	0.8	0.8	0.8
	20	54.5	54.5	54.2	0.6	0.6	0.7	0.8	0.8	0.8
	30	54.5	54.5	54.3	0.6	0.6	0.6	0.8	0.8	0.8
	nx
80	20	37.5	37.5	37.3	0.6	0.8	0.6	0.6	0.5	0.5
	40	54.2	54.2	53.9	0.5	0.5	0.5	0.8	0.8	0.8
	60	65.0	65.0	64.7	0.1	0.1	0.1	0.5	0.5	0.5
100	30	41.9	42.0	41.8	0.8	0.9	0.8	0.6	0.6	0.6
	50	55.4	55.4	55.1	0.7	0.7	0.8	0.9	0.9	0.9
	70	65.7	65.7	65.4	0.4	0.4	0.4	0.9	0.9	0.9
	$\lambda$
80	0.25	82.1	82.0	81.9	0.6	0.5	0.5	0.9	0.8	0.8
	0.50	42.5	42.5	42.0	0.4	0.7	0.5	0.6	0.6	0.6
	0.75	32.2	32.2	32.0	0.2	0.2	0.2	0.4	0.4	0.4
100	0.25	85.5	85.6	85.3	0.9	1.0	1.0	1.2	1.2	1.2
	0.50	43.8	43.8	43.5	0.7	0.7	0.7	0.7	0.7	0.7
	0.75	33.7	33.7	33.5	0.3	0.3	0.3	0.5	0.5	0.5
Total mean		53.3	53.3	53.1	0.6	0.6	0.6	0.7	0.7	0.7

summarizes the impacts of n_x, m, and λ.

Regarding the impact of machine number (m), the number of x-orders, and λ on the three heuristics and six SAH algorithms, the trends of the RPD for the nine heuristics/algorithms are quite similar to those of the small-sized number of orders, except that the RPDs for the six SAH algorithms are now very close, and the difference in RPD between any of the two SAHs is no greater than 0.1 (Table 5).

Moreover, from Table 5 and Figure 5, it can easily be seen that the better group is (SAHM, SAHm, SAHmean, SAHMb, SAHmb, SAHmeanb), with lower the RPDs at (0.6, 0.6, 0.6, 0.7, 0.7, 0.7), and that the worse group is (Mpi, mpi, meanpi), with RPDs at (53.3, 53.3, 53.1).

To explore the statistical significance of the differences, a Friedman test was used to differentiate the performance of the nine heuristics and SAH algorithms. Based on the sum of ranks of the RPD calculated for each of the 54 blocks of 100 test problem instances, the Freidman test showed that the p-value was less than 0.0001 (with a chi-square value at 331.1 and eight degrees of freedom). These test results verify that the RPD samples do not follow the same distribution.

For the nine heuristics/algorithms, Table 4 (column 4) shows the sums of the rank of the RPD across the 54 blocks. The WNMT test shows that (SAHM, SAHm, SAHmean, SAHmb) is the best group with lower rank sums at (157.0, 138.0, 140.0, 223.0), and the group (Mpi, mpi, meanpi) is the worst group with rank sums at (453.0, 461.0, 382.0). This test confirms that the two groups of performances are statistically different at α = 0.05. The SAHm is the best overall performer for a large-sized number of orders.

Overall, the performance of the proposed simulated annealing hyper-heuristic algorithms (SAHM, SAHm, and SAmean) appear to be experimentally validated for the considered problem.

5. Conclusions

A two-agent customer order scheduling with scenario-dependent component processing times and ready (release) times is considered in this study. The contributions are as follows: We proposed a lower bound and established six dominance properties in order to infuse them into a branch-and-bound algorithm to optimally solve this intractable problem. To obtain approximate solutions, we then developed three heuristics, Mpi, mpi, and meanpi, and six simulated annealing hyper-heuristic algorithms, SAHM, SAHm, SAHmean, SAHMb, SAHmb, and SAHmeanb. The computer simulation’s experimental results showed that SAHM, SAHm, and SAHmean performed better than the other SHA algorithms and the three heuristics. Overall, although the performance of SAHmb for large-sized problem instances is not statistically and significantly different from that of the SAHM, SAHm, and SAHmean, the results suggest that SAHM, SAHm, and SAHmean can be used to solve the problem because they are very efficient and can quickly find solutions that are very close to the optimal solutions.

For a future study, we may consider performing a sensitivity analysis of the proposed algorithms. Table 2 shows that B&B can only solve up to n = 10 problem instances due to the weak ability of its properties or a lower bound. Therefore, another future study might create more powerful properties to increase the efficacy of the B&B.