Optimization of Sustainable Supply Chain Network for Perishable Products

Abstract

In today’s perishable products industry, the importance of sustainability as a critical consideration has significantly increased. This study focuses on the design of a sustainable perishable product supply chain network (SPPSCN), considering the factors of economics cost, environmental impacts, and social responsibility. The proposed model is a comprehensive production–location–inventory problem optimization framework that addresses multiple objectives, echelons, products, and periods. To solve this complex problem, we introduce three hybrid metaheuristic algorithms: bat algorithm (BA), shuffled frog leaping algorithm (SFLA), and cuckoo search (CS) algorithm, all hybrid with variable neighbourhood search (VNS). Sensitivity to input parameters is accounted for using the Taguchi method to tune these parameters. Additionally, we evaluate and compare these approaches among themselves and benchmark their results against a reference method, a hybrid genetic algorithm (GA) with VNS. The quality of the Pareto frontier is evaluated by six metrics for test problems. The results highlight the superior performance of the bat algorithm with variable neighbourhood search. Furthermore, a sensitivity analysis is conducted to evaluate the impact of key model parameters on the optimal objectives. It is observed that an increase in demand has a nearly linear effect on the corresponding objectives. Moreover, the impact of extending raw material shelf life and product shelf life on these objectives is limited to a certain range. Beyond a certain threshold, the influence becomes insignificant.

1. Introduction

Perishable products refer to goods whose quality or quantity gradually decays or deteriorates with time (e.g., pharmaceutical products, dairy products, fruits, vegetables, flowers, and blood products) [1,2]. These products have attracted increasing attention from academics and industries in recent years for their fundamental importance in many industries and trade environments and close relations to human lives [3,4,5].

Globally, about US $1 trillion worth of food is either wasted or lost annually, accounting for nearly one-third of the food produced. Europe can incur losses estimated at 30% due to fruit and vegetable spoilage and damage during production, storage, transportation, and selling periods. The situation in the US is almost the same, as approximately 60% of food waste is caused by inefficiencies in food supply chains [6,7]. These statistics call for the need to improve perishable product supply chains. Furthermore, strong economic turmoil and fierce market competition drive companies to design cost-effective, efficient, and responsive supply chain networks (SCNs) [8]. Perishable product producers and traders are under tremendous pressure to find efficient coordination solutions among demand, production, and supply.

Despite the extensive coverage and comprehensiveness of existing optimization-based supply chain management (SCM) literature, many of these are not directly applicable to perishable products [6]. Designers must make a wide range of SCM decisions at different levels (strategic, tactical, and operational) while considering demand uncertainties and longevity fluctuations in products.

Over the past two decades, researchers have used a sequential method to address strategic and tactical decisions in supply chains. This method separates long-term strategic decisions, such as facility location, from short-term tactical decisions like production and inventory decisions. However, research has demonstrated that this sequential approach results in suboptimal supply chain design and management. Consequently, the joint location–inventory problem was proposed [8,9,10,11,12]. This study extends the traditional location–inventory optimization framework by introducing a combined production–location–inventory problem for perishable products’ supply chain networks. It incorporates demand uncertainty and the limited shelf life of raw materials and products, thereby making the problem more realistic.

Sustainability encompasses three pillars: economic, environmental, and social effects [13,14]. Perishable product supply chains, due to the need for stringent temperature and humidity control during storage and transportation, contribute substantially to carbon emissions [15,16]. Approximately 80% of greenhouse gases and 36% of carbon emissions are direct results of supply chain activities [17,18]. Employing additional workers and managers leads to increased total costs for companies. Sustainable supply chains should maximize the availability of job opportunities to enhance social benefits [19]. These factors motivate SCM policy assessment by including economic, environmental, and social responsibility aspects and redesigning supply chains [20,21]. SC managers are provided with effective decision support tools through supply chain mathematical modelling, which integrates carbon emissions and social responsibility considerations. Therefore, we consider sustainability factors, such as economic, environmental, and social impacts, for perishable product supply chains.

Section 2 reviews the existing literature, while Section 3 introduces problem definitions and mathematical models for the SPPSCN. Solution technique is discussed in Section 4, followed by computational experiments in Section 5. Finally, Section 6 presents the conclusions.

2. Literature Review

Our literature review is divided into two main parts. The first part covers studies that have utilized integrated decision models to design supply chain designs of perishable products. The second part focuses on the challenges and problems associated with the design of SPPSCN.

2.1. Integrated Decision Models of Perishable Products Supply Chain Design

Previous research paid attention to inventory control in perishable product supply chains [22,23,24,25]. Meanwhile, current studies have become increasingly interested in integrated decision making in perishable product supply chains due to problem complexity. Research on the location–inventory problem for perishable products’ supply chain networks is relatively new. The location–inventory problem concerns integrated decisions, including facility location and inventory [9,11,26]. Dai et al. [27] put forward two heuristics to cope with the SCN location–inventory model of perishable products with fuzzy capacity and emission constraints. Based on a real case, Hosseini-Motlagh et al. [28] proposed a bi-objective two-stage stochastic programming model to address the design of a blood SCN. The model incorporates location-allocation and inventory management decisions. Combining production and inventory decisions in the supply chain designs of perishable products is another research line [29,30,31,32]. Dolgui et al. [33] provided a production–inventory–distribution model in a multistage perishable product supply chain. They used a modified genetic algorithm to solve the problem involved in the model. As inventory decisions play crucial roles in perishable products’ supply chain designs, studies that only account for location and production decisions are scarce. Carrizosa et al. [34] embedded a production policy in a facility location–allocation problem in which a single perishable product is considered. Through an empirical analysis, they revealed that the production policy and other factors can affect SCND.

Some scholars have recently paid attention to integrated location, inventory, and production decisions in perishable products’ supply chain designs. Motivated by the real world, Aazami and Saidi-Mehrabad [35] developed a novel production–distribution problem involving multiple periods for perishable products’ supply chain design with three levels. To solve this complex problem, they employed a hierarchical heuristic technique derived from the Benders decomposition algorithm and genetic algorithm. Hamdan and Diabat [36] introduced a two-stage stochastic programming model, incorporating three objectives. The model involved optimal location, production, and inventory decisions for red blood cells and accounted for blood type substitution, perishability, and demand and supply uncertainty. The problem was solved using the epsilon (ε)-constraint method.

On the basis of the above discussions, no model which integrates production, location, and inventory decisions in perishable products’ supply chain designs has been developed in a systematic way. We thus propose a comprehensive model, which considers these three elements, to extend existing research.

2.2. Sustainable Perishable Products’ Supply Chain Design

SPPSCN design mainly considers economy and environment decisions. The commonly used methods in investigating carbon emissions are divided into three categories. First, carbon emissions are regarded as certain constraints [27]. Second, carbon emission consideration is transformed into a part of a cost term in an economic objective [37]. Third, some researchers involve carbon emissions to construct models with two to three objectives to handle problems [38]. Our study adopts the third category.

Perishable product supply chains gain increasing interest because the sustainable development concept is gradually being supported by society [15,39,40]. Liu, Zhu, Xu, Lu, and Fan [16] explored SPPSCN problems in emerging markets. They established a functional model that integrates location, inventory, and routing factors to minimize economic costs and carbon emissions while maximizing product freshness. Biuki, Kazemi, and Alinezhad [3] introduced a two-phase approach for designing an SPPSCN. The proposed methodology employed the genetic algorithm and particle swarm optimization to effectively solve the model.

We now discuss the most similar studies to our work, including similar characteristics and major differences. Goodarzian, Taleizadeh, Ghasemi, and Abraham [13] proposed a mathematical model for the sustainable medical supply chain network (SCN), addressing the production–distribution–inventory–allocation–location problem across multiple products and periods. Sustainability has three pillars, which were regarded as one of their research contributions. They suggested three hybrid metaheuristic algorithms, ant colony optimization, fish swarm algorithm, and firefly algorithm, to effectively solve the proposed model. These algorithms were hybridized through variable neighbourhood search. Our work is different from Goodarzian, Taleizadeh, Ghasemi, and Abraham [13] in four ways. First, we include raw material procurement decisions in the supply chain and account for raw material procurement costs and related inventories in manufacturers. Second, we consider the first-in-first-out (FIFO) policy and deterioration cost for raw materials and products to alleviate waste. Third, we model three sustainable development aspects into corresponding objective functions and convert environmental impacts into fraction costs for operations. Fourth, due to the complexity of the problem, we propose three new hybrid metaheuristics, bat algorithm, shuffle frog-jumping algorithm, and cuckoo search, each of which is hybridized with the variable neighbourhood.

The contributions and innovations of this study can be summarized as follows:

• In this study, a novel multi-echelon, multi-objective, multi-period optimization framework for the production–location–inventory problem has been developed, tailored specifically for the perishable products industry. A new mixed-integer linear programming model, which incorporates demand uncertainties and longevity fluctuations in both raw materials and finished products, has been designed, setting it apart from traditional SCM models.
• The proposed model considers three pillars of sustainability: total cost, environmental impact, and social impact, aiming to construct a sustainable supply chain network for perishable products.
• Three swarm intelligence metaheuristics have been developed: the bat algorithm (BA), the shuffled frog leaping algorithm (SFLA), and the cuckoo search (CS) algorithm. Each of these metaheuristics has been hybridized with variable neighbourhood search (VNS) to solve this proposed model.
• The Taguchi method is applied to train the parameters, and the effectiveness of the proposed algorithms is evaluated using six assessment metrics: NPS, MID, SNS, DM, DEA, and POD metrics. Additionally, sensitivity analysis is conducted to further validate the suggested model.

3. Problem Description

The perishable supply chain network consists of four echelons: suppliers, manufacturers, distribution centres (DCs), and retailers. First, one must select the locations in which to operate DCs from candidate facilities. In each period, manufacturers procure raw materials from suppliers, process the different raw materials into different finished products, and deliver the finished products immediately to the operating DCs. That means there is no inventory for products in the manufacturer’s location. Following that, the products would be distributed to retailers by operating DCs. In this study, various decisions are taken into consideration, including the determination of the optimal locations of DCs, the allocation of raw materials from suppliers to manufacturers, the flow of products from manufacturers to retailers, the processing of products, and the management of inventory for both raw materials in manufacturers and products in DCs. Manufacturers at different capacity levels are also considered. In addition, considering the shelf lives of raw materials and products according to the possibility of perishing is one of the significant factors in this study.

The main goal of this paper is to optimize the supply chain network model for perishable products, considering the three pillars of sustainability. First, the aim is to minimize the total costs associated with established facilities, production, inventory holding, raw material procurement, deterioration, and transportation. Second, efforts are made to reduce the environmental impacts caused by established facilities, production, inventory, and transportation activities. Lastly, the focus is on maximizing social impacts by creating job opportunities in manufacturers and DCs to the greatest extent possible.

In our research, we regard total cost, environmental impacts, and social impacts as distinct objectives to optimization.

To make the problem realistic, our study supports the following assumptions:

• Raw materials and products are perishable, and they both have limited shelf lives.
• Each product requires several types and coefficient consumption levels of raw materials.
• The production capacity of each manufacturer is limited.
• Candidate DC locations and quantities are known.
• Retailer demands are uncertain and follow specific probability distributions.

This section may be divided by subheadings. These should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

Mathematical Model

The parameters and variables used in this paper are described in

Table 1. Description of parameters and variables.

Definition	Description
Sets and index
$i$	Supplier index i $\in I$
$j$	Manufacturer index j $\in J$
$k$	DC index k $\in K$
$l$	Retailer index l $\in L$
$r$	Raw material in r $\in R$
$p$	Perishable product index p $\in P$
$t$	Period index t $\in T$
Parameters
${rc}_{ir}$	The unit price of raw material $r$ provided by supplier $i$ (CNY)
${tc}_{ijr}^1$	The transportation cost of raw material $r$ from supplier $i$ to manufacturer $j$
$f_j^1$	The operating cost for manufacturer $j$ (CNY)
${hc}_{jr}^1$	The inventory cost of raw material $r$ in manufacturer $j$ (CNY)
${\alpha }_{rp}$	The coefficient of raw material $r$ consumption for perishable product $p$ (CNY)
${pc}_{jp}$	The production cost of perishable product $p$ in manufacturer $j$ (CNY)
${\tau }_r^1$	The shelf life of raw material $r$ (days)
${lc}_{jr}^1$	The deterioration cost of raw material $r$ in manufacturer $j$ (CNY)
${tc}_{jkp}^2$	The transportation cost of perishable product $p$ from manufacturer $j$ to DC $k$ (CNY)
$f_k^2$	The operating cost of DC $k$ (CNY)
${hc}_{kp}^2$	The inventory cost of perishable product $p$ in DC $k$ (CNY))
${\tau }_p^2$	The shelf life of perishable product $p$ (days)
${lc}_{kp}^2$	The deterioration cost of perishable product $p$ in DC $k$ (CNY)
${tc}_{klp}^3$	The transportation cost of perishable product $p$ from DC $k$ to retailer $l$ (CNY)
$d_{lpt}$	The demand of perishable product $p$ in retailer l at period $t$ (kg)
${et}_{ijr}^1$	The environmental impact of transporting raw material $r$ from supplier $i$ to manufacturer j (kg${\mathrm{c}\mathrm{o}}_2$)
${lc}_{kp}^2$	The deterioration cost of perishable product $p$ in DC $k$ (CNY)
${ef}_j^1$	The environmental impact of operating manufacturer $j$ (kg${\mathrm{c}\mathrm{o}}_2$)
${eh}_{jr}^1$	The environmental impact of storing raw material $r$ in manufacturer $j$ (kg${\mathrm{c}\mathrm{o}}_2$)
${ep}_{jr}$	The environmental impact of producing perishable product $p$ in manufacturer $j$ (kg${\mathrm{c}\mathrm{o}}_2$)
${ef}_k^2$	The environmental impact of operating DC $k$ in manufacturer $j$ (kg${\mathrm{c}\mathrm{o}}_2$)
${eh}_{kp}^2$	The environmental impact of storing perishable product $p$ in DC $k$ (kg${\mathrm{c}\mathrm{o}}_2$)
${et}_{jkp}^2$	The environmental impact of transporting perishable product $p$ from manufacturer $j$ to DC $k$ (kg${\mathrm{c}\mathrm{o}}_2$)
${et}_{klp}^3$	The environmental impact of transporting perishable product $p$ from DC $k$ to retailer l (kg${\mathrm{c}\mathrm{o}}_2$)
$a_j^1$	The number of permanent job opportunities created by manufacturer $j$
$a_{jp}^2$	The number of variable job opportunities created by manufacturer $j$ for producing perishable product $p$
$b_k^2$	The number of permanent job opportunities created by DC $k$
$b_{kp}^2$	The number of variable jobs opportunities created by DC $k$ for storing perishable p
$b_{kp}^2$	The number of variable jobs opportunities created by DC $k$ for storing perishable p
$m_{jpt}^{cap}$	The production capacity of manufacturer $j$ during period $t$
Decision variables
$x_k$	1 if DC $k$ is operated, 0 otherwise
$y_{ijrt}^1$	1 if supplier $i$ is assigned to manufacturer $j$ during period $t$, 0 otherwise
$y_{jkpt}^2$	1 if manufacturer $j$ is allocated to DC $k$ during period $t$, 0 otherwise
$y_{klpt}^3$	1 if DC $k$ is assigned to retailer $l$ during period $t$, 0 otherwise
$q_{ijrt}^1$	The allocation quantity of raw material $r$ from supplier $i$ to manufacturer $j$ during period $t$
$q_{jrt}^2$	The total inventory of raw material $r$ at manufacturer $j$ during period $t$
$q_{jpt}^3$	The producing quantity of raw material $r$ by manufacturer $j$ during period $t$
$q_{jrt}^4$	The cumulative deterioration of raw material $r$ at manufacturer $j$ during period $t$
$w_{jkpt}^1$	The quantity of perishable product $p$ purchased from manufacturer $j$ to DC $k$ during period $t$
$w_{kpt}^2$	The inventory quantity of perishable product $p$ at DC $k$ during period $t$
$w_{kpt}^3$	The cumulative deterioration of perishable product $p$ at DC $k$ during period $t$

.

Regarding the presented notations, the proposed multi-objective model (1)–(17) is formulated as follows:

$$\begin{array}{ll}Min{f}_1& ={\displaystyle \sum _{i,j,r,t}}{rc}_{ir}{q}_{ijrt}^1+{\displaystyle \sum _{i,j,r,t}}{tc}_{ijr}^1{q}_{ijrt}^1+{\displaystyle \sum _j}{f}_j^1+{\displaystyle \sum _{j,r,t}}{hc}_{jr}^1{q}_{jrt}^2+{\displaystyle \sum _{j,p,t}}{pc}_{jp}{q}_{jpt}^3\\ & +{\displaystyle \sum _{j,r,t}}{lc}_{jr}^1{q}_{jrt}^4+{\displaystyle \sum _{j,k,p,t}}{tc}_{jkp}^2{w}_{jkpt}^1+{\displaystyle \sum _k}{f}_k^2{x}_k+{\displaystyle \sum _{k,p,t}}{hc}_{kp}^2{w}_{kpt}^2\\ & +{\displaystyle \sum _{k,p,t}}{lc}_{kp}^2{w}_{kpt}^3+{\displaystyle \sum _{k,l,p,t}}{tc}_{klp}^3{d}_{lpt}{y}_{klpt}^3\end{array}$$

(1)

$$\begin{array}{ll}Min{f}_2& ={\displaystyle \sum _{i,j,r,t}}{et}_{ijr}^1{q}_{ijrt}^1+{\displaystyle \sum _j}{ef}_j^1+{\displaystyle \sum _{j,r,t}}{eh}_{jr}^1{q}_{jrt}^2+{\displaystyle \sum _{j,p,t}}{ep}_{jp}{q}_{jpt}^3\\ & +{\displaystyle \sum _{j,k,p,t}}{et}_{jkp}^2{w}_{jkpt}^1+{\displaystyle \sum _k}{ef}_k^2{x}_k+{\displaystyle \sum _{k,p,t}}{eh}_{kp}^2{w}_{kpt}^2+{\displaystyle \sum _{k,l,p,t}}{et}_{klp}^3{d}_{lpt}{y}_{klpt}^3\end{array}$$

(2)

$$Max{f}_3=\sum _{j,p,t}\left(a_j^1+a_{jp}^2{q}_{jpt}^3\right)+\sum _{k,p,t}\left(b_k^1{x}_k+b_{kp}^2{d}_{lpt}{y}_{klpt}^3\right)$$

(3)

Subject to the following:

$$\sum _i{y}_{ijrt}^1=1\hspace{1em}\forall j,r,t$$

(4)

$$\sum _{i,t}{q}_{ijrt}^1=\sum _{p,t}{\alpha }_{rp}{q}_{jpt}^3+\sum _t{q}_{jrt}^2+\sum _t{q}_{jrt}^4\hspace{1em}\forall j,r$$

(5)

$$q_{jrt}^2=\sum _i{q}_{ijrt}^1-\sum _p{\alpha }_{rp}{q}_{jpt}^3-q_{jrt}^4\hspace{1em}\forall j,r,t=1$$

(6a)

$$q_{jrt}^2=q_{jr\left[t-1\right]}^2+\sum _i{q}_{ijrt}^1-\sum _p{\alpha }_{rp}{q}_{jpt}^3-q_{jrt}^4\hspace{1em}\forall j,r,t\ge 2$$

(6b)

$$q_{jpt}^3\le m_{jpt}^{cap}\hspace{1em}\forall j,p,t$$

(7)

$$q_{jpt}^3=\sum _k{w}_{jkpt}^1\hspace{1em}\forall j,p,t$$

(8)

$$q_{jrt}^4=\sum _i{q}_{ijr\left(t-{\tau }_r^1\right)}^1-\sum _{t-{\tau }_r^1\le \sigma \le t}{\alpha }_{rp}{q}_{jp\sigma }^3\hspace{1em}\forall j,r,t,t>{\tau }_r^1$$

(9a)

$$q_{jrt}^4=0\hspace{1em}\forall j,r,t,t\le {\tau }_r^1$$

(9b)

$$\sum _j{y}_{jkpt}^2\le x_k\hspace{1em}\forall k,p,t$$

(10)

$$\sum _{j,t}{w}_{jkpt}^1=\sum _t{w}_{kpt}^2+\sum _t{w}_{kpt}^3+\sum _{l,t}{d}_{lpt}{y}_{klpt}^3\hspace{1em}\forall k,p$$

(11)

$$w_{kpt}^2=\sum _j{w}_{jkpt}^1-w_{kpt}^3-\sum _t{d}_{lpt}{y}_{klpt}^3\hspace{1em}\forall k,p,t=1$$

(12a)

$$w_{kpt}^2=w_{kp\left[t-1\right]}^2+\sum _j{w}_{jkpt}^1-w_{kpt}^3-\sum _t{d}_{lpt}{y}_{klpt}^3\hspace{1em}\forall k,p,t\ge 2$$

(12b)

$$\sum _k{y}_{klpt}^3=1\hspace{1em}\forall l,p,t$$

(13)

$$y_{klpt}^3\le x_k\hspace{1em}\forall k,l,p,t$$

(14)

$$w_{kpt}^3=\sum _j{w}_{jkp\left(t-{\tau }_p^2\right)}^1-\sum _{l,t-{\tau }_p^2\le \sigma \le t}{d}_{lpt}{y}_{klpt}^3\hspace{1em}\forall k,p,t,t>{\tau }_p^2$$

(15a)

$$w_{kpt}^3=0\hspace{1em}\forall k,p,t,t>{\tau }_p^2$$

(15b)

$$x_k,y_{ijrt}^1,y_{jkpt}^2,y_{klpt}^3\in \left\{\mathrm{0,1}\right\}$$

(16)

$$q_{ijrt}^1,q_{jpt}^2,q_{jpt}^3,q_{jpt}^4,w_{jkpt}^1,w_{kpt}^2,w_{kpt}^3\ge 0$$

(17)

Objective function (1) decreases total system costs, which comprise raw material procurement, manufacturer and DC operation, inventory, production, deterioration, and transportation costs. Objective function (2) minimizes SCN environmental effects associated with operation, production, inventory, and transportation activities. Objective function (3) maximizes job creation to promote sustainable development social factors. Constraint (4) indicates that a raw material type can be assigned to only one supplier. Constraint (5) calculates the total number of raw materials procured by manufacturers in the planning horizon. Constraint (6) demonstrates the inventory quantities of raw materials in manufacturers. Constraint (7) specifies manufacturer production capacity. Constraint (8) calculates production orders in each period. Constraint (9) reveals the deterioration number of raw materials in manufacturers. Constraint (10) imposes that manufacturers can only serve operating DCs. Constraint (11) calculates the total number of perishable products purchased by DCs in the planning horizon. Constraint (12) represents the inventory quantity of perishable products in DCs. Constraints (13) and (14) state that a retailer must be served by only one operating DC. Constraint (15) calculates the number of deteriorated products incurred in DCs. Finally, Constraints (16) and (17) represent decision variable sorting and range.

4. Solution Technique

The proposed mathematical model, as mentioned earlier, is initially designed to handle multiple objectives. To address the problem effectively, we adopt a multi-objective solution technique inspired by Goodarzian, Taleizadeh, Ghasemi, and Abraham [13], which converts the model into a single-objective formulation. Subsequently, we propose three hybrid meta-heuristic algorithms to tackle the problem. To ensure a fair comparison, all algorithms incorporate the variable neighbourhood search (VNS) algorithm as a local search component. These algorithms include HBA-VNS (a hybrid of bat algorithm (BA) and VNS), HSFLA-VNS (a hybrid of shuffled frog leaping algorithm (SFLA) and VNS), and HCS-VNS (a hybrid of cuckoo search (CS) and VNS). Additionally, we evaluate and compare these approaches among themselves and benchmark their results against a reference method, a hybrid genetic algorithm with VNS (HGV-VNS). In the following subsections, we provide an overview of the framework used by the presented algorithms. First, we discuss the solution representation and the hybrid of VNS as a vital component within the algorithm’s structure. Subsequently, we describe the presented algorithms in detail.

4.1. Solution Representation

To empower the hybrid meta-heuristic algorithms that have been developed to tackle discrete problems, the Random-Key (RK) technique is utilized for solution representation. This technique serves to encode the original solution, thereby equipping continuous algorithms with the ability to effectively address discrete problems [41]. Within the context of the developed problem, the sub-solutions are categorized into two distinct types:

(1)

The first solution scheme involves the method for determining the operating DCs. To initiate the process, a matrix comprising |J| elements is generated using a uniform distribution U (0,1). Subsequently, the first Max J units with the highest values are selected as the operating DCs. For instance, the solution {0.51, 0.39, 0.89, 0.21, 0.43} with $J_{max}$ = 3 illustrates the decoded solution {1, 0, 1, 0, 1}.

(2)

The second solution scheme involves the allocation relationship, which aims to allocate suppliers to manufacturers. As shown in Figure 1, a vector, max $\left(1\times \left(i_{max}+\left|j\right|-1\right)\right)$, is generated with a uniform distribution in U (0, 1). This vector serves as a permutation matrix, indicating the degree of each value. Values greater than $\left|i_{max}\right|$ represent manufacturer’s symptoms, while values less than $\left|i_{max}\right|$ correspond to suppliers included in the permutation matrix. The manufacturer’s symptoms act as dividers, distinguishing different suppliers. Thus, the sequence of suppliers among these dividers determines the allocation points for a manufacturer.

4.2. Variable Neighbourhood Search Algorithm

The VNS algorithm is a popular local search method utilized for optimizing different problems. It employs a neighbourhood structure (NS) and systematically modifies it to achieve a near-optimal solution [13,41,42]. The VNS algorithm has demonstrated remarkable performance and efficiency in solving various optimization problems in supply chain management such as location–distribution [43], production–distribution–inventory–allocation–location [13], and network design [44]. The VNS algorithm follows an initialization procedure and employs a set of neighbourhood structures (NSs), denoted as $N_{ns}$($ns$ = 1, 2, …, ${ns}_{max}$), where ${ns}_{max}$ represents the maximum number of NSs with a predetermined sequence. Algorithm 1 outlines the pseudocode for the VNS algorithm.

Algorithm 1: Pseudocode of VNS
1	Initialization
	Set parameters (maximum iteration (${Mi}_{VNS}),$ neighborhood structures $N_{ns}\left(ns=\mathrm{1,2},3,\dots ,{ns}_{max}\right)$)
	Generate or receive a population of size nPop and evaluate individual
	Select the best solution from the evaluated individuals as the initial solution $S^{*}$
2	for k to Vmax_it do:
3	for i to ${ns}_{max}$ do
4	S$\leftarrow S^{*}$
5	Randomly select a solution $S^{\prime }$ from $N_{ns}$// shaking
6	Apply a neighbourhood search $N_{ns}$ to $S^{\prime }$ to obtain the improved solution $S^{″}$
7	Evaluate $S^{″}$
8	If f$\left(S^{″}\right)<$ f$\left(S^{*}\right)$
9	$S^{*}\leftarrow S^{″}$ and break
10	Else if f$\left(S^{″}\right)=$ f$\left(S^{*}\right)$
11	Select one of them randomly
12	End if
13	End for
14	End for
15	Report  $S^{*}$
16	End

4.3. HBA-VNS Algorithm

The BA is an innovative swarm intelligence technique specifically designed to mimic the echolocation system of bats. First introduced by Yang [45], the BA is a stochastic algorithm that operates on a population-based approach. Over the past few years, the BA has demonstrated successful applications in various combinatorial optimization domains, including inventory problems [46], SCND [47], vehicle-routing problems [48], and p-facility maximum location [49]. During the optimization process, the update of a bat $i$, positioned at $x_i^t,$ with a velocity $v_i^t,$ at time $t$, is governed by Equations (18)–(20).

$$f_i=f_{min}+\left(f_{max}-f_{min}\right)\times \alpha$$

(18)

$$v_i^t=v_i^{t-1}+\left(x_i^t-x^{\mathrm{*}}\right)\times f_i$$

(19)

$$x_i^t=x_i^{t-1}+v_i^t$$

(20)

where $f_i$ denotes the pulse frequency emitted by bat i which is uniformly distributed U ($f_{min}$, $f_{max}$). $\alpha$ is a randomly generated number between 0 and 1. $x^{*}$ represents the current global optimal position at the time t.

Typically, during the search process, the loudness of a bat decreases while the pulse increases gradually when it finds a prey, as shown in Equations (21) and (22).

$$A_i^t=\beta \times A_i^{t-1}$$

(21)

$$r_i^t=r_i^0\times \left[1-exp\left(-\gamma t\right)\right]$$

(22)

where $\beta$ and $\gamma$ are positive constants. For any $0<\beta <1$ and $\mathsf{\gamma }>0$, it is shown that when $\mathrm{t}\to +\mathrm{\infty }$, the loudness $A_i^t\to 0$ and the wave length $r_i^t\to r_i^0.$ VNS is a widely utilized method known for its effective local optimization capabilities, often integrated into hybrid meta-heuristic algorithms. The VNS algorithm utilizes local modifications to explore the solution space, iteratively moving from the current solution to a new one in order to identify a locally optimal solution. Algorithm 2 presents the pseudocode for the HBA-VNS algorithm.

Algorithm 2: The pseudocode of HBA-VNS
1	Define the initial parameter Popsize (population size) and MaxIter (maximum iteration)
2	Generate initial bats population $x_i$ and $v_i$
3	Set the initial values of the pulse rate $r_i^0$ and the loudness $A_i^0$
4	While The algorithm does not reach the terminal criterion (MaxIter) do
5	for Each bat i in the population do
6	Generate a new solution by applying Equations (18)–(20).
7	if $\left(rand>r_i\right)$
8	Select a solution from the best individuals in the current population
9	Generate a local solution around the selected best solution
10	end if
11	Generate a new solution by flying randomly
12	if $\left(rand<A_i \mathrm{a}\mathrm{n}\mathrm{d} f\left(x_i\right)<f\left(x^{*}\right)\right)$
13	Accept the new solution
14	Increase parameters $r_i$ and reduce $A_i$ according to Equations (21) and (22).
15	end if
16	end for
17	Rank the current population and return to the optimal solution $x^{*}$
18	Pick up the optimal solution $x^{*}$ as initial solution, apply VNS to local search/perform Algorithm 1
19	end while

4.4. HSFLA-VNS Algorithm

The SFLA, a nature-inspired algorithm, emulates the foraging behaviour of frogs and was first introduced by Eusuff and Lansey [50] for optimizing water distribution network designs. The SFLA has demonstrated successful applications in solving diverse problems, such as the assembly line sequencing problem [51], flexible job shop scheduling [52], vehicle routing problem [53,54], and logistic scheduling problem [55]. However, traditional SFLA exhibits some limitations, such as inefficient local search and premature convergence. To optimize the efficiency of the local search procedure in SFLA, we propose the hybrid of VNS to tackle our specific problem. For more detailed information about SFLA, researchers are advised to refer to the work of [53]. Algorithm 3 presents the pseudocode for the HSFLA-VNS algorithm.

Algorithm 3: The pseudocode of HSFLA-VNS
1	Set parameters ($MaxIter$ (global iterations number), $L_{ite}$ (local iterations number), Popsize (population size), m (memeplexes)
2	Initialise frogs randomly and assess them
3	While the algorithm does not reach the terminal criterion (MaxIter) do
4	for j = 1 to m do
5	for k = 1 to $L_{ite}$ do
6	Sort the frogs in the memeplexes j
7	Update $x_w,$ generate $x_w^{new}$
8	if $f\left(x_w^{new}\right)<f\left(x_w\right)$
9	$x_w$ $\leftarrow x_w^{new}$, continue
10	$x_g$ replaces $x_b$ to update $x_w$
11	if $f\left(x_w^{new}\right)<f\left(x_w\right)$
12	$x_w$ $\leftarrow x_w^{new}$, continue
13	$x_w$ is replaced by random solution
14	end if
15	end for
16	end for
17	Combine all frogs in each memeplex into a new population
18	Rank the new population and return to the optimal solution $x_s^{*}$
19	Pick up the optimal solution ${\mathit{x}}_{\mathit{s}}^{\mathit{*}}$ as initial solution, apply VNS to local search//perform Algorithm 1
20	end while
21	Return to the global best solution $x_g^{*}.$

4.5. HCS-VNS Algorithm

The CS algorithm, introduced by Yang and Deb [56], is a meta-heuristic intelligence algorithm that draws inspiration from the breeding behaviour of cuckoos in nature. Its core principles encompass the parasitic breeding style of cuckoos and the Lévy flight mechanism. CS has demonstrated successful applications in various domains, including the flow shop scheduling problem [57,58], production routing problem [59], and production–inventory–routing problem [37]. For a more comprehensive understanding of the CS, researchers are encouraged to consult the works of [56,59]. To further enhance the efficiency of the local search procedure in CS, we have integrated it with VNS to address our specific problem. Algorithm 4 outlines the pseudocode for the HCS-VNS algorithm.

Algorithm 4: The pseudocode of HCS-VNS
1	Set parameters Popsize (host nets), MaxIter (maximum iterations number), $P_{\alpha }$ (host detection probability)
2	Generate a population of Popsize host nets $\left(i\in \left[1, N\right]\right)$
3	Evaluate fitness value, $obj\left(x_i\right)$ of every nest i
4	Record initial global best solution as $x_b$
5	While the algorithm does not reach the terminal criterion (MaxIter) do
6	Lévy flight, update all host nests position
7	Evaluate its objective value $obj\left(x_i\right)$
8	Choose a nest $l\in \left[1, N\right]$ randomly
9	if $\left(obj\left(x_i\right)<obj\left(x_l\right)\right)$
10	Replace solution l with solution i
11	end if
12	A fraction $\left(P_{\alpha }\right)$ of worst nets are abandoned
13	New nests/solution are built using Lévy flight
14	Preserve the optimal solutions (or nests with quality solutions)
15	Evaluate the solutions and identify the current optimal solution as $x_b^{*}$
16	Pick up the optimal solution ${\mathit{x}}_{\mathit{b}}^{\mathit{*}}$ as initial solution, apply VNS to local search//perform Algorithm 1
17	end while
18	Return to the global best solution $x_g^{*}$
19	Output:  $x_g^{*}$

5. Computational Experiments

To evaluate the efficacy of the proposed metaheuristic algorithm, we first introduce the data generation process. Subsequently, a detailed explanation of the evaluation metrics and techniques used for parameter tuning are provided. After that, the obtained results by the algorithms, HBA-VNS, HSFLA-VNS, HCS-VNS, and the benchmark algorithm HGA-VNS, are analysed and compared based on different evaluation metrics. In addition, the proposed algorithms are evaluated according to the obtained Pareto frontier. Finally, to validate the SPPSCN model, a set of sensitivity analyses was employed. The implementation of the proposed algorithms was conducted using MATLAB^® 2020a software on a PC with 2.4 GHz.

5.1. Instances

In this paper, a numerical experiment was devised to compare the suggested algorithms with HGA-VNS. A search of the literature revealed a lack of standard benchmark dataset for our research problem; randomly generated test instances are required to represent the practical problems. The data generated are presented in

Table 2. Parameter generation scheme.

Parameter	Range	Parameter	Range
${rc}_{ir}$	Uniform (80, 90)	$d_{klp}$	Uniform (30, 50)
${tc}_{ijr}^1$	Uniform (20, 50)	${et}_{ijr}^1$	Uniform (3.0, 7.5)
$f_j^1$	105 Uniform (3.0, 3.6)	${ef}_j^1$	104 Uniform (2.4, 3.6)
${hc}_{jr}^1$	Uniform (35, 45)	${eh}_{jr}^1$	Uniform (1.5, 4.5)
${\alpha }_{rp}$	Uniform (0.1, 0.3)	${ep}_{jr}$	Uniform (7.5, 20)
${pc}_{jp}$	Uniform (90, 110)	${ef}_k^2$	104 Uniform (1.2, 1.8)
${\tau }_r^1$	Uniform (5, 7)	${eh}_{kp}^2$	Uniform (1.5, 4.0)
${lc}_{jr}^1$	Uniform (12, 50)	${et}_{jkp}^2$	Uniform (3.0, 9.0)
${tc}_{jkp}^2$	Uniform (25, 50)	${et}_{klp}^3$	Uniform (3.0, 9.0)
$f_k^2$	105 Uniform (1.6, 2.0)	$a_j^1$	Uniform int (100, 170)
${hc}_{kp}^2$	Uniform (30, 40)	$a_{jp}^2$	Uniform int (1.5, 1.6)
${\tau }_p^2$	Uniform (5, 7)	$b_k^2$	Uniform int (70, 90)
${lc}_{kp}^2$	Uniform (50, 70)	$b_{kp}^2$	Uniform int (1.1, 1.2)
${tc}_{klp}^3$	Uniform (25, 50)	$m_{jpt}^{cap}$	Uniform (800, 1200)

. Additionally, 10 instances with different sizes are illustrated in

Table 3. Instances size.

No.	|I|–|J|–|K|–|L|–|R|–|P|–|T|	No.	|I|–|J|–|K|–|L|–|R|–|P|–|T|
1	2–1–3–8–2–2–2	6	7–2–7–18–6–6–6
2	3–1–4–10–2–2–2	7	8–3–7–20–6–6–6
3	4–1–5–12–4–4–4	8	9–3–8–25–8–8–8
4	5–2–5–12–4–4–4	9	11–3–9–27–8–8–10
5	6–2–6–15–6–6–6	10	14–5–12–50–10–10–10

, which comprises suppliers (I), manufacturers (J), the candidate DCs (K), retailers (L), raw materials (R), perishable product (P), and the periods (T).

5.2. The Taguchi Method for Parameter Setting

The setting of parameters or metaheuristic algorithms is crucial as it directly impacts their performance and efficiency [60,61]. In this study, the Taguchi method was employed to tune these parameters. A key concept in the Taguchi method is the signal-to-noise ratio, which measures the ratio between the signal and the noise, serving as an important performance indicator for the algorithm. Maximizing the signal-to-noise ratio during the parameter tuning process helps identify the optimal parameter combination, ensuring optimal algorithm performance and desired outcomes. Equation (23) presents the formulation of the S/N value.

$$\mathrm{S}/\mathrm{N} \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}=-10\times \mathrm{l}\mathrm{o}\mathrm{g}\left(\sum _{i=1}^n{Y}_i^2/n\right)$$

(23)

where $Y_i$ represents the response value for the $i$th orthogonal array, and n denotes the total number of orthogonal arrays.

By employing orthogonal arrays, the Taguchi method effectively reduces the total number of tests required to control the proposed algorithms within a reasonable timeframe. In this study, the L27 orthogonal array was chosen for the all algorithms, HBA-VNS, HSFLA-VNS, HCS-VNS, and HGA-VNS, using Minitab 19 Software. The levels of the factors (corresponding to the algorithm’s parameters) are presented in

Table 4. Factors and their levels.

Algorithm	Factors	Levels
HBA-VNS	Popsize: the population size of bats	$\mathrm{A}1$:60 A2:90 A3:120
	MaxIter: the maximum iterations of BA	$\mathrm{B}1$:200 B2:300 B3:400
	$A_i^0:$ the initial loudness of bat i	$C$1:0.8 C2:1.0 C3:1.2
	$r_i^0:$ the initial wavelength of bat i	$\mathrm{D}$1:0.2 D2:0.4 D3:0.6
	$\beta :$ the constant value for adjusting the loudness	$\mathrm{E}$1:0.85 E2:0.90 E3:0.95
	$\gamma :$ the constant value for adjusting the wave length	$\mathrm{F}$1: 0.85 F2:0.90 F3:0.95
	${Mi}_{VNS}$: the maximum number of iterations for VNS	G1:100 G2:150 G3:200
	$N_l:$ the number of neighborhood structure	$\mathrm{H}$1:3 H2:4 H3:5
HSFLA-VNS	Popsize: the population size of frogs	$\mathrm{A}1$:60 A2:90 A3:120
	$MaxIter:$ the maximum iterations of SFLA	$\mathrm{B}1$:200 B2:300 B3:400
	$L_{iter}:$ the evolutionary times within each memeplex	$C$1:20 C2:25 C3:30
	$m:$ the number of memeplexes	$\mathrm{D}$1:6 D2:10 D3:15
	${Mi}_{VNS}$: the maximum number of iterations for VNS	E1:100 E2:150 E3:200
	$N_l:$ the number of neighborhood structure	F1:3 F2:4 F3:5
HCS-VNS	Popsize: the population size of cuckoos	$\mathrm{A}1$:60 A2:90 A3:120
	$MaxIter:$ the maximum iterations of CS	$\mathrm{B}1$:200 B2:300 B3:400
	$P_{\alpha }:$ the host detection probability	$C$1:0.2 C2:0.3 C3:0.4
	${Mi}_{VNS}$: the maximum number of iterations for VNS	D1:100 D2:150 D3:200
	$N_l:$ the number of neighborhood structure	E1:3 E2:4 E3:5
HGA-VNS	Popsize: the population size of GA	$\mathrm{A}1$:60 A2:90 A3:120
	$MaxIter:$ the maximum iterations of GA	$\mathrm{B}1$:200 B2:300 B3:400
	$P_c:$ the crossover probability	$C$1:0.7 C2:0.8 C3:0.9
	$P_m:$ the mutation probability	D1:0.1 D2:0.2 D3:0.3
	${Mi}_{VNS}$: the maximum number of iterations for VNS	E1:100 E2:150 E3:200
	$N_l:$ the number of neighborhood structure	F1:3 F2:4 F3:5

, where three levels are specified for the algorithm’s factors.

When comparing results obtained from scientific experiments, challenges may arise due to differences in the test problem or objective function size. To address this issue, researchers have devised a method to convert experimental data into relative percentage deviations (RPD). Equation (24) is commonly employed to calculate the RPD value of experimental data.

$$\mathrm{R}\mathrm{P}\mathrm{D}=\frac{{Alg}_{sol}-{Best}_{sol}}{{Best}_{sol}}\times 100\%$$

(24)

The terms ${Alg}_{sol}$ and ${Best}_{sol}$ denote the objective function values obtained for each iteration of the algorithm and the best solution discovered by the algorithm, respectively. The objective function values are then converted to the RPDs for each experiment, and the mean RPD is computed.

The S/N ratio output is analysed by Minitab software to identify the optimal levels for each algorithm, as illustrated in Figure 2. Higher values of the S/N ratio indicate more robust algorithms, and these levels are considered optimal. As depicted in Figure 2, the optimal levels for the HBA-VNS algorithm are established as follows: A3, B3, C2, D3, E2, F2, G1, and H2. For the HSFLA-VNS algorithm, the optimal levels are identified as follows for the respective parameters: A3, B2, C1, D1, E1, and F2. Similarly, for the HCS-VNS algorithm, the optimal levels are set as follows for the respective parameters: A3, B2, C2, D1, and E3. Lastly, for the HGA-VNS algorithm, the optimal parameter levels are determined to be A3, B3, C1, D3, E3, and F3.

5.3. Comparison of Results

Based on the relevant model parameter data extracted from Table 2 and Table 3, alongside the optimal parameter configurations obtained through the Taguchi experimental method, a comparative analysis was conducted on the proposed algorithms: HBA-VNS, HSFLA-VNS, HCS-VNS, and HGA-VNS. Figure 3 showcases the Pareto fronts generated by these algorithms for P3, P6, and P9, capturing various metrics such as total cost, environmental impact, and social impact. The findings underscore the superior efficiency of HBA-VNS, closely followed by HCS-VNS. The remaining two algorithms exhibit comparable performance, with marginal differentials between them.

5.4. Assessment Metrics

Moreover, different efficiency assessment metrics are explained to facilitate a comparative evaluation of the quality of the Pareto-optimal solutions attained by the algorithms. Our study utilizes six assessment metrics: (a) number of Pareto solutions (NPS), (b) mean ideal distance (MID), (c) the spread of non-dominance solution (SNS), (d) diversification metric (DM), (e) data envelopment analysis (DEA), and (f) percentage of domination (POD). These assessment metrics have been employed in previous studies by Govindan, Jafarian, and Nourbakhsh [41] and Kamran et al. [62].

The metrics are derived from the attained Pareto fronts across all test scenarios, offering a comprehensive evaluation of each method’s performance. The outcomes of these assessment metrics are detailed in

Table A1. The obtained metrics for algorithms performance (NPS, MID, SNS, DM, DEA and POD).

Assessment Metrics	Instances	HBA-VNS	HSFLA-VNS	HCS-VNS	HGA-VNS	Assessment Metrics	Instances	HBA-VNS	HSFLA-VNS	HCS-VNS	HGA-VNS
NPS	P1	11	7	10	9	MID	P1	15,256	25,843	13,470	16,982
	P2	12	7	12	8		P2	12,785	20,835	17,482	20,567
	P3	14	9	11	9		P3	20,551	22,009	21,344	22,915
	P4	11	10	9	7		P4	19,504	24,955	14,949	24,744
	P5	10	9	11	9		P5	10,293	13,933	15,007	24,187
	P6	12	9	11	10		P6	14,778	22,922	14,083	13,051
	P7	10	8	7	6		P7	11,766	13,145	14,778	22,530
	P8	9	6	8	9		P8	17,936	24,488	17,845	24,172
	P9	11	8	9	6		P9	17,579	20,181	17,769	23,385
	P10	14	10	11	6		P10	15,174	13,042	12,709	18,411
SNS	P1	4748	3124	2536	3299	DM	P1	871	281	1159	516
	P2	4584	4910	4596	3202		P2	999	716	1157	213
	P3	4816	2474	2886	3546		P3	1356	469	717	304
	P4	4831	2361	2889	4837		P4	917	294	1297	1094
	P5	4871	2123	4023	2650		P5	1283	731	933	259
	P6	4871	2284	3276	3155		P6	779	222	1139	507
	P7	4989	3241	2880	2483		P7	1134	944	1018	292
	P8	4990	2477	3130	2323		P8	772	567	1246	755
	P9	3496	2043	3392	2997		P9	1155	887	980	631
	P10	4626	4351	4442	4993		P10	922	812	690	892
DEA	P1	0.43	0.53	0.44	0.28	POD (%)	P1	30	25	29	15
	P2	0.5	0.23	0.3	0.27		P2	25	15	33	20
	P3	0.49	0.44	0.51	0.34		P3	27	19	22	11
	P4	0.57	0.24	0.3	0.35		P4	34	33	30	20
	P5	0.56	0.53	0.22	0.26		P5	35	12	21	16
	P6	0.60	0.26	0.58	0.33		P6	30	10	16	29
	P7	0.53	0.37	0.45	0.21		P7	31	19	22	11
	P8	0.47	0.39	0.28	0.6		P8	24	16	22	10
	P9	0.23	0.23	0.52	0.58		P9	28	10	23	19
	P10	0.48	0.24	0.6	0.58		P10	36	10	30	18

of Appendix A. To ensure comparability, the results presented in Table A1 are standardized into the relative deviation index (RDI) [41], as defined by Equation (25).

$$\mathrm{R}\mathrm{D}\mathrm{I}=\frac{\left|{Alg}_{sol}-{Best}_{sol}\right|}{{Max}_{sol}-{Min}_{sol}}\times 100\%$$

(25)

where ${Alg}_{sol}$ represents the objective value, ${Max}_{sol}$ and ${Min}_{sol}$ are the maximum and minimum values for each performance metric, respectively, and ${Best}_{sol}$ denotes the best solution in the algorithm (either ${Max}_{sol}$ or ${Min}_{sol}$ depending on the metric’s nature). It should be noted that a lower RDI value indicates better performance.

A statistical analysis of algorithm effectiveness was performed by computing a 95% confidence interval for the performance metrics of all algorithms. Figure 4 displays the means plot and least significant difference (LSD) intervals for all algorithms. Figure 4a,f show statistically significant differences among the methods concerning the NPS and POD metrics, with HBA-VNS outperforming the others. HCS-VNS follows, while HSFLA-VNS and HGA-VNS show no significant differences between them, indicating HBA-VNS obtains more non-dominated solutions. Figure 4b,d reveal no statistically significant differences between HBA-VNS and HCS-VNS in the MID and DM metrics but significant differences compared to HSFLA-VNS and HGA-VNS. These results suggest that HBA-VNS and HCS-VNS provide solution sets that are more uniformly distributed along the Pareto front, exhibit better convergence, and show superior diversity. Figure 4c indicates a statistically significant difference between HBA-VNS and the other three algorithms in the SNS metric. The non-dominated solution set obtained by HBA-VNS is more extensive and closer to the Pareto front, while the other three methods show no significant differences among themselves. Finally, Figure 4e shows no statistically significant differences among the algorithms, indicating similar efficiency in obtaining all non-dominated solutions. Overall, HBA-VNS demonstrates the best performance, followed by HCS-VNS, with HSFLA-VNS and HGA-VNS performing similarly.

5.5. Analysis and Discussion

To better understand the behaviour of the SPPSCN model and evaluate the practical benefits of the proposed approach, several sensitivity analyses were conducted on key model parameters. In this regard, a test problem, P9, was chosen, and the HBA-VNS algorithm was identified as the most efficient meta-heuristic algorithm for this study. Various variations of the extended SPPSCN model, including demand, raw material shelf life, and product shelf life, were examined. The detailed results of the sensitivity analysis can be found in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

Figure 5 demonstrates the direct relationships between demand change and total costs, environmental impacts, and social impacts. However, these factors exhibit varying fluctuation rates. Specifically, when demand increases by 40%, total costs, environmental impact, and social impact rates increase by 40%, 39%, and 33%, respectively.

Figure 6 illustrates how extending the shelf life of raw materials directly impacts costs, the environment, and society. However, these reductions and variations vary across time intervals. During the transition from one to two periods, there is a notable decrease in total costs (21%), environmental impacts (12%), and social impacts (3%). This period demonstrates a relatively higher rate of decrease compared to the approximately 10%, 8%, and 0% reductions observed between two and eight periods. Notably, between eight and ten periods, all factors remain constant.

To explore the underlying factors, we performed a sensitivity analysis on raw material shelf life across all components in the objective functions. Figure 7 demonstrates that transitioning from one to two periods significantly affects deterioration cost. Procurement, transportation, and inventory costs are also sensitive to shelf life changes, albeit to a lesser extent. A shorter shelf life increases waste, requiring more purchases, transportation, and inventory management, leading to heightened environmental impacts. However, job creation by manufacturers and DCs remains unaffected, indicating limited social impacts. Between two and eight periods, shelf life primarily influences raw material inventory performance. Beyond eight periods, all factors remain constant.

Similar to the previous analysis on raw material shelf life, we examined the impact of product shelf life on objective functions. Figure 8 reveals that extending product shelf life reduces total cost and environmental impact within the supply chain system. However, job quantity may experience a slight decrease. Reductions and variations vary across intervals. Transitioning from one to two periods yields significant decreases in total cost (43%), environmental impact (40%), and social impact (25%) rates. Between two and ten periods, the reductions follow a slope of approximately 6%, 5%, and 0%.

Figure 9 reports the sensitivity analysis of product shelf life on various objective functions. From one to two periods, deterioration cost exhibits high sensitivity to product shelf life, followed by procurement, production, transportation, operation, and inventory costs. A longer shelf life reduces product waste, resulting in decreased raw material purchases, production and inventory operations, DC utilization, and transportation within the network, thereby lowering environmental impacts. Social impacts are mainly attributed to reduced production and the limited opening of new DCs. Between two and ten periods, shelf life primarily affects inventory performance, while other factors show no sensitivity to product shelf life.

Through a sensitive analysis conducted on instance P9, it is observed that an increase in demand has a nearly linear effect on the corresponding objectives. Conversely, the impact of extending raw material shelf life and product shelf life on these objectives is limited to a certain range. Beyond a certain threshold, the influence becomes insignificant.

5.6. Insights

• Theoretically, this study considers various aspects of sustainability, including economic, environmental, and social impacts. Based on the shelf life of raw materials and finished products and considering the uncertainties of different parameters in the model, a multi-objective optimization model was constructed.
• To solve the proposed model, three heuristic methods based on metaheuristic algorithms were extended for the first time. Six evaluation metrics—NPS, MID, SNS, DM, DEA, and POD—were used to validate the proposed hybrid algorithms. The results of these metrics indicate that the HBA-VNS algorithm demonstrates high performance and efficiency.
• The established joint optimization model enables supply chain decision-makers to make comprehensive decisions under uncertain conditions, thereby determining the optimal logistics within the network.
• Focusing on the shelf life of raw materials and finished products can reduce total costs and environmental impacts within a certain range, but it may also reduce social benefits. Decision-makers can balance the costs with environmental and social benefits of the supply chain network by designing reasonable range values.

6. Conclusions and Future Works

This paper has presented a comprehensive model for a perishable supply chain network, encompassing the three pillars of sustainability. It has considered the unique characteristics of perishable products holistically. The developed model addresses the multi-echelon, multi-product, and multi-period production–location–inventory problem in the SPPSCN. To solve this model, three hybrid metaheuristic algorithms (HBA-VNS, HSFLA-VNS, and HCS-VNS) have been proposed, each combining bat algorithm (BA), shuffled frog leaping algorithm (SFLA), and cuckoo search (CS) with variable neighbourhood search (VNS). The Taguchi method has been utilized to tune the parameters, considering the sensitivity of metaheuristic methods. The quality of the obtained Pareto frontier has been assessed by comparing the proposed approaches with a benchmark algorithm, HGA-VNS, using six metrics. The numerical experiments have demonstrated the superior performance of HBA-VNS. Moreover, the sensitivity analysis on instance P9 indicates that an increase in demand has an almost linear impact on the corresponding objectives. In contrast, extending the shelf life of raw materials and products only affects these objectives within a specific range. Beyond this point, the effect becomes negligible.

Future research can explore alternative uncertainty programming techniques, including robust optimization and fuzzy set, to enhance the modelling of uncertainties. The incorporation of resilient concepts, vehicle routing planning, diverse transportation modes, and cold chain logistics offers an intriguing direction for problem solving in SCNs. Moreover, scholars can investigate the integration of artificial intelligence algorithms into the framework of metaheuristic algorithms to improve the solution approach for the proposed problem. Lastly, leveraging the potential of Blockchain and Internet of Things technology can be explored in future studies to enhance the performance of the SPPSCN.