1 Introduction

Studying parallel machines are of high importance both theoretically and practically. In terms of theory, it is an extension of the single machine and a particular instance of the flexible flow shop. In terms of practice, it is vitally important for the sake of common parallel resources in the real world, but very difficult to find an optimum solution for the problem. Likewise, decomposition methods used in multi-stage systems sometimes utilize parallel machine techniques.

Sometimes the parallel machine scheduling (PMS) is considered as a two-step procedure. First, assigning jobs to available machines. Second, when allocating jobs to machines has been done, the sequence of the jobs needs to be allocated to each machine [1]. Since we cannot consistently find an absolute solution for parallel machines and for most of the criteria, especially those based on tardiness - that don’t have a linear relationship with completion time - the problem is non-polynomial (NP). Based on the division problem, Lenstra et al. [2] proved that a parallel machine problem of minimizing the total tardiness is a binary non-polynomial problem even for two machines. Koulamas [3], based on the single machine model in Du and Leung [4], indicated that the parallel machine problem is at least a binary non-polynomial. About the minimization of weighted tardiness, Lawler et al. [5] proved that when jobs have different weights, the problem would be NP-hard.

One of the deterministic methods for solving such a category of problems is dynamic programming [6], although the extent of the application of this method is very limited due to the diversity and multiplicity of options whether in job allocation or in determining their sequence on each machine. Branch and bound approach has also been used in many studies to solve parallel machine problems, for example, Elmaghraby and Park’s study [7] for jobs with identical processing time and due date, Barnes and Brennans’s study [8] for mostly 20 jobs and 4 machines, and Yalaoui and Chu’s study [9] for minimizing total tardiness.

Moreover, some papers investigate solutions to parallel machine problems in particular situations, for example, Lawler’s study [10] with respect to transportation model and jobs with an identical processing time, Root’s study [11] with identical due date and Pritsker et al’s study [12] in which binary linear programming is used. Ranjbar et al. [13] also used a branch and bound approach for PMS with stochastic processing time. Figielska [14] proposed a problem of scheduling pre-emptive jobs on parallel unrelated machines in the presence of renewable resource constraints and sequence-dependent setup costs. The problem was solved using the column generation method and ant colony optimization (ACO) algorithm. Chen [15] developed a column generation based branch and bound method to solve simultaneous job scheduling and resource allocation problems.

Reviewing parallel machine flexible resource scheduling (PMFRS) problems and an unspecified version of it (UPMFRS), Edis and Oguz [16] gave some extensions to the model of dynamic PMFRS and proposed integer programming (IP) models for static and dynamic UPMFRS problems with the aim of minimizing the makespan. Lee et al. [17] proposed a uniform parallel machine problem in which the objective is to jointly find an optimal assignment of operators to machines and an optimal schedule to minimize the makespan. The problem was solved using a genetic algorithm (GA) and particle swarm optimization (PSO).

Due to the complexity of PMS problems and the size of real-world problems, deterministic methods cannot be effective. Thus, heuristic algorithms are often used. These methods are often to determine a sequence based on a priority criterion and then to assign the jobs to machines based on the time of availability and also they are single-objective. Various methods are proposed in this field, among which montage’s ratio method [18], cost over time [19], apparent urgency [20], traffic congestion [21], DS [22], hybrid simulated annealing [3], PDEC [9], high priority job first [23] and modified due date [24] can be indicated. Besides, in terms of the multi-objective algorithms we can count recently introduced algorithms by [25] in the field of inventory control, and [26, 27] in the field of transportation science.

Correa et al. [28] suggested a strong model for unrelated PMS. Tao [29] proposed an online algorithm for unrelated PMS with a minimization of weighted total completion time. Ouazenea et al. [30] proposed a mixed-integer mathematical model for Job loading balance on parallel machines. Also, Hsu et al. [31] studied rate-modifying operations with the aim of minimization of the total completion time on unrelated PMS.

Regarding the complexity of real-world problems, it is necessity to study scheduling problems with multi-objectives. Some scholars take the advantages of recently multi-objective algorithms extended to make contributions for solving scheduling problems to reach better solutions. Therefore, in this study, we extend a multi-objective version of the VPL algorithm using crowding distance in comparison with the original ones which is archive based.

In this study, we concentrate on identical PMS considering splitting jobs together with sequence-dependent setup waste and time. The remaining structure of the paper is shown as follows. Section 2 gives a brief literature review on PMS and then determines a theoretical gap in this research area. Section 3 is devoted to the mathematical modeling of the problem. Section 4 comes to extensively explain multi-objective Pareto-based meta-heuristic algorithms and solution representation of our approach. Sections 5 and 6 introduce our proposed algorithm. Section 7 investigates the performance of our proposed algorithm using different metrics and a set of numerical examples. Finally, Section 8 will draw the conclusion and future researches.

2 Literature review

The allocation of assets to a certain number of jobs in a predetermined time period defines scheduling while optimizing one or more objectives [1]. The scientific approach to operation planning originated from the industrial revolution and Henry Laurence Gantt’s studies [32], and since then it has been widely studied. In this section, a review of PMS and the approached environment in this paper (i.e. machine scheduling with spilling jobs) is outlined. PMS problems are among the most complex machine scheduling problems. Based on the operating velocity of each machine for various jobs, there are three parallel machine scheduling (PMS) environments according to [33]: identical machines run at the same speed; heterogeneous machines run at various velocities but its operating rate/velocity is compatible for every machine operating various jobs. In the sense of unrelated machines, defining a generalization of a heterogeneous environment as an unrelated machine collection can consist of a non-identical machine collection. Some studies have been done on PMS with and without splitting jobs as follows.

Serafini [34] studied a real problem of scheduling looms in the textile industry. Two environments of uniform machines and unrelated machines have been considered in this study, while the first one was modeled as network flow and the latter one was modeled using linear programming. The objective is defined as minimizing maximum tardiness and weighted tardiness. Xing and Zhang [35] considered the PMS problem with splitting job property and independent times of setup. Some simple polynomial solvable cases and a heuristic maximum completion time estimation procedure and set-up-time list scheduling (ML) including worst-case analysis with the aim of minimization makespan were presented as well.

Sarıçiçek and Çelik [36] presented a mathematical model using mixed-integer variables for a parallel machine scheduling problem with job split, with the aim of minimizing total tardiness. They applied two metaheuristics (i.e. Tabu search and simulated annealing) to solve three groups of test problems. Results suggest that the performance of the simulated annealing algorithm is far better than the tabu search algorithm in terms of optimal solution deviation and computational time. Park et al. [37] addressed the same problem considering sequence-dependent major/minor times of setting up, with the aim of minimizing total tardiness. They developed three heuristics: the slack-based heuristic, the dynamic scheduling window-based heuristic, and the estimated latest starting time-based heuristic. The results illustrate that the former and the latter heuristics have the same performance and they have supremacy in terms of performance than the dynamic scheduling window-based heuristic.

Shim and Kim [38] proposed a branch and bound algorithm to schedule jobs composed of some unit jobs considering independent setup times on parallel machines. The objective has been stated to be the minimization of total tardiness. Kim et al. [39] also considered the same objective. They proposed a heuristic algorithm which is consists of two-phase to find a solution for the problem. In phase one, an initial sequence is suggested for the PMS problem. In phase two, the initial scheduling plan is rescheduled on machines considering the jobs and sub-jobs.

Yalaoui and Chu [40] studied the scheduling problem of splitting jobs on parallel machines considering sequence-dependent setup times with the aim of minimization of the maximum makespan. They provided a heuristic method to solve the problem. Firstly the problem is simplified into a subject called single machine scheduling and is then metamorphosed into a well-known traveling salesman problem and solved utilizing Little’s method. Using the results of the first part, the initial solution is improved in a stepwise manner taking into account the setup times and job splitting property. Zhu and Heady [41] used a mixed-integer programming approach for minimizing job earliness and tardiness in unrelated PMS problem with sequence-dependent setup time. They stated that the proposed model is beneficial to researchers for testing the performance of earliness and tardiness heuristics.

The unrelated PMS problem with the aim of minimization total tardiness has also been studied by Shim and Kim [42]. Using several dominance properties and lower bounds, they suggested a branch and bound algorithm to solve the problem. Logendran and Subur [43] Considered the PMS problem with the aim of minimization of total weighted tardy jobs in the field of dynamic job releases and machine availability. They proposed a mixed-integer linear programming model containing constraints including hard operations to guarantee that Just In Time (JIT) production system is running. Four methods are implemented to generate initial solutions that are inputs for performing a Tabu search algorithm for reaching an optimal solution. Logendran et al. [44] investigated research with the same method but without splitting job property and sequence dependent setup times. Six various Tabu search-based algorithms and four various primary solutions finding mechanisms have been extended for generating a finest scheduling plan that minimizes weighted tardiness. Based on computational results, a fixed Tabu list size and a short-term memory algorithm were suggested for small instances whereas, for medium and large instances minimum iteration and long-term memory algorithm were used.

Liaw et al. [45] addressed the issue of unrelated PMS considering the minimization of the total weighted tardiness as an objective function. A branch and bound algorithm is proposed to find a solution using supremacy rules to reduce hopeless solutions. The results show that the proposed algorithm is appropriate for up to 18 jobs and 4 machines. Wang et al. [46] investigated the unrelated PMS problem considering splitting job property with the objective of minimization of makespan. Differential evolution is applied as the solution method and a new method for crossover and mutation are presented for global search procedure according to the job splitting constraints. Experimental results show the efficiency and feasibility of the proposed hybrid differential evolution.

Chen and Wu [47] studied the unrelated PMS problem with auxiliary constraints and total tardiness minimization objective. In this study it is assumed that each job requires one operation and has a due date, and to switch from one type of job to another type, setup time for die is needed. To solve this problem they proposed a heuristic-based threshold-accepting methods, Tabu lists, and improvement procedures. The experiment results show that the proposed algorithm can obtain optimal solutions for small size problems and performs much better than SA and apparent-tardiness cost-with-setup for problems of larger sizes.

Vallada and Ruiz [48] developed a GA with a new crossover operator and a very fast local search for unrelated PMS problems with sequence-dependent setup times. They suggested a MIP model with the objective of minimization of makespan. The experiment results show that the proposed method significantly outperforms the best methods known from the literature. An unrelated PMS problem with sequence and machine-dependent setup times and weighted jobs was studied by [49]. They applied a branch and bound algorithm in which the upper bound is obtained from the solution provided by GRASP. Two MIP models are solved using CPLEX and the results are compared with the proposed algorithm. Outcomes prove the superiority of the performance of the algorithm on instances with about 30 jobs. Lin et al. [50] compared the performance of different heuristic and metaheuristic methods for a variety of unrelated PMS environments with three performance measures. The objective function is minimization of makespan, total weighted completion time and total weighted tardiness, and the test is made using least significant difference method.

Fanjul-Peyro and Ruiz [51] studied two situations of unrelated PMS problem under makespan minimization: one situation in which not all available parallel machines should be utilized, and another situation with no obligation to do all the jobs. MIP models are proposed for both situations and it is shown that the latter situation could be solved with commercial solvers efficiently. The PMS with splitting jobs and learning effects to minimize total completion time in a labor-intensive industry was presented by Wang et al. al. [52], solved by branch and bound algorithm and greedy search in small and large scale, respectively. A Tabu-enhanced iterated Pareto greedy algorithm was proposed by Lin et al. [53] to solve the unrelated PMS of completion times. Lara et al. [54] studied identical PMS problems with a different problem. Correa et al. [28] present a scheduling problem with splitting jobs which each part needs to be set up to minimize the weighted sum release date to minimize total tardiness. In this study preemption and splitting are not allowed. Shahvari and Logendran [55] developed a Tabu search heuristic for unrelated PMS problem to minimize the total weighted tardiness and total weighted completion time.

Yin et al. [56] investigated the PMS of deteriorating jobs in a disruptive environment where some of the machines are out of reach at a specific time. To minimize the total completion time, they applied different algorithms in polynomial time. Mensendiek et al. [57] studied identical PMS problems with a fixed due date which is prevalent in the industry. This problem is solved by two exact and approximate methods, branch and bound and Tabu search and genetic algorithms; respectively. Lee and Kim [58] studied two identical PMS problems in the case of machines are not uninterrupted available to minimize total tardiness and solved it by branch and bound algorithm.

Mora and Mosheiov [59] studied scheduling problems of batch processing machines with controllable processing times both on single and identical PMS. In scheduling problems with controllable processing times, the processing times are controlled by additional resources. The aim of this research is to study two exemplars of the presented model. They minimized the total flow time and the compression cost, and in the next step minimized the total flowtime subject to an upper bound on the maximum compression. Lia et al. [60] addressed the unrelated parallel batch processing machine problems. In this study, jobs are not identical and each machine is able to process several jobs as a batch. The size of the jobs does not exceed the capacity of each machine. The objective of the problem is the minimization of the makespan. The model is solved by two groups of heuristics. In the first group, at the outset, jobs scheduled in each batch, and then each batch is assigned to each machine. By contrast, in the second group, jobs are assigned to each machine firstly and then scheduling is done in each batch.

Zhoua et al. [61] studied the parallel machine subject with the aim of minimization of the makespan of arbitrary job size. In this model, machines have different capacity and processing speed. They presented an effective differential evolution-based hybrid algorithm for the large scale problems, in which batches are generated and assigned to machines. Li et al. [62] studied the integrated production scheduling and delivery on identical parallel batch processing machines with the aim of maximization of total profit company. Multiple jobs can be loaded on each machine as a batch provided that the size of the jobs does not exceed to machine capacity. The company may earn a profit if the jobs are done on due dates. They demonstrated that, in the event that the sizes of the jobs are identical, the problem could be solved in polynomial time. Otherwise, it is known as NP-hard problem. Thus, they proposed and compared heuristics to solve NP-hard cases.

Identical PMS is studied by [63] in the field of JIT (Just-In-Time) environment in the presence of WIP (Work In Process) and preemptive jobs. Two objectives are considered in this paper: the former is the minimization of total weighted earliness and tardiness, holding cost of all jobs which are waiting to be processed as WIP costs, and the latter is related to the number of interrupted jobs. Moreover, this study utilizes two multi-objective algorithms known as the NSGAII and the NRGA.

In a research done by [64], the unrelated PMS with sequence-dependent setup time and splitting jobs is combined with the heterogeneous vehicle routing problem to tackle a real case study in the field of the metal packaging industry. The objectives of the problems are minimization of setup costs and distribution costs. Related to solving the problem they introduced a two-stage iterative heuristic. Two mathematical models are represented by [65] for identical PMS with specified operating times and tool requirements with the objective of makespan minimization. The problem is solved by an adaptive large neighborhood search metaheuristic they introduced. A memetic differential evolutionary algorithm is represented in [66] for minimizing energy-efficient and the makespan in the field of PMS.

A study done by [67] a heuristic called SLMR is presented for solving the scheduling problems of splitting jobs on parallel machines with learning effects and the vital-few law with the objective of minimization of makespan. The uniform PMS problem with splitting job and setup resources in the field of dedicated machines was studied by [68] in which the setup times are sequence-independent with the minimization of makespan. To solve the problem they introduced a heuristic which is enough good for a real case study in the field of manufacturing of automotive pistons in Korea. The uniform PMS problem considering the green approach to minimize emissions of pollutions with the objective of minimization of makespan simultaneously was studied by [69]. To solve the problem they improved H1 and H2 heuristics. A memetic differential evolution algorithm was developed by [66] for solving unrelated PMS problems with the aim of minimization makespan and total energy consumption consequently.

The PMS problem considering job splitting and hazardous wastes in the field of fuzzy environment is formulated for the injection molding industry by [70] with the aim of minimization of the total tardiness and the total hazardous waste. For reaching a good solution for this industry, a genetic algorithm is implemented.

Identical parallel batch processing machine with the two objectives of minimizing makespan and maximum tardiness has been studied [71]. To find out the solution, a multi-objective ant colony optimization approach called Pareto-based ant colony system (PACS) is used. Wang and Leung [72] addressed parallel batch processing with identical processing time jobs on machines with different capacities to minimize the makespan. In this model, the size of the jobs does not exceed the capacity of each machine. Finally, a polynomial-time approximation algorithm was proposed and showed to be a very good algorithm in practice. Lopes and de Carvalho [73] considered the unrelated PMS problem with setup times associated with the sequence and availability dates for machines and release date for jobs with the aim of minimization of whole weighted tardiness. They developed a branch and price method to solve the problem and a new column generation acceleration method was developed in this paper. Results indicate the reasonable computational time for large size problems. In addition to above studies, it is worth mentioning that some excellent algorithms have been presented for solving combinatorial optimization problems such as [74,75,76,77], and some of them are specialized in the field of scheduling, for instance [78,79,80,81], and which all of them give us valuable insight to extend our algorithm.

The literature we have reviewed are selected from globally well-known scholarly data centers, for instance, Elsevier, Springer, etc. We search for articles published between 1996 and 2020. Since we aim to investigate recent researches in the field of PMS, we will be focused on published papers between the years 2010 and 2020. We concentrate on two aspects in the seeking procedure: mathematical modeling and solving methodologies. In terms of mathematical modeling, the type of jobs is considered (splitting or not) and the makespan, total weighted completion time, lateness, total weighted tardiness, total number of tardy jobs, maximum tardiness, maximum weighted tardiness, total tardiness, maximum makespan, earliness and tardiness, total completion time, the weighted sum of completion times, total flow time and the completion cost, total profit, and total production costs form the objective functions. In terms of solving methodologies, a variety of techniques for instance meta-heuristics methods (e.g. GA, PSO, TS, SA) and exact methods (e.g. constraint programming, branch and bond) have been used. Due to the aforementioned review, the research gap is identified as follows: first and foremost, the objective of minimization of setup wastes together with minimization of total tardiness in the form of a bi-objective problem is not studied so far. Furthermore, another restriction that increases the complexity of the problem is the limited number of molds for each job rather than a machine; in which we consider it as a constraint in our model. It is worth noting that another main contribution is presenting a Multi-Objective version of the Volleyball Premier League (MOVPL) algorithm using crowding distance borrowed from NSGA-II while the original MOVPL is archive based. Consequently, this modification of MOVPL makes it a more proper and accurate algorithm for solving identical unrelated PMS with divisible jobs with the aim of reducing total tardiness and hazardous waste simultaneously. This supremacy will prove in the following sections.

3 Problem description and mathematical model

In the plastic injection industries, the granule shaped raw plastic materials need to be melted and subsequently injected with pressure to the specified molds to making goods. In the continuing, plastic goods have been chilled to become hard, then the mold will open to deliver the goods. This procedure is performed successively. Each one of the goods needs a specific mold to be produced, and each mold could be set up on a number of specified machines. For those goods with more than one mold, we can take advantage of this opportunity by setting up molds on a number of machines concurrently to reduce the total tardiness and the total amount of wastes per machine.

One of the most common situations in the plastic injection molding industry that is an extension of the classical parallel machines model is investigated in this study. A sequence dependent wastes and mold resources constraints are added to the PMS model with splitting jobs, in addition to the sequence dependent setup times. These modifications boost the complexity of the problem. Therefore, this complexity need to be considered. Firstly, when setting up a job on a machine, it is vitally important which kind of job has been done before, since the color and material type have direct effect on the waste produced. For instance, suppose we have a sequence of two jobs on a machine where the first one has black color and the second one has white color. When a product with black color has done, the cylinder is full of black color and we need to waste material to clean the cylinder for setting up the white color product. Therefore, this process takes time and waste. Secondly, mold is the important source of manufacturing in addition to machines. In the molding industry, commonly we have more than one mold for a specific product. We can take the advantages of this resource for reducing the completion time of each job by allocating one specified job to more than one machine at same time considering the number of mold resources of the job. So, the mold resource constraint plays a special role in the molding industry as a machine and it makes the PMS problem more complex. Therefore, we need to consider these situations to generate a good sequence of jobs on each machine.

This research suggests an MIP model for identical PMS problem. Jobs can be split on machines, but they have dependent waste and it is worth noting that there are precedence constraints amongst jobs in terms of their material and color. To clear the waste dependent problem on parallel machines we assume 4 jobs with one machine. Let job 1 be the first one to process, then all remaining jobs could be processed after job one but with different amount of setup waste which we call waste dependent. Due to the nature of injection molding industry, it takes times to shift from one job to another job, and consequently it makes different waste. So, the job with lower setup waste should be allocated after job one. The resource constraints and sequence dependent wastes were added regarding the essence of the injection industry. Minimizing wastes output and the total tardiness constitute the objective function. All machines are available throughout the time horizon. A job is defined as a single operation job which has a volume, process time, setup time and waste dependent, due date and also a limited source for proceeding called mould. Some jobs has more than one mould for proceeding which means these moulds can setup on identical parallel machines and run simultaneously to deliver the total amount of a job more quickly. Therefore, a job can be split and each part can be processed on at least two machines at a time. Jobs can be done as various lots, independently and asynchronously, but the completion time is calculated based on the biggest part of each job. Jobs are classified into different groups regarding their features. For each new group, all machines need to be installed and installation time is sequence-dependent.

For reaching the optimum solution for the splitting jobs on the parallel machine scheduling with sequence-dependent wastes, a MIL mathematical programming model is presented. The model including N jobs. Processing time of any given job i on machine K illustrates by P(i) is. Due date of every job represents by d(i). The parameters, decision variables and mathematical model are given as follows.

Indices:

M  Set of machines indexed by k = 1, …, m

N  Set of jobs indexed by i and j = 1, …, n

Parameters:

MOi  Total number of machines which can be used to perform job i .

Pi  Processing time of job i.

Qi  Amount of job i.

Si, j  Setup time in case job j coms after job i.

di  Due date of job i.

cik  completion time of job i on the machine k.

c0k  completion time of first job on the machine k.

Wi, j  Sequence dependent setup defective output of job j if it coms after job i.

ak  Amount of defective output of machine k.

Decision variables:

X0, j, k 1 if job j is the first job to process on machine k, 0 otherwise.

yi, k  1 if job i is processed on machine k, 0 otherwise.

Xi, 0, k  1 if job i is the last job to process on machine k, 0 otherwise.

qi, k  Amount of job i on machine k.

Xi, j, k  1 if job j comes after job i on machine k, 0 otherwise.

Zi  Tardiness for job i.

The mathematical model consists of the following equations.

$$ F(1):\mathit{\operatorname{Min}}{\sum}_i{Z}_i $$
(1)
$$ F(2):\mathit{\operatorname{Min}}\ {\sum}_k{\sum}_i{\sum}_j{X}_{i,j,k}\ {W}_{i,j} $$
(2)
$$ {\sum}_k{q}_{ik}={Q}_i\kern4.75em i=1\dots N $$
(3)
$$ {y}_{ik}\ge \frac{q_{ik}}{Q_i}\kern1.75em ,i=1\dots N;\kern0.5em k=1\dots M $$
(4)
$$ {y}_{ik}\le {q}_{ik}\kern1.25em ,\kern0.75em i=1\dots N;\kern0.5em k=1\dots M $$
(5)
$$ {\sum}_j{X}_{ijk}={y}_{ik}\kern1.5em i=1\dots N;\kern0.5em k=1\dots M $$
(6)
$$ {\sum}_i{X}_{ijk}={y}_{ik}\kern2.5em i=1\dots N;\kern0.5em k=1\dots M\kern1.25em $$
(7)
$$ {\sum}_{i=1}^N{y}_{ik}\le 1\kern2em ,\kern1em k=1\dots M $$
(8)
$$ {X}_{ij k}\times {S}_{ij}+{C}_{ik}+{q}_{jk}\times {P}_j\le {C}_{jk}\kern2.5em i=0\dots N;\kern0.5em j=1\dots N;k=1\dots M $$
(9)
$$ {C}_{ik}-{d}_i\le {Z}_i\kern2em i=1\dots N;\kern0.75em k=1\dots M $$
(10)
$$ {C}_{ok}=0\kern1.75em k=1\dots M $$
(11)
$$ {X}_{ijk}\in \left\{0,1\right\}\kern2em i,j=0\dots N;\kern0.5em k=1\dots M $$
(12)
$$ {q}_{ik}\ge 0;\kern0.75em {C}_{ik}\ge 0;\kern0.75em {Z}_i\ge 0;\kern0.5em {y}_{ik}\ge 0\kern1.25em {X}_{iik}=0;\kern0.5em i=1\dots N;\kern1.25em k=1\dots M $$
(13)
$$ {\sum}_i{\sum}_j{X}_{ij k}\times {W}_{ij}\le {a}_k\kern1.25em k=1\dots M $$
(14)
$$ {\sum}_k{y}_{ik}\le {Mo}_i,\kern4.5em i=0\dots N $$
(15)

The objective function is presented by Eqs. (1) and (2) for minimizing total tardiness and total wastes, respectively. Processing all parts of each job is guaranteed by means of Eq. (3). Equations (4) and (5) guarantee the assigning all parts of each job to machines. Equations (6)–(8) represent that each part of the job will not process more than once. Equation (9) calculates the completion time. Calculating tardiness for each job is done using Eq. (10). Equations (11)–(13) define variables of the model. Waste of each machine is shown in Eq. (14). Finally, Eq. (15) is used for defining the resource limitation of molds.

4 Optimization by volleyball premier league (VPL) algorithm

VPL algorithm is a newly proposed meta-heuristic algorithm founded on some realistic metaphors of a volleyball tournament [82]. To show validity of the proposed algorithm, a comprehensive analyse based on 23 test instances is organized in which nine eminent metaheuristic comprising GA, Harmony Search algorithm (HS), Fire Fly Algorithm (FA), Soccer League Competition (SLC), Differential Evolution (DE), Sin Cosine Algorithm (SCA), PSO, League Championship Algorithm (LCA) and Artificial Bee Colony (ABC) were compared to the suggested VPL. The results of the VPL show its superiority in comparison to its rivals. To show the main structure of this algorithm, we have to consider the coach, formation of team and substitute bench. These parts play important role in designing the algorithm which has heightened the contribution of VPL in exploration and exploitation. The solution representation of this algorithm contains two parts: firstly, active part contains formation (Ƒ) used for evaluating fitness function, and secondly, passive part stores the information of substitution (Ѕ) which is used for specific strategy in this part. In this paper, we use ℤƑ and ℤЅ for evaluating the active and passive parts, respectively. Typically, there are more than two teams in a division, in which all the teams must play to each other in a specific period. To reach this goal, an organized schedule is proposed by authority of division about the time and place of each game that are exactly determined. In VPL, an specefic procedure, Single Round Robin (SRR) method, is used to generate league schedule. Logically, winner and loser teams are determined after the match. In VPL, the loser team follows some specific process, i.e. knowledge sharing strategy, repositioning strategy and finally substitution strategy to enhance its wretched performance, and accordingly, the winner team tracks the other process to keep its encouraging performance. In the next stage of knowledge sharing strategy, the valuable information is shared among all team members to reach a comprehensive insight throughout the competition, and then, some players are swapped with ones who are positioned in substitutions bench with respect to repositioning strategy in which the information of Ƒ and Ѕ property is exchanged to each other. Another process used in this algorithm which is followed by winner teams is winner strategy, in which the winner teams in competition reconstruct their position according to the position of the best team in the league. The most important stage of this algorithm is implemented as learning strategy in which all position of teams are updated for Ƒ and Ѕ property based on top 3 teams participating in the current tournament. In the season transfer strategy, players can transfer to other teams. And finally, like a regular sport, the best teams of lower division are replaced to the worst teams of upper division by using promotion and relegation process. So that, the worst teams are removed from the league, and the new solutions generated randomly are added to the league as promoted teams from the lower division, respectively.

5 Multi-objective Pareto-based Metaheuristic algorithm

In this section, the basic assumptions of constrained multi-objective optimization are considered. A common continuous and unconstrained multi-objective optimization problem can be written as follows [83].

$$ {\min}_{\mathrm{x}\upepsilon \mathcal{F}}={\left({\mathrm{f}}_1\left(\mathrm{u}\right),{\mathrm{f}}_2\left(\mathrm{u}\right),\dots, {\mathrm{f}}_{\mathrm{m}}\left(\mathrm{u}\right)\right)}^{\mathrm{T}}, $$
(16)

where \( \mathcal{F} \) is defined as the decision space, u = (u1, u2, …, un) is denoted as a decision vector defined in \( \mathcal{F} \), and m is the number of objective functions.

The following mapping function defines the objective function bounded in the decision space:

$$ \mathrm{F}:\mathcal{F}\to {\mathrm{R}}^{\mathrm{m}}, $$
(17)

where Rm is denoted as the objective space of the problem. In any multi-objective problem, it is inevitable to sacrifice one goal while the decision maker is trying to improve the other goal.

Definition 1. (Pareto-dominance): Let u = (u1, u2, …, un) and v = (v1, v2, …, vn) are defined as decision vectors, u is said to dominate vector v (shown as u ≻ v) if and only if

$$ \left(\forall \mathrm{i}\in \left\{1,2,\dots, \mathrm{n}\right\}:{\mathrm{f}}_{\mathrm{i}}\left({\mathrm{u}}_{\mathrm{i}}\right)\le {\mathrm{f}}_{\mathrm{i}}\left({\mathrm{v}}_{\mathrm{i}}\right)\right)\wedge \Big(\exists \mathrm{j}\left\{1,2,\dots, \mathrm{n}\right\}:{\mathrm{f}}_{\mathrm{i}}\left({\mathrm{u}}_{\mathrm{i}}\right)<{\mathrm{f}}_{\mathrm{i}}\left({\mathrm{v}}_{\mathrm{i}}\right) $$
(18)

Definition 2. (Pareto-optimal): A decision vector u is defined as a Pareto optimal if and only if

$$ \neg \exists \mathrm{v}\in \mathcal{F}:\mathrm{v}\succ \mathrm{u} $$
(19)

Definition 3. (Pareto-optimal set): The set PS includes all Pareto optimal decision vectors, defined as

$$ \mathrm{PS}=\left\{\mathrm{u}|\neg \exists \mathrm{v}\in \mathcal{F}:\mathrm{v}\succ \mathrm{u}\right\}` $$
(20)

Definition 4. (Pareto-optimal front): The set PF comprises the values of all the objective functions point out the Pareto-optimal decision vectors in PS.

$$ \mathrm{PF}=\left\{\mathrm{F}\left(\mathrm{u}\right)={\left({\mathrm{f}}_1\left(\mathrm{u}\right),{\mathrm{f}}_2\left(\mathrm{u}\right),\dots, {\mathrm{f}}_{\mathrm{m}}\left(\mathrm{u}\right)\right)}^{\mathrm{T}}|\mathrm{u}\in \mathrm{PS}\right\} $$
(21)

5.1 Multi-objective volleyball premier league

This part devoted to elucidating the main structure of Multi-Objective Volleyball Premier League (MOVPL). To date various methods have been developed and introduced to cope with multi-objective function in MOEAs. To manage the problem, in this algorithm we use crowding distance concept used in NSGA-II. Generally speaking, the execution of MOVPL follows the main structure of basic VPL algorithm which performs the crowding distance concept (Table 1).

Table 1. Psoudeo code for MOVPL.
Full size table

As can be seen in the above Psoudeo code, the MOVPL begins with defining parameters and then letting the set PF equal null, the next step is setting first population as Pt. Evaluating Pt and generating fronts based on non-domination sorting method and crowding distance form the next step. After that set the non-dominated solutions as set of PF.

Afterward, the main loop will be started, and the loop is repeated until the specified number of iteration is done. At the beginning in the main loop, set the Rt = Pt and then do the following steps. In the first step in the main loop, winner and loser teams should be determined for Rt and follow the knowledge sharing and winner strategy process. In the following step of learning phase, the season transfer process and repositioning strategy must be done for Rt. The Qt will be generated by the non-dominated sorting method. Consequently, according to each level of non-dominated sorting method, the solution ranks are estimated where the former stage consists of resolutions which have the top rates. Subsequently, the crowding distance among individuals is calculated for every stage with the purpose of using the selection procedure. For calculating the crowding distance, total resolutions for the first stage are denoted as Dl = 0, afterwards per any objective function, the set is sorted in the worst-to-best order. Therefore, the vector of sorted indices is presented as follows:

$$ {I}_m= sort\left({f}_m>\right) $$
(22)

Let m represent the number of objectives, crowding distance of lth resolution per any level is calculated using a linear distance criterion mention below:

$$ {D}_{I_l^m}={D}_{I_l^m}+\frac{f_m^{\left({I}_{l+1}^m\right)}-{f}_m^{\left({I}_{l-1}^m\right)}}{f_m^{max}-{f}_m^{min}} $$
(23)

Where \( {f}_m^{max} \) and \( {f}_m^{min} \) denote the maximum and minimum values of the m-th objective function, respectively. Moreover, lowest and highest amounts of function are assigned by unbounded amounts of distance for each objective function.

An operator called binary tournament is used for selection. Calculating the rating of the crowding distance is required for each member, therefore the first two members are chosen in the population. For the sake of selecting a member, the crowding distance and the rating for each solution are considered. Afterwards the solution is chosen based on lower rating. In case of having an equal rank, the one with the larger crowding distance will be chosen. Then, the selection, crossover and mutation operators are used to generate new crowd of offspring (Qt) with a size of N. The algorithm introduces Rt as an integration of parents (Pt) and offspring, therefore ∣Rt ∣  = 2N. Next, Rt the above method is used for sorting, and N best members are chosen for the subsequent descendant (Pt + 1). The number of iteration of this procedure relies on the stopping condition. In the last step of MOVPL, the best solutions for the multi-objective optimization are calculated and illustrated in Pareto fronts set. Solution representation relies on the type of optimization problem.

This section defines the solution representation for solving the problem using in MOEAs algorithm. This special case of the scheduling problem is to deal with job sequencing and limitation of machine availability regarding related molds. Therefore, we suggest a three-part encoding scheme, namely, α = [x1; x2; x3] where the size of all three parts is same. The following formula is used to determine the size of each part.

$$ n={\sum}_{i=1}^N Mo(i), $$
(24)

where Mo(i) represents the total number of eligible injection machines to carry out the operations Oi, and n defines the size of each section for any given solutions. On sequent, the size of any solution is equal to 3 × n. Principally, n represents the length of any part of the solution representation. The outline of each solution is illustrated in Figs. 1, 2 and 3

Fig. 1
figure 1

General scheme of a solution

Full size image
Fig. 2
figure 2

The output of the Taguchi method

Full size image
Fig. 3
figure 3

Graphical comparisons of all MOEAs in terms of different metrics

Full size image
$$ \alpha =\left\{{x}_1=\left\{{x}_{11},{x}_{12},\dots {x}_{1n}\right\}\kern0.5em {x}_2=\left\{{x}_{21},{x}_{22},\dots, {x}_{2n}\right\}\kern0.5em {x}_3=\left\{{x}_{31},{x}_{32}\dots {x}_{3n}\right\}\right\} $$
$$ A=\left\{{X}_1=\left\{{X}_{11},{X}_{12},\dots {X}_{1n}\right\}\kern0.5em {X}_2=\left\{{X}_{21},{X}_{22},\dots, {X}_{2n}\right\}\kern0.5em {X}_3=\left\{{X}_{31},{X}_{32}\dots {X}_{3n}\right\}\right\} $$

Based on Fig. 3, set α would be changed to set Α throughout solution procedure. Consequently, set Α has three parts as the same as set α. The former section of the set α illustrates weight of any given solution to determine the final solution. Moreover, to transform x1 to X1, a method has been considered below:

No table of figures entries found. Psoudeo code for calculating X1

figure a

For illustrating this encoding related to the former section of the solution, in the following the solution representation for 4 machines and 5 jobs is given for set x, Mo and Q:

$$ x=\left\{\begin{array}{c}{x}_1=\left\{0.82\ 0.98\ 0.73\ 0.34\ 0.58\ 0.11\ 0.91\ 0.88\ 0.81\ 0.26\ 0.59\ 0.02\right\}\\ {}{x}_2=\left\{0.43\ 0.31\ 0.16\ 0.18\ 0.42\ 0.09\ 0.60\ 0.47\ 0.70\ 0.71\ 0.64\ 0.03\right\}\\ {}{x}_3=\left\{0.07\ 0.32\ 0.53\ 0.65\ 0.41\ 0.82\ 0.72\ 0.97\ 0.53\ 0.33\ 0.11\ 0.61\right\}\end{array}\right. $$
$$ Mo=\left[2\kern0.5em 2\kern0.5em 3\kern0.5em 3\kern0.5em 2\right] $$
$$ Q=\left[10000\kern0.5em 12000\kern0.5em 14000\kern0.5em 8000\kern0.5em 25000\right] $$

In this example, and using Eq. (24), the size of solution is n = 12. So, due to the aforementioned method, the former section of X is presented in the following:

$$ {X}_1=\left[4562\kern0.5em 5438\kern0.5em 8158\kern0.5em 3842\kern0.5em 5117\kern0.5em 944\kern0.5em 7939\kern0.5em 3594\kern0.5em 3341\kern0.5em 1065\kern0.5em 24088\kern0.5em 912\right] $$

The latter section related to machine assignment. For determining X2, the following function has been considered.

No table of figures entries found. Psoudeo code for calculating X2

figure b

Using the aforementioned method, the following vector will be generated as the second part of set X. The first operation is assigned to machine 1, the second to machine 2, and so on.

$$ {X}_2=\left[2\kern1.25em 1\kern1.25em 1\kern1.25em 1\kern1.25em 2\kern1.25em 1\kern1.25em 2\kern1.25em 2\kern1.25em 3\kern1.25em 3\kern1.25em 2\kern1.25em 1\right] $$

The final part of the solution representation is related to the sequencing of jobs. On that account, X3 is computed according to sorting ranks of the second part of set X members in ascending order. So, X3 is given as follows:

$$ {X}_3=\left[1\kern0.5em 11\kern0.5em 2\kern0.5em 10\kern0.5em 5\kern0.5em 3\kern0.5em 9\kern0.5em 12\kern0.5em 4\kern0.5em 7\kern0.5em 6\kern0.5em 8\right] $$

6 Experiments and output analysis

The following five standard metrics are represented with the aim of comparing the performances of the proposed Pareto-based multi-objective algorithms.

  • Number of Pareto Solutions (NPS): this metric measures the number of the Pareto solutions in Pareto optimal front [84]:

$$ NPS=\mid PF\mid $$
(25)

  • The Spacing (SP): this metric [85] is used to quantify and measure the convergence of the algorithm:

$$ SP=\sqrt{\left(\frac{1}{\left| PF\right|}\right)\sum \limits_{i=1}^n{\left(\overline{d}-{d}_i\right)}^2} $$
(26)

The following metric measures the standard deviation of the distances between members, where\( {d}_i={\mathit{\min}}_{\overrightarrow{y}\in PF}{\sum}_{k=1}^K\left(|{o}_k\left(\overrightarrow{x}\right)-{o}_k\left(\overrightarrow{y}\right)|\right) \) represents the minimum distance of resolution \( \overrightarrow{x} \) from the rest of resolutions and \( \overline{d} \) represents the median amount for di.

  • CPU Time (CT): this metric represents the computational (CPU) time of running the algorithm [84].

  • Maximum Spread (MS): Eq. (27) illustrates the spread metric of Pareto solutions [84]:

$$ MS=\sqrt{\sum \limits_{i=1}^I{\left(\min {f}_i-\max {f}_i\right)}^2} $$
(27)

Where minfi and maxfi for each fitness function are the lowest and highest values which are calculated for all non-dominated members by various manners [86].

  • Error Ratio (ER): the ER is computed using the following formula:

$$ ER\frac{\sum_{h=1}^n{e}_h}{n} $$
(28)

Where eh is equal to one if a solution is not related to Pareto front otherwise zero, and n denotes the population size.

Ten test problems are utilized for experimentation. In addition, each problem is run thirty times with the aim of reducing uncertainties. Afterward, the last stage is done by calculating the averages of thirty runs. The Taguchi method [87] is utilized to detect the optimum level for any parameters in various algorithms. To calculating the amount of response variations, the signal to noise ratio (S/N ) is used [88]:

$$ S/N=-10\times \mathit{\log}\left(S\left({Y}^2/N\ \right)\right), $$
(29)

Where N denotes the number of orthogonal arrays and Y idicates the response. To accomplish the procedure, identifying the levels of parameters needs to be done as in Table 2.

Table 2 Spans and stages of MOEA parameters
Full size table

To make comparison among different algorithms, Table 2. presents the multi-objective metrics amounts on the 10 test problems. It should be noted that Matlab software (Version 8.3.0.532, R2014a) is used to code the proposed meta-heuristic algorithms, and the programs are executed on a 2.50 GHz core i5 CPU with 4.00 GB of RAM under iOS operating system.

Orthogonal arrays are used in Taguchi method with the aim of studying all factors concurrently. Afterwards, the L9 design is used for NSGA-II, MOEA/D, MOGWO, and MOALO at the same time, L27 is performed for MOPSO, SPEA-II, and MOVPL up to now. This approach has been executed in various researches, for instance [89, 90]. Accordingly, the schemes of consequence for S/N ratios of all four compaed algorithms are illustrated in Fig. 2. Note that WOGWO and MOALO are two relatively new bio-inpired metaheuristics, where MOGWO stands for multi-objective grey wolf optimization [91], and MOALO stands for multi-objective ant lion optimization [92].

According to Fig. 2, the best levels of each algorithm are determined. To show compatibility of obtained results, we design 10 problems in this way. Results of different metrics for these problems are shown in Table 3.

Table 3 MOEA metrics for proposed Pareto-based metaheuristics on 10 problems
Full size table

Additionally, graphical comparisons of 10 problems solved using the proposed algorithms are presented in Table 2 for five multi-objective metrics, i.e. CPU time consumed (CP), number of pareto solutions (NPS), spacing (SP), maximum spread (MS), and error ratio (ER).

As can be seen in Table 2 and Figs. 4 and 5, in terms of computation time (i.e. CP metric) our proposed algorithm consumes significantly less time for almost all problems in comparison with the other six metaheuristic algorithms. Only for problems 3 and 6 the NSGA-II algorithm takes less time in comparison with the MOVPL (166 vs. 175 and 77 vs. 82, respectively). In terms of the NPS metric, just the MOEA/D algorithm gains a little bit better Pareto front solutions for problems 1, 3, 4, 6 and 9 compared to the MOVPL. In spite of the fact that the MOGWO algorithm defeats the MOVPL algorithm just for problem 2 in terms of the SP metric, in all other cases, the situations are vise versa. In the MS metric, the MOVPL algorithm defeats all algorithms except the NSGA-II algorithm for problem 3. For the last ER metric, the MOVPL wins all algorithms for all problems, except the MOEA/D algorithm for problems 3, 4, 9 and 10. Therefore, it can be inferred that the MOVPL is superior to the other six well-known metaheuristic algorithms based on the five aforementioned metrics.

Fig. 4
figure 4

Pareto fronts for 10 test problems

Full size image
Fig. 5
figure 5

Solutions in the final Pareto front of MOEAs

Full size image

The Pareto fronts of 10 problems obtained using proposed algorithms are illustrated in Figs. 4 and 5. As can be seen, each algorithm has a different Pareto front for its test function. The Pareto front of the tardiness objective function (F1) and the defective objective function (F2) are illustrated on all 10 problems. The solutions obtained by MOVPL and MOGWO indicate that both F1 and F2 are equally preferred, while F1 objective function for problems 2–8 is less important and dominated by F2 objective function from solutions obtained by MOPSO, MOALO and MOEA/D. By contrast, the F2 objective function is dominated by these problems by other algorithms. It is worth noting that for all problems the MOVPL equally prefers both F1 and F2 and gains minimal solutions in comparison with other metaheuristic algorithms.

In the above figure, for the purpose of comparing all algorithms, all obtained Pareto fronts are presented simply. To gain the best solution, it needs to cumulate Pareto front solutions offered by all. Let \( \mathcal{L} \) is a set containing every individual of PF offered by all MOEAs, defined as follows:

$$ \mathcal{L}=\left\{x\in {PF}_i\right|i\in \left\{ MOEAs\right\}\Big\} $$
(30)

Because there is no general agreement between researchers about which MOEA has supreme performance in comparison with other algorithms, the Pareto fronts discover process has been accomplished for set \( \mathcal{L} \) to specify ultimate PF for all test problems. Based on the feasibility of solutions and in comparison to the other algorithms, one (or maybe two) algorithm(s) will be better solution technique. Fig. 5 illustrates the digit of solutions in the ultimate Pareto front based on all MOEAs for total test problems.

As can be seen in Fig. 5, MOVPL has the best performance for all problems. One thing which is equally important that in problems 4 and 5, the NSGA-II, the MOPSO, and the MOALO have better performance among other problems; though the competition is between MOVPL and the NSGA-II, and the former has superiority on the latter. It can be also argued that in problem 5 the NSGA-II loses its performance and the MOVPL shows a huge superiority on the MOPSO and MOALO algorithms.

The suggested MOEAs are actuarially compared from various metrics with no major difference. Hence, the ANOVA experiments (analysis of variance) are guaranteed at at 95% confidence level. The outcomes of ANOVA test are shown in Table 4.

Table 4 Results of ANOVA test at a 95% of confidence interval
Full size table

The seven algorithms have been ranked according to the five comparison metrics. As shown in Table 5, in terms of CP, SP, MS and ER, the MOVPL is able to defeat the other rival algorithms. The NPS is the only metric in which the MOEA/D has pushed down the MOVPL to the second rank. Overall, the MOVPL has won four out of five metrics against the other six well-known multi-objective evolutionary algorithms.

Table 5 Rankings of the MOEAs according to the five metrics
Full size table

7 Conclusions

This research focused on the parallel machine scheduling with splitting jobs on a set of identical machines minimizing wastes and total tardiness. The mixed integer programming model is formulated using data from a real plastic injection industry. The nature of the problem in this industry is that two conflicting objectives, minimizing total tardiness and minimizing total wastes, should be taken into account. To make the model more realistic, it has been developed under a multi-objective framework, making the model difficult to solve in optimality using classical solution methods.

Regarding the fact that the problem at hand is NP-hard, metaheuristic algorithms have been applied. Since the mathematical model of the problem is bi-objective within the multi-objective framework, MOPSO, NSGA-II, MOEA/D MOGWO, MOALO and SPEA2 are compared with the proposed new algorithm called MOVPL. The MOVPL uses the crowding distance concept applied in NSGA-II as an extension of the Volleyball Premier League (VPL). For comparing the efficiency of all rival algorithms, the five most common standard metrics in the field of multi-objective algorithms (i.e. NPS, SP, CT, MS and ER) are used. As can be seen from the computational results, the proposed MOVPL algorithm has reached the first rank in four standard metrics except for the Number of Pareto Solutions (NPS) in which the MOVPL has been ranked the second position.

The originality of the paper is twofold. Firstly, we build a mathematical model for a challenging multi-objective scheduling problem in real-world plastic injection molding industry. Secondly, based on the original VPL algorithm which was proposed for the optimization on single objective and was evaluated on some simple test functions, in this paper we redesign the algorithm to enable global optimization on multiple objectives and evaluate it on the NP-hard real scheduling problem against some state-of-the-art algorithms.

As our future work, we will continue to refine the proposed MOVPL algorithm by taking into account new constraints to the problem such as preemption, maintenance times, and uncertainty in the volume of jobs.