Population-Based Metaheuristic Algorithms for a Hybrid Batch-Continuous Production Scheduling Problem in a Distributed Pharmaceutical Supply Chain

1. Introduction

Pharmaceutical medicines for treating or preventing diseases are divided into two main categories depending upon their manufacturing requirements: small molecules chemically synthesized and biologics derived from biological sources [1]. The demand on the pharmaceutical medicines is continuously increasing due to the impact of the COVID-19 pandemic, growing healthcare, etc. By 2027, spending on global medicine is predicted to reach $1.9 trillion [2]. One of the key concerns of the pharmaceutical manufacturing industry is establishing efficient production plans and scheduling to quickly respond to the demand fluctuations during pandemic periods [1]. To agilely handle the supply chain risk caused by the demand fluctuation,

In pharmaceutical supply chain, the production companies have been shifted their manufacturing plants from centralized to decentralized. Due to the globalization of markets, technological advancements, and other various factors in the pharmaceutical supply chain, it is common to operate multi-production sites in various regions effectively, and integrating these operations is very important [3,4,5]. In a distributed manufacturing (DM) strategy, production orders are decentralized from one large plant to multiple smaller sites to enhance responsiveness and system reliability while reducing production costs and shutdown risks. The COVID-19 pandemic has accelerated the shift toward decentralized pharmaceutical manufacturing systems [6,7]. According to Marques et al [8], DM strategy is an effective way for developing countries, where epidemics occur frequently, to provide flexibility in response to sudden demand for medicines.

In pharmaceutical manufacturing process, a hybrid batch-continuous manufacturing process has recently gained attention to pharmaceutical manufacturing plants as an alternative manufacturing process that enables flexible responses to uncertain pharmaceutical demand [8,9,10,11,12]. Pharmaceutical batch manufacturing simultaneously feeds inputs and produces intermediates or final products at once. In contrast, continuous manufacturing steadily inputs materials and continuously produces products during the operation [9,10]. The appropriate manufacturing varies depending upon the characteristics of the pharmaceutical process [8]. For example, batch manufacturing is effective for long processing times and multi-stage processes such as reaction or refining. Continuous manufacturing is suitable for processes with short processing times and high exothermic reactions, such as crystallization or tablet processing [9]. Therefore, many pharmaceutical manufacturing plants require to equip the hybrid batch–continuous manufacturing process that reflects the characteristics of pharmaceutical processes to enhance both productivity and quality.

In this paper, we study an integrated scheduling problem with heterogeneous manufacturing plants and one distribution center (DC). In the two-stage pharmaceutical supply chain, hybrid batch-continuous flowshop manufacturing processes and a direct-shopping delivery are used in response to the flexible response to uncertain demand. The integrated manufacturing and delivery scheduling problem in the two-stage pharmaceutical supply chain is called PS-HF&DS.

In each plant, an unrelated permutation flow shop layout consisting of a hybrid batch-continuous production line is configured to produce various pharmaceuticals. Each pharmaceutical order is split and produced in a distributed manner across the plants. The pharmaceutical products are directly shipped to DC after production. To minimize the makespan, we should (i) split each pharmaceutical order into sub-orders, (ii) assign the sub-orders to plants, and (iii) determine the permutation sequence of production and direct shipment. We formulate PS-HF&DS as a mathematical model. We propose genetic algorithm (GA) and particle swarm optimization (PSO), and two solution structures, each with a unique decoding process. As a result, we develop four variants, comprising two distinct solution structures for each metaheuristic to solve this model. We conduct numerical experiments to verify and evaluate the proposed metaheuristics’ performance.

The remainder of our study is organized as follows: In Section 2, we provide the literature survey. In Section 3, we describe the addressed problem and the mathematical formulation. In Section 4, we describe the proposed metaheuristics. In Section 5 and Section 6, we show the results of the numerical experiment and sensitivity analysis, respectively. In Section 7, we summarize the conclusion.

2. Literature Survey

The addressed problem is related to the pharmaceutical scheduling and distributed permutation flowshop scheduling (DPFS) problems. We review relevant studies focusing on production systems, manufacturing strategy, scheduling constraints, methodologies, and objective functions. Table 1 summarizes and classifies the reviewed literature.

A real-life pharmaceutical batch scheduling problem formed as a flowshop with multiple processors (FSMP) was addressed in [13]. FSMP considered cleaning, changeover, etc., to minimize total completion time. A hybrid genetic algorithm was proposed and demonstrated robust performance through comparative experiments. The pharmaceutical scheduling problems in a multi-stage batch plant were presented in [14,15]. A MIP-based decomposed heuristic was proposed to solve the scheduling problem and verified by using the benchmark instances of a pharmaceutical facility [14]. To minimize the makespan, they proposed a constraint programming model in [15]. The planning and scheduling of multi-product in the multi-stage continuous plant problem was addressed in [16]. To deal with the addressed problem, three novel MILP models were proposed. A campaign planning and scheduling problem in biopharmaceutical manufacturing was tackled in [11]. The tackled problem considered batch and continuous processes, sequence-dependent changeovers, etc. A mathematical model was developed, and managerial insights were presented by comparison with the previously developed models.

A distributed permutation flowshop scheduling (DPFS) problem, first discussed in [17], has attracted significant attention in the scheduling study. In the last decade, the DPFS problem was extended with variants such as blocking and buffer constraints. The DPFS problem, considering a blocking constraint, was studied in [18,19]. A variant of the iterated greedy algorithm, PBIG, was proposed to minimize total flowtime [18]. To minimize makespan, an evolution strategy (ES) and two mathematical models were proposed in [19]. DPFS with no-wait constraint was addressed in [20]. To minimize the makespan, the general variable neighborhood search (GVNS) algorithm was proposed, and the performance of the proposed GVNS was verified through computational experiments. A distributed permutation heterogeneous flowshop scheduling with order acceptance constraint was studied in [21]. To maximize profit, they proposed GA and PSO. In addition to these studies, DPFS is extended with variants such as setup times constraint [22,23,24], order constraint [25], integration with transportation [26,27], heterogeneous distributed plants [24,28], distributed hybrid flowshop framework [29,30], and distributed assembly flowshop framework [31,32,33]. The variants of DPFS have been studied with various methodologies: mathematical model [19,31], iterated greedy algorithm [18,28], variable neighborhood search [20], evolutionary algorithm [19], fruit fly optimization[23], etc. Additionally, variants of the DPFS problem are reviewed and classified in [5,34].

We study the scheduling model simultaneously considering distributed manufacturing, hybrid batch-continuous manufacturing, order splitting, direct shipping policy, and sequence-dependent changeovers to minimize the makespan in the two-stage pharmaceutical supply chain. Despite the practical advantages and necessity of the proposed integrated scheduling in the pharmaceutical industry, to the best of our knowledge, there is no study for the scheduling model related in the pharmaceutical supply chain due to a problem-complexity. Therefore, the contributions of this study are as follows:

(1)

We firstly address PS-HF&DS and formulate it as a MILP model to minimize the makespan.

(2)

To efficiently solve the large-sized model, we independently apply GA and PSO to two distinct solution representations. As a result, we propose four population-based metaheuristics variants.

(3)

We present managerial insights derived from a sensitivity analysis evaluating the performance of the metaheuristics depending upon the maximum algorithm computational time.

3. Mathematical Model

We formulate the addressed problem as a mixed-integer-linear-programming (MILP) model. We present the description and an illustrative example of the addressed problem in Section 3.1. In Section 3.2, we present the MILP formulation. The sets, parameters, and variable are defined as follows:

Set & Parameter
$DF$	Distributed plant set
$TB$	Batch task set
$TC$	Continuous task set
$T$	Total task set
$I$	Ordered pharmaceutical product set
$L$	Sub-order set related to pharmaceutical order, $L=I$
$ID$	Pharmaceutical products and dummy product set, $\left\{0\right\}\cup I$
$LD$	Sub-orders related to pharmaceutical order and dummy set, $LD=ID$
$W_i$	Order amount of pharmaceutical order$i\in I$
${F1}_i$	Type of pharmaceutical order$i\in I$
${F2}_l$	Type of production sub-orders$l\in L$
${PR}_{ti}^f$	Production rate of task $t\in T$for pharmaceutical order $i\in I$in distributed plant $f\in DF$
$Y_{ti}^f$	Production yield of task $t\in T$for pharmaceutical order $i\in I$in distributed plant $f\in DF$
${CO}_{ij}^f$	Change-over time between pharmaceutical orders $i,j\in ID$in distributed plant $f\in DF$
${DT}_f$	Delivery time between distributed plant $f\in DF$and DC
$M$	Large number
Decision Variable
$z_{fl}$	Equal to 1, if sub-order $l\in L$in distributed plant $f\in DF$is assigned; Otherwise, 0;
${pq}_{tl}^f$	Production amount of sub-order $l\in L$at task $t\in T$in distributed plant $f\in DF$
$q_{fl}$	Production amount of sub-order $l\in L$in distributed plant $f\in DF$
$x_{ij}^f$	Equal to 1, if sub-order $j\in LD$is immediately preceded sub-order $i\in LD$in distributed plant $d\in DF$; Otherwise, 0;
${um}_l^f$	Manufacturing sequence of sub-order $l\in LD$in distributed plant $f\in DF$
${st}_{tl}^f$	Start time of sub-order$l\in LD$at task $t\in T$in distributed plant $f\in DF$
${pt}_{tl}^f$	Productions time of sub-order $l\in LD$at task $t\in T$in distributed plant $f\in DF$
${ct}_{tl}^f$	Completion time of sub-order $l\in LD$at task $t\in T$in distributed plant $f\in DF$
${cm}_l^f$	Manufacturing completion time of sub-order $l\in L$in distributed plant $f\in DF$
${dc}_l^f$	Delivery completion time of sub-order $l\in L$in distributed plant $f\in DF$
$C_i$	Completion of pharmaceutical order $i\in I$
$C_{max}$	Makespan

3.1. Problem Description

We assume that DC has the same number of vehicles as the number of plants to pick up the produced pharmaceutical products from the plants. The detailed descriptions of the addressed PS-HF&DS are as follows:

(1)

Each pharmaceutical order is split into sub-orders, and the production amount of each sub-order is determined in the order splitting phase. The sum of the production amount of the sub-orders must be equal to the order amount of the related order.

(2)

The sub-orders are assigned to plants in the sub-order assigning phase. The processing time of each sub-order differs at each task depending upon the distributed plant, and it is calculated based on [16].

(3)

In the permutation sequencing of production and direct shipment phase, we determine not only the permutation sequence but also the start time and the completion time of each sub-order at each task in the distributed plant. Between two consecutive continuous-task stages, for each sub-order, the next task starts after the previous task starts, and the next task completes after the previous task completes [16].

(4)

The DC picks up the pharmaceutical products according to a direct shipping policy. Since vehicle fleet utilization and routing are not within the scope of this study, the products are directly delivered without allowing order-splitting and combining. The completion of the pharmaceutical order is defined as the latest arrival time at the DC among all related sub-orders.

We illustrate the PS-HF&DS in Figure 1 with three pharmaceutical orders and two distributed plants layouts as one batch task and two consecutive continuous task stages. Figure 1 (a) represents the result of the order splitting and assigning, where bar charts represent the production amount of each sub-order in every plant as a ratio. For bar charts of Figure 1 (a), the Y-axis represents the sub-orders in each plant. Index $O_i^f$ means the sub-order for the pharmaceutical order $i$ produced in the plant $f$. The X-axis represents the amount produced by each plant as a percentage of the total order amount. We show the Gantt charts for plants and DC in Figure 1 (b). The permutation manufacturing sequences of the sub-orders correspond to $\left[O_1^1\to O_2^1\to O_3^1\right]$ and $\left[O_1^2\to O_2^2\right]$ in the D-Plants 1 and 2, respectively. For each distributed plant, since stages C2 and C3 are two consecutive continuous task stages, each sub-order begins task C3 before task C2 is finished. However, each sub-order conducts task C2 only after task B1 is finished. In the shipping stage, each sub-order is directly delivered to the DC after production. The completion time of each order is defined as the latest arrival time of the sub-order related to the order.

3.2. Mixed-Integer-Linear-Programming Model

The objective function and constraints of MILP model can be formulated as follows:

$Minimize{C}_{max}$		(1)
Subject to
$q_{fl}\le M\times z_{fl}$	$\forall f\in DF,l\in L$	(2)
$W_i\le \sum _{f\in DF}{q}_{fl}$	$\forall i\in I,l\in L:{F1}_i={F2}_l$	(3)
$x_{ii}^f=0$	$\forall f\in DF,i\in LD$	(4)
$x_{0i}^f+\sum _{j\in L:i\ne j}{x}_{ji}^f=z_{fl}$	$\forall f\in DF,i\in L$	(5)
$x_{i0}^f+\sum _{j\in L:i\ne j}{x}_{ij}^f=z_{fl}$	$\forall f\in DF,i\in L$	(6)
$\sum _{i\in L}{z}_{fi}\le M\times \sum _{j\in L}{x}_{0j}^f$	$\forall f\in DF$	(7)
$\sum _{i\in L}{z}_{fi}\le M\times \sum _{j\in L}{x}_{j0}^f$	$\forall f\in DF$	(8)
${um}_0^f=0$	$\forall f\in DF$	(9)
${um}_i^f\le M\times z_i^f$	$\forall f\in DF,i\in L$	(10)
${um}_i^f\le \sum _{j\in L}{z}_j^f$	$\forall f\in DF,i\in L$	(11)
${um}_i^f+x_{ij}^f\le {um}_j^f+M\times \left(1-x_{ij}^f\right)$	$\forall f\in DF,i\in LD,j\in L:i\ne j$	(12)
${pq}_{tl}^f={pt}_{tl}^f\times {PR}_{ti}^f$	$\forall f\in DF,i\in I,l\in L:{F1}_i={F2}_l,t\in T$	(13)
${pq}_{tl}^f\times Y_{tl}^f={pq}_{t+1,l}^f$	$\forall f\in DF,i\in I,l\in L:{F1}_i={F2}_l,t\in T:T\ne \left|T\right|$	(14)
${pq}_{tl}^f\times Y_{ti}^f=q_{fl}$	$\forall f\in DF,i\in I,l\in L:{F1}_i={F2}_l,t\in T:T=\left|T\right|$	(15)
${st}_{tl}^f\le M\times z_{fl}$	$\forall f\in DF,l\in L,t\in T$	(16)
${ct}_{tl}^f\le M\times z_{fl}$	$\forall f\in DF,l\in L,t\in T$	(17)
${st}_{tl}^f+{pt}_{tl}^f={ct}_{tl}^f$	$\forall f\in DF,l\in L,t\in T$	(18)
${ct}_{ti}^f+{CO}_{ij}^f\le {st}_{tj}^f+M\times \left(1-x_{ij}^f\right)$	$\forall f\in DF,i,j\in L:i\ne j,t\in T$	(19)
${CO}_{ti}^f\le {st}_{ti}^f+M\times \left(1-x_{0i}^f\right)$	$\forall f\in DF,i\in L:i\ne j,t\in T$	(20)
${ct}_{tl}^f\le {st}_{t+1,l}^f+M\times \left(1-z_{fl}\right)$	$\forall f\in DF,l\in L,t\in TB$	(21)
${st}_{il}^f\le {st}_{jl}^f+M\times \left(1-z_{fl}\right)$	$\forall f\in DF,l\in L:i,j\in TC:i<j$	(22)
${ct}_{il}^f\le {ct}_{jl}^f+M\times \left(1-z_{fl}\right)$	$\forall f\in DF,l\in L:i,j\in TC:i<j$	(23)
${ct}_{tl}^f\le {cm}_l^f+M\times \left(1-z_{fl}\right)$	$\forall f\in DF,l\in L,t\in T:t=\left|T\right|$	(24)
${cm}_l^f+{DT}_f\le {dc}_l^f+M\times \left(1-z_{fl}\right)$	$\forall f\in DF,l\in L$	(25)
${dc}_i^f+2\times {DT}_f\le {dc}_j^f+M\times \left(1-x_{ij}^f\right)$	$\forall f\in DF,i,j\in L:i\ne j$	(26)
${cm}_l^f+{DT}_f\le {dc}_l^f+M\times \left(1-x_{0l}^f\right)$	$\forall f\in DF,l\in L$	(27)
${dc}_l^f\le C_i$	$\forall f\in DF,i\in I,l\in L:{F1}_i={F2}_l$	(28)
$C_i\le C_{max}$	$\forall i\in I$	(29)
$z_{fl},x_{ij}^f\in \left\{0,1\right\}$	$\forall f\in DF,l\in L,j\in LD,j\in LD$	(30)
${pq}_{tl}^f,q_{fl},{cm}_l^f,{dc}_l^f\ge 0$	$\forall f\in DF,l\in L,t\in T$	(31)
${st}_{tl}^f,{pt}_{tl}^f,{ct}_{tl}^f\ge 0$	$\forall l\in LD,t\in T,f\in DF$	(32)
$C_i\ge 0$	$\forall i\in I$	(33)

Constraint (1) represents the objective function of the proposed PS-HF&DS to minimize the makespan. Constraints (2) and (3) determine the production amount of each sub-order in every plant. The permutation manufacturing sequence in each distributed plant is determined by constraints (4) to (12). The production sequence of the sub-orders $i$ and $j$ in the distributed plant $f$ are determined by in Constraints (5) and (6). Constraints (7) - (8) ensure that if there is any sub-order $i$ assigned to the distributed plant $f$, the first and last manufacturing sequence order must exist in the distributed plant $f$. Constraints (9) – (12) prevent sub-tours in distributed plant $f$. Constraints (13) – (24) represent the permutation sequencing phase. Constraints (13) – (15) calculate the production time of the sub-order $l$ in each task and its production amount based on [35]. Constraints (16) – (18) represent the relationship among start, production, and completion times of sub-order $l$. Constraints (19) – (20) calculate start and completion times of sub-orders $i$ and $j$ within a task. Constraint (21) calculates start and completion times between consecutive batch tasks $t$ and $t+1$. Constraints (22) and (23) calculate start and completion times of sub-order $l$ between continuous tasks $t$ and $t+1$. The manufacturing completion time of sub-order $l$ in distributed plant $f$ is calculated by constraint (24). Constraints (25) – (29) represent the shipping stage. Constraints (25) and (26) calculate the delivery completion time of sub-order $l$ in distributed plant $f$. For each plant $f$, the first manufactured sub-order $l$ is delivered and calculated the delivery completion time by Constraint (27). The completion time of pharmaceutical sub-order $i$ and makespan are calculated by constraint (28) and (29), respectively. Constraints (30) – (33) are binary and non-negative constraints for decision variables.

4. Metaheuristics

Since DPFS problem has been proven NP-hard problem [25], the large-sized proposed problem is not able to be solved within a limited time by the MILP model. Hence, we propose genetic algorithm (GA) and particle swarm optimization (PSO) to solve and find the optimal solutions of the reasonable computational time. GA and PSO have demonstrated good performance in related studies [36,37,38,39,40]. We introduce two distinct solution structures with unique decoding processes. The first structure determines order splitting, assignment, and permutation sequence through solution structure in the decoding process, and we refer to a metaheuristic using this solution structure as *_OFP. The second structure one also determines order splitting and permutation sequence by solution structure. However, the plant assignment is determined via a completion time-based assignment heuristic (CAH) in the decoding process. We refer to a metaheuristic using this solution structure as *_OP-CAH. Consequently, we develop four variants, GA_OFP, GA_OP-CAH, PSO_OFP, and PSO_OP-CAH, and these four variants are classified in Figure 2.

4.1. Solution Structure and Decoding PROCESS for *_OFP

*_OFP uses three one-dimensional arrays as an encoded solution. Arrays 1 and 3 are represented as decimals generated by following a uniform distribution $U\left(0,1\right)$ with $\left|DF\right|\times \left|I\right|$ length. Array 2 is represented as an integer generated by following a discrete uniform distribution $unif\left\{1,\left|DF\right|\right\}$ with $\left|DF\right|\times \left|I\right|$ length. An example of an encoded solution is illustrated in Figure 3. Figure 3 (a) shows an encoded solution using the data from Figure 1. Arrays 1, 2, and 3 correspond to [0.210, 0.030, 0.042, 0.0148, 0.563, 0.317], [1,2,2,1,1,1], and [0.071, 0.384, 0.859, 0.413, 0.802, 0.273], respectively.

Array 1 splits each pharmaceutical order into $\left|DF\right|$ virtual sub-orders. It also calculates the production amount of each virtual sub-order by normalizing the decimal elements of Array 1 related to order $i$. Each virtual sub-order is assigned to a plant based on the integer elements of Array 2. The actual sub-orders and production amount of the sub-order related to order $i$ are determined for each plant via decoding Arrays 1 and 2. Figure 3 (b) illustrates the process of the order splitting using Array1. In Figure 3 (b), Array 1 is decoded as the production percentage array. Array 1 [0.210, 0.030, 0.042, 0.148, 0.563, 0.317] becomes $\left[{\displaystyle \frac{0.210}{0.210+0.030}},{\displaystyle \frac{0.030}{0.210+0.030}},{\displaystyle \frac{0.148}{0.148+0.042}},{\displaystyle \frac{0.042}{0.148+0.042}},{\displaystyle \frac{0.563}{0.563+0.317}},{\displaystyle \frac{0.317}{0.563+0.317}}\right]$ as [87.5 (%), 12.5 (%), 77.8 (%), 22.2 (%), 64.0 (%), 36.0 (%)]. Since the elements 0.210 and 0.030 are related to pharmaceutical order 1, these elements are normalized as ${\displaystyle \frac{0.210}{0.210+0.030}}$ and ${\displaystyle \frac{0.030}{0.210+0.030}}$, respectively. Figure 3 (c) illustrates assigning each virtual sub-order into plant $f$ using Array 2. For pharmaceutical order $3$, since every virtual sub-order is assigned to the distributed plant 1, the production amount of the actual sub-order related to order 3 in the distributed plant 1, $q_{13}$, is represented as 100.0 (%).

Array 3 indicates the production and direct shipping sequence of the sub-order $O_i^f$ by sorting the elements of Array 3 in increasing order. Each sub-order $O_i^f$ is produced at the earliest possible start time ${st}_{tl}^f$ in each task $t$ in distributed plant $f$. Figure 3 (d) illustrates the production and direct shipping sequence of each sub-order within a distributed plant corresponding to the Gantt chart in Figure 1 (a).

4.2. Solution Structure and Decoding Process for *_OP-CAH

*_OP-CAH uses two one-dimensional arrays as an encoded solution. Arrays 1 and 2 are represented as decimals generated following a uniform distribution $U\left(0,1\right)$ with $\left|DF\right|\times \left|I\right|$ length. The production amount of each virtual sub-order is determined by decoding Array 1 as illustrated in Figure 3 (b). Array 2 indicates the assigning sequence of the virtual sub-orders by sorting its elements in increasing order. We determine the sub-order assignment by using a completion time-based assignment heuristic (CAH). In CAH, the current dispatching sub-order is virtually assigned to all plants and then assigned into the plant with the minimum delivery completion time. The pseudocode of CAH is as follows:

4.3. Genetic Algorithm

GA utilizes a set of encoded solutions called chromosomes as a population. It consists of several processes: initialization, decoding, evaluation, selection, and crossover and mutation. A population is initially created using distribution $U\left(0,1\right)$ and $unif\left\{1,\left|DF\right|\right\}$ in the initialization process. In the decoding process, to obtain $C_{max}$, each chromosome is decoded by following Section 4.1 and Section 4.2. The fitness of each chromosome $i$, $F\left(i\right)$ is calculated by Eq. (34) in the evaluation process.

$$F\left(i\right)={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$C_{max}^i$}\right.}$$

(34)

, where $C_{max}^i$ denotes $C_{max}$ obtained by chromosome $i$. In the selection process, we adopt a selection mechanism, the roulette wheel method.. For the crossover operation, a one-cut crossover is conducted. For the mutation operation, Arrays 1 and 2 are tweaked by uniform and swap mutations, respectively. Figure 4 indicates the flow charts of the proposed GA.

.4.4. Particle Swarm Optimization

PSO utilizes a set of encoded solutions called particles as a swarm. It consists of several processes: initialization, decoding, evaluation, best position update, and velocity and position update. The initialization, decoding, and evaluation processes of PSO are identical with GA. In the best position update process, the personal and global best position vectors are updated by following Eqs.(35) - (36)

$$P_p^t\leftarrow \left\{\begin{array}{cc}{X}_p^t,& F\left(X_p^t\right)>F\left(P_p^t\right)\\ P_p^t,& otherwise\end{array}\right.$$

(35)

$$G^t\leftarrow \left\{\begin{array}{cc}{X}_p^t,& F\left(X_p^t\right)>F\left(G^t\right)\\ G^t,& otherwise\end{array}\right.$$

(36)

, where $X_p$, $P_p$, and $G$ represent position and personal best position vectors of particle $p$ and global best position vector in iteration $t$, respectively. Each array of the position and velocity vectors are independently updated by following Eqs. (37) – (39) in the velocity and position update processes, respectively. The vectors represented as $U\left(0,1\right)$ are updated by following Eqs. (37) - (38)

$$V_p^{t+1}\leftarrow w\times V_p^t+c_1\times U\left(\mathrm{0,1}\right)\times \left(P_p^t-X_p^t\right)+c_2\times U\left(\mathrm{0,1}\right)\times \left(G^t-X_p^t\right)$$

(37)

$$X_p^{t+1}\leftarrow X_p^t+V_p^{t+1}$$

(38)

, where $w$, $c_1$, and$c_2$ represent weight, cognitive coefficient and social coefficient, respectively. $V_p^{t+1}$ represents velocity vector in iteration $t$. The position vectors represented as $unif\left\{1,n_{DF}\right\}$ are updated by following Eq. (39) referring pan et al. [41]

$$X_p^{t+1}\leftarrow c_2⨂F_3\left(c_1⨂F_2\left(w⨂F_1\left(X_p^t\right),P_p^t\right),G^t\right)$$

(39)

, where $F_1$, $F_2$, and $F_3$ are mutation and crossover operators [41]. The pseudocode of Eq. (39) is as follows: , where $mutation$ and $crossover$ mean swap mutation and one-cut crossover, respectively.

Figure 5 indicates the flow charts of the proposed PSO.

5. Experiments

We conduct numerical experiments to validate the proposed metaheuristics’ performance. The complexity of PS-HF&DS is affected by the number of virtual sub-orders, which is the interaction between the number of ordered pharmaceutical products and distributed plants, $\left|DF\right|\times \left|I\right|$. Therefore, depending upon the size of $\left|DF\right|$ and $\left|I\right|$, we generate two types of instances, which are small-sized and large-sized. $\left|DF\right|$ is $\left\{2,3\right\}$, and $\left|I\right|$ is $\left\{6,8,10\right\}$ in the small-sized instances. $\left|DF\right|$ is $\left\{4,5,6\right\}$, and $\left|I\right|$ is $\left\{40,50,60\right\}$ in the large-sized instances.

To reflect the features of the pharmaceutical supply chain, we generate the experimental parameters by referring to other related studies. For the manufacturing process layouts, we adopt three flowshop layouts as $\left\{\left(C,C\right),\left(B,C,C\right),\left(B,C,C,B\right)\right\}$ presented in [11], in which $B$and $C$ represent batch and continuous tasks, respectively. Since no study simultaneously considers sequence-dependent change-over time, batch and continuous tasks, and distributed manufacturing, we randomly generate change-over time, production rate, and production yield by referring to pharmaceutical scheduling studies [11,14,16]. In addition, we determine weight of pharmaceutical order, $W_i$, by discrete uniform distribution $unif\left\{\left|DF\right|,3\times \left|DF\right|\right\}$ and the delivery time between distributed plant and DC as Manhattan distance.

For every instance, we replicate the proposed metaheuristics 30 times with the maximum computational time as $\left|DF\right|\times \left|I\right|\times 0.5$. The performance measures of small and large size instances are absolute and relative percentage deviation ($ARPD$) as follows:

$$A{RPD}_i\left(\mathrm{\%}\right)={\displaystyle \frac{best-{alg}_i}{best}}\times 100,$$

(40)

, where ${alg}_i$ and $best$ are solution obtained by metaheuristic $i$ and the best solution, respectively. Within the limited time, if the MILP model finds the optimal solutions of the small-sized instances, $best$ means optimal solution. MILP model and metaheuristics are implemented by IBM ILOG CPLEX 22.1.1 and Python IDE, respectively. We use experimental environment as an Intel Core i7-13700KF 3.40GHz CPU.

5.1. Parameter Calibration

We calibrate metaheuristic parameters by using the two-phase calibration method [42,43]. For GA_OFP and GA_OP-CAH, the popsize $\left(P\right)$ and crossover and mutation rates $\left(CR,MR\right)$ are calibrated. For PSO_OFP and PSO_OP-CAH, the swarm size $\left(S\right)$, weight $\left(w\right)$, and cognitive and social coefficients $\left(c_1,c_2\right)$ are calibrated. Every metaheuristic metaheuristic use the stopping criterion as a maximum CPU time,$\left|DF\right|\times \left|I\right|\times 0.5$ seconds. Especially, we set the parameters $P$ and $S$ as $\left|DF\right|\times \left|I\right|\times x$, and we determine $x$ in the two-phase calibration. To calibrate the parameters, three experimental instances are randomly generated and run five replications.

Table 2 shows the result of the parameter calibration with the best parameters of each metaheuristic in bold. For each parameter, if the best value is the largest or smallest candidate in the first calibration phase, we conduct the second phase to finely tune the parameter. For the GA_OFP, $P$, $CR$, and $MR$ are 1.00, 0.70, and 0.15. For the GA_OP-CAH, $P$, $CR$, and $MR$ are 1.00, 0.40, and 0.15. For PSO_OFP, $S$, $w$, $c_1$, and $c_2$ are 0.75, 0.45, 0.20, and 0.30. For PSO_OP-CAH, $S$, $w$, $c_1$, and $c_2$ are 1.00, 0.40, 0.15, and 0.50. We set the termination condition of every metaheuristic as the maximum computational time, $\left|DF\right|\times \left|I\right|\times 0.5$ seconds.

5.2. Experiment Results

In Table 3, the results of the small-sized experiments are shown with the median of and computational time (CPU). We validate the proposed four metaheuristics with the MILP model described in Section 3.2. We limit the maximum computational time of the MILP model to 1,800 seconds. If the MILP model does not solve the instance within 1,800 seconds, we mark the OBJ and CPU time as the best feasible solution obtained by CPLEX IBM ILOG and 1,800++, respectively. Otherwise, we mark the OBJ as the optimal solution with a star symbol. The averages of median of GA_OFP, PSO_OFP, GA_OP-CAH, and PSO_OP-CAH are 3.43, 4.72, 1.94, and 3.15, respectively. Every proposed metaheuristic shows the average of median within 5.0%.

In Table 4, the results of the large-sized experiments are shown with the median of . We compare the proposed metaheuristics by using the . The averages of median of GA_OFP, PSO_OFP, GA_OP-CAH, and PSO_OP-CAH are 25.87, 16.93, 3.16, and 2.93, respectively. *_OP-CAHs show a lower average of median than *_OFPs. Within the OP-CAH and *_OFPs, PSO also shows a lower average of median than GA.

We conduct the ANOVA test with a 0.95 confidence level and present the results as the interval plot in Figure 6. There are no overlaps with the intervals of GA_OFP and PSO_OFP. For the intervals of GA_OP-CAH and PSO_OP-CAH, these do not overlap with the intervals of GA_OFP and PSO_OFP. However, the intervals of GA_OP-CAH and PSO_OP-CAH overlap. It means that *_OFP and *_OP-CAH are statistically significant, and GA and PSO show no statistical significance within the *_OP-CAH. To analyze the robustness of the proposed metaheuristics depending upon the changes in the experimental parameters, we present the interval plots with a 0.95 confidence level for the flowshop layout, $\left|DF\right|$, and $\left|I\right|$ in Figure 7. In Graphs (a) to (c), GA_OFP and PSO_OFP exhibit significant fluctuation in the mean value and confidence interval width depending upon the experimental parameter changes. Especially, the performance of the *_OFP decreases with a decrease in the number of batch tasks and increases in and $\left|I\right|$. In contrast, GA_OP-CAH and PSO_OP-CAH show robust mean values and tight confidence intervals depending upon the experimental parameter changes. Since *_OP-CAH determines the sub-order assignment with the CAH, unlike *_OFP, *_OP-CAH shows robust performance despite the experimental parameter changes.

6. Sensitivity Analysis

We conduct sensitivity analyses to observe the difference in performance gap between *_OFP and *_OP-CAH, and suggest managerial insights for operational strategies to a decision maker. In this section, we analyze the performance gap depending on the maximum CPU time for identical large-sized instances.

To observe the difference in performance gap between *_OFP and *_OP-CAH, we propose a sensitivity analysis by using the identical instances and the experimental environment in the large-sized experiments. We repeat PSO_OFP and PSO_OP-CAH 10 times for each instance with maximum CPU time as $\left|DF\right|\times \left|I\right|\times t$ seconds by using a time parameter $t$. Figure 8 presents an interval plot of the sensitivity analysis by using $ARPD$. PSO_OP-CAH shows the interval plots with no significant difference in the mean value and confidence intervals as the $t$ values increase. In contrast, PSO_OFP shows the interval plots with exponentially decreasing mean value and confidence intervals as increasing the $t$ values increase. Especially, PSO_OFP shows non-overlap with other interval plots and the lowest mean $ARPD$ when the $t$ value exceeds 4.0.

The performance difference in the $ARPD$ between *_OFP and *_OP-CAH is highly related to the maximum CPU time. *_OFP globally searches for the order splitting, assignment, and permutation sequence by the encoded solution, unlike *_OP-CAH. This implies that *_OFP explores a more spacious solution space than *_OP-CAH. This also implies that *_OP-CAH more effectively exploits the solution space within the limited maximum CPU time than *_OFP.

We additionally analyze the relationship with the makespan and distributed plant utilization by using the results of the sensitivity analysis in this section. The average plant utilization rate at $t$ (${AUR}^t$) of each distributed plant is calculated as follows:

$${UR}_f^t\left(\mathrm{\%}\right)={\displaystyle \frac{\underset{f\in \mathit{DF}}{\max }{cm}_f^t-{cm}_f^t}{\underset{f\in \mathit{DF}}{\max }{cm}_f^t}}\times 100$$

(41)

$${AUR}^t\left(\mathrm{\%}\right)=\overline{{UR}_f^t}$$

(42)

, where ${cm}_f^t$ and ${UR}_f^t$ mean manufacturing completion time and utilization rate of distributed plant $f$ at time parameter $t$, respectively. The X-axis means the $ARPD$ value (RPD), and the Y-axis means the $AUR$ value (AUR). In Figure 9, there are four scatter plots, in which time parameters $t$ are 0.5, 2.0, 3.0, and 5.0, respectively.

For PSO_OP-CAH, as shown in scatter plots (a) – (d) of Figure 9, the cluster becomes denser as the $t$ value increases. However, the cluster does not dynamically move in the direction of smaller RPD and AUR. PSO_OFP forms a cluster with a lower density than PSO_OP-CAH when the $t$ value is 0.5, as shown in the scatter plot (a) of Figure 9. It also forms a denser cluster as the t value increases, and the cluster dynamically moves in the direction of smaller RPD and AUR, as shown in the scatter plots (a)-(d) of Figure 9.

PSO_OP-CAH assigns the sub-order to the plant by focusing on minimizing the completion time of sub-orders via CAH. PSO_OP-CAH shows no noticeable difference in AUR and RPD despite the t-value increasing, unlike PSO_OFP. It implies the relationship between the makespan and plant utilization in the distributed manufacturing environment. Thus, the decision maker has to schedule the distributed plants by considering the plant utilization to minimize the makespan.

From the results of the sensitivity analysis, we propose operational strategies to a decision maker depending upon the decision-making time. If the decision maker has sufficient time for decision-making, we recommend determining the main decision contents of the PS-HF&DS by using the *_OFP to globally search the solution space. However, if the decision maker should decide within a limited time, we recommend using the *_OP-CAH to rapidly converge on a good solution.

7. Conclusion

In this paper, we address PS-HF&DS considering multiple unrelated plants, a single distinct permutation flowshop-line, hybrid batch-continuous manufacturing, sequence-dependent changeover time, and direct shipping policy within the pharmaceutical supply chain. We formulate the PS-HF&DS as the MILP model and implement *_OFP and *_OP-CAH by applying GA and PSO to effectively and efficiently solve the model. From the numerical experiments, we demonstrate that *_OP-CAH more effectively finds the near-optimal solutions than *_OFPs within a limited computational time. We additionally conduct the sensitivity analysis of the computational time to analyze the performance gap depending upon the time by comparing PSO_OFP with PSO_OP-CAH. As the computational time increases, *_OFP shows exponentially enhanced performances. Based on the sensitivity analysis result, we propose the operational strategies related to the decision-making time and plant utilization.

Two future research should extend our problem structure. First, we should extend PS-HF&DS with a hybrid flowshop layout with a dynamic order environment. We also should have an extended study of PS-HF&DS related to the routing problem with multiple vehicle-types in the delivery-stage