A bi-criterion sequence-dependent scheduling problem with order deliveries

Introduction

In today’s manufacturing landscape, challenges such as intense competition, a proliferation of product varieties, heightened customization demands, and shortened product life cycles have become more pressing than ever. Consequently, there is greater focus on the need to optimize manufacturing systems. Among the well-established scheduling problems, the customer order scheduling problem (COSP) has always held significant importance in the realm of manufacturing. The primary challenge of the COSP lies in determining the sequence of jobs to meet the diverse demands of customers ordering various products to be processed on a single machine.

The COSP holds particular significance in make-to-order and make-to-assembly production systems, wherein a single machine is responsible for creating various types of products. In this context, a customer may place an order for one or multiple such products, with the shipment dispatched only after all specified items have been manufactured and packaged together. The concept of COSP was initially introduced by Julien & Magazine (1990). Their work introduced a dynamic programming model and delved into the fundamental characteristics of the optimal solution for the broader case. Over the past decades, COSP on both single and multiple machines has attracted numerous researchers. For instance, Erel & Ghosh (2007) addressed COSP with the objective of minimizing the total order lead time and established for the first time that the problem is NP-hard. Su, Chen & Chen (2013) considered COSP on parallel machines dispatched in batches to minimize maximum lateness and developed three heuristics based on scheduling rules. Framinan & Perez-Gonzalez (2018) proposed a mixed-integer linear programming formulation, a constructive heuristic, and two matheuristics to minimize the total tardiness of orders. Wu et al. (2019) addressed COSP with learning effects on multiple machines to minimize total tardiness, subsequently developing four heuristics, three metaheuristics, and a branch-and-bound to solve it.

In practical scenarios, a setup time is required whenever a machine transitions from processing one product to another. Optimizing setup times is an important issue for several reasons, particularly in manufacturing and production environments. For example, shorter setup times mean less downtime between production runs. Reduced setup times enhance a company’s ability to respond quickly to changes in customer demand or shifts in the market. Longer setup times often involve increased labor costs, as more time and manpower are required for the setup process. Optimizing setup times is a strategic initiative that can positively impact various aspects of a business, including productivity, flexibility, costs, and overall competitiveness. It aligns with modern manufacturing principles such as lean manufacturing and just-in-time production, which aim to eliminate waste and enhance efficiency throughout the production process. Allahverdi & Soroush (2008) emphasized that streamlining setup times or reducing costs contributes to the timely delivery of dependable products or services. For a comprehensive review of setup time applications, readers are directed to the in-depth survey papers authored by Allahverdi (2015), Yang & Liao (1999), Allahverdi, Gupta & Aldowaisan (1999), Allahverdi et al. (2008), Cheng, Gupta & Wang (2000), and other notable works.

Despite the acknowledged significance of setup time in job scheduling, the literature on the COSP with setup time is limited. Hazır, Günalay & Erel (2008) conducted a comparative study involving four different metaheuristics to minimize average customer order flow time. Recently, de Athayde Prata, Rodrigues & Framinan (2021a) proposed two mixed-integer linear programming models and a fixed variable list algorithm for solving the COSP with sequence-dependent setup time (SDST). Furthermore, they introduced a discrete differential evolution algorithm for the same problem in de Athayde Prata, Rodrigues & Framinan (2021b).

To align with the reality of customer order scheduling, this study considers the classification of jobs from different orders into distinct classes. Every order consists of at least one job from each job class, necessitating a setup time whenever the machine transitions between classes. In the domain of single-machine scenarios with multiple objectives and various job classes, Liao (1993) introduced a branch and bound (B&B) algorithm to systematically investigate all viable solutions. Gupta, Ho & van der Veen (1997) were pioneers in investigating COSP with multiple classes and setup times, introducing a bi-criteria problem. They aimed to minimize both makespan and the total carrying costs of customer orders, with carrying costs computed based on the time difference between the completion times of the first and last job in each customer’s order. Lin, Yin & Liu (2013) addressed a sequence-dependent scheduling problem that incorporates order delivery, presenting a binary integer program and a dynamic programming algorithm. They also introduced tabu search, iterated local search, and genetic algorithm approaches for obtaining approximate solutions. In related studies on customer order scheduling with setup times, Erel & Ghosh (2007) examined customer orders with varying quantities of products from different product families processed on a continuously available machine in any sequence. Liu (2009), Liu (2010) and Hsu & Liu (2009) explored various scenarios involving multiple jobs within a job shop setting. Recently, Li et al. (2023) expanded upon the work laid down by Gupta, Ho & van der Veen (1997) by introducing an innovative approach that incorporates four heuristics, three locally enhanced search methods, and a branch and bound algorithm. Additionally, they applied a water wave optimality algorithm, offering four distinct wavelength variants. This collaborative effort is aimed at minimizing a mixed criterion including the order ranges, order tardiness, and total job completion times. In a parallel effort, Lin et al. (2023) presented four heuristics complemented by a local search method, four theoretical simulated annealing techniques, a cloudy theoretical simulated annealing hyperheuristic, and a B&B method to solve the problem. Their collective endeavors focused on minimizing a mixed criterion including the total order ranges and total order tardiness. Meanwhile, Gupta et al. (2023) introduced six two-phase heuristics, six variants of water-wave optimization algorithms, a mixed-integer linear programming formulation, and a B&B algorithm to solve their problem. This concerted effort is directed at producing approximate solutions that effectively minimize a mixed criterion of the cumulative holding costs for all orders and makespan for all jobs.

Setup times present a significant concern for manufacturing companies. Lengthy setup times can jeopardize the on-time delivery of customer orders, leading to potential losses in customer satisfaction and incurring direct or indirect costs. Effectively incorporating and managing the setup time factor in scheduling decisions can yield substantial benefits (Kim & Bobrowski, 1997; Allahverdi, Gupta & Aldowaisan, 1999; Allahverdi et al., 2008; van Donk & van Doorne, 2016; Zhao et al., 2018; Ying et al., 2023).

Motivated by the scarcity of research that address sequence-dependent bi-criteria scheduling problems that incorporate order deliveries, we introduce a novel problem in this domain. Our aim is to pinpoint a schedule that optimally balances the makespan and the weighted completion time of designated orders through a linear combination. Expanding on these insights, this study addresses a sequence-dependent bi-criteria scheduling challenge involving order deliveries. The primary aim is to devise an optimized schedule that minimizes a mixed criterion of the overall job makespan and the aggregated weighted completion time of all orders.

The main contributions of this study are summarized as follows: (a) we address a sequence-dependent bi-criteria scheduling challenge involving order deliveries; (b) we propose three lemmas and a lower bound in a branch-and-bound algorithm for finding an optimal solution; (c) we propose a simple heuristic according to the smallest first value of d_uv and their three improvements by three local searching methods and (d) we build three population-based genetic algorithms. The subsequent sections of this study are structured as follows. ‘Problem Statement’ outlines the proposed problem formulation. In ‘Heuristic methods and heuristic-based genetic algorithm’ , we present three dominant properties, a proposed lower bound, a heuristic, and its three improved schemes (i.e., pairwise interchange, extraction and forward-shifted reinsertion, and extraction and backward-shifted reinsertion). Additionally, we propose three variants of the genetic algorithm (GA) based on these three improved heuristics. ‘Tuning the Related Parameter of the Three Heuristic-Based Genetic Algorithms’ focuses on testing the appropriate parameter values in the GA. ‘Computational Simulations and Discussions’ provides a comprehensive report of all test results. The final section offers conclusions and outlines future research directions.

Problem statement

The formal framework of our proposed study is outlined as follows. Let us consider a set Ω = O₁, …, O_m comprising m orders designated for processing on a single machine. The machine cannot break down while performing all jobs. Each order O_i consists of n_i jobs, each assigned a weight w_i (w_i ∈ (0, 1)). The total number of jobs across all orders amounts to n, which defines the set N = {J₁, …, J_n}.

A critical assumption In this model (see examples in Lin, Yin & Liu (2013)) is that, irrespective of their respective orders, if job J_v immediately follows job J_u in the schedule, then there is a sequence-dependent setup time d_uv > 0 required. Additionally, we define d_0Ju to signify the machine status before the initial scheduling of job J_u, elucidating the setup required.

Notably, in this model, we assume that the processing times for all jobs in set N are negligible. Given a complete schedule σ encompassing all jobs, we introduce the notation C_Op(σ) to represent the completion time for order O_p, which is the time at which the last job of O_p is finished. Simultaneously, the makespan, denoted as C_max(σ), signifies the completion time of the last job in set N.

The primary objective of this study is to determine an optimal schedule that minimizes a linear combination of the sum of makespan and the weighted completion times of all m orders. In mathematical terms, we aim to find a solution to the following optimization problem:

Minimize $\mathrm{g}\left(\sigma \right)=\theta C_{max}\left(\sigma \right)+\left(1-\theta \right){\sum }_{p=1}^m{w}_p{C}_{Op}\left(\sigma \right).$

Here, θ represents a weighting parameter that allows for the adjustment of the trade-off between makespan and weighted completion times.

To simplify matters, we use the terms “node,” “job,” and “city” interchangeably to describe the fundamental elements of our scheduling problem. Additionally, we introduce two notations, $C_u\left(\sigma \right)$ and $C_v\left(\sigma \right)$ to represent the completion times of jobs u and v in schedule σ, respectively. $C_v\left({\sigma }^{\prime }\right)$ and $C_u\left({\sigma }^{\prime }\right)$ represent the completion times of jobs v and u in schedule σ′, respectively, where σ = (π, J_u, J_v, π^c) and ${\sigma }^{\prime }=\left(\pi ,J_v,J_u,{\pi }^c\right)$. Here, π and π^c are two distinct sequences within set N.

To enhance the efficiency of our B&B method, we establish two properties that assist in narrowing down the nodes that satisfy the conditions outlined in these properties. Before presenting the properties, we designate job J_L as the final job in sequence π.

Property 1: As jobs J_u ∈ O_p, J_v ∈ O_q, p ≠ q, if J_u is the last assigned job of O_p, and J_v is the last assigned job of O_q, $d_{Lu}+\frac{w_q}{w_p+w_q}{d}_{uv}<d_{Lv}+\frac{w_p}{w_p+w_q}{d}_{vu}$ and d_Lu + d_uv + max_x∈π^c{d_vx} < d_Lv + d_vu + min_x∈π^c{d_ux}, then σ dominates σ′.

Proof: In following, one will claim that (i) $C_{max}\left({\sigma }^{\prime }\right)-C_{max}\left(\sigma \right)>0$, and (ii) ${\sum }_{p=1}^m{w}_p{C}_{Op}\left({\sigma }^{\prime }\right)-{\sum }_{p=1}^m{w}_p{C}_{Op}\left(\sigma \right)>0$. Then, $\mathrm{g}\left({\sigma }^{\prime }\right)-\mathrm{g}\left(\sigma \right)>0$.

Let t_π be the completion time of the last job in π, and t_p1 (t_q1) be the completion time of the first job of O_p (O_q). Note that t_p1 < t_π (t_q1 < t_π). Moreover, for any possibly arranged job sequence of all jobs in π^c, let J_s be the first assigned job in π^c. Note that J_s ≠ J_u, J_s ≠ J_v. Then, one has the following results, (i) and (ii).

(i) $C_{max}\left({\sigma }^{\prime }\right)-C_{max}\left(\sigma \right)=\left[d_{Lv}+d_{vu}+d_{us}\right]-\left[d_{Lu}+d_{uv}+d_{vs}\right]>0.$ The last inequality follows from the given condition d_Lu + d_uv + max_x∈π^c{d_vx} < d_Lv + d_vu + min_x∈π^c{d_ux}.

(ii) For the orders O_p and O_q, $\left[w_q{C}_{Oq}\left({\sigma }^{\prime }\right)+w_p{C}_{Op}\left({\sigma }^{\prime }\right)\right]-\left[w_p{C}_{Op}\left(\sigma \right)+w_q{C}_{Oq}\left(\sigma \right)\right]=\left[w_q\left(t_{\pi }+d_{Lv}-t_{q1}\right)+w_p\left(t_{\pi }+d_{Lv}+d_{vu}-t_{p1}\right)\right]-\left[w_p\left(\left(t_{\pi }+d_{Lu}-t_{p1}\right)\right)+w_q\left(t_{\pi }+d_{Lu}+d_{uv}-t_{q1}\right)\right]=w_q\left[d_{Lv}-d_{Lu}-d_{uv}\right]+w_p\left[d_{Lv}+d_{vu}-d_{Lu}\right]=\left(w_q+w_p\right)\left[d_{Lv}-d_{Lu}+\frac{w_p}{w_q+w_p}{d}_{vu}-\frac{w_q}{w_q+w_p}{d}_{uv}\right]>0$.

For these orders (not including O_p and O_q), which are the last jobs of the order assigned in π^c, the sum of all $w_p{C}_{Op}\left({\sigma }^{\prime }\right)-w_p{C}_{Op}\left(\sigma \right)$ is ${\sum }_p^{}{w}_p{C}_{Op}\left({\sigma }^{\prime }\right)-{\sum }_p^{}{w}_p{C}_{Op}\left(\sigma \right)={\sum }_p^{}{w}_p\left[\left(d_{Lv}+d_{vu}+d_{us}\right)-\left(d_{Lu}+d_{uv}+d_{vs}\right)\right]>0$. Therefore, ${\sum }_{p=1}^m{w}_p{C}_{Op}\left({\sigma }^{\prime }\right)-{\sum }_{p=1}^m{w}_p{C}_{Op}\left(\sigma \right)>0$. Combining the results of (i) and (ii), $\mathrm{g}\left({\sigma }^{\prime }\right)>\mathrm{g}\left(\sigma \right)$ follows. □

Property 2: As jobs J_u ∈ O_p, J_v ∈ O_q, p ≠ q, if J_u is not the last assigned job of O_p, and J_v is also not the last assigned job of O_q, d_Lu + d_uv < d_Lv + d_vu and max_x∈π^cd_vx < min_y∈π^cd_uy, then σ dominates σ′.

Proof: It needs to be shown that $\mathrm{g}\left(\sigma \right)<\mathrm{g}\left({\sigma }^{\prime }\right)$.

For any possibly arranged job sequence of all jobs in π^c, one can be shown that $C_{max}\left(\sigma \right)<C_{max}\left({\sigma }^{\prime }\right)$, and ${\sum }_{p=1}^m{w}_p{C}_{Op}\left(\sigma \right)-{\sum }_{p=1}^m{w}_p{C}_{Op}\left({\sigma }^{\prime }\right)<0$.

Let J_s be the first assigned job in π^c. Note that J_s ≠ J_u, J_s ≠ J_v. Then, $C_{max}\left(\sigma \right)-C_{max}\left({\sigma }^{\prime }\right)=\left[d_{Lu}+d_{uv}+d_{vs}\right]-\left[d_{Lv}+d_{vu}+d_{us}\right]=\left[\left(d_{Lu}+d_{uv}\right)-\left(d_{Lv}+d_{vu}\right)\right]+\left(d_{vs}-d_{us}\right)$. Applying two given inequalities, one has $C_{max}\left(\sigma \right)-C_{max}\left({\sigma }^{\prime }\right)$<0. By the same argument, ${\sum }_{p=1}^m{w}_p{C}_{Op}\left(\sigma \right)-{\sum }_{p=1}^m{w}_p{C}_{Op}\left({\sigma }^{\prime }\right)={\sum }_p^{}{w}_p\left[\left(d_{Lu}+d_{uv}\right)-\left(d_{Lv}+d_{vu}\right)+\left(d_{vs}-d_{us}\right)\right]$<0 because 0 <w_p < 1. The last sum is over all orders on which the last job of the order is in π^c. Thus, $\mathrm{g}\left(\sigma \right)<\mathrm{g}\left({\sigma }^{\prime }\right).\square$

Property 3: Consider a scheduled σ = (PS, US), where PS and US denote the scheduled part with k jobs and the unscheduled part with (n-k) jobs, and PS∪US = N. If $d_{j_{\left[k-1\right]}{J}_{\left[k\right]}}>d_{J_{\left[k-1\right]}{J}_v}$for all J_v ∈ US, and J_v and $J_{\left[k\right]}$ belong to the same order, then σ = (PS, US) can be eliminated.

In the ensuing discussion, our objective is to establish a lower bound for the schedule denoted as σ = (π, π^c), where π represents the scheduled portion comprising k jobs, and π^c signifies the set of n-k unscheduled jobs. Let t_π denote the completion time of the final job, J_L, in π, and ∑_i∈πw_iC_Oi stand for the aggregate weighted completion time of the scheduled orders in π.

In light of the sequence-dependent aspect, we incorporate J_LL into the set π^c, resulting in the creation of a new set denoted as π^c∗ = π^c∪{J_L}, now comprising (n − k + 1) jobs. To approximate the completion times in π^c, we embark on the following procedure to extract the (n − k) smallest sequence-dependent setup times within π^c∗. These recorded values are denoted as D(1), D(2), …, D(n − k).

Notably, the values of D(1), …, D(n − k) exhibit a non-decreasing trend. Consequently, we can estimate the lower bound for each $C_{\left[k+i\right]}\left(\sigma \right)$, with i ranging from 1 to n − k, by employing the following formula: (1)${\displaystyle C_{\left[k+i\right]}\left(\sigma \right)=t_{\pi }+\sum _{v=1}^i{D}_v,i=1,\dots ,n-k,}$

Additionally, we maintain a record of the frequency number f_i corresponding to each order in π^c, where i spans from 1 to q, with 1 ≤ f₁ ≤ f₂ ≤ … ≤ f_q and q<m. Utilizing Eq. (1), we derive completion times for the outstanding orders in π^c as $C_{\left[k+f_1\right]}\left(\sigma \right)$, $C_{\left[k+f_1+f_2\right]}\left(\sigma \right)$, …, $C_{\left[k+f_1+f_2+\dots +f_q\right]}\left(\sigma \right)$.

Furthermore, we arrange the weights $w_{\left(1\right)}\le w_{\left(2\right)}\le \dots \le w_{\left(q\right)}$ in non-decreasing order within the set {w_v ∈ π^c, v =1 , …, q}. Applying the lemma established by Hardy, Littlewood & Polya (1967), we arrive at a lower bound for the estimated weighted completion times of the unscheduled q orders, represented as ${\sum }_{i=1}^q{w}_{\left(q-i+1\right)}{C}_{\left[k+{\sum }_{v=1}^i{f}_v\right]}\left(\sigma \right)$. Hence, the lower bound is derived as follows:

$\theta C_{max}\left(\sigma \right)+\left(1-\theta \right){\sum }_{p=1}^m{w}_p{C}_p\left(\sigma \right)\ge \theta \left(t_{\pi }+{\sum }_{v=1}^{n-k}{D}_v\right)+\left(1-\theta \right)\left({\sum }_{v=1}^b{w}_v{C}_{\left[v\right]}\left(\sigma \right)+{\sum }_{i=1}^q{w}_{\left(q-i+1\right)}{C}_{\left[k+{\sum }_{v=1}^i{f}_v\right]}\left(\sigma \right)\right)$, where b +q =m.

Tuning the related parameter of three heuristic-based Genetic Algorithms

In this section, we conduct several experimental tests to explore the values of related parameters of three heuristic-based genetic algorithms (GAs). To determine the parameters used in three GAs as chromosome count nsize, evolution generations gsize and mutation rate p, we consider an experimental design as the number of jobs at n =12, the number of orders at m =2, and the weight of objective function at θ = 0.25. One hundred instances are tested for this case. The initial setup times d_0Ju and sequence-dependent setup time d_uv are generated from uniformly distributed U[1, 20], while the weights of orders w_i are generated from uniformly distributed U[1, 10]. The optimal solutions (Opt_i) and the near-optimal solutions (GAH_i) are produced by performing the B&B method and three heuristic-based GAs (DPIGA, DFORGA, or DBKGA). The measurement criterion is based on the average error percentage (AEP) or AEP = $\left\{{\sum }_{i=1}^{100}\left(GA{H}_i-Op{t}_i\right)/Op{t}_i\right\}/100$. The parameter tuning process is based on a single-factor experimental method applied to a set of generated problem instances.

To evaluate the value of nsize, gsize was fixed initially at 50 and p was fixed at 0.02. We executed DPIGA, DFORGA, and DBKGA several times and changed the values of nsize over the interval [10, 50] with an increment of four. As shown in Fig. 1 (top), as the nsize value increases, the trends of the AEPs rapidly decrease and become flat at point 50 obtained from DBKGA. Therefore, in subsequent simulation experiments, the nsize is adopted at 50.

With the previously determined nsize at 50 and fixed mutation rate at 0.02, we performed three GAs to determine the values of gsize over [10, 200] with an increment of 10. Figure 1 (middle) shows that the AEP of DBKGA has a minimum point at 180. Thus, gsize is set at 180 for all three proposed GAs in later experiments.

After investigating the parameters nsize and gsize, we established nsize =50 and gsize =180. We conducted three GAs to assess the mutation rate p, ranging from 0.01 to 0.2 with a consistent increment of 0.01. As depicted in Fig. 1 (bottom), AEP demonstrates a gradual convergence, while DBKGA attains the minimal error at a mutation rate of 0.2. Accordingly, we adopt a mutation rate of 0.2 as the benchmark for all three GAs in the subsequent stages.

In a similar way, we conducted another one hundred random instances at n =240 to explore the values of three parameters of nsize, gsize, and p for three GAs (DPIGA, DFORGA, and DBKGA). Based on previous studies, to capture parameters that outperform three simple heuristic algorithms (DPI, DFOR, and DBK), we only performed DFORGA to change the values of one parameter but fixed other parameters each time and compared the objective function of DFORGA with those from DPI, DFOR, and DBK. In this case, we performed three simple heuristic algorithms (DPI, DFOR, and DBK) and DFORGA to solve the data set and computed each average of their objective functions. The final results are summarized and reported in Table 1. It can be observed in Table 1 that nsize =100, gsize =1500, and p =0.2 are appropriate values for using DPIGA, DFORGA, and DBKGA.

Table 1:

The summary of changing related parameters in DFORGA.

	DPI	DFOR	DBK	DFORGA	DFORGA	DFORGA	DFORGA	DFORGA
(nsize, gsize, p)				(100, 1500, 0.2)	(150, 1500, 0.2)	(200, 1500, 0.2)	(250, 1500, 0.2)	(300, 1500, 0.2)
mean	7076.30	7144.24	7072.63	7069.48	7076.26	7101.08	7125.53	7128.32
(nsize, gsize, p)				(100,1200, 0.2)	(100, 1300, 0.2)	(100, 1400, 0.2)	(100, 1500, 0.2)	(100, 1600, 0.2)
mean	7076.30	7144.24	7072.63	7079.52	7077.23	7072.63	7069.48	7071.43
(nsize, gsize, p)				(100,1500, 0.16)	(100,1500, 0.17)	(100, 1500, 0.18)	(100, 1500, 0.19)	(100, 1500, 0.2)
mean	7076.30	7144.24	7072.63	7085.795	7072.375	7072.69	7084.76	7069.48

DOI: 10.7717/peerjcs.1763/table-1

Based on the above tests, the population size (nsize), generation size (gsize), and mutation rate (p) are set at (50, 180, 0.2) for the small n case and at (100, 1500, 0.2) for big n case in all three GAs (DPIGA, DFORGA, and DBKGA).

Computational simulations and discussions

This section conducts some tested designs to assess the performances of the proposed B&B, three heuristic methods (DPI, DFOR, and DBK), and three heuristic-based GAs (DPIGA, DFORGA, and DBKGA). The solution methodologies for the proposed problem are summarized as follows.

All the related parameters are designed as follows. The initial setup times d_0Ju and sequence-dependent setup time d_uv are generated from uniformly distributed U[1, 20], and the weights of orders w_i are generated from uniformly distributed U[1, 10]. The experimental designs include the weight of the objective function at θ = 0.25, 0.5, and 0.75, the number of jobs at n =12, 16, 20, and 24 for the small-job case, and n =60, 120, 180 and 240 for the large-number-of-job case. It should be emphasized that n =m*nc signifies that the total of n jobs can be apportioned among m orders, with each order containing nc jobs. The details are summarized in Table 2. The experiment investigates 24 cases for both small and large job numbers. In each case, 100 instances are generated to assess the performance of all proposed methods. In total, 2,400 instance problems are examined for both small and large job sets. The findings are summarized in Tables 3, 4, and 5.

When dealing with a case with a small number of jobs, we assess the effectiveness of the three heuristic algorithms (DPI, DFOR, and DBK) and the DPIGA, DFORGA, and DBKGA algorithms by examining their average error percentage and maximum error across various values of n, m, and θ. The AEP is computed as follows: AEP = $\left\{{\sum }_1^{100}\left(Hso{l}_i-Op{t}_i\right)/Op{t}_i\right\}/100$, where Opt_i is derived from the B&B method, while Hsol_i isproduced using the six proposed approximate techniques. Additionally, we furnish the mean and maximum values of the explored nodes, along with their respective CPU times (measured in seconds), to evaluate the performance of B&B for each scenario. It is worth noting that an instance problem is considered feasible if the B&B process requires fewer than 10⁸ nodes.

With regard to the performance of B&B, columns 3 and 5 of Table 3 show that B&B consumes more nodes or CPU time at smaller order numbers (m =2) compared to larger order numbers (m =3 or 4). This implies that the properties and lower bounds become more useful as the number of orders increases. The effect of θ in columns 3 and 5 of Table 3 does not show a clear trend. Additionally, the number of nodes increases significantly as the number of jobs (n) increases. This phenomenon aligns with the characteristics of the NP-hard problem.

In regard to the obtained results for the three heuristics (DPI, DFOR, and DBK) and three genetic algorithms (DPIGA, DFORGA, and DBKGA) applied to cases with small n, Table 4 and Fig. 2 present their respective average error percentages. To summarize, the mean AEP values are 0.4733, 0.5246, and 0.4551 for DPI, DFOR, and DBK and 0.3412, 0.3504, and 0.3358 for DPIGA, DFORGA, and DBKGA, respectively. The DBKGA algorithm performed better than all other approaches. Table 4 provides the mean AEP values corresponding to different values of parameters θ, m, and n. Upon inspection of Table 4, it is evident that the AEPs for the three heuristics and three GAs experience a marginal rise as n increases from 12 to 24. In contrast, these AEPs exhibit a comparatively stable trend when subjected to changes in the values of m or θ.

According to the criterion of average relative percentage deviation (RPD), RPD is calculated as follows: RPD = $\left\{{\sum }_1^{100}\left(Hso{l}_i-H_{min}\right)/H_{min}\right\}/100,$ where Hsol_i represents the output obtained using the six proposed approximate methods, and H_min is the minimum value among these six solutions. Table 5 and Fig. 3 present an evaluation of the performance of the heuristics and genetic algorithms in the context of larger n values. The data in Table 5 reveal a consistent trend: regardless of the value of n, DBKGA consistently achieves the lowest average RPD among the three heuristics (DPI, DFOR, and DBK) and the two GAs (DPIGA, DFORGA). Specifically, the average RPD values are 0.0383, 0.0548, and 0.0341 for DPI, DFOR, and DBK, respectively, and 0.0136, 0.0201, and 0.0099 for DPIGA, DFORGA, and DBKGA, respectively.

Additionally, box plots are presented in Fig. 3 to visually depict the mean RPD for the heuristics and GAs. These plots clearly show that the overall range and interquartile ranges of the mean RPDs are broader for DPI, DFOR, and DBK than for DPIGA, DFORGA, and DBKGA.As indicated in Table 5, as n increases from 60 to 240, the average RPD values exhibit a declining trend.

Regarding computational times, Fig. 4 displays violin plots that visualize the distribution of average CPU time for both the three heuristic algorithms and the three GAs under the conditions of large n. These plots clearly indicate that the overall range and interquartile ranges of average CPU times for DPIGA, DFORGA, and DBKGA are broader compared to DPI, DFOR, and DBK. Specifically, the average CPU time values are 0.0196, 0.0195, and 0.0197 for DPI, DFOR, and DBK, respectively, while they are 2.5596, 2.6028, and 2.5983 for DPIGA, DFORGA, and DBKGA, respectively. We do not provide CPU time reports for these six methods in the case of small n because they are less than one second of CPU time.

Below, we delve into the statistical significance of distinctions among the three heuristic approaches and the three GA algorithms. Initially, we applied an analysis of variance (ANOVA) technique utilizing a linear model on the AEP data for small ‘n’ or RPD data for large ‘n’ within SAS 9.4. The ANOVA, shown in Table 6 and 7, revealed that the model was significant for both small ‘n’ and large ‘n’. Additionally, the algorithm and job-size (n) factors exhibited significance for both small ‘n’ with an F test value of 456.67 and five degrees of freedom and for RPD with an F test value of 50.13 and five degrees of freedom. It is also noticeable that the factors of weight (θ) for the objective function and number of orders (m) are significant in Table 6 for small ‘n’, whereas they did not demonstrate significance for the large ‘n’ case shown in Table 7. The rationale behind this observation is that the effect of ‘n’ outweighs the values of weight (θ) for the objective function or the number of orders (m) in the context of the large ‘n’ case.

The normality assumption of the linear model was validated for both the AEP and RPD data. This was evidenced by the Kolmogorov–Smirnov normality test, which produced p values exceeding 0.05 (with statistic values of 0.053336 for small ‘n’ and 0.058442 for large ‘n’). Additionally, the AEP and RPD Tukey grouping for means of algorithms were also tested. Figure 5 shows that the algorithms can be ranked into five groups. The five groups range from the worst to the best as DFOR, DPI, DBK, DFORGA, and (DPIGA, DBKGA) for the AEPs, shown in Fig. 5 (left) as well as the four groups from the worst to the best as DFOR, (DPI, DBK), (DFORGA, DPIGA), and (DPIGA, DBKGA) for the RPDs shown in Fig. 5 (right).

The three heuristics (DPI, DFOR, and DBK) rely on predefined rules or strategies to find solutions quickly. While effective in some cases, they may get stuck in local optima and fail to explore the entire solution space. On the other hand, the three GAs (DPIGA, DFORGA, and DBKGA) employ population-based searches and a genetic-inspired selection, crossover, and mutation process, which allows them to explore a broader range of the solution space, potentially discovering better solutions. The advantages of genetic algorithms over heuristics in setup time optimization often stem from their ability to handle complex relationships, explore diverse solution spaces, and converge toward globally optimal solutions by using evolutionary mechanisms. Overall, it can be confirmed that the DBKGA is the best approach with the smallest value of AEP or RPD among all six proposed algorithms.

Conclusions and suggestions

In practical settings, a machine incurs a setup time whenever it switches from processing one product to another. Scheduling jobs while taking into account setup times and associated costs has attracted significant attention in both the manufacturing and service sectors, driving extensive research endeavors. While prior studies on customer order scheduling often focused on jobs to be processed across multiple machines, they frequently overlooked the critical factor of setup time. This study addresses a sequence-dependent bi-criterion scheduling problem, integrating considerations for order delivery. The primary aim is to pinpoint an optimal solution that minimizes the combined makespan and weighted sum of completion times for all specified orders.

To tackle this formidable challenge, we propose a branch-and-bound method incorporating two dominance properties and a lower bound. This approach is finely tuned to yield optimal solutions in scenarios involving a limited number of jobs, exemplified by cases with n =12, 16, 20, and 24. Demonstrating notable efficacy, this method excels for up to n =24 jobs, accommodating diverse combinations of order quantities and objective weights. For both small n (=12, 16, 20, and 24) and large n (= 60, 120, 180, and 240) tested instances, we introduce a tailored heuristic approach. This heuristic is complemented by three strategic techniques: pairwise interchange, extraction and forward-shifted reinsertion, and extraction and backward-shifted reinsertion (DPI, DFOR, and DBK). Additionally, we employ three heuristic-based genetic algorithms (DPIGA, DFORGA, and DBKGA) to obtain reliable approximate solutions.

Notably, it is evident that DBKGA performed significantly better than the five other proposed heuristics for both small and large test cases. A potential direction for future research could involve expanding the model to a sequence-dependent flow shop setting that incorporates order deliveries. Furthermore, it may be worthwhile to investigate extending this framework to minimize the total tardiness of all fulfilled orders. As another future topic, we can also consider that the machine can break down or extend this model to the multiple-machine setting.

Supplemental Information