A hybrid multi-objective gray wolf optimization algorithm for a fuzzy blocking flow shop scheduling problem

Abstract

Blocking flow shop scheduling problems have important applications in manufacturing. Because of the imprecise and vague temporal parameters in real-world production, this article formulates a fuzzy blocking flow shop scheduling problem with fuzzy processing time and fuzzy due date in order to minimize the fuzzy makespan and maximize the average agreement index. To solve this combinational optimization problem, a hybrid multi-objective gray wolf optimization algorithm is proposed. The hybrid multi-objective gray wolf optimization utilizes the largest position value rule for solution representation, employs a dynamic maintenance strategy to maintain an archive, and develops a thorough mechanism for leader selection. In the hybrid multi-objective gray wolf optimization, a novel heuristic process is designed to generate initial solutions with a certain quality, and a local search strategy is embedded to improve the exploitation capability. The performance of the hybrid multi-objective gray wolf optimization is tested on the production instances of panel block assembly in shipbuilding. Computational comparisons of the hybrid multi-objective gray wolf optimization with two other well-known multi-objective evolutionary algorithms demonstrate the feasibility and effectiveness of the hybrid multi-objective gray wolf optimization in generating optimal solutions to the bi-criterion fuzzy blocking flow shop scheduling problem.

Introduction

In the classical flow shop scheduling problem, an assumption is made that there are intermediate buffers of infinite capacity between two consecutive machines where jobs can be stored to wait for subsequent operations. However, some production environments, such as concrete block production, ¹ steel manufacturing, ² and panel block assembly in shipbuilding, lack intermediate buffers between machines, mainly because of technical requirements, process characteristics, or some other constraints. The flow shop scheduling problem without intermediate buffers is known as the blocking flow shop scheduling problem. It considers that a set of $n$ jobs are to be processed on $m$ machines in the same order, from machine 1 to machine $m$ ; moreover, a job which has completed processing on a machine must remain on and block the machine until the downstream machine is available for processing. This problem has received increasing attention in the literature, and the majority of the studies on this topic consider minimizing the makespan, that is, the maximum completion time of all jobs. The blocking flow shop scheduling problem with a makespan criterion is usually denoted as $F_m|\mathrm{blocking}|C_{max}$ . The problem with more than two machines $(m>2)$ has been proved to be non-deterministic polynomial-time (NP)-hard in the strong sense. ³

For more than a decade, various heuristics have been proposed for the $F_m|\mathrm{blocking}|C_{max}$ problem. For example, McCormick et al. ⁴ developed the profile fitting (PF) heuristic, which attempts to arrange jobs for a minimum sum of the idle and blocking time on machines. Based on the makespan properties proposed by Ronconi and Armentano, ⁵ Ronconi ⁶ presented the MinMax (MM) heuristic and constructed the MME and PFE heuristics by combining MM and PF with the enumeration procedure of the Nawaz–Enscore–Ham ⁷ (NEH) heuristic. The MME and PFE heuristics were subsequently shown to be superior to the NEH, MM, and PF heuristics. A growing number of meta-heuristics have also been developed for this type of problem, such as the genetic algorithm (GA), ⁸ tabu search approaches (TS and TS + M), ⁹ discrete particle swarm optimization (DPSO) algorithm, ¹⁰ and discrete self-organizing migrating algorithm (DSOMA). ¹¹ Recently, Shao et al. ¹² proposed an estimation of distribution algorithm with path relinking technique (P-EDA) and compared it with other meta-heuristics, such as the hybrid discrete differential evolution (HDDE) algorithm, ¹³ iterated greedy (IG) algorithm, ¹⁴ revised artificial immune system (RAIS) algorithm, ¹⁵ memetic algorithm (MA), ¹⁶ and modified fruit fly optimization (MFFO) algorithm. ¹⁷

In most previous works, including those mentioned above, scheduling problems were considered in deterministic environments where the parameters, such as the processing time and due date, were taken as crisp values. However, the assumption of precise temporal parameters is frequently violated in practice. During production, the processing time is often affected by uncertainty arising from both machine and human factors. The due date is expected to be met; however, certain earliness and tardiness limits can be tolerated, and extended limits will have lower values. Thus, modeling scheduling problems with imprecise and vague temporal parameters is more reasonable in real-world production. Fortunately, based on the concept of fuzzy sets, ¹⁸ uncertain processing time and due date can be represented as fuzzy parameters to model a fuzzy scheduling problem.^19,20 For instance, Wu ²¹ formulated a fuzzy flow shop scheduling problem in order to minimize the weighted sum of fuzzy earliness and fuzzy tardiness. Lai and Wu ²² considered a fuzzy flow shop scheduling problem with the fuzzy makespan objective function. Because simultaneously optimizing the makespan and due date satisfaction is usually involved in decision-making, Sakawa and Kubota ²³ proposed a GA for a multi-objective fuzzy job shop scheduling problem that maximized the minimum agreement index and average agreement index and minimized the fuzzy makespan. Lei ²⁴ developed a Pareto archive particle swarm optimization (PAPSO) algorithm for a multi-objective fuzzy job shop scheduling problem that maximized the minimum agreement index and minimized the fuzzy makespan and mean fuzzy completion time. Behnamian and Ghomi ²⁵ proposed a bi-level algorithm for a bi-objective fuzzy hybrid flow shop scheduling problem that minimized the makespan and the sum of earliness and tardiness. Wang et al. ²⁶ designed an MA for a multi-objective fuzzy flexible job shop scheduling problem that simultaneously optimized the fuzzy makespan, minimum agreement index, and average agreement index. In consideration of the vague temporal parameters as well as the simultaneous optimization of the makespan and due date satisfaction, this article formulates a bi-criterion fuzzy blocking flow shop scheduling problem with fuzzy processing time and fuzzy due date in order to minimize the fuzzy makespan and maximize the average agreement index.

The gray wolf optimizer (GWO) developed by Mirjalili et al. ²⁷ is a new population-based and easy-to-implement meta-heuristic that has been successfully applied to various engineering problems, such as unmanned combat aerial vehicle path planning ²⁸ and economic load dispatch. ²⁹ Importantly, the GWO was employed to solve a flow shop scheduling problem and exhibited high competitiveness. ³⁰ Recently, Mirjalili et al. ³¹ proposed a multi-objective GWO for continuous multi-objective problems, and Lu et al. ³² designed an effective multi-objective discrete GWO for a scheduling problem in welding production. In this article, a hybrid multi-objective gray wolf optimization (HMOGWO) algorithm is proposed to solve the bi-criterion fuzzy blocking flow shop scheduling problem.

The remainder of this article is organized as follows. Section “Fuzzy blocking flow shop scheduling problem” describes the bi-criterion fuzzy blocking flow shop scheduling problem. Section “Operations on fuzzy numbers” introduces the operations on fuzzy numbers that are required to formulate the scheduling problem. Section “HMOGWO for the fuzzy blocking flow shop scheduling problem” elaborates the proposed algorithm for solving the scheduling problem. Computational analysis based on real-time production data is reported in section “Computational analysis,” followed by the conclusions in section “Conclusion and future studies.”

Fuzzy blocking flow shop scheduling problem

Problem description

The fuzzy blocking flow shop scheduling problem can be defined as follows. A set of jobs are to be processed sequentially on machine 1, machine 2, and so on until the final machine. All machines are continuously available, and all jobs are ready for processing at time zero. At any time, each machine can process at most one job, and each job can be processed on at most one machine. The processing of a job on a machine cannot be interrupted. A job that has completed processing on a machine cannot leave the machine until the downstream machine is free. The fuzzy processing time and fuzzy due date are represented by fuzzy numbers.

Problem formulation

The notations of the fuzzy blocking flow shop scheduling problem are defined as follows:

• $n$ is the number of jobs;

• $m$ is the number of machines;

• $i$ is the index of jobs, $i\in \{1,2,\dots ,n\}$ ;

• $j$ is the index of machines, $j\in \{1,2,\dots ,m\}$ ;

• ${\tilde{p}}_{i,j}$ is the fuzzy processing time of job $i$ on machine $j$ ;

• ${\tilde{d}}_i$ is the fuzzy due date of job $i$ ;

• ${\tilde{D}}_{i,j}$ is the fuzzy departure time of job $i$ on machine $j$ ;

• ${\tilde{C}}_i$ is the fuzzy completion time of job $i$ ;

• $A{I}_i$ is the agreement index of job $i$ .

In this study, the fuzzy processing time is a triangular fuzzy number denoted by ${\tilde{p}}_{i,j}=(p_{i,j}^O,p_{i,j},p_{i,j}^P)$ , which includes the following three parameters: the optimistic value $(p_{i,j}^O)$ , the most plausible value $(p_{i,j})$ , and the pessimistic value $(p_{i,j}^P)$ . The membership function of the triangular fuzzy processing time is shown in Figure 1(a). The fuzzy due date is considered a trapezoidal fuzzy number denoted by ${\tilde{d}}_i=(d_i^L,d_i^{E_1},d_i^{E_2},d_i^U)$ , where $d_i^L$ and $d_i^U$ are the lower and upper bounds of the fuzzy due date, respectively, and $d_i^{E_1}$ and $d_i^{E_2}$ represent the expected due date interval $(d_i^{E_1},d_i^{E_2})$ . The membership function of the trapezoidal fuzzy due date, which is shown in Figure 1(b), represents the degree of satisfaction with respect to the completion time.

Figure 1.

(a) Membership function of the triangular fuzzy processing time; (b) membership function of the trapezoidal fuzzy due date; and (c) agreement index (AI).

Let $\mathbf{\pi }\mathbf{=}\mathbf{[}\mathbf{\pi }\mathbf{(}1\mathbf{)}\mathbf{,}\mathbf{\pi }\mathbf{(}2\mathbf{)}\mathbf{,}\dots \mathbf{,}\mathbf{\pi }\mathbf{(}\mathbf{n}\mathbf{)}\mathbf{]}$ denote a processing sequence of $n$ jobs. Suppose that job $i$ is allocated at the kth dimension of $\mathbf{\pi }$ . Then, the fuzzy completion time of each job can be calculated using the following formulas

$${\tilde{D}}_{\pi (1),0}=(0,0,0)$$

(1)

$${\tilde{D}}_{\pi (1),j}={\tilde{D}}_{\pi (1),j-1}+{\tilde{p}}_{\pi (1),j},j\in \{1,2,\dots ,m\}$$

(2)

$${\tilde{D}}_{\pi (k),0}={\tilde{D}}_{\pi (k-1),1},k\in \{2,3,\dots ,n\}$$

(3)

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{D}}_{\pi (k),j}=max({\tilde{D}}_{\pi (k),j-1}+{\tilde{p}}_{\pi (k),j},{\tilde{D}}_{\pi (k-1),j+1}),}\\ k\in \{2,3,\dots ,n\},j\in \{1,2,\dots ,m-1\}\hfill \end{array}$$

(4)

$${\tilde{D}}_{\pi (k),m}={\tilde{D}}_{\pi (k),m-1}+{\tilde{p}}_{\pi (k),m},k\in \{2,3,\dots ,n\}$$

(5)

$${\tilde{C}}_i={\tilde{D}}_{\pi (k),m}$$

(6)

where ${\tilde{D}}_{\pi (k),0}$ denotes the starting time of job $i$ on the first machine, and ${\tilde{D}}_{\pi (k),j}$ represents the fuzzy departure time of job $i$ on machine $j$ . The fuzzy completion time, which is denoted by ${\tilde{C}}_i=(C_i^O,C_i,C_i^P)$ , includes three parameters: the optimistic value $(C_i^O)$ , the most plausible value $(C_i)$ , and the pessimistic value $(C_i^P)$ .

The completion time is always expected to meet the due date. The agreement index $(\mathrm{AI})$ ²³ of job $i$ , which is defined in equation (7) and is shown in Figure 1(c), indicates the degree of compliance between ${\tilde{C}}_i$ and ${\tilde{d}}_i$ .

$$A{I}_i=\frac{\mathrm{area}({\tilde{C}}_i\bigcap {\tilde{d}}_i)}{\mathrm{area}({\tilde{C}}_i)}$$

(7)

This article formulates the fuzzy scheduling problem with the bi-objective of minimizing the fuzzy makespan and maximizing the average agreement index $(\overline{\mathrm{AI}})$

$$\mathrm{Minimize}\mathrm{makespan}=\underset{i=1,2,\dots ,n}{max}{\tilde{C}}_i$$

(8)

$$\mathrm{Maximize}\overline{\mathrm{AI}}=\frac{1}{n}\sum _{i=1}^{n}A{I}_i$$

(9)

Operations on fuzzy numbers

Several operations, including the addition operation (+), max operation $(\vee )$ , and ranking method for two or more fuzzy numbers, are essential for the formulation of the scheduling problem.

For two triangular fuzzy numbers $\tilde{A}=(a^1,a^2,a^3)$ and $\tilde{B}=(b^1,b^2,b^3)$

$$\tilde{A}+\tilde{B}=(a^1+b^1,a^2+b^2,a^3+b^3)$$

(10)

With respect to the max operation, an approximation proposed by Lei ³³ is applied in this article. According to Lei’s criterion, the approximate max is either $\tilde{A}$ or $\tilde{B}$ , as given below

$$\mathrm{If}\tilde{A}>\tilde{B},\mathrm{then}max(\tilde{A},\tilde{B})=\tilde{A}\vee \tilde{B}=\tilde{A}$$

(11)

Equations (4), (8) and (11) show that a ranking method for fuzzy numbers is required. This article uses the following three criteria to rank triangular fuzzy numbers: ²³

• Criterion 1.

The greatest associate ordinary number

$$C_1(\tilde{A})=\frac{a^1+2a^2+a^3}{4}$$

(12)

is used as the first criterion to rank the triangular fuzzy numbers.

• Criterion 2.

If $C_1$ does not rank the fuzzy numbers, then the best maximal presumption

$$C_2(\tilde{A})=a^2$$

(13)

is selected as the second criterion.

• Criterion 3.

If $C_1$ and $C_2$ do not rank the fuzzy numbers, then the difference of the spreads

$$C_3(\tilde{A})=a^3-a^1$$

(14)

is utilized as the third criterion.

HMOGWO for the fuzzy blocking flow shop scheduling problem

GWO

The GWO ²⁷ is inspired by the social hierarchy and hunting behavior of gray wolves. To mimic the social hierarchy of gray wolves when developing the GWO, the fittest solution is regarded as the alpha $(\alpha )$ wolf. In addition, the second and third best solutions are expressed as the beta $(\beta )$ and delta $(\delta )$ wolves, respectively. The remaining candidate solutions are deemed as the omega $(\omega )$ wolves. The search (hunt) process is guided by $\alpha$ , $\beta$ , and $\delta$ , which are concomitantly followed by the $\omega$ wolves. Gray wolves habitually encircle prey during the hunt. The equations that simulate this encircling behavior are presented as follows

$$\mathit{D}=|\mathit{C}·{\mathit{X}}_p(t)-\mathit{X}(t)|$$

(15)

$$\mathit{X}(t+1)={\mathit{X}}_p(t)-\mathit{A}·\mathit{D}$$

(16)

where $t$ denotes the current iteration, $X_p$ represents the position vector of the prey, and $\mathit{X}$ represents the position vector of a gray wolf. $\mathit{A}$ and $\mathit{C}$ are coefficient vectors that are calculated from the following equations

$$\mathit{A}=2·\mathit{a}·{\mathit{r}}_1-\mathit{a}$$

(17)

$$\mathit{C}=2·{\mathit{r}}_2$$

(18)

where the elements of $\mathit{a}$ are linearly decreased from 2 to 0 during the search process, and ${\mathit{r}}_1$ and ${\mathit{r}}_2$ are random vectors in the interval [0, 1]. As the position of the optimum (prey) in an abstract search space is not known, the GWO supposes that the $\alpha$ , $\beta$ , and $\delta$ wolves have better knowledge of the location of the prey. Thus, in the GWO, the potential solutions are required to update their positions according to the positions of $\alpha$ , $\beta$ , and $\delta$ , using the formulas given below

$$\begin{array}{c}\hfill {\mathit{D}}_{\alpha }=|{\mathit{C}}_1·{\mathit{X}}_{\alpha }(t)-\mathit{X}(t)|\hfill \\ \multicolumn{1}{c}{{\mathit{D}}_{\beta }=|{\mathit{C}}_2·{\mathit{X}}_{\beta }(t)-\mathit{X}(t)|}\\ \multicolumn{1}{c}{{\mathit{D}}_{\delta }=|{\mathit{C}}_3·{\mathit{X}}_{\delta }(t)-\mathit{X}(t)|}\end{array}$$

(19)

$$\begin{array}{c}\hfill {\mathit{X}}_1={\mathit{X}}_{\alpha }(t)-{\mathit{A}}_1·{\mathit{D}}_{\alpha }\hfill \\ \multicolumn{1}{c}{{\mathit{X}}_2={\mathit{X}}_{\beta }(t)-{\mathit{A}}_2·{\mathit{D}}_{\beta }}\\ \multicolumn{1}{c}{{\mathit{X}}_3={\mathit{X}}_{\delta }(t)-{\mathit{A}}_3·{\mathit{D}}_{\delta }}\end{array}$$

(20)

$$\mathit{X}(t+1)=\frac{{\mathit{X}}_1+{\mathit{X}}_2+{\mathit{X}}_3}{3}$$

(21)

The GWO starts with a random population of gray wolves and is easy to implement because of its simple parameters and simple mechanism.

HMOGWO

This article proposes a HMOGWO algorithm for the fuzzy blocking flow shop scheduling problem. As a multi-objective algorithm, the HMOGWO is developed based on the concept of Pareto optimality. It is described in detail below.

Solution representation

$\mathit{X}=[x(1),x(2),\dots ,x(n)]$ , the position of a gray wolf in the GWO, is an n-dimensional real number vector. In the HMOGWO, the largest position value (LPV) rule ³⁴ is utilized to convert the continuous position of a potential solution into a discrete n-job permutation $\mathbf{\pi }\mathbf{=}\mathbf{[}\mathbf{\pi }\mathbf{(}1\mathbf{)}\mathbf{,}\mathbf{\pi }\mathbf{(}2\mathbf{)}\mathbf{,}\dots \mathbf{,}\mathbf{\pi }\mathbf{(}\mathbf{n}\mathbf{)}\mathbf{]}$ . A simple example is illustrated in Figure 2. Each dimension is represented as a typical job, and the permutation $\mathbf{\pi }$ is generated by sorting the position values of $\mathit{X}$ in non-increasing order. ${\mathit{X}}^o=[x^o(1),x^o(2),\dots ,x^o(n)]$ is the vector of the non-increasingly ordered position values.

Figure 2.

Example of the mapping between the position and the job permutation.

Population initialization

The PFE heuristic ⁶ is an effective heuristic for a blocking flow shop scheduling problem with makespan criterion. To produce initial solutions with a certain quality for the bi-objective scheduling problem considered in this article, a variant of the PFE heuristic denoted by PFE-Pareto and a simple yet effective heuristic termed AAI-Pareto are designed. The PFE-Pareto is used to produce non-dominated solution(s) that are biased toward minimizing the makespan, whereas the AAI-Pareto is for the non-dominated solution(s) that are biased toward maximizing the average agreement index $(\overline{\mathrm{AI}})$ . The procedures of the PFE-Pareto and AAI-Pareto are presented in Figures 3 and 4, respectively. In these figures, $n$ denotes the number of jobs, $m$ represents the number of machines, ${\tilde{p}}_{i,j}$ denotes the fuzzy processing time of job $i$ on machine $j$ , and ${\tilde{D}}_{i,j}$ denotes the fuzzy departure time of job $i$ on machine $j$ .

Figure 3.

Pseudocode of the PFE-Pareto heuristic.

Figure 4.

Pseudocode of the AAI-Pareto heuristic.

Forming a set of non-dominated individuals by gathering all the solutions obtained by these heuristics is an indispensable step in the heuristic process. Subsequently, the member(s) of the set with the largest crowding distance ³⁵ is (are) selected as the heuristic initial solution(s). To ensure the diversity of the initial population, the rest of the individuals are randomly generated, and the position of a solution, $\mathbf{X}=[x(1),x(2),\dots ,x(n)]$ , is produced using the following equation

$$x(k)=x_{max}(k)-(x_{max}(k)-x_{min}(k))\times r,k=1,2,\dots ,n$$

(22)

where $x_{min}(k)=0$ and $x_{max}(k)=10$ for each dimension $k$ , and $r$ is a uniform random number in the interval [0, 1]. Because a heuristic initial solution is a job permutation $\mathbf{\pi }\mathbf{=}\mathbf{[}\mathbf{\pi }\mathbf{(}1\mathbf{)}\mathbf{,}\mathbf{\pi }\mathbf{(}2\mathbf{)}\mathbf{,}\dots \mathbf{,}\mathbf{\pi }\mathbf{(}\mathbf{n}\mathbf{)}\mathbf{]}$ , it needs to be converted into a position vector. Fortunately, this conversion can be easily accomplished using the following equation

$$\begin{array}{c}\multicolumn{1}{c}{x(\pi (k))=x_{max}-\frac{x_{max}-x_{min}}{n}}\\ \times (k-1+r),k=1,2,\dots ,n\hfill \end{array}$$

(23)

where $x_{max}$ and $x_{min}$ are the minimum and maximum values in the original position vector obtained by equation (22), respectively.

Archive and its maintenance

In the HMOGWO, an archive is used to store the non-dominated solutions that are produced during the search process. The predetermined maximum archive size $(S_{\mathrm{AR}})$ limits the number of non-dominated solutions in the archive. The archive is maintained at each iteration, and the procedure of archive maintenance is presented as follows

1. Add every non-dominated solution obtained in the current iteration to the archive and then remove the dominated member(s) from the archive.

2. If the number of the non-dominated individuals in the archive does not exceed the $S_{\mathrm{AR}}$ , then keep the current archive; else, go to (3).

3. Apply a dynamic crowding distance-based maintenance strategy to the archive, as provided below:

$$\mathrm{fitnes}{s}_1=(1-\exp (-\mathrm{CD}){)}^{-\lambda }$$

(24)

where $\mathrm{CD}$ represents the crowding distance of an archive member and $\lambda$ , the deletion selection pressure coefficient, is a positive constant that takes the value of 4 in this article. It can be inferred that the solution with the smallest crowding distance is preferably removed. Since deleting a member will change the crowding distances of its neighbors, only one non-dominated solution is removed at each turn and then the crowding distances are recalculated. The dynamic maintenance strategy not only improves the distribution of non-dominated solutions but also prevents optimal solutions from being removed when removing several individuals located in a crowded region at a time.

Leader selection

The leaders ( $\alpha$ , $\beta$ , and $\delta$ wolves) in the HMOGWO tend to be selected from the archive. The leaders are selected using a fitness value-based roulette wheel, as defined in equation (25)

$$\mathrm{fitnes}{s}_2=(1-\exp (-\mathrm{CD}){)}^{\gamma }$$

(25)

where $\mathrm{CD}$ denotes the crowding distance of an archive member and $\gamma$ , the leader selection pressure coefficient, is a positive constant that is set to 2 in this article. This equation implies that a member possessing a larger crowding distance has a higher probability of becoming a leader. In the leader selection procedure, $\alpha$ is first selected, followed by $\beta$ and $\delta$ . The member selected as a leader will not participate in the next selection. Importantly, three leaders might be selected in certain special cases. If there are fewer than three members in the archive, the other leader(s) are selected from the second non-domination level using the abovementioned roulette wheel. Similarly, if there is only one solution in the second level, and it is $\beta$ , then $\delta$ would be selected from the third level.

Local search strategy

Various studies have shown that incorporating a local search into an evolutionary optimization algorithm improves the exploitation capability of the algorithm. In addition, the insert neighborhood structure is recognized as being superior for yielding neighboring solutions for flow shop scheduling problems. ³⁶ Thus, two local search procedures based on the insert neighborhood structure are incorporated into the HMOGWO. To perform the local search strategy, each gray wolf (potential solution) is assumed to possess four attributes: current position, best previous position, and their corresponding job permutations. A potential solution whose current job permutation does not dominate its best previous job permutation is treated as unimproved. The local search procedure I shown in Figure 5 is implemented on each unimproved potential solution with a certain probability $p_{\mathrm{LS}}^I$ to produce an improved permutation that dominates the best previous one. The local search procedure II is applied to all archive members, that is, the candidate selected leaders, with a certain probability $p_{\mathrm{LS}}^{\mathrm{II}}$ . As illustrated in Figure 6, the local search procedure II can yield one or more improved permutations for each archive member.

Figure 5.

Pseudocode of local search procedure I.

Figure 6.

Pseudocode of local search procedure II.

After the local search strategy is implemented, the position of an improved solution should be adjusted accordingly. Let ${\mathit{X}}^o=[x^o(1),x^o(2),\dots ,x^o(n)]$ be the vector of non-increasingly ordered position values, $\mathbf{\pi }\prime \mathbf{=}\mathbf{[}\mathbf{\pi }\prime \mathbf{(}1\mathbf{)}\mathbf{,}\mathbf{\pi }\prime \mathbf{(}2\mathbf{)}\mathbf{,}\dots \mathbf{,}\mathbf{\pi }\prime \mathbf{(}\mathbf{n}\mathbf{)}\mathbf{]}$ be an improved permutation, and $\mathit{X}\prime =[x\prime (1),x\prime (2),\dots ,x\prime (n)]$ be its corresponding position vector. As illustrated by the example in Figure 7(a), the adjustment can be easily achieved using the following formula

$$x\prime (\pi \prime (k))=x^o(k),k=1,2,\dots ,n$$

(26)

Note that the above adjustment can be used only if $x^o(k)\ne x^o(d)$ , $\forall k,d\in 1,2,\dots ,n\wedge k\ne d$ . However, during the evolution of the HMOGWO, a position vector with some identical values might be yielded. If the adjustment is applied to such a position vector, the permutation implied by $\mathit{X}\prime$ may differ from $\mathit{\pi }\prime$ . To overcome this drawback, a simple revision process ³⁴ for ${\mathit{X}}^o$ is adopted to modify some of its position values. As shown by the example in Figure 7(b), the revision process can be easily performed using the procedure presented in Figure 8.

Figure 7.

(a) Example of the adjustment of $\mathit{X}\prime$ in the local search; (b) example of the revision process for ${\mathit{X}}^o$ with some identical values ( $\epsilon =0.01$ and $n=6$ ).

Figure 8.

Pseudocode of the revision process for ${\mathit{X}}^o$ .

Outline of the proposed HMOGWO

In summary, the HMOGWO utilizes the LPV rule for solution representation, introduces a heuristic process for population initialization, applies a dynamic maintenance strategy for archiving, develops a thorough mechanism for leader selection, and incorporates a local search strategy based on the insert neighborhood structure. The pseudocode of the HMOGWO is provided in Figure 9, where $\mathrm{Ps}$ denotes population scale, $\mathrm{nFEs}$ represents the current number of function evaluations, $\mathrm{nFE}{s}_{max}$ represents the maximum number of function evaluations, and $r$ is a uniform random number in the interval [0, 1].

Computational analysis

In this study, the performance of the HMOGWO is tested on real-time production data of panel block assembly in shipbuilding. A typical assembly line for panel blocks consists of seven main processes: baseplate splicing, baseplate welding, longitudinal assembly, longitudinal welding, girder and floor assembly, girder and floor welding, and checking and carting. Each process is implemented at its corresponding station. The panel blocks are usually of large size. Accordingly, no intermediate buffer between two consecutive stations exists, mainly because of a lack of space. In addition, the processing time of each process is often affected by uncertainty due to both machine and human factors, and the due date takes the form of an interval that is related to the satisfaction degree of the demand side. Therefore, the scheduling problem of panel block construction should be considered as a fuzzy blocking flow shop scheduling problem. The real-time data of the fuzzy processing time and the fuzzy due date of panel blocks to be constructed originate from a large shipyard in Shanghai, China. An example of the data is shown in Table 1. The most plausible value $(p_{i,j})$ of the fuzzy processing time $({\tilde{p}}_{i,j}=(p_{i,j}^O,p_{i,j},p_{i,j}^P))$ is established as the mean value of the historical processing times of the same or very similar panel blocks. The optimistic value $(p_{i,j}^O)$ and the pessimistic value $(p_{i,j}^P)$ are randomly obtained from $[{\delta }_{11}{P}_{i,j},{\delta }_{12}{P}_{i,j}]$ and $[{\delta }_{21}{P}_{i,j},{\delta }_{22}{P}_{i,j}]$ . ¹⁹ In this study ${\delta }_{11}$ , ${\delta }_{12}$ , ${\delta }_{21}$ , and ${\delta }_{22}$ are set to 0.85, 0.90, 1.10 and 1.25, respectively, according to historical data and the advice of experienced engineers. The fuzzy due date $({\tilde{d}}_i)$ of each panel block is provided by the hull block assembly shop, which is the demand side. The dataset used for performance evaluation comprises 40 instances grouped into four sets by $n$ (the number of jobs), where $n=[20,30,40,50]$ . In the performance evaluation experiment, the proposed HMOGWO is compared with two well-known population-based multi-objective evolutionary algorithms, the non-dominated sorting genetic algorithm II (NSGA-II) ³⁵ and the multi-objective particle swarm optimization (MOPSO) (constriction factor version ³⁷ ). This experiment is implemented in MATLAB R2015b on a computer with an Intel Core i5 2.50 GHz CPU and 4 GB RAM.

Table 1.

Example of fuzzy processing time and fuzzy due date.

Processing time/due date (min.)	Panel block Pb1	Panel block Pb2
${\tilde{p}}_{i,1}$	121 139 174	87 97 122
${\tilde{p}}_{i,2}$	206 240 271	150 172 191
${\tilde{p}}_{i,3}$	118 131 161	87 97 110
${\tilde{p}}_{i,4}$	124 146 182	103 114 143
${\tilde{p}}_{i,5}$	143 166 199	106 120 135
${\tilde{p}}_{i,6}$	204 240 300	139 164 200
${\tilde{p}}_{i,7}$	111 125 145	72 80 94
${\tilde{d}}_i$	1000 1300 2000 2500	2000 2200 2700 3000

Performance metrics

To evaluate the performance of multi-objective optimization algorithms, inverted generational distance (IGD), ³⁸ spread, ³⁵ and coverage ³⁹ are used as the metrics.

IGD

IGD measures the distance between each point composing the optimal Pareto front $(P{F}^{*})$ and the obtained front and is defined as

$$\mathrm{IGD}=\frac{\sqrt{\sum _{s=1}^N{d}_s^2}}{N}$$

(27)

where $N$ is the number of points in the $P{F}^{*}$ , and $d_s$ denotes the Euclidean distance between each point in the $P{F}^{*}$ and its nearest neighbor in the obtained front. Because the $P{F}^{*}$ of the scheduling problem is not known, a reference front established by gathering all obtained fronts of all the compared algorithms is used as a replacement. The $\mathrm{IGD}$ indicator is used to evaluate convergence; a front with a lower $\mathrm{IGD}$ value is desirable.

Spread $(\Delta )$

$\Delta$ measures the extent of spread achieved in the obtained front. This indicator is defined as follows

$$\Delta =\frac{d_f+d_l+\sum _{s=1}^{N-1}|d_c-\overline{d}|}{d_f+d_l+(n-1)\overline{d}}$$

(28)

where $N$ is the number of points in the obtained front, $d_c$ represents the Euclidean distance between consecutive points in the obtained front, $\overline{d}$ is the average of these distances, and $d_f$ and $d_l$ are the Euclidean distances of boundary points of the obtained front to the extreme points of the $P{F}^{*}$ . This indicator is 0 for an ideal distribution.

When calculating the $\mathrm{IGD}$ and $\Delta$ indicators, the objective values are normalized into values in the interval [1, 2]. Here, the makespan, which is assumed to be a triangular fuzzy number, is defuzzified using equation (12). Lower $\mathrm{IGD}$ and $\Delta$ indicators correspond with better obtained fronts.

Coverage $(C)$

$C$ is a binary indicator. Let $S_A$ and $S_B$ represent two sets of approximate Pareto optimal solutions that are generated by algorithm A and algorithm B. $C(S_A,S_B)$ , which is defined in equation (29), measures the fraction of the members of $S_B$ that are covered by members of $S_A$ , reflecting the dominance and equivalence relation between the two sets

$$C(S_A,S_B)=\frac{|\{X^B\in S_B|\exists X^A\in S_A,X^A{X}^B}{|S_B|}$$

(29)

The $C$ indicator maps the ordered pair $(S_A,S_B)$ to the interval [0, 1]. $C(S_A,S_B)=1$ indicates that all of the solutions in $S_B$ are covered by individuals in $S_A$ , whereas $C(S_A,S_B)=0$ implies that none of the solutions in $S_B$ are dominated or equaled by the members of $S_A$ .

Parameter tuning

The NSGA-II has three inherent parameters, namely, population scale $(\mathrm{Ps})$ , crossover probability $(p_c)$ , and mutation probability $(p_m)$ . $p_c$ and $p_m$ determine how often the crossover and the mutation operators will be performed, respectively. The MOPSO has five parameters: population scale $(\mathrm{Ps})$ , the maximum archive size $(S_{\mathrm{AR}})$ , constriction factor $(\chi )$ , and two learning factors ( $c_1$ and $c_2$ ), where $S_{\mathrm{AR}}$ is set to one-fourth of $\mathrm{Ps}$ , $\chi =2/|2-\phi -\sqrt{{\phi }^2-4\phi }|$ and $\phi =c_1+c_2$ . $\chi$ is used to alleviate the swarm convergence issue. $c_1$ and $c_2$ determine the influences of personal best and global best positions, respectively. Typically, $c_1$ and $c_2$ are set to 2.05, and $\chi$ is thus 0.7298. ³⁷ The HMOGWO has four parameters: population scale $(\mathrm{Ps})$ , the maximum archive size $(S_{\mathrm{AR}})$ and two local search probabilities ( $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ ), where $S_{\mathrm{AR}}$ is set to one-fourth of $\mathrm{Ps}$ . As population-based multi-objective evolutionary algorithms, the NSGA-II, MOPSO, and HMOGWO have a common parameter, that is, population scale $(\mathrm{Ps})$ . It is set to be 40, 60, 80, and 80 for the instances with 20, 30, 40, and 50 jobs, respectively. The maximum number of function evaluations $(\mathrm{nFE}{s}_{max})$ is used as the stopping criterion, and all the algorithms are executed with $\mathrm{nFE}{s}_{max}=1000\times n$ for a n-job instance.

The performance of an algorithm is significantly affected by the values of parameters. The parameters $p_c$ and $p_m$ of the NSGA-II and $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ of the HMOGWO need to be tuned. Full factorial experiments are conducted to tune these parameters. The considered levels of these parameters are as follows: $p_c=[0.7,0.8,0.9]$ and $p_m=[0.1,0.2,0.3]$ ; $p_{\mathrm{LS}}^I=[0,0.1,0.2,0.3]$ and $p_{\mathrm{LS}}^{\mathrm{II}}=[0,0.7,0.8,0.9,1.0]$ . $p_{\mathrm{LS}}^I=0$ and $p_{\mathrm{LS}}^{\mathrm{II}}=0$ are for investigating the effects of the absences of the two parameters. Different combinations of parameters result in different configurations of an algorithm, each of which is independently run 30 times against all instances of the dataset for parameter tuning. Importantly, the dataset for parameter tuning is different from the dataset for performance evaluation. The former dataset consists of 20 instances divided into four groups by $n$ (the number of jobs), where $n=[20,30,40,50]$ . The $\mathrm{IGD}$ and $\Delta$ indicators are employed to evaluate each parameter combination. The parameter combinations that exhibit better performance are presented in Table 2. Figure 10 depicts the estimated marginal means of the $\mathrm{IGD}$ and $\Delta$ indicators for $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ of the HMOGWO over 20 instances. It can be seen that $p_{\mathrm{LS}}^{\mathrm{II}}=[0.7,0.8,0.9,1.0]$ is always better than $p_{\mathrm{LS}}^{\mathrm{II}}=0$ for both the $\mathrm{IGD}$ and $\Delta$ indicators and that $p_{\mathrm{LS}}^I=[0.1,0.2]$ yields a better $\mathrm{IGD}$ indicator than $p_{\mathrm{LS}}^I=0$ . This observation indicates that the introduction of $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ in the HMOGWO is necessary and rewarding. Considering both the $\mathrm{IGD}$ and $\Delta$ indicators, $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ are established as 0.1 and 0.9, respectively (Figure 10).

Table 2.

Main parameter settings of NSGA-II, MOPSO, and HMOGWO.

$n$	$\mathrm{nFE}{s}_{max}$	NSGA-II			MOPSO					HMOGWO
		$\mathrm{Ps}$	$p_c$	$p_m$	$\mathrm{Ps}$	$S_{\mathrm{AR}}$	$\chi$	$c_1$	$c_2$	$\mathrm{Ps}$	$S_{\mathrm{AR}}$	$p_{\mathrm{LS}}^I$	$p_{\mathrm{LS}}^{\mathrm{II}}$
20	20,000	40	0.8	0.1	40	10	0.7298	2.05	2.05	40	10	0.1	0.9
30	30,000	60	0.8	0.1	60	15	0.7298	2.05	2.05	60	15	0.1	0.9
40	40,000	80	0.8	0.1	80	20	0.7298	2.05	2.05	80	20	0.1	0.9
50	50,000	80	0.8	0.1	80	20	0.7298	2.05	2.05	80	20	0.1	0.9

NSGA-II: non-dominated sorting genetic algorithm II; MOPSO: multi-objective particle swarm optimization; HMOGWO: hybrid multi-objective gray wolf optimization.

Figure 9.

Pseudocode of the HMOGWO.

Figure 10.

Estimated marginal means of the $\mathrm{IGD}$ and $\Delta$ indicators for $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ : (a) IGD and (b) Δ.

Experimental evaluation of HMOGWO

In the performance evaluation experiment, the NSGA-II, MOPSO, and HMOGWO are configured by the parameters listed in Table 2. Each algorithm is independently run 30 times against all 40 instances of the dataset for performance evaluation. The means of the $\mathrm{IGD}$ and $\Delta$ indicators of each algorithm for each instance are computed from 30 replicates. Table 3 summarizes the mean values of the $\mathrm{IGD}$ and $\Delta$ indicators for the 10 instances in each group and indicates that on average, the HMOGWO outperforms the NSGA-II and MOPSO regarding the two indicators. To compute the $C$ indicator of two algorithms for each instance, a two-step procedure is proposed. First, form three sets that are, respectively, denoted as $S_{\mathrm{NSGA}-\mathrm{II}}$ , $S_{\mathrm{MOPSO}}$ , and $S_{\mathrm{HMOGWO}}$ by gathering all the non-dominated solutions produced by the corresponding algorithm in 30 runs. Second, compute $C(S_{\mathrm{HMOGWO}},S_{\mathrm{NSGA}-\mathrm{II}})$ , $C(S_{\mathrm{NSGA}-\mathrm{II}},S_{\mathrm{HMOGWO}})$ , $C(S_{\mathrm{HMOGWO}},S_{\mathrm{MOPSO}})$ , and $C(S_{\mathrm{MOPSO}},S_{\mathrm{HMOGWO}})$ using equation (29). Table 4 lists the mean values of the $C$ indicator for the 10 instances in each group. It can be observed that for each group, few optimal solutions of the HMOGWO are covered by those of the NSGA-II or MOPSO, whereas most or almost all optimal solutions of both the NSGA-II and MOPSO are covered by the individuals of the HMOGWO. This observation demonstrates that the HMOGWO is superior concerning the $C$ indicator compared to the other two algorithms. Figure 11 illustrates the distributions of the solutions in $S_{\mathrm{NSGA}-\mathrm{II}}$ , $S_{\mathrm{MOPSO}}$ , and $S_{\mathrm{HMOGWO}}$ in objective space for four instances with different number of jobs, that is, instance 20–10, 30–6, 40–3, and 50–1. The instance label, for example, 20–10 represents the 10th instance in the group of 20 jobs. The coverage of the optimal solutions obtained by different algorithms can be seen from this figure.

Table 3.

Results of the $\mathrm{IGD}$ and $\Delta$ indicators for the instances grouped by $n$ .

$n$	$\mathrm{IGD}$ , Mean (SD)			$\Delta$ , Mean (SD)
	NSGA-II	MOPSO	HMOGWO	NSGA-II	MOPSO	HMOGWO
20	0.066 (0.015)	0.078 (0.016)	0.020 (0.006)	0.714 (0.045)	0.824 (0.037)	0.473 (0.045)
30	0.094 (0.037)	0.082 (0.028)	0.018 (0.008)	0.857 (0.041)	0.895 (0.038)	0.618 (0.079)
40	0.105 (0.035)	0.087 (0.032)	0.013 (0.003)	0.939 (0.107)	0.958 (0.085)	0.644 (0.111)
50	0.109 (0.012)	0.084 (0.007)	0.012 (0.002)	0.941 (0.013)	0.946 (0.013)	0.641 (0.074)

IGD: inverted generational distance; NSGA-II: non-dominated sorting genetic algorithm II; MOPSO: multi-objective particle swarm optimization; HMOGWO: hybrid multi-objective gray wolf optimization; Bold signifies best values.

Table 4.

Results of the $C$ indicator for the instances grouped by $n$ .

$n$	HMOGWO versus NSGA-II, Mean (SD)		HMOGWO versus MOPSO, Mean (SD)
	$C(S_{\mathrm{HMOGWO}},S_{\mathrm{NSGA}-\mathrm{II}})$	$C(S_{\mathrm{NSGA}-\mathrm{II}},S_{\mathrm{HMOGWO}})$	$C(S_{\mathrm{HMOGWO}},S_{\mathrm{MOPSO}})$	$C(S_{\mathrm{MOPSO}},S_{\mathrm{HMOGWO}})$
20	0.817 (0.143)	0.116 (0.116)	0.808 (0.127)	0.118 (0.093)
30	1.000 (0.000)	0.000 (0.000)	0.961 (0.083)	0.010 (0.021)
40	0.992 (0.026)	0.006 (0.019)	0.971 (0.048)	0.008 (0.013)
50	0.986 (0.045)	0.003 (0.009)	0.958 (0.060)	0.006 (0.013)

HMOGWO: hybrid multi-objective gray wolf optimization; NSGA-II: non-dominated sorting genetic algorithm II; MOPSO: multi-objective particle swarm optimization; Bold signifies best values.

Figure 11.

Distributions of the gathered solutions of each algorithm for four instances: (a) 20–10, (b) 30–6, (c) 40–3, and (d) 50–1.

To test if the outperformance of the HMOGWO is statistically significant, for each group, pairwise comparisons of HMOGWO versus NSGA-II and HMOGWO versus MOPSO are performed by applying the Wilcoxon signed rank test to the $\mathrm{IGD}$ , $\Delta$ , and $C$ indicators for the involved instances. The Wilcoxon test which is a non-parametric test and does not require the assumption of normal distributions has been widely used to analyze the obtained results of evolutionary algorithms.^32,40,41 The function wilcox.test(HMOGWO, NSGA-II or MOPSO, paired = TRUE, alternative = “less”) implemented in R version 3.2.4 is applied to the $\mathrm{IGD}$ and $\Delta$ indicators, and wilcox.test(HMOGWO, NSGA-II or MOPSO, paired = TRUE, alternative = “greater”) is used for the $C$ indicator. The statistical results in Table 5 suggest that for each group, the HMOGWO provides statistically better $\mathrm{IGD}$ , $\Delta$ , and $C$ indicators than the other algorithms. Moreover, Tables 3–5 show that the HMOGWO outperforms the other two algorithms and exhibits robust superiority for different problem sizes. In summary, the results of the $\mathrm{IGD}$ , $\Delta$ , and $C$ indicators demonstrate that the HMOGWO is competitive and effective in solving the fuzzy blocking flow shop scheduling problem in terms of the convergence, spread, and coverage of the optimal solutions.

Table 5.

Wilcoxon signed rank test for the $\mathrm{IGD}$ , $\Delta$ , and $C$ indicators grouped by $n$ .

$n$	$\mathrm{IGD}$ , p-value		$\Delta$ , p-value		$C$ , p-value
	HMOGWO versus NSGA-II	HMOGWO versus MOPSO	HMOGWO versus NSGA-II	HMOGWO versus MOPSO	HMOGWO versus NSGA-II	HMOGWO versus MOPSO
20	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04
30	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04
40	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04
50	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04	9.8e–04

IGD: inverted generational distance; HMOGWO: hybrid multi-objective gray wolf optimization; NSGA-II: non-dominated sorting genetic algorithm II; MOPSO: multi-objective particle swarm optimization.

Effects of the heuristic process and local search strategy

The heuristic process and local search strategy are two important components of the HMOGWO. An experiment is conducted to investigate the effects of them. For simplicity, hereafter $\mathrm{HP}$ denotes heuristic process and $\mathrm{LS}$ represents local search strategy. In this experiment, the levels of $\mathrm{HP}$ and $\mathrm{LS}$ are as follows: $\mathrm{HP}=[\mathrm{without},\mathrm{with}]$ and $\mathrm{LS}=[\mathrm{without},\mathrm{with}]$ . Therefore, four algorithms denoted by MOGWO-B, MOGWO-HP, MOGWO-LS, and HMOGWO (i.e. MOGWO-HP-LS) can be assembled, where MOGWO-B denotes basic MOGWO. The parameters, namely, the maximum number of function evaluations $(\mathrm{nFE}{s}_{max})$ , population scale $(\mathrm{Ps})$ , the maximum archive size $(S_{\mathrm{AR}})$ and two local search probabilities ( $p_{\mathrm{LS}}^I$ and $p_{\mathrm{LS}}^{\mathrm{II}}$ ) follow the values presented in Table 2. Each algorithm is independently run 30 times for each instance of the dataset for performance evaluation. The $\mathrm{IGD}$ and $\Delta$ indicators of each algorithm are computed for performance evaluation. The estimated marginal means of the $\mathrm{IGD}$ and $\Delta$ indicators for $\mathrm{HP}$ and $\mathrm{LS}$ over 40 instances are shown in Figure 12. Clearly, $\mathrm{HP}=\mathrm{with}$ is always better than $\mathrm{HP}=\mathrm{without}$ , and $\mathrm{LS}=\mathrm{with}$ is always better than $\mathrm{LS}=\mathrm{without}$ , demonstrating that the heuristic process and the local search strategy are effective and contribute to improving the performance of the proposed algorithm.

Figure 12.

Estimated marginal means of the $\mathrm{IGD}$ and $\Delta$ indicators for $\mathrm{HP}$ and $\mathrm{LS}$ : (a) IGD and (b) Δ.

Conclusion and future studies

This article formulates a fuzzy blocking flow shop scheduling problem with fuzzy processing time and fuzzy due date and considers simultaneously minimizing the fuzzy makespan and maximizing the average agreement index. The HMOGWO algorithm is proposed to solve this combinational optimization problem. It utilizes the LPV rule to convert the continuous position of potential solutions into discrete job permutations, employs a dynamic maintenance strategy to maintain the archive, and develops a thorough mechanism to select the leaders. In addition, the HMOGWO introduces a novel heuristic process to initialize population and generate a portion of initial solutions with a certain quality, and incorporates a local search strategy to improve its exploitation capability. Computational experiments are conducted on instances of panel block assembly in shipbuilding. Computational comparisons of the proposed HMOGWO with the NSGA-II and MOPSO demonstrate the effectiveness of the HMOGWO in solving the bi-criterion fuzzy blocking flow shop scheduling problem in terms of the convergence, spread, and coverage of the optimal solutions. Moreover, the heuristic process and the local search strategy are verified as being effective and improving the performance of the HMOGWO.

Future studies focus on the following aspects: (1) extending the HMOGWO to more complex fuzzy blocking flow shop scheduling problem considering other restrictions, such as sequence-dependent setup times; (2) investigating the hybridization of the HMOGWO with other algorithms to achieve even better results; and (3) studying the fuzzy blocking flow shop scheduling problem from the view of green manufacturing.