Next Article in Journal
The Efficacy of Micronutrient Fertilizers on the Yield Formulation and Quality of Wheat Grains
Previous Article in Journal
The Current Distribution of Carex alatauensis in the Qinghai–Tibet Plateau Estimated by MaxEnt
Previous Article in Special Issue
Research on Control Strategy of Light and CO2 in Blueberry Greenhouse Based on Coordinated Optimization Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Sustainable Phosphorus Fertilizer Supply Chain Management to Improve Crop Yield and P Use Efficiency Using an Ensemble Heuristic–Metaheuristic Optimization Algorithm

by
Mohammad Shokouhifar
1,*,
Mahnaz Sohrabi
2,
Motahareh Rabbani
3,
Seyyed Mohammad Hadji Molana
3 and
Frank Werner
4,*
1
Department of Electrical and Computer Engineering, Shahid Beheshti University, Tehran 1983969411, Iran
2
Faculty of Industrial and Mechanical Engineering, Qazvin Branch, Islamic Azad University, Qazvin 3471993116, Iran
3
Department of Industrial Engineering, Science and Research Branch, Islamic Azad University, Tehran 1477893855, Iran
4
Faculty of Mathematics, Otto-Von-Guericke-University, 39016 Magdeburg, Germany
*
Authors to whom correspondence should be addressed.
Agronomy 2023, 13(2), 565; https://doi.org/10.3390/agronomy13020565
Submission received: 27 December 2022 / Revised: 5 February 2023 / Accepted: 14 February 2023 / Published: 16 February 2023

Abstract

:
Phosphorus (P) is the most important substance in inorganic fertilizers used in the agriculture industry. In this study, a multi-product and multi-objective model is presented considering economic and environmental concerns to design a renewable and sustainable P-fertilizer supply chain management (PFSCM) strategy. To handle the complexities of the model, an ensemble heuristic–metaheuristic algorithm utilizing the heuristic information available in the model, the whale optimization algorithm, and a variable neighborhood search (named H-WOA-VNS) is proposed. First, a problem-dependent heuristic is designed to generate a set of near-optimal feasible solutions. These solutions are fed into a population-based whale optimization algorithm which benefits from exploration and exploitation strategies. Finally, the single-solution variable neighborhood search is applied to further improve the quality of the solution using local search operators. The objective function of the algorithm is formulated as a weighted average function to minimize total economic cost while increasing crop yield and P use efficiency. The experimental results for a real case study of the P-fertilizer supply chain confirm the effectiveness of the proposed method in improving the crop yield and P use efficiency by 33% and 27.8%, respectively. The results demonstrate that the proposed H-WOA-VNS algorithm outperforms the Heuristic, WOA, and VNS models in reducing the total objective function value of the PFSCM model by 9.8%, 2.9%, and 4%, respectively.

1. Introduction

Phosphorus (P) is an indispensable nutrient that is essential to global food security and plays a vital role in crop growth and soil productivity. Nevertheless, this essential substance has several environmental effects [1]. Based on the simulation results of a recent work in 2020 [2], the production of phosphate rock (PR) needs to double by 2050 compared to the present levels in order to match the total P requirements. The limitation of P is often viewed as a “bottleneck” in the agricultural industry.
Nowadays, investigating supply chain management (SCM) is increasingly encouraged in mining, manufacturing, and agriculture. Optimizing SCM can be described as a strategic management decision considering other decisions in the supply chain. In today’s economy, sustainability is one of the crucial issues [3]. To design a renewable and sustainable P-fertilizer supply chain management (PFSCM) various challenges including PR mining, assessing P-loss throughout the supply chain, and environmental protection issues should be taken into account [4]. Mathematical modeling and optimizing of existing PFSCM strategies are a neglected concern in the context of operational research so far.
Many techniques have been developed to solve complex issues of the SCM problems. Traditional optimization methods inherently suffer from high computation complexity and local optimal stagnation. These methods mainly rely on fitness function’s derivatives to determine the direction to search in to achieve the final solution. Moreover, they are very sensitive to the initial estimate, which may lead to converge into a locally optimal solution. As a result, these techniques are not effective enough to efficiently solve SCM problems. Over the past years, many metaheuristics have been introduced to solve real-world SCM optimization problems. However, their performance is not guaranteed for various SCM problems as they are random search methods. Therefore, there is a need to utilize practical algorithms utilizing the heuristic information available in the problem [5]. Accordingly, hybrid heuristic–metaheuristic techniques have recently become more favored for solving complex optimization problems [6,7,8,9].
Contrary to the previous studies, this paper aims to address a combined three-stage heuristic–metaheuristic approach with global and local search strategies to find a high-speed, high-precision, and high-solution-quality method to solve the sustainable PFSCM problem. In this regard, a problem-dependent heuristic is proposed to generate a set of feasible solutions as the initial population of the population-based metaheuristic, i.e., a whale optimization algorithm (WOA), followed by a single-solution metaheuristic, i.e., a variable neighborhood search (VNS). The key contributions in this paper are as follows:
  • Designing a renewable and sustainable closed-loop supply chain network for the chemical P-fertilizer industry;
  • Utilizing a P recycling policy to reduce economic costs and improve sustainability;
  • Proposing a hybrid three-stage ensemble heuristic–metaheuristic algorithm (named H-WOA-VNS) utilizing heuristic information and global/local search metaheuristics based on WOA and VNS with multiple local search operators for the first time to solve supply chain management problems;
  • Developing a backward knowledge-based heuristic to provide the metaheuristic phases with a set of near-optimal solutions as the start point of searching;
  • Applying VNS on the global best solution found by the WOA;
  • Alleviating the effects of the PFSCM on the environment by promoting P use efficiency and crop yield improvement.
In the rest of this paper, Section 2 studies the literature, including fertilizers supply chains and solution methods. Section 3 is dedicated to sustainable PFSCM modeling. In Section 4, the H-WOA-VNS algorithm is presented. Section 5 provides the case study and discusses the results. Section 6 concludes this paper with future perspectives.

2. Literature Review

In recent years, along with the increasing development in various industries, academic research on optimizing and finding practical solutions to solve the SCM problems also prospered. According to the existing literature, the scope of this paper is related to two main flows of SCM: research related to P-fertilizer supply chains and studies devoted to solution methods, which are described in the following.

2.1. P-Fertilizer Supply Chain

The P-fertilizer industry is associated with the mining and processing of raw materials such as sulfur ores, phosphate ores, potassium salts, etc.; manufacturing fertilizers; and distributing them [10]. Among the essential nutrients for plant growth, P is the pillar of global food security. Despite this fact, PR is a finite and non-renewable mineral resource [11]. Because of the highly dissipative nature of P, improving P use efficiency (PUE) is essential to minimize its effects on aquatic systems and biodiversity [12].
Chemical P-fertilizers include low, medium, and high PR processing. The most widely used low to medium PR processing P-fertilizers are single super phosphate (SSP) and triple super phosphate (TSP), while high PR processing P-fertilizers include mono-ammonium phosphate (MAP), di-ammonium phosphate (DAP), etc. [13]. The excessive application of such P-fertilizers harms the environment and ecosystems by emitting greenhouse gasses (GHG) to the atmosphere, which cannot be neglected in the agricultural industry [14].
Optimizing the P flow from suppliers (PR mines) to consumers (farms) is considered to be the most practical approach to improving the PUE. However, appropriate strategies are lacking to maximize PUE by linking soil–crop systems and fertilizer types with the P flow from a PFSCM perspective. The currently available PR reserves are only available in a few countries. The abundance of currently known PR reserves can avoid supply bottlenecks in the short and medium terms. Over the past years, various policies have been employed to enhance PUE in agriculture industry [15]. Feeding roots rather than soil has been used to exploit soil legacy P to achieve higher PUE and crop yield [16]. Recently, concerns about agricultural sustainability have focused on developing techniques that are effective for farmers but do not have adverse impacts on the environment [17].
The impact of P-fertilizers on PUE was analyzed from a supply chain point of view in different studies [18,19]. The effectiveness of P-fertilizers and the production of achenes considering different levels of P is a strategy to allow its cultivation in soils with different P concentrations without compromising the crop yield or risking the environmental issues associated with the applications of P-fertilizers. The authors in [20] have examined the effects of different levels of P-fertilizers on sunflower regarding the PUE in two agricultural crops. The low availability of P is an important issue limiting sustainable crop production. The performance of wheat production with and without the application of DAP fertilizers was evaluated in [21]. Their results showed that the physiological parameters of the plants were enhanced by utilizing the recommended rate of DAP fertilizers. Moreover, the biochemical properties of wheat grain, including fat, ash, fiber, and protein, were enhanced by 1.24%, 2.26%, 1.47%, and 11.2%, respectively.
To maintain food sustainability, it is essential to ensure the continued supply of PR and find new strategic options that can respond to supply chain problems. To help communities, research on the interactions and trade-offs among the different levels of the supply chain are warranted to aid decision makers. However, improving the P flow among different echelons of the P supply chain and P-fertilizer distribution in such a way that can develop a sustainable PFSCM strategy is still an unclear area in the literature. Promoting PFSCM efficiencies, such as in mining, production, crop yield, distribution, inventory flow, etc., is a hot subject in this field of research. To the best of our knowledge, a comprehensive study on optimizing a closed-loop renewable P supply chain network by considering economic and environmental aspects has yet to be reported elsewhere.

2.2. Solution Methods

There are different methods for the optimization of various SCM problems. Generally, the existing solution search methods can be classified into exact, heuristic, and random (metaheuristic) algorithms. Exact (complete or exhaustive search) techniques are guaranteed to obtain optimal solutions for finite-size combinatorial optimization problems within a finite running time or otherwise prove that no feasible solution exists using a reasonable run-time [22]. The main disadvantage of exact methods is that their computational time grows exponentially with the problem’s size. In this regard, exact methods cannot be applied to solve real-size problems such as sustainable SCM [23]. Moreover, a high-performance exact algorithm for a specific problem is often challenging to extend to another problem as its formulation would be changed [24,25]. In these cases, the optimal solution is sacrificed by finding near-optimal solutions using heuristic or metaheuristic algorithms in a reasonable time [26].
The majority of real-world SCM problems are complex and large in size, so they cannot be optimally solved by exact search methods within an acceptable and reasonable time. Heuristics are an alternative way to find acceptable solutions with a reasonable time. They may also offer an optimal solution in some cases [27]. However, one of their weaknesses is falling into local optima traps and not crossing them. Therefore, metaheuristics have been proposed to effectively deal with such limitations, as the application of these methods in complicated problems is growing extensively.
In recent years, a great deal of effort has been invested into the field of developing metaheuristic algorithms to solve medium- and large-size SCM optimization problems. Their use yields a computationally efficient and convergent procedure for such problems [28]. These methods can be categorized into evolutionary, swarm intelligence, physics-based, and human behavior-based algorithms. The main advantage of such techniques is their utilization of the ‘trial-and-error’ principle in their search process for an optimal solution. Over the past years, metaheuristics have been successfully performed to solve complex problems such as SCM problems [29,30,31,32,33].

2.3. Our Contribution to Existing Methods

The literature proves that hybridizing metaheuristics with other soft computing methods is a suitable way to benefit from the advantages of the basic algorithms [34,35,36,37]. Generally, metaheuristics outperform heuristics in terms of converging to a solution with higher quality. However, they require more computational resources and longer running times than heuristic models. As an appealing solution, heuristic-empowered metaheuristics can achieve a better trade-off between complexity and efficiency using heuristic information in different phases of the metaheuristic algorithm (e.g., initial population generation and population updating) by exploiting the problem-dependent heuristic information available in the problem model [5]. Accordingly, this study focuses on optimizing a sustainable PFSCM by applying a hybrid three-stage algorithm based on heuristic information and two popular metaheuristics (population-based WOA and solution-based VNS) to obtain a high-quality solution while reducing the running time. In this regard, a problem-dependent heuristic is performed to generate a set of near-optimal feasible initial solutions for the first metaheuristic phase (i.e., the WOA phase); then, the second metaheuristic phase (i.e., the VNS phase) is applied to further enhance the global best solution found by WOA through multiple local search operators.

3. Sustainable PFSCM Model

3.1. Supply Chain Network

This paper designs a six-echelon closed-loop PFSCM model comprising P-mines (iI), suppliers of raw materials (jJ), fertilizer manufacturers (mM), distribution centers (dD), farms (fF), and a recycling center, as shown in Figure 1. The model is formulated in T time periods (tT), i.e., months. The raw materials (rR) include sulfuric acid (SA), phosphoric acid (PA), and ammonium (A), where PA is made from P and SA. Moreover, four types of products (pP) including low-to-medium-grade P-fertilizers (SSP and TSP) and high-grade P-fertilizers (MAP and DAP), are considered. In every month t, each P-mine i or supplier j may supply raw materials to one or more manufacturers up to its capacity. Each manufacturer m may be sourced by one, two, or more P-mines and suppliers. Each distribution center d can purchase the required fertilizers from various manufacturers until its total demand is satisfied. Then, each distributor delivers the received SSP, TSP, MAP, and DAP fertilizers to the corresponding farms (i.e., demand nodes). Every farm f has a demand DPft of the total P-uptake by different fertilizers. Moreover, the P-uptake from each fertilizer p should be within [LPfpt,UPfpt]. At the end of every time period t, P-leaching from the different farms is collected and transferred to the recycling center to recycle them and provide phosphate, which can be used as a P supplier in the next time period t+1.
The following assumptions are considered in the proposed PFSCM model:
  • P is the main resource of producing P-fertilizers.
  • There are three raw materials: SA, PA, and A.
  • There are four P-fertilizer products: SSP, TSP, MAP, and DAP.
  • To produce a unit of SSP, U 1 P = 0.64 units of P and U 31 S = 0.37 units of SA are used.
  • To produce a unit of TSP, U 2 P = 0.4 units of P and U 22 S = 0.34 units of PA are used.
  • To produce a unit of MAP, U 13 S = 0.12 units of A and U 23 S = 0.51 units of PA are used.
  • To produce a DAP, U 14 S = 0.23 units of A and U 24 S = 0.47 units of PA are consumed.
  • PA is considered as a raw material ignoring the production procedure.
  • Each manufacturer has a limited capacity to store P-fertilizers.
  • Distribution centers can store a part of the received P-fertilizers.
  • Lead time of P-recycling is assumed to be one month: the collected P-leaching from farms in each month would be recycled and can be used in the next month.
  • Each mine has a limited capacity to provide P each month.
  • Each supplier has a specific capacity to provide raw material each month.
  • Each manufacturer may produce one, two, three, or four P-fertilizer types, each with a limited capacity.
  • The demand of each farm in terms of the total P-uptake and minimum/maximum required level of each type of fertilizer is assumed to be known for every month.

3.2. Notations

The notations of the PFSCM model are provided in Table 1. To solve the model via H-WOA-VNS and empower it to satisfy the problem constraints, we distinguish two types of decision variables (DVs): direct DVs and indirect DVs. The direct DVs are those which are encoded into feasible solutions and optimized using H-WOA-VNS, while the indirect DVs are decoded based on the model parameters and the direct DVs.

3.3. Problem Formulation

In the following, the PFSCM model is formulated as a multi-objective problem to minimize the total economic cost, while maximizing the crop yield and PUE. By using mixed integer linear programming (MILP), a mathematical model of the PFSCM model is provided.

3.3.1. Economic Objective Function

The total economic cost comprises the purchasing cost of P from P-mines (CPP), the purchasing cost of raw materials from different suppliers (CPS), the production cost of the manufacturers (CPR), the recycling cost in the recycling center (CRC), the inventory cost (CIH), and the transportation cost (CTR), which are addressed in Equations (1)–(6), respectively:
C P P = i B C i P t m Y i m t P X i m t P
C P S = j r B C j r S t m Y j m r t S X j m r t S
C P R = m p P C m p t P m p t
C R C = t R C f p μ p Y f p t F X f p t F
C I H = t m p I C m M I P m p t + t d p I C d D I D d p t + t I R R I R t
C T R = t i m T C i m P M Y i m t P X i m t P + t m T C m R M Y m t R X m t R + t j m r T C j m S M Y j m r t S X j m r t S + t m d p T C m d M D Y m d p t M X m d p t M + t d f p T C d f D F Y d f p t D X d f p t D + t f p T C f F R Y f p t F X f p t F
By aggregating the six economic costs as mentioned above, the total economic objective function can be expressed as:
minimize   Z E C = C P P + C P S + C P R + C R C + C I H + C T R = i B C i P t m Y i m t P X i m t P + j r B C j r S t m Y j m r t S X j m r t S + m p P C m p t P m p t + t R C f p μ p Y f p t F X f p t F + t m p I C m M I P m p t + t d p I C d D I D d p t + t I R R I R t + ( t i m T C i m P M Y i m t P X i m t P + t m T C m R M Y m t R X m t R + t j m r T C j m S M Y j m r t S X j m r t S + t m d p T C m d M D Y m d p t M X m d p t M + t d f p T C d f D F Y d f p t D X d f p t D + t f p T C f F R Y f p t F X f p t F )

3.3.2. Environmental Objective Functions

In the proposed model, two environmental objectives are used to increase the crop yield (ZCY) and PUE (ZPUE), which are formulated in Equations (8) and (9), respectively. As P is a non-sustainable resource, it is of utmost importance to reduce the total P-loss, which is described as the total losses from the exploited PR to that of uptake by the plants:
maximize   Z C Y = f p β f p t d Y d f p t D X d f p t D
maximize   Z P U E = p 1 λ p f α f p d t Y d f p t D X d f p t D + γ f p d t Y d f p t D X d f p t D

3.3.3. Multi-Objective Function

To solve the multi-objective problem, the objective functions ZEC, ZCY, and ZPUE are aggregated into a single-objective formula by means of a weighted averaging method. The economic objective ZEC is calculated in USD to be minimized, while environmental objectives are to be maximized. To be able to combine the different objective functions, all objective functions are normalized into Z E C ¯ , Z C Y ¯ , and Z P U E ¯ to be minimized within [0,1] considering the minimum and maximum expected values for each objective function. As a result, the total objective function can be expressed as:
minimize   O F = w E C Z E C ¯ + w C Y 1 Z C Y ¯ + w P U E 1 Z P U E ¯ × P F
where wEC, wCY, and wPUE (wEC + wCY + wPUE = 1) are constant parameters which determine the relative impact of ZEC, ZCY, and ZPUE, in the objective function. Moreover, PF is a penalty function calculated as PF = 2NumP, where NumP is the number of constraints in Equations (11)–(34) which have not been satisfied.
Subject to:
p a f p d Y d f p t D X d f p t D B f d D P f t               f , t
L P f p t d Y d f p t D X d f p t D B f d U P f p t               f , p , t
d Y d f p t D X d f p t D B f d = R f p t B S p               f , p , t
X d f p t D B f d               d , f , p , t
d X d f p t D = 1               f , p , t
I D d p 0 = 0               d , p
I D d p t 1 + m Y m d p t M X m d p t M f Y d f p t D X d f p t D B f d               d , p , t
I D d p t = I D d p t 1 + m Y m d p t M X m d p t M f Y d f p t D X d f p t D B f d               d , p , t
p I D d p t W d D               d , t
P m p t C a p m p M                 m , p , t
i Y i m t P X i m t P + Y m t R X m t R p U p P P m p t               m , t
j Y j m r t S X j m r t S p U r p S P m p t               m , r , t
I P m p 0 = I m p 0               m , p
I P m p t 1 + P m p t d Y m d p t M X m d p t M               m , p , t
I P m p t = I P m p t 1 + P m p t d Y m d p t M X m d p t M               m , p , t
p I P m p t W m M               m , t
m Y j m r t S X j m r t S C a p j r t S             j , r , t
m Y i m t P X i m t P C a p i t P               i , t
f p Y f p t F X f p t F d f p γ f p μ p Y d f p z D X d f p z D B f d               t
m Y m t R X m t R z = 1 t 1 f p μ p Y f p t F X f p t F z = 1 t 1 m Y m z R X m z R               t
I R 0 = 0
m Y m t R X m t R I R t 1               t
I R t = I R t 1 m Y m t R X m t R + f p μ p Y f p t F X f p t F               t
I R t W R               t
Constraints (11) and (12) are related to farms. Constraint (11) expresses the P-uptake demand of farm f at each time period t. In other words, the monthly demand of each farm must be met by the P-uptake from all P-fertilizers received from the different distribution centers. Constraint (12) expresses the boundaries of the required amounts of fertilizers by each demand node. In fact, the amounts of fertilizers purchased from all distributors to farm f should be between the minimum and maximum amounts of the demand for each fertilizer p.
Constraints (13–19) describe the constraints of the distribution centers. Constraint (13) shows that the amount of fertilizer p delivered to farm f at time t is equal to the summation of the delivered fertilizer p from different distributors. Constraint (14) indicates whether farm f is supported by distributor d to purchase fertilizer p or not. Constraint (15) ensures that farm f is sourced for fertilizer p by one distributor at each period t. Constraint (16) expresses that the initial inventory of distributor d for fertilizer p is zero. Constraint (17) illustrates that the total inventory and purchased product items of fertilizer p through all manufacturers by distribution center d at time t must be at least equal to under-support farms’ demands. Constraint (18) calculates the total inventory of fertilizer p in distribution center d at the end of time t. Moreover, constraint (19) ensures that the total inventory of all types of fertilizers held by distributor d at time t should not exceed the corresponding warehouse capacity.
Constraints (20–26) are related to the manufacturers. Constraint (20) indicates that the produced fertilizer p by manufacturer j at time t must not exceed a specified capacity. Constraint (21) indicates that the total amounts of carried PR from the mines and carried P from recycling center to manufacturer m at time t must at least fulfill that manufacturer P requirement to satisfy its demands. Constraint (22) addresses that total raw material r transferred from all suppliers to manufacturer m at time t must at least fulfill the required raw materials for the production of all types of P-fertilizers. Constraint (23) describes the initial inventories of the different manufacturers for each P-fertilizer. Constraint (24) ensures that the total inventory of fertilizer p of manufacturer m at time t satisfies the requirements of its under-support distribution centers. Constraint (25) updates the inventory level of the manufacturer m for fertilizer p at the end of time t. Constraint (26) indicates that the total inventory of all fertilizers at time t do not exceed the capacity of manufacturer m.
Constraints (27) and (28) address the supplying capacity of the P and raw materials. Constraint (27) expresses that the amount of raw material r carried from the supplier j to manufacturer m at time t must not exceed the production capacity of supplier j. Constraint (28) indicates that the amount of carried P from the P-mine i to manufacturer m at time t does not exceed the capacity of P-mine i.
Constraints (29–34) are related to the recycling center. Constraint (29) calculates the total collected P obtained by P leach of all fertilizers in all farms to the recycling center at period t. Constraint (30) expresses the total amount of manufacturers’ requirements sourced through utilizing recycled P by the recycling center at period t. Constraint (31) denotes that the initial inventory of the recycling center is zero. Constraint (32) limits the inventory of recycled P at the recycling center at the end of every period. Constraint (33) expresses the updated inventory of the recycling center at time t. Constraint (34) ensures that the total inventory held at the recycling center at time t must not be higher than its warehouse capacity.

4. Solution Method: H-WOA-VNS

Generally, SCM is an NP-hard combinatorial optimization problem [38,39,40], and thus, implementing exact methods with an exhaustive search strategy cannot be applicable due to the required computational resources for real-size SCMs. In this case, heuristic and/or metaheuristic algorithms can be adopted. Heuristics are problem-dependent techniques which are specifically designed for solving a problem utilizing the available information in the problem model [41]. These algorithms are rapid and easy to use; however, they do not effectively investigate the whole search space. In contrast, metaheuristic algorithms are higher-level, iterative-based random search methods which can obtain a better solution quality by allowing for a greater running time [42].
Based on the number of solutions encountered in each iteration, metaheuristics can be categorized into population-based (global search) and solution-based (local search) methods [43]. Population-based metaheuristics investigate the search space in parallel (i.e., more exploration), while solution-based metaheuristics attempt to locally improve the quality of the current solution using local search operators (i.e., more exploitation). To achieve a better trade-off between solution quality and speed, we propose a combined three-stages heuristic–metaheuristic approach with balanced exploration–exploitation strategies to efficiently solve the PFSCM model. In the first stage, a heuristic algorithm is designed to generate feasible solutions while satisfying all constraints. These solutions are fed into the population-based metaheuristic based on WOA; then, a solution-based metaheuristic based on VNS is performed to further improve the best solution of WOA through local search operators.
The overall flowchart of the proposed combined H-WOA-VNS algorithm is shown in Figure 2. Because a knowledge-based heuristic algorithm is used, the search process of the metaheuristic algorithms starts from a set of near-optimal feasible solutions rather than starting from random solutions, and thus, fewer iterations are required to effectively search among the search space. Moreover, because the exploitation-specific local search VNS algorithm is applied, the global best solution of the WOA can be further improved via the different local search operators stored in the VNS. For more insight into the proposed H-WOA-VNS algorithm, its pseudo-code is provided in Algorithm 1.
Algorithm 1. Proposed H-WOA-VNS algorithm for solving the PFSCM problem.
Input:
   Model parameters of the PFSCM model
Output:
   GBestSOL: optimized solution for the PFSCM model
Heuristic Phase:
1.  for (s ≤ PopSizeWOA)
2.     Generate a feasible solution SOL(s) using the heuristic algorithm
3.     Add SOL(s) into a population set, s={1,2,…,PopSizeWOA}
4.  end for
WOA Phase:
1.  Considering the heuristic solutions as initial population of WOA
2.  for (s ≤ PopSizeWOA)
3.     Evaluate the quality of SOL(s) using OF according to Equation (10)
4.  end for
5.  for1 (n ≤ MaxIterWOA)
6.     for2 (s ≤ PopSizeWOA)
7.       Update the values of p, l, and a , and A
8.       if1 (p < 0.5)
9.          if2 (| A | ≥ 1)
10.            Update solution s via search for prey as Equation (45)
11.         else if2
12.            Update solution s via encircling prey using Equation (46)
13.         end if2
14.       else if1
15.         Update solution s via bubble-net attacking using Equation (47)
16.       end if1
17.       Revise solution s, if it goes beyond the search space
18.       Evaluate the quality of SOL(s) using OF using Equation (10)
19.     Updating of the global best solution: GBestSOL
20.  end for1
VNS Phase:
1.  Considering GBestSOL as initial solution of VNS: SOLcurrent
2.  for1 (n ≤ MaxIterVNS)
3.     for2 (l ≤ NumLSVNS)
4.       Generate SOLnew in vicinity of SOLcurrent via local search operator l
5.       Evaluate the quality of SOLnew using OF using Equation (10)
6.       if MoFnew<MoFcurrent
7.         Replace SOLcurrent with SOLnew
8.         Replace MoFcurrent with MoFnew
9.         Break for2
10.       end if
11.     end for2
12.     Updating of the global best solution: GBestSOL
13.  end for1
Output: Return the GBestSOL as the final optimized PFSCM model

4.1. Solution Representation

As mentioned above, we use direct and indirect DVs to encode and decode feasible solutions, respectively. Accordingly, a solution is represented via the direct DVs, while it is evaluated by decoding it and extracting the indirect DVs.

4.1.1. Solution Encoding

An example for the encoding of the direct DVs to solve the PFSCM model is shown in Figure 3. A solution, SOL, can be encoded via the direct DVs as multiple integer and permutation matrices as follows:
  • SOL.P ( P m p t ) is an integer matrix of dimension P×M×T, which determines the amount of fertilizer p which is produced in manufacturer m at time t (per ton).
  • SOL.R ( R f p t ) is an integer matrix of dimension P×F×T to determine the amount of each fertilizer p delivered from the corresponding distribution center to each farm f at every time period t (per BSp).
  • SOL.PLM ( P L M ) is a permutation of M manufacturers to determine the priority of the manufacturers to order from different P-mines and suppliers.
  • SOL.SLP ( S L P m ) is a matrix of dimension M×(I + 1) comprising M permutation vectors of I + 1 P suppliers, i.e., each row m is a permutation of I mines plus the recycling center, which determines the selection list of the different P-mines for the manufacturer m.
  • SOL.SLS ( S L S m ) is a matrix of dimension M×J, where each row m specifies selection list of the different suppliers to supply raw materials for the manufacturer m.
  • SOL.PLD ( P L D ) is a permutation of D distribution centers which determines priority of them in ordering fertilizers from the different manufacturers.
  • SOL.SLM ( S L M d ) is a matrix of dimension D×M, where each row d is a permutation of M manufacturers to determine the selection list of the different manufacturers for the distributor d.

4.1.2. Solution Decoding

To evaluate a feasible solution, SOL, its direct DVs must be decoded into the indirect DVs, including S P m r p t , S m r t , R P t ,   T D d p t ,   X i m t P , X m t R , X j m r t S , X m d p t M , X d f p t D , X f p t F , Y i m t P , Y m t R , Y j m r t S , Y m d p t M , Y d f p t D , Y f p t F , I P m p t , I D d p t , and I R t , based on the direct DVs and other model parameters. To achieve this purpose, a backward mechanism from the demand nodes to the suppliers is applied to determine the value of the indirect DVs, as follows:
  • Decoding of farms-distribution centers DVs: in the first step, X d f p t D and Y d f p t D are obtained according to the requested demands stored in R f p t as:
    X d f p t D = 1                   i f   B f d = 1   a n d   R f p t > 0 0                   e l s e                                                                                               p , d , f , t
    Y d f p t D = R f p t B S p                 i f   B f d = 1 0                                           e l s e                                           p , d , f , t
  • Decoding of farms-recycling center DVs: at the end of each time period, after delivering fertilizers Y d f p t D to the farm f and the irrigation procedure, the amount of P-leaching for fertilizer p can be estimated by multiplying the total recyclable fertilizers to the P-leaching coefficient of fertilizer p, according to Equation (37). Then, the total amount of recycled phosphate at the end of time t is achieved according to Equation (39). It should be emphasized that the lead time of P recycling is one month, and consequently, the recycled P at time t can be transferred to the manufacturers in times t+1, t+2, and so on. To this end, the values of X f p t F , Y f p t F , and R P t , are extracted:
    Y f p t F = d γ f p μ p Y d f p t D X d f p t D B f d                 p , f , t
    X f p t F = 1                     i f   Y f p t F > 0   0                     e l s e                                             p , f , t
    R P t = f p μ p Y f p t F X f p t F = f p d γ f p Y d f p t D X d f p t D B f d                     t
  • Decoding of distribution centers–manufacturers DVs: Based on the received demands by the farms, the total demand of each distribution center d for each fertilizer p at every time period t is calculated according to Equation (40). Then, the different distributors are given one by one according to their priorities in P L D . For each distribution center d, different manufacturers according to their selection list in S L M d are evaluated one by one to deliver fertilizers to the distributor d until satisfying its total demand T D d p t . As a result, the values of T D d p t , X m d p t M , and Y m d p t M are extracted:
    T D d p t = f Y d f p t D X d f p t D B f d                     p , d , t
  • Decoding of manufacturers–raw material suppliers DVs: Every time t, the demand of each raw material r for each manufacturer m for the production of fertilizer p can be calculated by multiplying U r p S to the amount of the produced fertilizer p according to Equation (41). Accordingly, the total demand of raw material r for each manufacturer m at every time t is obtained by Equation (42). Then, to satisfy the demand of the manufacturers, they are evaluated one by one based on their priorities in P L M . The manufacturer m may be sourced by different suppliers according to their selection list in S L S m until satisfying the total demand of the manufacturer m for raw materials. As a result, S m r p t S , S T m r t S , X j m r t S , and Y j m r t S , are obtained:
    S m r p t S = U r p S P m p t                         r , p , m , t
    S T m r t S = p S m r p t S                         r , m , t
  • Decoding of manufacturers–phosphate suppliers DVs: The phosphate demand of each manufacturer m to produce the fertilizer p at every time period t can be expressed as in Equation (43), and then, the total phosphate demand for each manufacturer m at every time t is calculated according to Equation (44). Similar to the supplier selection, to satisfy the P demand of the manufacturer m, different P-mines as well as the recycling center (I + 1 phosphate suppliers) are evaluated one by one based on their selection list in S L P m until satisfying the phosphate demand of the manufacturer m. To this end, S m p t P , S T m t P , X i m t P , X m t R , Y i m t P , and Y m t R are achieved:
    S m r p t S = U r p S P m p t                         r , p , m , t
    S T m r t S = p S m r p t S                         r , m , t
  • Decoding of inventories: At the end of each time period t, the inventory of distribution centers, manufacturers, and recycling center, are updated according to Equations (18), (25) and (33), respectively.

4.2. Initialization Using a Heuristic Algorithm

The solution is generated via the proposed heuristic algorithm using a backward flow from the demand nodes (farms) to the suppliers, as follows:
  • At first, SOL.R is constructed using the P-uptake demand and minimum/maximum requested fertilizers by the farms. More specifically, the demand of farm f for the fertilizer p at every time period t, R f p t , is randomly generated within [LPfpt,UPfpt] to satisfy Equation (12), while the total P-uptake from the different fertilizers fulfills the total demand given in Equation (11).
  • The total demand of each distribution center for all fertilizers at all time periods is calculated by the summation of the demands of the corresponding farms. Then, the priority list of the different distribution centers (SOL.PLD) is obtained by sorting them from the highest demand to the lowest demand.
  • For each distribution center d, different manufacturers are sorted according to the total cost per unit of fertilizers (including purchasing and transportation costs), from the lowest to the highest cost. This procedure is repeated for all distribution centers to construct SOL.SLM.
  • The amount of fertilizer p produced by manufacturer m at time t is considered to be a random value within [0.5 × C a p m p M , C a p m p M ] to construct SOL.P.
  • The priority list of the manufacturers (SOL.PLM) is determined by sorting them from the most significant to the least significant according to their total production at all time periods.
  • For each manufacturer m, the different raw material suppliers are sorted according to the total cost of purchasing and transportation from the lowest to the highest cost, and eventually, SOL.SLS is obtained. The same procedure is repeated for the I + 1 phosphate suppliers (P-mines plus the recycling center) to construct SOL.SLP.

4.3. Global Search Using WOA

WOA is a swarm intelligence algorithm introduced by Mirjajili and Lewis [44], which mimics the hunting behavior of whales. In the proposed algorithm, the feasible solutions generated by the heuristic algorithm are considered as the initial population of the WOA. Then, at every iteration of the WOA, the quality of the current solutions is evaluated by the OF using Equation (10), and eventually, the entire population is updated using search for prey, encircling prey, and bubble-net attack. To update the position of each whale w, a uniform random number p in [0,1] is generated. If p ≥ 0.5, the solution is updated using the bubble-net attack. Otherwise, a vector A is generated as A = 2 a · r a , where the components of a are linearly decreased from 2 to 0 during the course of the iterations, and r is a vector whose elements are randomly generated within [0,1]. If | A | > 1, the search for prey is applied; otherwise, the solution is subjected to be updated via the encircling prey. In the following, these updating operators are illustrated.

4.3.1. Search for Prey

If | A |≥1, the whale is moved toward a random solution, which emphasizes more exploration (global search). Based on the search for prey, the whale w is updated as:
SOL w t + 1 = SOL r a n d t A · 2 r · SOL r a n d t SOL w t

4.3.2. Encircling Prey

In each iteration, it is assumed that the best solution found so far (GBestSOL) is in the vicinity of the global best solution (i.e., an optimal solution). Accordingly, if | A |<1, the whale w tries to update its position toward the GBestSOL, as follows:
SOL w t + 1 = GBestSOL A · 2 r · GBestSOL SOL w t

4.3.3. Bubble-Net Attack

The bubble-net attacking behavior can be modelled as a spiral equation between the current position of whale w and GBestSOL, which simulates the helix-shaped movement of the whales. It can be formulated as:
SOL w t + 1 = GBestSOL SOL w t · e b l cos ( 2 π l ) + GBestSOL
where l is a random number within [−1,1], and b is a constant which defines the logarithmic spiral shape.

4.4. Local Search Using VNS

After termination of the WOA phase, its final solution is used as the initial solution of VNS, i.e., SOLcurrent = GBestSOL. In each iteration of VNS, a new solution, namely SOLnew, is constructed in the vicinity of the current solution, SOLcurrent, using multiple local search operators. If the quality of the new solution has been enhanced, the current solution is replaced with the new one; otherwise, it is not changed. In the case of integer matrices SOL.P and SOL.R, either Integer Swap (I-Swap) or Integer Exchange (I-Exchange) may be applied, as illustrated in Figure 4. Moreover, one of the permutation operators, including Exchange, Relocate, OrOpt, TwoOpt, and Reverse, may be performed on the permutation matrices. These operators can be seen in Figure 5.

5. Experimental Results

The proposed PFSCM model and H-WOA-VNS algorithm have been coded in MATLAB R2020b. All experiments have been conducted on a personal computer with 2.59 GHz i7 CPU and 16 GB RAM running on windows 10.

5.1. Case Study

To validate the H-WOA-VNS algorithm for solving the PFSCM model, a real case study based on the collected data from Iran has been used. Significant lands available for cultivation in Iran have poor soil fertility. Due to the availability of natural gas reserves in Iran, there are many facilities for the production of the required raw materials for the P-fertilizers. However, Iran has few P-mines, which results in a shortage of the required P to produce P-fertilizers.
The distribution centers are located in 32 provinces in Iran. The annual demands of different distribution centers in terms of the total P-uptake, as well as the lower and upper bounds of the required P-fertilizers, are reported in Table 2. Four types of P-fertilizers in Iranian agriculture are used, including SSP, TSP, MAP, and DAP. The composition of the main raw materials (P, A, SA, and PA) required to produce the four P-fertilizers (SSP, TSP, MAP, and DAP) gathered from experts in Iranian manufacturers are provided in Table 3. Moreover, the number of suppliers of raw material, as well as their supplying capacities, are reported in Table 4, where the P suppliers include two internal P-mines in Iran, two external suppliers in Iraq and Syria, and one recycling center in Iran. According to Table 2, Table 3 and Table 4, the required A, SA, and PA, can be satisfied via internal suppliers. However, the main problem in the Iranian P-fertilizers industry is that only about 30–40% of the need for P can be met through the domestic sources (including two internal P suppliers and one recycling center). Therefore, the remaining requirement for P may be imported from the two external suppliers, i.e., other countries. However, there are more than 120 manufacturers of P-fertilizers in Iran, with different production capacities. In this paper, we have chosen the 21 biggest P-fertilizers manufacturers in Iran with more than 2000 tons of production (monthly). The production capacity of the manufacturers is addressed in Table 5.

5.2. Settings

To set the parameters of the H-WOA-VNS algorithm, different options have been evaluated, and the best values and operators have been set for the final simulations. The parameters of H-WOA-VNS are provided in Table 6. The number of iterations and population size in WOA are considered as 500 and 70, respectively. To achieve the same number of objective function evaluations (NFE) in both the WOA and VNS phases, the number of iterations in VNS has been set 10 times larger than that of in the WOA phase, i.e., to obtain MaxIterVNS × NumLSVNS = MaxIterWOA × PopSizeWOA. The weights of the three objective functions ZEC, ZCY, and ZPUE, for the reported results in Section 5.3 and Section 5.4 have been considered as 0.5, 0.3, and 0.2, respectively. However, these weights will be changed in Section 5.5 to evaluate their effects on the different objective functions through a sensitivity analysis. The PFSCM model has been discussed and solved in T = 12 time periods (months), i.e., one year. The other parameters of the model have been set according to the case study data, as summarized in Table 7.

5.3. Results

5.3.1. Optimization Results

Because of the random-based nature of the H-WOA-VNS algorithm (as well as other metaheuristics), it was performed in 10 successive runs for solving the PFSCM model. Considering the weights of the multi-objective function in Equation (10) as wEC = 0.5, wCY = 0.3, and wPUE = 0.2, the obtained results for 10 runs are reported in Table 8, including the economic cost (ZEC), the average crop yield (ZCY), the average PUE (ZPUE), the penalty function (PF), and the total objective function (OF). It can be seen in Table 8 that the difference between the worst and best solutions is 0.57% (OFworst − OFbest = 0.0028). Moreover, the standard deviation (STD%) of the obtained OF in 10 runs is 0.19%. The results indicate that there is no significant difference between the results of H-WOA-VNS in different runs, which shows a high reliability in a single run.
The convergence of the H-WOA-VNS algorithm (on average over 10 runs) can be seen in Figure 6 and Figure 7. Figure 6 illustrates the best OF versus iterations in WOA and VNS phases. To gain more insights about the convergence of the algorithm, Figure 7 merges the two phases into a single diagram, representing the best OF versus NFE. The H-WOA-VNS algorithm starts by generating an initial population for WOA via the Heuristic algorithm, wherein the best solution has obtained an OF = 0.5491. In the beginning of the WOA phase, the algorithm’s convergence is sharp, and the best OF decreases from about 0.55 to 0.51 in only 10% of the early iterations (50 out of all 500 iterations). This convergence speed is not only because of the nature of population-based metaheuristics, but also due to a set of near-optimal heuristic solutions being generated in different positions of the whole search space. As WOA progresses, its convergence speed gradually decreases, and finally, it converges to OF = 0.5014. Then, by calling VNS to improve the global best solution of WOA through multiple local-search operators, a sudden shock in the convergence speed can be observed at early iterations of the VNS phase. It efficiently helps the H-WOA-VNS algorithm to further enhance the global best solution found by WOA, resulting in a 0.006 reduction in the OF and obtaining the final solution with OF = 0.4954.

5.3.2. PFSCM Results

In the following, the results of the best solution among 10 runs (Run 4 in Table 8) with OF = 0.4941, are reported. The total economic cost 446.38 M$ includes six sub-costs, as summarized in Table 9. The total productions of the P-fertilizers in all time periods for each manufacturer are provided in Table 10. According to Table 10, a total of 346,001 tons SSP, 334,133 tons TSP, 229,976 tons MAP, and 185,799 tons DAP, have been produced by all manufacturers. The required raw materials supplied by the P-mines and raw material suppliers to produce the mentioned P-fertilizers are provided in Table 11. Considering the total P-fertilizers in the PFSCM for all time periods (sum of the initial inventories and the total production in all manufacturers), there are 416,001 tons SSP, 404,133 tons TSP, 279,976 tons MAP, and 255,799 tons DAP, which can be delivered to the distributors or held in the manufacturers’ warehouses. The total delivered P-fertilizers to each distributor can be seen in Table 12. Among the total P-fertilizers in all time periods, 409,332 tons SSP, 394,231 tons TSP, 274,317 tons MAP, and 248,570 tons DAP have been delivered to all distribution centers, while 6,669 tons SSP, 9,902 tons TSP, 5,659 tons MAP, and 7,229 tons DAP have been held in the manufacturers’ warehouses at the end of the 12th time period. Moreover, Table 13 summarizes initial inventories, total productions, total P-fertilizers, total delivered P-fertilizers, and final inventories.

5.4. Validation

In this section, the H-WOA-VNS algorithm is validated by comparing its results on some test problems against the exact method (Section 5.4.1). Then, the performance of H-WOA-VNS is justified by comparing it with its three components including the Heuristic, WOA, and VNS methods and with three existing metaheuristics (Section 5.4.2).

5.4.1. Validation of H-WOA-VNS against the Exact Method

To validate the optimality of the H-WOA-VNS algorithm, its results on different PFSCM test problems as well as on the real case study are compared with an exact search. As mentioned above, the PFSCM (like other combinatorial SCMs), is an NP-hard problem, and thus, an exact method with an exhaustive search strategy is not applicable in terms of computational time complexity for the real-world SCM problems. However, to justify the H-WOA-VNS algorithm against the exact search method, we applied both techniques on five synthetic datasets (SDs) with small- and medium-sizes. The details of the synthetic and real PFSCM test problems are summarized in Table 14.
The obtained results of performing the exact method as well as the H-WOA-VNS algorithm for solving the different PFSCM problems in 10 successive runs are summarized in Table 15. According to the obtained results, the running time of the exact algorithm exponentially grows with the size of the problem, and thus, it cannot solve medium- and large-size test problems (SD 5 and Case Study) in a reasonable time and faced with a low memory error. However, running time of the H-WOA-VNS algorithm increases almost linearly with the problem size, as the NFE of the algorithm is not changed and only the size of feasible solutions is increased. In terms of the optimization results, the obtained OF using H-WOA-VNS has a deviation of 0% to 0.13% from the optimal solution in four small-size datasets, which ensures high-quality near-optimal solutions are achieved for the larger test problems such as the Case Study.

5.4.2. Validation of H-WOA-VNS against Other Heuristics and Metaheuristics

To find the effect of the different stages of the H-WOA-VNS, it was compared with its individual components (Heuristic, WOA, and VNS) when applied separately to solve the PFSCM model. For a fair comparison between the different methods, the number of iterations in WOA and VNS have been set to twice of those in the combined H-WOA-VNS algorithm, i.e., MaxIterWOA = 1000 and MaxIterVNS = 10,000, to obtain the same NFE = 70,000 for all algorithms. Because of the random-based nature of the algorithms, each method was applied in 10 runs. The convergence of the different techniques in terms of the best OF versus NFE is shown in Figure 8. Although WOA and VNS have a good convergence in the early 10% NFEs, they often become trapped in local optima points. To compare the proposed combined H-WOA-VNS algorithm (comprising both population-based and solution-based metaheuristics) against the existing methods in the literature for SCM, a population-based metaheuristic using a genetic algorithm (GA) [45], a solution-based metaheuristic using simulated annealing (SA) [46], and a combined metaheuristic based on GA and SA (GLGASA) [39] have also been used for the optimization of PFSCM considering our Case Study dataset. The OF obtained by different techniques are provided in Table 16. The results demonstrate the effectiveness of the H-WOA-VNS algorithm against other methods, in terms of the best and average results over 10 runs. Furthermore, the proposed algorithm obtains less STD%, which shows a higher reliability in a single run.

5.5. Discussion

Generally, decision makers may consider different preferences in designing supply chains, depending on the importance of different objectives. In the design of the PFSCM model in this paper, we have set the weights of the multi-objective function in Equation (10) as wEC = 0.5, wCY = 0.3, and wPUE = 0.2. By changing the relative impacts of the objective functions between 0 and 1 (while their summation remains fixed equal to 1), we can find the effect of the different weights on the sub-objectives and total objective function. Table 17 reports some samples of the sensitivity analysis, where the first row reports the default weights. For each sample, the model has been solved in 10 successive runs, and the average results are reported in Table 17. Based on the obtained results, the first environmental objective function ZCY has a huge conflict with the economic objective function. It aims to improve crop yield, which needs more high PR processing P-fertilizers (MAP and DAP) than the low to medium PR processing P-fertilizers (SSP and TSP). It consequently increases the total economic cost, as high PR processing P-fertilizers have more purchasing and production costs. However, the second environmental objective ZPUE has less conflict with the total economic cost. To evaluate the effect of changing the weights of the different objectives in Equation (10), we have changed the two environmental weights wCY and wPUE from 0 to 0.5 in steps of 0.1. Eventually, the weight of the economic objective function is considered as wEC = 1 − (wCY + wPUE). By applying the H-WOA-VNS algorithm to solve the model for each combination, the normalized economic cost ( Z E C ¯ ) versus wCY and wPUE can be illustrated in Figure 9. The figure confirms more conflict of the wCY on the total economic costs.
To give an insight into the improved effectiveness of the proposed H-WOA-VNS algorithm over its components (Heuristic, WOA, and VNS) and the existing metaheuristics (GA, SA, and GLGASA), the worst, mean, and best OFs obtained by the different methods over 10 successive runs are compared in Figure 10. Moreover, improvement % of H-WOA-VNS against other techniques is provided in Table 18. The results demonstrate the superiority of the proposed algorithm against the compared methods. It not only achieves the best results among all techniques, but also, more importantly, it has a much smaller deviation than the other methods. The STD% of the H-WOA-VNS algorithm over 10 runs is only 0.19%, while it varies from 1.45% to 2.33% for the other methods.

5.6. Managerial Insights

In recent years, efficiently managing soil resources is a strategic approach taken by decision makers to maintain global food security. In this regard, the optimization of soil nutrition supply chain networks to increase the fertility of the agricultural productions is of the utmost importance and is critical for administrators in agricultural industries. Generally, phosphorus, the main raw material used in the production of P-fertilizers, is a non-substitutable nutrient. Therefore, it is necessary to design a renewable and sustainable supply chain for the management of P-fertilizer production to deliver the required P-fertilizers to crop production systems and farms. Although the mathematical modeling and optimization of the PFSCM network is of utmost importance to enhance crop yield and PUE in achieving global food security, it has not received enough attention thus far by researchers from an operations research point of view. To address these gaps, we have examined the problem of PFSCM network design by considering the phosphorus renewability and sustainability issues in terms of economic and environmental concerns. The proposed model can help decision makers and agricultural administrators to develop proper strategies and make better decisions to achieve a better trade-off between the different objectives.

6. Conclusions

This study has introduced a new model for designing a renewable and sustainable supply chain in the P-fertilizer industry, and the model has been validated using a real dataset in Iran. To solve the established model, a combined heuristic–metaheuristic algorithm with global and local search strategies, named H-WOA-VNS, has been introduced. In the first stage, the model performs a problem-dependent heuristic to generate an initial population in order to guide the metaheuristic algorithm to start from a set of near-optimal solutions rather than starting from fully random solutions. At the second stage, a population-based metaheuristic with exploration and exploitation strategies is used to globally investigate the search space, and finally, an exploitation-oriented metaheuristic with multiple local search operators is used to further improve the quality of the global best solution found by the population-based metaheuristic algorithm.
Experimental results for five synthetic datasets and one real case study have shown the validity and effectiveness of the proposed PFSCM model and H-WOA-VNS solution algorithm. A comparison of the obtained results with other heuristic and metaheuristic methods demonstrated the effectiveness of the H-WOA-VNS algorithm to solve the sustainable PFSCM model, which can also be applied to other SCM models. Due to a lack of attention in the literature to the PFSCM from an operations research point of view, future research is required by focusing on mathematical modeling and optimizing techniques for sustainable PFSCM by utilizing other modeling approaches and solutions. Moreover, the uncertainties occurring in the system (parameters, objectives, and constraints), as well as the disruptions, are also important issues which could be under consideration in future works. In this paper, controllable parameters of the H-WOA-VNS algorithm have been set by trial-and-error. As a future work, design of experiment (DOE) techniques, such as the Taguchi method, can be considered to tune the parameters of the algorithms.

Author Contributions

Conceptualization, M.S. (Mohammad Shokouhifar), M.S. (Mahnaz Sohrabi), M.R. and F.W.; methodology, M.S. (Mohammad Shokouhifar), M.S. (Mahnaz Sohrabi) and M.R.; software, M.S. (Mohammad Shokouhifar); validation, M.S. (Mohammad Shokouhifar), M.S. (Mahnaz Sohrabi) and F.W.; formal analysis, M.S. (Mahnaz Sohrabi), S.M.H.M. and F.W.; investigation, F.W.; resources, M.R.; data curation, M.R. and S.M.H.M.; writing—original draft preparation, M.S. (Mohammad Shokouhifar), M.S. (Mahnaz Sohrabi), M.R. and S.M.H.M.; writing—review and editing, M.S. (Mohammad Shokouhifar), M.S. (Mahnaz Sohrabi) and F.W.; visualization, M.S. (Mohammad Shokouhifar) and M.S. (Mahnaz Sohrabi); supervision, S.M.H.M. and F.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used in the study are available from the authors and can be shared upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Huang, H.; Li, J.; Li, B.; Zhang, D.; Zhao, N.; Tang, S. Comparison of different K-struvite crystallization processes for simultaneous potassium and phosphate recovery from source-separated urine. Sci. Total Environ. 2019, 651, 787–795. [Google Scholar] [CrossRef]
  2. Nedelciu, C.E.; Ragnarsdottir, K.V.; Schlyter, P.; Stjernquist, I. Global phosphorus supply chain dynamics: Assessing regional impact to 2050. Glob. Food Secur. 2020, 26, 100426. [Google Scholar] [CrossRef]
  3. Barbosa, M.W. Uncovering research streams on agri-food supply chain management: A bibliometric study. Glob. Food Secur. 2021, 28, 100517. [Google Scholar] [CrossRef]
  4. Kang, S.; Kim, G.; Roh, J.; Jeon, E.C. Ammonia Emissions from NPK Fertilizer Production Plants: Emission Characteristics and Emission Factor Estimation. Int. J. Environ. Res. Public Health 2022, 19, 6703. [Google Scholar] [CrossRef]
  5. Shokouhifar, M. FH-ACO: Fuzzy heuristic-based ant colony optimization for joint virtual network function placement and routing. Appl. Soft Comput. 2021, 107, 107401. [Google Scholar] [CrossRef]
  6. Al-Qaness, M.A.; Ewees, A.A.; Fan, H.; Abd Elaziz, M. Optimized forecasting method for weekly influenza confirmed cases. Int. J. Environ. Res. Public Health 2020, 17, 3510. [Google Scholar] [CrossRef]
  7. Khairunissa, M.; Lee, H. Hybrid Metaheuristic-Based Spatial Modeling and Analysis of Logistics Distribution Center. ISPRS Int. J. Geo-Inf. 2022, 11, 5. [Google Scholar] [CrossRef]
  8. Shokouhifar, M.; Ranjbarimesan, M. Multivariate time-series blood donation/demand forecasting for resilient supply chain management during COVID-19 pandemic. Clean. Logist. Supply Chain 2022, 5, 100078. [Google Scholar] [CrossRef]
  9. Shehadeh, H.A.; Mustafa, H.M.J.; Tubishat, M. A Hybrid Genetic Algorithm and Sperm Swarm Optimization (HGASSO) for Multimodal Functions. Int. J. Appl. Metaheuristic Comput. (IJAMC) 2022, 13, 1–33. [Google Scholar] [CrossRef]
  10. Ilinova, A.; Dmitrieva, D.; Kraslawski, A. Influence of COVID-19 pandemic on fertilizer companies: The role of competitive advantages. Resour. Policy 2021, 71, 102019. [Google Scholar] [CrossRef]
  11. Rabbani, M.; Molana, S.M.H.; Sajadi, S.M.; Davoodi, M.H. Sustainable fertilizer supply chain network design using evolutionary-based resilient robust stochastic programming. Comput. Ind. Eng. 2022, 174, 108770. [Google Scholar] [CrossRef]
  12. Scholz, R.W.; Ulrich, A.E.; Eilittä, M.; Roy, A. Sustainable use of phosphorus: A finite resource. Sci. Total Environ. 2013, 461, 799–803. [Google Scholar] [CrossRef]
  13. Gong, H.; Meng, F.; Wang, G.; Hartmann, T.E.; Feng, G.; Wu, J.; Jiao, X.; Zhang, F. Toward the sustainable use of mineral phosphorus fertilizers for crop production in China: From primary resource demand to final agricultural use. Sci. Total Environ. 2022, 804, 150183. [Google Scholar] [CrossRef]
  14. Godde, C.M.; Mason-D’Croz, D.; Mayberry, D.E.; Thornton, P.K.; Herrero, M. Impacts of climate change on the livestock food supply chain; a review of the evidence. Glob. Food Secur. 2021, 28, 100488. [Google Scholar] [CrossRef]
  15. Yu, X.; Keitel, C.; Dijkstra, F.A. Global analysis of phosphorus fertilizer use efficiency in cereal crops. Glob. Food Secur. 2021, 29, 100545. [Google Scholar] [CrossRef]
  16. Withers, P.J.A.; Sylvester-Bradley, R.; Jones, D.L.; Healey, J.R.; Talboys, P.J. Feed the Crop not the Soil: Rethinking Phosphorus Management in the Food Chain; ACS Publications: Washington, DC, USA, 2014. [Google Scholar]
  17. Gaffney, J.; Bing, J.; Sawyer, J.; Byrne, P.; Cassman, K.; Ciampitti, I.; Warner, D. Science-based intensive agriculture: Sustainability, food security, and the role of technology. Glob. Food Secur. 2019, 23, 236–244. [Google Scholar] [CrossRef]
  18. Chowdhury, R.B.; Moore, G.A.; Weatherley, A.J.; Arora, M. A novel substance flow analysis model for analysing multi-year phosphorus flow at the regional scale. Sci. Total Environ. 2016, 572, 1269–1280. [Google Scholar] [CrossRef]
  19. Gurevitch, J.; Koricheva, J.; Nakagawa, S.; Stewart, G. Meta-analysis and the science of research synthesis. Nature 2018, 555, 175–182. [Google Scholar] [CrossRef]
  20. De Oliveira, A.K.S.; Soares, E.B.; dos Santos, M.G.; Lins, H.A.; de Freitas Souza, M.; dos Santos Coêlho, E.; Silveira, L.M.; Mendonça, V.; Barros Júnior, A.P.A.; de Araújo de Araújo Rangel Lopes, W. Efficiency of Phosphorus Use in Sunflower. Agronomy 2022, 12, 1558. [Google Scholar] [CrossRef]
  21. Tabbassum, R.; Naveed, M.; Mehboob, I.; Babar, M.H.; Holatko, J.; Akhtar, N.; Rafique, M.; Kucerik, J.; Brtnicky, M.; Kintl, A. Comparative Response of Fermented and Non-Fermented Animal Manure Combined with Split Dose of Phosphate Fertilizer Enhances Agronomic Performance and Wheat Productivity through Enhanced P Use Efficiency. Agronomy 2022, 12, 2335. [Google Scholar] [CrossRef]
  22. Fakhravar, H. Combining heuristics and Exact Algorithms: A Review. arXiv 2022, arXiv:2202.02799. [Google Scholar]
  23. Tirkolaee, E.B.; Goli, A.; Ghasemi, P.; Goodarzian, F. Designing a sustainable closed-loop supply chain network of face masks during the COVID-19 pandemic: Pareto-based algorithms. J. Clean. Prod. 2022, 333, 130056. [Google Scholar] [CrossRef]
  24. Sohrabi, M.; Zandieh, M.; Afshar-nadjafi, B. A simple empirical inventory model for managing with demand differentiation. Comput. Appl. Math. 2021, 40, 1–38. [Google Scholar] [CrossRef]
  25. Sohrabi, M.; Zandieh, M.; Nadjafi, B.A. Dynamic demand-centered process-oriented data model for inventory management of hemovigilance systems. Healthc. Inform. Res. 2021, 27, 73–81. [Google Scholar] [CrossRef]
  26. Shokouhifar, M. Swarm intelligence RFID network planning using multi-antenna readers for asset tracking in hospital environments. Comput. Netw. 2021, 198, 108427. [Google Scholar] [CrossRef]
  27. Rajabi Moshtaghi, H.; Toloie Eshlaghy, A.; Motadel, M.R. A comprehensive review on meta-heuristic algorithms and their classification with novel approach. J. Appl. Res. Ind. Eng. 2021, 8, 63–89. [Google Scholar]
  28. Farag, M.A.; El-Shorbagy, M.A.; Mousa, A.A.; El-Desoky, I.M. A new hybrid metaheuristic algorithm for multiobjective optimization problems. Int. J. Comput. Intell. Syst. 2020, 13, 920. [Google Scholar] [CrossRef]
  29. Chouhan, V.K.; Khan, S.H.; Hajiaghaei-Keshteli, M. Metaheuristic approaches to design and address multi-echelon sugarcane closed-loop supply chain network. Soft Comput. 2021, 25, 11377–11404. [Google Scholar] [CrossRef]
  30. Gholizadeh, H.; Goh, M.; Fazlollahtabar, H.; Mamashli, Z. Modelling uncertainty in sustainable-green integrated reverse logistics network using metaheuristics optimization. Comput. Ind. Eng. 2022, 163, 107828. [Google Scholar] [CrossRef]
  31. Karmakar, S.; Kundu, A.; John, B. Optimizing a Supply Chain Network Using Metaheuristic for Pre and Post Pandemic Scenario. In Proceedings of the 2021 IEEE International Conference on Industrial Engineering and Engineering Management (IEEM), Singapore, 13–16 December 2021; pp. 41–45. [Google Scholar]
  32. Nath, M.P.; Priyadarshini, S.B.B.; Mishra, D. Supply chain management (SCM): Employing various big data and metaheuristic strategies. In Advances in Machine Learning for Big Data Analysis; Springer: Berlin/Heidelberg, Germany, 2022; pp. 145–165. [Google Scholar]
  33. Nahavandi, B.; Homayounfar, M.; Daneshvar, A.; Shokouhifar, M. Hierarchical structure modelling in uncertain emergency location-routing problem using combined genetic algorithm and simulated annealing. Int. J. Comput. Appl. Technol. 2022, 68, 150–163. [Google Scholar] [CrossRef]
  34. Wu, W.; Zhou, W.; Lin, Y.; Xie, Y.; Jin, W. A hybrid metaheuristic algorithm for location inventory routing problem with time windows and fuel consumption. Expert Syst. Appl. 2021, 166, 114034. [Google Scholar] [CrossRef]
  35. Shokouhifar, M.; Pilevari, N. Combined adaptive neuro-fuzzy inference system and genetic algorithm for e-learning resilience assessment during COVID-19 pandemic. Concurr. Comput. Pract. Exp. 2022, 34, e6791. [Google Scholar] [CrossRef]
  36. Liu, Y.; Chen, C. Improved RFM Model for Customer Segmentation Using Hybrid Meta-heuristic Algorithm in Medical IoT Applications. Int. J. Artif. Intell. Tools 2022, 31, 2250009. [Google Scholar] [CrossRef]
  37. Ghasemi Darehnaei, Z.; Shokouhifar, M.; Yazdanjouei, H.; Rastegar Fatemi, S.M.J. SI-EDTL: Swarm intelligence ensemble deep transfer learning for multiple vehicle detection in UAV images. Concurr. Comput. Pract. Exp. 2022, 34, e6726. [Google Scholar] [CrossRef]
  38. Shokouhifar, M.; Sabbaghi, M.M.; Pilevari, N. Inventory management in blood supply chain considering fuzzy supply/demand uncertainties and lateral transshipment. Transfus. Apher. Sci. 2021, 60, 103103. [Google Scholar] [CrossRef]
  39. Naderi, R.; Shafiei Nikabadi, M.; Alem Tabriz, A.; Pishvaee, M.S. Supply chain sustainability improvement using exergy analysis. Comput. Ind. Eng. 2021, 154, 107142. [Google Scholar] [CrossRef]
  40. Sohrabi, M.; Zandieh, M.; Shokouhifar, M. Sustainable inventory management in blood banks considering health equity using a combined metaheuristic-based robust fuzzy stochastic programming. Socio. -Econ. Plan. Sci. 2022. In press. [Google Scholar] [CrossRef]
  41. Sörensen, K. Metaheuristics—The metaphor exposed. Int. Trans. Oper. Res. 2015, 22, 3–18. [Google Scholar] [CrossRef]
  42. Esmaeili, H.; Hakami, V.; Bidgoli, B.M.; Shokouhifar, M. Application-specific clustering in wireless sensor networks using combined fuzzy firefly algorithm and random forest. Expert Syst. Appl. 2022, 210, 118365. [Google Scholar] [CrossRef]
  43. Shokouhifar, M.; Jalali, A. Optimized sugeno fuzzy clustering algorithm for wireless sensor networks. Eng. Appl. Artif. Intell. 2017, 60, 16–25. [Google Scholar] [CrossRef]
  44. Mirjalili, S.; Lewis, A. The whale optimization algorithm. Adv. Eng. Softw. 2016, 95, 51–67. [Google Scholar] [CrossRef]
  45. Wang, K.M.; Wang, K.J.; Chen, C.C. Capacitated production planning by parallel genetic algorithm for a multi-echelon and multi-site TFT-LCD panel manufacturing supply chain. Appl. Soft Comput. 2022, 127, 109371. [Google Scholar] [CrossRef]
  46. Fathollahi-Fard, A.M.; Govindan, K.; Hajiaghaei-Keshteli, M.; Ahmadi, A. A green home health care supply chain: New modified simulated annealing algorithms. J. Clean. Prod. 2019, 240, 118200. [Google Scholar] [CrossRef]
Figure 1. The PFSCM model.
Figure 1. The PFSCM model.
Agronomy 13 00565 g001
Figure 2. Overall flowchart of the proposed H-WOA-VNS algorithm.
Figure 2. Overall flowchart of the proposed H-WOA-VNS algorithm.
Agronomy 13 00565 g002
Figure 3. Encoding of the direct DVs into a feasible solution.
Figure 3. Encoding of the direct DVs into a feasible solution.
Agronomy 13 00565 g003
Figure 4. Integer local search operators: (a) I-Swap, and (b) I-Exchange. The variables selected for change in SOLcurrent and the changed variables in SOLnew are shown in blue and green, respectively.
Figure 4. Integer local search operators: (a) I-Swap, and (b) I-Exchange. The variables selected for change in SOLcurrent and the changed variables in SOLnew are shown in blue and green, respectively.
Agronomy 13 00565 g004
Figure 5. Permutation operators: (a) Exchange, (b) Relocate, (c) OrOpt, (d) TwoOpt, and (e) Reverse. The selected variables for change are shown in colors.
Figure 5. Permutation operators: (a) Exchange, (b) Relocate, (c) OrOpt, (d) TwoOpt, and (e) Reverse. The selected variables for change are shown in colors.
Agronomy 13 00565 g005
Figure 6. Best OF versus iteration in WOA and VNS phases: WOA phase (left), VNS phase (right).
Figure 6. Best OF versus iteration in WOA and VNS phases: WOA phase (left), VNS phase (right).
Agronomy 13 00565 g006
Figure 7. Convergence of the H-WOA-VNS algorithm in terms of the best OF versus NFE.
Figure 7. Convergence of the H-WOA-VNS algorithm in terms of the best OF versus NFE.
Agronomy 13 00565 g007
Figure 8. The best OF versus NFE obtained by WOA, VNS, and H-WOA-VNS.
Figure 8. The best OF versus NFE obtained by WOA, VNS, and H-WOA-VNS.
Agronomy 13 00565 g008
Figure 9. Sensitivity analysis of the environmental weights on the normalized economic cost: the lower the value for the weights of environmental objectives (either wCY or wPUE) is, the higher the value for the economic objective (wEC) is considered, and thus, a lower economic cost (ZEC) is obtained.
Figure 9. Sensitivity analysis of the environmental weights on the normalized economic cost: the lower the value for the weights of environmental objectives (either wCY or wPUE) is, the higher the value for the economic objective (wEC) is considered, and thus, a lower economic cost (ZEC) is obtained.
Agronomy 13 00565 g009
Figure 10. Comparison of the worst, mean, and best results, obtained by different methods.
Figure 10. Comparison of the worst, mean, and best results, obtained by different methods.
Agronomy 13 00565 g010
Table 1. List of indices, sets, parameters, and decision variables.
Table 1. List of indices, sets, parameters, and decision variables.
Sets and IndicesDefinition
iISet of P-mines (suppliers of phosphorus)
jJSet of suppliers of raw materials
mMSet of manufacturers
dDSet of distributors
fFSet of farms (demand nodes)
tTSet of months (time periods)
rRSet of raw materials: A (r = 1), PA (r = 2), and SA (r = 3)
pPSet of products: SSP (p = 1), TSP (p = 2), MAP (p = 3), and DAP (p = 4)
ParametersDefinition
B f d 1 if farm f is supported by distribution center d; 0 otherwise
B S p Batch size of ordering fertilizer p by farms (ton)
C a p i t P Capacity of supplying P by mine i in time t (ton)
C a p j r t S Capacity of supplier j to supply raw material r in time t (ton)
C a p m p M Capacity of manufacturer m to produce fertilizer p in each time (ton)
W m M Warehouse capacity of manufacturer m (ton)
W d D Warehouse capacity of distribution center d (ton)
W R Warehouse capacity of recycling center (ton)
U p P Required P rock to produce unit fertilizer p (%)
U r p S Required raw material r for the production of unit fertilizer p (%)
a f p Phosphorus uptake by farm f from unit fertilizer p (%)
β f p Crop yield increase in farm f from unit fertilizer p (%)
γ f p Recyclable phosphorus in farm f from unit fertilizer p (%)
λ p Total phosphorus loss to produce a unit of fertilizer p (%)
μ p Amount of recyclable P per unit P-leaching of fertilizer p (%)
D P f t Phosphorus demand by farm f in time period t (ton)
L P f p t Minimum amount of fertilizer p which should be delivered to farm f in time t (ton)
U P f p t Maximum amount of fertilizer p which should be delivered to farm f in time t (ton)
T C i m P M Transportation cost of carrying phosphorus from mine i to manufacturer m (USD/ton)
T C m R M Transportation cost of carrying phosphorus from recycling center to manufacturer m (USD/ton)
T C j m S M Transportation cost from supplier j to manufacturer m for carrying raw materials (USD/ton)
T C m d M D Transportation cost of carrying fertilizers from manufacturer m to distribution center d (USD/ton)
T C d f D F Transportation cost of carrying fertilizers from distribution center d to farm f (USD/ton)
T C f F R Transportation cost of carrying phosphorus leaching from farm f to recycling center (USD/ton)
B C i P Purchasing cost of phosphate from mine i (USD/ton)
B C j r S Purchasing cost from supplier j for raw material r (USD/ton)
P C m p Production cost in manufacturer m for unit fertilizer p (USD/ton)
R C Recycling cost of unit P in recycling center (USD/ton)
I m p 0 Initial inventory of fertilizer p in manufacturer m (ton)
I C m M Inventory cost in manufacturer m for each month (USD/ton/month)
I C d D Inventory cost in distribution center d for each month (USD/ton/month)
I C R Inventory cost in recycling center for each month (USD/ton/month)
Direct DVsDefinition
P m p t Amount of fertilizer p produced in time t by manufacturer j (ton)
R f p t Amount of fertilizer p delivered in time t to farm f (BSp)
P L M Priority list of manufacturers determining their order in P-mine and supplier assignment
S L P m Selection list of P-mines (plus recycling center) for manufacturer m
S L S m Selection list of suppliers for manufacturer m
P L D Priority list of distributors determining the order of distributors in manufacturer assignment
S L M d Selection list of manufacturers for distributor d
Indirect DVsDefinition
S P m r p t Required amount of raw material r in manufacturer m to produce fertilizer p in time t (ton)
S m r t Required raw material r in manufacturer m to produce all fertilizers in time t (ton)
R P t Total P recycled at time t (ton)
T D d p t Total demand of fertilizer p in distribution center d in time t (ton)
X i m t P 1 if mine i supplies P for manufacturer m in time t; 0 otherwise
X m t R 1 if the recycling center supplies P to manufacturer m in time t; 0 otherwise
X j m r t S 1 if supplier j supplies raw material r to manufacturer m in time t; 0 otherwise
X m d p t M 1 if fertilizer p is delivered from manufacturer m to distribution center d in time t; 0 otherwise
X d f p t D 1 if fertilizer p is delivered from distribution center d to farm f in time t; 0 otherwise
X f p t F 1 if P-leaching of fertilizer p is delivered from farm f to the recycling center in time t; 0 otherwise
Y i m t P Amount of P delivered in time t from mine i to manufacturer m (ton)
Y m t R Amount of P transferred in time t from the recycling center to manufacturer m (ton)
Y j m r t S Amount of raw material r transferred in time t from supplier j to manufacturer m (ton)
Y m d p t M Amount of fertilizer p transferred in time t from manufacturer m to distribution center d (ton)
Y d f p t D Amount of fertilizer p transferred in time t from distribution center d to farm f (ton)
Y f p t F Amount of P-leaching of fertilizer p transferred in time t from farm f to the recycling center (ton)
I P m p t Inventory of manufacturer m for fertilizer p at time t (ton)
I D d p t Inventory of distribution center d for fertilizer p at time t (ton)
I R t Inventory of recycled P in the recycling center at time t (ton)
Table 2. Demands of the distribution centers in terms of the total P-uptake and lower/upper bound on the required fertilizers (ton/year).
Table 2. Demands of the distribution centers in terms of the total P-uptake and lower/upper bound on the required fertilizers (ton/year).
DistributorP-UptakeTSP MAP DAP
LowerUpperLowerUpperLowerUpper
125,820912035,82010,66041,580557025,380
223,130998033,66012,86041,400720023,040
322,18010,27035,280720052,020576019,080
417,210691030,780778031,860384013,680
58610326014,220336013,50020208460
610,260346015,300547022,14025908640
7281013404860144046806702700
87260154012,600422013,32016306480
95140192079202780864014404320
105810259084601820990015405940
1141,33014,50070,02017,86066,24011,62042,300
1210,760317020,160490019,980259011,160
1352,81020,74097,20026,40092,34013,63047,520
1412,180653024,120480020,880288010,980
154350173054002300756012503420
169810470016,020566019,800202010,440
1737,27014,59068,40014,78059,940701032,760
1815,000566022,500557028,080374011,700
19418019206120202081005803780
2019,200768032,760883033,120432019,080
2115,410643025,380758028,440240012,960
2224,67014,30046,98012,77046,980547026,100
234530202073801340792012504860
2425,05011,14040,680931052,020547022,680
2516,460691024,840768031,860336015,300
2620,360912034,560662030,780442017,820
2711,090451019,980499021,060240010,800
289610403013,500451017,10024008640
296830355010,440394010,98021107740
3019,480710033,300893036,000422018,540
31356013405580125059408603600
3213,840643021,240710021,600259013,860
Sum506,010208,490845,460226,730905,760118,850473,760
Table 3. Raw materials composition (%) for the different P-fertilizers.
Table 3. Raw materials composition (%) for the different P-fertilizers.
P-FertilizerPASAAP
SSP064064
TSP3440040
MAP510120
DAP470230
Table 4. Capacity of the suppliers (ton/year).
Table 4. Capacity of the suppliers (ton/year).
Raw MaterialMinimum Capacity Maximum CapacityTotal Capacity
PA12,000240,000870,000
SA15,000620,0003,420,000
A18,0001,550,0006,200,000
P80,000250,000830,000
Table 5. Monthly capacity of the manufacturers for the P-fertilizer production (ton).
Table 5. Monthly capacity of the manufacturers for the P-fertilizer production (ton).
ManufacturerSSPTSPMAPDAP
115,0005000910011,860
2666013,34011,00014,500
32500250025603400
43340334000
52500166024803180
62500250000
73340166025203460
84160000
983405000816011,320
10016,66013,38023,860
1116,660021,32029,360
12016,66000
1310,000208056409080
142400450000
15834016,66000
16416084032203840
174160000
187500250000
198340834000
206000600063608860
213340150000
Sum119,240110,74085,740122,720
Table 6. Parameters of the H-WOA-VNS algorithm.
Table 6. Parameters of the H-WOA-VNS algorithm.
PhaseParameterValue/Description
HeuristicMaximum Iterations (MaxIterWOA)70 (=PopSizeWOA)
WOAPopulation Size (PopSizeWOA)500
Population updating mechanisms70
Maximum Iterations (MaxIterVNS)Encircling prey, Search for prey, Bubble-net attack
VNSNumber of Local Search operators (NumLSVNS)5000
Local search operators7
Weight of economic cost (wEC)I-Swap, I-Exchange, Exchange, Relocate, OrOpt, TwoOpt, and Reverse
OF weightsWeight of crop yield (wCY)0.5
Maximum Iterations (MaxIterWOA)0.3
Weight of PUE (wPUE)0.2
Table 7. Setting of the PFSCM model parameters.
Table 7. Setting of the PFSCM model parameters.
(a) Purchasing Cost of the Raw Materials (USD/ton)
Raw MaterialPASAAP
Purchasing cost 310~37060~7020~2580~95
(b) Production cost of P-fertilizers (USD/ton)
FertilizerSSPTSPMAPDAP
Production cost 110~130120~150210~250230~280
(c) P-fertilizers coefficients (%)
ParameterSSPTSPMAPDAP
a fp 0.26~0.320.29~0.350.42~0.470.46~0.52
β fp 0.21~0.280.24~0.30.38~0.460.45~0.51
γ fp 0.23~0.270.21~0.270.12~0.170.16~0.2
λ p 0.530.380.640.62
μ p 0.090.10.070.06
(d) Other parameters
ParameterValueParameterValue
BS p 1 (ton) TC im PM 1.5~2 (USD/truck/km)
Truck size25 (ton) TC m RM 1.5~2 (USD/truck /km)
RC 20 (USD/ton) TC jm SM 2~3 (USD/truck /km)
IC m M 1.5 (USD/ton/month) TC md MD 2~2.5 (USD/truck /km)
IC d D 2.5 (USD/ton/month) TC df DF 2~2.5 (USD/truck /km)
IC R 2 (USD/ton/month) TC f FR 2.5~3 (USD/truck /km)
Table 8. Results of the H-WOA-VNS algorithm in 10 successive runs to solve the sustainable PFSCM model for the Case Study.
Table 8. Results of the H-WOA-VNS algorithm in 10 successive runs to solve the sustainable PFSCM model for the Case Study.
RunTotal Cost (M$) Z E C ¯ Z C Y ¯ Z P U E ¯ PFOF
1451.810.30120.32910.278900.4961
2454.390.30290.330.27800.4969
3451.870.30120.33130.277800.4957
4446.380.29760.33020.278100.4941
5449.970.30.33170.277600.495
6450.40.30030.33030.27800.4954
7448.040.29870.33030.277700.4947
8453.470.30230.32980.277900.4966
9449.10.29940.32950.278600.4951
10446.550.29770.32930.278500.4944
Worst454.390.30290.32910.277600.4969
Best446.380.29760.33170.278900.4941
Mean450.20.30010.33010.278100.4954
STD%0.610.610.250.1500.19
Table 9. Consumed raw materials for the production of the different fertilizers (ton).
Table 9. Consumed raw materials for the production of the different fertilizers (ton).
FunctionCPPCPSCPRCRCCIHCTRTotal (ZEC)
Cost (M$)12.45118.68178.738.4325.23102.86446.38
Table 10. Total production of the P-fertilizers by each manufacturer (ton).
Table 10. Total production of the P-fertilizers by each manufacturer (ton).
ManufacturerSSPTSPMAPDAP
162,20116,81423,75513,767
220,88034,68424,41129,242
334436301761010,408
410,60313,51300
59711748579303633
66847691500
712,862402042877042
811,579000
936,14219,01628,0274825
10040,03027,37429,133
1155,695065,21852,210
12029,81800
1324,000736520,8056362
14452117,70400
1510,99954,91400
1610,2592594516110,004
1714,697000
1816,603830100
1912,71135,75100
2012,26424,22015,39819,173
219984468800
Total346,001334,133229,976185,799
Table 11. Consumed raw materials (ton) for the production of the different fertilizers.
Table 11. Consumed raw materials (ton) for the production of the different fertilizers.
P-FertilizerPASAAP
SSP0128,0200221,441
TSP113,60500133,653
MAP117,288027,5970
DAP87,326042,7340
Total318,219128,02070,331355,090
Table 12. Total delivered P-fertilizers (ton) and P-uptake (ton) in each distribution center.
Table 12. Total delivered P-fertilizers (ton) and P-uptake (ton) in each distribution center.
DistributorDelivered P-FertilizersP-Uptake
SSPTSPMAPDAPDemandSatisfied
122,37118,38413,86912,44225,82025,990
218,54019,41111,70711,30923,13023,536
317,26518,40312,01610,49622,18022,531
413,53612,66210,395823517,21017,542
5675870884386444386108806
6779983655305557810,26010,572
7247921931453150328102957
8553662863956356372607503
9400642922750243451405216
10499944023140271158105888
1133,18937,55419,86319,20041,33042,001
12833581056165513310,76010,825
1340,30438,31331,62527,52652,81054,315
14879190926695646512,18012,234
15335434892418223643504490
16770477574956502798109891
1731,44728,59821,20917,34437,27038,173
1813,16111,2128016717415,00015,278
19321633752488196341804301
2015,69514,54010,099954319,20019,373
2113,02211,9817935730315,41015,518
2221,68318,11812,70911,91924,67024,872
23374934602403227245304615
2419,97720,23512,62912,22925,05025,165
2513,56012,52610,850777516,46017,480
2615,49116,72311,036963020,36020,515
27876580256110562811,09011,164
28821165725299486196109732
29598658143638342468307265
3016,46314,592949910,40719,48019,795
31318326411846171235603616
3210,75710,0237852708513,84014,017
Total409,332394,231274,317248,570506,010515,176
Table 13. Distribution of the different P-fertilizers (ton) in the system.
Table 13. Distribution of the different P-fertilizers (ton) in the system.
SSPTSPMAPDAP
Initial inventories70,00070,00050,00070,000
Total productions346,001334,133229,976185,799
Total P-fertilizers416,001404,133279,976255,799
Total satisfied demands409,332394,231274,317248,570
Final inventories6,6699,9025,6597,229
Table 14. PFSCM test examples.
Table 14. PFSCM test examples.
DatasetRPI + 1JMDFT
SD 121122223
SD 221133343
SD 3222535106
SD 43327510106
SD 5443105103012
Case Study44552213256512
Table 15. Comparison of the average results of H-WOA-VNS over 10 runs against exact search.
Table 15. Comparison of the average results of H-WOA-VNS over 10 runs against exact search.
DatasetOptimal (Exact Search)H-WOA-VNSError (%)
OFTime (s)OFTime (s)
SD 10.47120.10.4712320.000
SD 20.3973.50.397480.000
SD 30.43512760.43521260.023
SD 40.461337,6500.46192350.13
SD 5N/AN/A0.4825846N/A
Case StudyN/AN/A0.49542725N/A
Table 16. Comparison of the H-WOA-VNS algorithm against other techniques, over 10 runs.
Table 16. Comparison of the H-WOA-VNS algorithm against other techniques, over 10 runs.
RunHeuristic
(Stage 1)
WOA
(Stage 2)
VNS
(Stage 3)
GA
[45]
SA
[46]
GLGASA
[39]
H-WOA-VNS (Proposed)
10.54950.50340.52830.52250.51270.51630.4961
20.54630.52670.50360.53440.55220.52220.4969
30.55710.5020.53280.51070.52410.52130.4957
40.54120.52040.50780.50630.52450.50970.4941
50.55340.51560.52750.52630.53220.5310.495
60.56510.49830.52240.51850.51670.52210.4954
70.53860.50660.50450.51220.5420.52360.4947
80.54380.5130.5210.53320.53330.50210.4966
90.54440.50780.50190.53270.51910.51110.4951
100.55120.50830.51170.52780.53730.5250.4944
Worst0.56510.52670.53280.53440.55220.5310.4969
Best0.53860.49830.50190.50630.51270.50210.4941
Mean0.54910.51020.51610.52250.52940.51840.4954
STD%1.451.722.241.942.331.660.19
Table 17. Effect of the weights of the total objective function on the different sub-objectives.
Table 17. Effect of the weights of the total objective function on the different sub-objectives.
w E C w C Y w P U E Z E C ¯ Z C Y ¯ Z P U E ¯
0.50.30.20.30010.33010.2781
0.30.50.20.4720.56440.2846
0.20.30.50.43460.40320.4328
0.50.500.41670.4750.2153
0.500.50.27650.23870.585
00.50.50.56230.6710.5012
Table 18. Improvement rate (%) of the H-WOA-VNS algorithm against other algorithms, in terms of the worst, best, mean, and STD% of the obtained OF over 10 runs.
Table 18. Improvement rate (%) of the H-WOA-VNS algorithm against other algorithms, in terms of the worst, best, mean, and STD% of the obtained OF over 10 runs.
MeasureGASAGLGASAHeuristicWOAVNS
Worst7.0210.016.4212.075.666.74
Best2.413.631.598.260.841.55
Mean5.196.424.449.782.94.01
STD%90.2191.8588.5586.988.9591.52
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Shokouhifar, M.; Sohrabi, M.; Rabbani, M.; Molana, S.M.H.; Werner, F. Sustainable Phosphorus Fertilizer Supply Chain Management to Improve Crop Yield and P Use Efficiency Using an Ensemble Heuristic–Metaheuristic Optimization Algorithm. Agronomy 2023, 13, 565. https://doi.org/10.3390/agronomy13020565

AMA Style

Shokouhifar M, Sohrabi M, Rabbani M, Molana SMH, Werner F. Sustainable Phosphorus Fertilizer Supply Chain Management to Improve Crop Yield and P Use Efficiency Using an Ensemble Heuristic–Metaheuristic Optimization Algorithm. Agronomy. 2023; 13(2):565. https://doi.org/10.3390/agronomy13020565

Chicago/Turabian Style

Shokouhifar, Mohammad, Mahnaz Sohrabi, Motahareh Rabbani, Seyyed Mohammad Hadji Molana, and Frank Werner. 2023. "Sustainable Phosphorus Fertilizer Supply Chain Management to Improve Crop Yield and P Use Efficiency Using an Ensemble Heuristic–Metaheuristic Optimization Algorithm" Agronomy 13, no. 2: 565. https://doi.org/10.3390/agronomy13020565

APA Style

Shokouhifar, M., Sohrabi, M., Rabbani, M., Molana, S. M. H., & Werner, F. (2023). Sustainable Phosphorus Fertilizer Supply Chain Management to Improve Crop Yield and P Use Efficiency Using an Ensemble Heuristic–Metaheuristic Optimization Algorithm. Agronomy, 13(2), 565. https://doi.org/10.3390/agronomy13020565

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop