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

Author(s): Mohammad Shokouhifar (corresponding author) [1,*]; Mahnaz Sohrabi [2]; Motahareh Rabbani [3]; Seyyed Mohammad Hadji Molana [3]; Frank Werner (corresponding author) [4,*]

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 (i?I), suppliers of raw materials (j?J), fertilizer manufacturers (m?M), distribution centers (d?D), farms (f?F), and a recycling center, as shown in Figure 1. The model is formulated in T time periods (t?T), i.e., months. The raw materials (r?R) include sulfuric acid (SA), phosphoric acid (PA), and ammonium (A), where PA is made from P and SA. Moreover, four types of products (p?P) 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[sub.1][sup.P]= 0.64 units of P and U[sub.31][sup.S]= 0.37 units of SA are used.

* To produce a unit of TSP, U[sub.2][sup.P]= 0.4 units of P and U[sub.22][sup.S]= 0.34 units of PA are used.

* To produce a unit of MAP, U[sub.13][sup.S]= 0.12 units of A and U[sub.23][sup.S]= 0.51 units of PA are used.

* To produce a DAP, U[sub.14][sup.S]= 0.23 units of A and U[sub.24][sup.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:(1)C[sub.PP]=?i(BC[sub.i][sup.P]?t?mY[sub.imt][sup.P]X[sub.imt][sup.P]) (2)C[sub.PS]=?j?r(BC[sub.jr][sup.S]?t?mY[sub.jmrt][sup.S]X[sub.jmrt][sup.S]) (3)C[sub.PR]=?m?pPC[sub.mp]?tP[sub.mpt] (4)C[sub.RC]=?t(RC?f?pµ[sub.p]Y[sub.fpt][sup.F]X[sub.fpt][sup.F]) (5)C[sub.IH]=?t?m?p(IC[sub.m][sup.M]IP[sub.mpt])+?t?d?p(IC[sub.d][sup.D]ID[sub.dpt])+?t(IR[sup.R]IR[sub.t]) (6)C[sub.TR]=?t?i?mTC[sub.im][sup.PM]Y[sub.imt][sup.P]X[sub.imt][sup.P]+?t?mTC[sub.m][sup.RM]Y[sub.mt][sup.R]X[sub.mt][sup.R]+?t?j?m?rTC[sub.jm][sup.SM]Y[sub.jmrt][sup.S]X[sub.jmrt][sup.S]+?t?m?d?pTC[sub.md][sup.MD]Y[sub.mdpt][sup.M]X[sub.mdpt][sup.M]+?t?d?f?pTC[sub.df][sup.DF]Y[sub.dfpt][sup.D]X[sub.dfpt][sup.D]+?t?f?pTC[sub.f][sup.FR]Y[sub.fpt][sup.F]X[sub.fpt][sup.F]

By aggregating the six economic costs as mentioned above, the total economic objective function can be expressed as:(7)minimizeZ[sub.EC]=C[sub.PP]+C[sub.PS]+C[sub.PR]+C[sub.RC]+C[sub.IH]+C[sub.TR]=(?i(BC[sub.i][sup.P]?t?mY[sub.imt][sup.P]X[sub.imt][sup.P]))+(?j?r(BC[sub.jr][sup.S]?t?mY[sub.jmrt][sup.S]X[sub.jmrt][sup.S]))+(?m?pPC[sub.mp]?tP[sub.mpt])+(?t(RC?f?pµ[sub.p]Y[sub.fpt][sup.F]X[sub.fpt][sup.F]))+(?t?m?p(IC[sub.m][sup.M]IP[sub.mpt])+?t?d?p(IC[sub.d][sup.D]ID[sub.dpt])+?t(IR[sup.R]IR[sub.t]))+(?t?i?mTC[sub.im][sup.PM]Y[sub.imt][sup.P]X[sub.imt][sup.P]+?t?mTC[sub.m][sup.RM]Y[sub.mt][sup.R]X[sub.mt][sup.R]+?t?j?m?rTC[sub.jm][sup.SM]Y[sub.jmrt][sup.S]X[sub.jmrt][sup.S]+?t?m?d?pTC[sub.md][sup.MD]Y[sub.mdpt][sup.M]X[sub.mdpt][sup.M]+?t?d?f?pTC[sub.df][sup.DF]Y[sub.dfpt][sup.D]X[sub.dfpt][sup.D]+?t?f?pTC[sub.f][sup.FR]Y[sub.fpt][sup.F]X[sub.fpt][sup.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:(8)maximizeZ[sub.CY]=?f?pß[sub.fp](?t?dY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]) (9)maximizeZ[sub.PUE]=?p(1-?[sub.p])(?f[a[sub.fp]?d?tY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]+?[sub.fp]?d?tY[sub.dfpt][sup.D]X[sub.dfpt][sup.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[sub.EC]¯, Z[sub.CY]¯, and Z[sub.PUE]¯ 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:(10)minimizeOF=[w[sub.EC]Z[sub.EC]¯+w[sub.CY](1-Z[sub.CY]¯)+w[sub.PUE](1-Z[sub.PUE]¯)]×PF 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:(11)?p(a[sub.fp]?dY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd])=DP[sub.ft]?f,?t (12)LP[sub.fpt]=?dY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]=UP[sub.fpt]?f,?p,?t (13)?dY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]=R[sub.fpt]BS[sub.p]?f,?p,?t (14)X[sub.dfpt][sup.D]=B[sub.fd]?d,?f,?p,?t (15)?dX[sub.dfpt][sup.D]=1?f,?p,?t (16)ID[sub.dp0]=0?d,?p (17)ID[sub.dpt-1]+?mY[sub.mdpt][sup.M]X[sub.mdpt][sup.M]=?fY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]?d,?p,?t (18)ID[sub.dpt]=ID[sub.dpt-1]+?mY[sub.mdpt][sup.M]X[sub.mdpt][sup.M]-?fY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]?d,?p,?t (19)?pID[sub.dpt]=W[sub.d][sup.D]?d,?t (20)P[sub.mpt]=Cap[sub.mp][sup.M]?m,?p,?t (21)?iY[sub.imt][sup.P]X[sub.imt][sup.P]+Y[sub.mt][sup.R]X[sub.mt][sup.R]=?pU[sub.p][sup.P]P[sub.mpt]?m,?t (22)?jY[sub.jmrt][sup.S]X[sub.jmrt][sup.S]=?pU[sub.rp][sup.S]P[sub.mpt]?m,?r,?t (23)IP[sub.mp0]=I[sub.mp][sup.0]?m,?p (24)IP[sub.mpt-1]+P[sub.mpt]=?dY[sub.mdpt][sup.M]X[sub.mdpt][sup.M]?m,?p,?t (25)IP[sub.mpt]=IP[sub.mpt-1]+P[sub.mpt]-?dY[sub.mdpt][sup.M]X[sub.mdpt][sup.M]?m,?p,?t (26)?pIP[sub.mpt]=W[sub.m][sup.M]?m,?t (27)?mY[sub.jmrt][sup.S]X[sub.jmrt][sup.S]=Cap[sub.jrt][sup.S]?j,?r,?t (28)?mY[sub.imt][sup.P]X[sub.imt][sup.P]=Cap[sub.it][sup.P]?i,?t (29)?f?pY[sub.fpt][sup.F]X[sub.fpt][sup.F]=?d?f?p?fp/µpY[sub.dfpz][sup.D]X[sub.dfpz][sup.D]B[sub.fd]?t (30)?mY[sub.mt][sup.R]X[sub.mt][sup.R]=?z=1t-1?f?pµ[sub.p]Y[sub.fpt][sup.F]X[sub.fpt][sup.F]-?z=1t-1?mY[sub.mz][sup.R]X[sub.mz][sup.R]?t (31)IR[sub.0]=0 (32)?mY[sub.mt][sup.R]X[sub.mt][sup.R]=IR[sub.t-1]?t (33)IR[sub.t]=IR[sub.t-1]-?mY[sub.mt][sup.R]X[sub.mt][sup.R]+?f?pµ[sub.p]Y[sub.fpt][sup.F]X[sub.fpt][sup.F]?t (34)IR[sub.t]=W[sup.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.

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[sub.mpt]) 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[sub.fpt]) 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 (PLM) is a permutation of M manufacturers to determine the priority of the manufacturers to order from different P-mines and suppliers.

* SOL.SLP (SLP[sub.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 (SLS[sub.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 (PLD) is a permutation of D distribution centers which determines priority of them in ordering fertilizers from the different manufacturers.

* SOL.SLM (SLM[sub.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 SP[sub.mrpt], S[sub.mrt], RP[sub.t],TD[sub.dpt],X[sub.imt][sup.P], X[sub.mt][sup.R], X[sub.jmrt][sup.S], X[sub.mdpt][sup.M], X[sub.dfpt][sup.D], X[sub.fpt][sup.F], Y[sub.imt][sup.P], Y[sub.mt][sup.R], Y[sub.jmrt][sup.S], Y[sub.mdpt][sup.M], Y[sub.dfpt][sup.D], Y[sub.fpt][sup.F], IP[sub.mpt], ID[sub.dpt], and IR[sub.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[sub.dfpt][sup.D] and Y[sub.dfpt][sup.D] are obtained according to the requested demands stored in R[sub.fpt] as:(35)X[sub.dfpt][sup.D]={1ifB[sub.fd]=1andR[sub.fpt]>00else?p,?d,?f,?t (36)Y[sub.dfpt][sup.D]={R[sub.fpt]BS[sub.p]ifB[sub.fd]=10else?p,?d,?f,?t

Decoding of farms-recycling center DVs: at the end of each time period, after delivering fertilizers Y[sub.dfpt][sup.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[sub.fpt][sup.F], Y[sub.fpt][sup.F], and RP[sub.t], are extracted:(37)Y[sub.fpt][sup.F]=?d?fp/µpY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]?p,?f,?t (38)X[sub.fpt][sup.F]={1ifY[sub.fpt][sup.F]>00else?p,?f,?t (39)RP[sub.t]=?f?pµ[sub.p]Y[sub.fpt][sup.F]X[sub.fpt][sup.F]=?f?p?d?[sub.fp]Y[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]?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 PLD. For each distribution center d, different manufacturers according to their selection list in SLM[sub.d] are evaluated one by one to deliver fertilizers to the distributor d until satisfying its total demand TD[sub.dpt]. As a result, the values of TD[sub.dpt], X[sub.mdpt][sup.M], and Y[sub.mdpt][sup.M] are extracted:(40)TD[sub.dpt]=?fY[sub.dfpt][sup.D]X[sub.dfpt][sup.D]B[sub.fd]?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[sub.rp][sup.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 PLM. The manufacturer m may be sourced by different suppliers according to their selection list in SLS[sub.m] until satisfying the total demand of the manufacturer m for raw materials. As a result, S[sub.mrpt][sup.S], ST[sub.mrt][sup.S], X[sub.jmrt][sup.S], and Y[sub.jmrt][sup.S], are obtained:(41)S[sub.mrpt][sup.S]=U[sub.rp][sup.S]P[sub.mpt]?r,?p,?m,?t (42)ST[sub.mrt][sup.S]=?pS[sub.mrpt][sup.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 SLP[sub.m] until satisfying the phosphate demand of the manufacturer m. To this end, S[sub.mpt][sup.P], ST[sub.mt][sup.P], X[sub.imt][sup.P], X[sub.mt][sup.R], Y[sub.imt][sup.P], and Y[sub.mt][sup.R] are achieved:(43)S[sub.mrpt][sup.S]=U[sub.rp][sup.S]P[sub.mpt]?r,?p,?m,?t (44)ST[sub.mrt][sup.S]=?pS[sub.mrpt][sup.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[sub.fpt], 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 × Cap[sub.mp][sup.M],Cap[sub.mp][sup.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?=2a?·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:(45)SOL[sub.w](t+1)=SOL[sub.rand](t)-A?·|2r?·SOL[sub.rand](t)-SOL[sub.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:(46)SOL[sub.w](t+1)=GBestSOL-A?·|2r?·GBestSOL-SOL[sub.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:(47)SOL[sub.w](t+1)=|GBestSOL-SOL[sub.w](t)|·e[sup.bl]cos(2pl)+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., SOL[sup.current] = GBestSOL. In each iteration of VNS, a new solution, namely SOL[sup.new], is constructed in the vicinity of the current solution, SOL[sup.current], 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 MaxIter[sub.VNS] × NumLS[sub.VNS] = MaxIter[sub.WOA] × PopSize[sub.WOA]. 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% (OF[sub.worst] - OF[sub.best] = 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., MaxIter[sub.WOA] = 1000 and MaxIter[sub.VNS] = 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[sub.EC]¯) 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.

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.

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.

References

1. H. Huang; J. Li; B. Li; D. Zhang; N. Zhao; S. Tang Comparison of different K-struvite crystallization processes for simultaneous potassium and phosphate recovery from source-separated urine., 2019, 651,pp. 787-795. DOI: https://doi.org/10.1016/j.scitotenv.2018.09.232.

2. C.E. Nedelciu; K.V. Ragnarsdottir; P. Schlyter; I. Stjernquist Global phosphorus supply chain dynamics: Assessing regional impact to 2050., 2020, 26,p. 100426. DOI: https://doi.org/10.1016/j.gfs.2020.100426.

3. M.W. Barbosa Uncovering research streams on agri-food supply chain management: A bibliometric study., 2021, 28,p. 100517. DOI: https://doi.org/10.1016/j.gfs.2021.100517.

4. S. Kang; G. Kim; J. Roh; E.C. Jeon Ammonia Emissions from NPK Fertilizer Production Plants: Emission Characteristics and Emission Factor Estimation., 2022, 19, 6703. DOI: https://doi.org/10.3390/ijerph19116703.

5. M. Shokouhifar FH-ACO: Fuzzy heuristic-based ant colony optimization for joint virtual network function placement and routing., 2021, 107,p. 107401. DOI: https://doi.org/10.1016/j.asoc.2021.107401.

6. M.A. Al-Qaness; A.A. Ewees; H. Fan; M. Abd Elaziz Optimized forecasting method for weekly influenza confirmed cases., 2020, 17, 3510. DOI: https://doi.org/10.3390/ijerph17103510.

7. M. Khairunissa; H. Lee Hybrid Metaheuristic-Based Spatial Modeling and Analysis of Logistics Distribution Center., 2022, 11, 5. DOI: https://doi.org/10.3390/ijgi11010005.

8. M. Shokouhifar; M. Ranjbarimesan Multivariate time-series blood donation/demand forecasting for resilient supply chain management during COVID-19 pandemic., 2022, 5,p. 100078. DOI: https://doi.org/10.1016/j.clscn.2022.100078.

9. H.A. Shehadeh; H.M.J. Mustafa; M. Tubishat A Hybrid Genetic Algorithm and Sperm Swarm Optimization (HGASSO) for Multimodal Functions., 2022, 13,pp. 1-33. DOI: https://doi.org/10.4018/IJAMC.292507.

10. A. Ilinova; D. Dmitrieva; A. Kraslawski Influence of COVID-19 pandemic on fertilizer companies: The role of competitive advantages., 2021, 71,p. 102019. DOI: https://doi.org/10.1016/j.resourpol.2021.102019.

11. M. Rabbani; S.M.H. Molana; S.M. Sajadi; M.H. Davoodi Sustainable fertilizer supply chain network design using evolutionary-based resilient robust stochastic programming., 2022, 174,p. 108770. DOI: https://doi.org/10.1016/j.cie.2022.108770.

12. R.W. Scholz; A.E. Ulrich; M. Eilittä; A. Roy Sustainable use of phosphorus: A finite resource., 2013, 461,pp. 799-803. DOI: https://doi.org/10.1016/j.scitotenv.2013.05.043.

13. H. Gong; F. Meng; G. Wang; T.E. Hartmann; G. Feng; J. Wu; X. Jiao; F. Zhang Toward the sustainable use of mineral phosphorus fertilizers for crop production in China: From primary resource demand to final agricultural use., 2022, 804,p. 150183. DOI: https://doi.org/10.1016/j.scitotenv.2021.150183.

14. C.M. Godde; D. Mason-D’Croz; D.E. Mayberry; P.K. Thornton; M. Herrero Impacts of climate change on the livestock food supply chain; a review of the evidence., 2021, 28,p. 100488. DOI: https://doi.org/10.1016/j.gfs.2020.100488.

15. X. Yu; C. Keitel; F.A. Dijkstra Global analysis of phosphorus fertilizer use efficiency in cereal crops., 2021, 29,p. 100545. DOI: https://doi.org/10.1016/j.gfs.2021.100545.

16. P.J.A. Withers; R. Sylvester-Bradley; D.L. Jones; J.R. Healey; P.J. Talboys, ACS Publications: Washington, DC, USA, 2014,

17. J. Gaffney; J. Bing; J. Sawyer; P. Byrne; K. Cassman; I. Ciampitti; D. Warner Science-based intensive agriculture: Sustainability, food security, and the role of technology., 2019, 23,pp. 236-244. DOI: https://doi.org/10.1016/j.gfs.2019.08.003.

18. R.B. Chowdhury; G.A. Moore; A.J. Weatherley; M. Arora A novel substance flow analysis model for analysing multi-year phosphorus flow at the regional scale., 2016, 572,pp. 1269-1280. DOI: https://doi.org/10.1016/j.scitotenv.2015.10.055.

19. J. Gurevitch; J. Koricheva; S. Nakagawa; G. Stewart Meta-analysis and the science of research synthesis., 2018, 555,pp. 175-182. DOI: https://doi.org/10.1038/nature25753.

20. A.K.S. De Oliveira; E.B. Soares; M.G. dos Santos; H.A. Lins; M. de Freitas Souza; E. dos Santos Coêlho; L.M. Silveira; V. Mendonça; A.P.A. Barros Júnior; W. de Araújo de Araújo Rangel Lopes Efficiency of Phosphorus Use in Sunflower., 2022, 12, 1558. DOI: https://doi.org/10.3390/agronomy12071558.

21. R. Tabbassum; M. Naveed; I. Mehboob; M.H. Babar; J. Holatko; N. Akhtar; M. Rafique; J. Kucerik; M. Brtnicky; A. Kintl 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., 2022, 12, 2335. DOI: https://doi.org/10.3390/agronomy12102335.

22. H. Fakhravar Combining heuristics and Exact Algorithms: A Review., 2022,

23. E.B. Tirkolaee; A. Goli; P. Ghasemi; F. Goodarzian Designing a sustainable closed-loop supply chain network of face masks during the COVID-19 pandemic: Pareto-based algorithms., 2022, 333,p. 130056. DOI: https://doi.org/10.1016/j.jclepro.2021.130056.

24. M. Sohrabi; M. Zandieh; B. Afshar-nadjafi A simple empirical inventory model for managing with demand differentiation., 2021, 40,pp. 1-38. DOI: https://doi.org/10.1007/s40314-021-01663-8.

25. M. Sohrabi; M. Zandieh; B.A. Nadjafi Dynamic demand-centered process-oriented data model for inventory management of hemovigilance systems., 2021, 27,pp. 73-81. DOI: https://doi.org/10.4258/hir.2021.27.1.73.

26. M. Shokouhifar Swarm intelligence RFID network planning using multi-antenna readers for asset tracking in hospital environments., 2021, 198,p. 108427. DOI: https://doi.org/10.1016/j.comnet.2021.108427.

27. H. Rajabi Moshtaghi; A. Toloie Eshlaghy; M.R. Motadel A comprehensive review on meta-heuristic algorithms and their classification with novel approach., 2021, 8,pp. 63-89.

28. M.A. Farag; M.A. El-Shorbagy; A.A. Mousa; I.M. El-Desoky A new hybrid metaheuristic algorithm for multiobjective optimization problems., 2020, 13,p. 920. DOI: https://doi.org/10.2991/ijcis.d.200618.001.

29. V.K. Chouhan; S.H. Khan; M. Hajiaghaei-Keshteli Metaheuristic approaches to design and address multi-echelon sugarcane closed-loop supply chain network., 2021, 25,pp. 11377-11404. DOI: https://doi.org/10.1007/s00500-021-05943-7.

30. H. Gholizadeh; M. Goh; H. Fazlollahtabar; Z. Mamashli Modelling uncertainty in sustainable-green integrated reverse logistics network using metaheuristics optimization., 2022, 163,p. 107828. DOI: https://doi.org/10.1016/j.cie.2021.107828.

31. S. Karmakar; A. Kundu; B. John Optimizing a Supply Chain Network Using Metaheuristic for Pre and Post Pandemic Scenario.,pp. 41-45.

32. M.P. Nath; S.B.B. Priyadarshini; D. Mishra Supply chain management (SCM): Employing various big data and metaheuristic strategies., Springer: Berlin/Heidelberg, Germany, 2022,pp. 145-165.

33. B. Nahavandi; M. Homayounfar; A. Daneshvar; M. Shokouhifar Hierarchical structure modelling in uncertain emergency location-routing problem using combined genetic algorithm and simulated annealing., 2022, 68,pp. 150-163. DOI: https://doi.org/10.1504/IJCAT.2022.123466.

34. W. Wu; W. Zhou; Y. Lin; Y. Xie; W. Jin A hybrid metaheuristic algorithm for location inventory routing problem with time windows and fuel consumption., 2021, 166,p. 114034. DOI: https://doi.org/10.1016/j.eswa.2020.114034.

35. M. Shokouhifar; N. Pilevari Combined adaptive neuro-fuzzy inference system and genetic algorithm for e-learning resilience assessment during COVID-19 pandemic., 2022, 34,p. e6791. DOI: https://doi.org/10.1002/cpe.6791.

36. Y. Liu; C. Chen Improved RFM Model for Customer Segmentation Using Hybrid Meta-heuristic Algorithm in Medical IoT Applications., 2022, 31,p. 2250009. DOI: https://doi.org/10.1142/S0218213022500099.

37. Z. Ghasemi Darehnaei; M. Shokouhifar; H. Yazdanjouei; S.M.J. Rastegar Fatemi SI-EDTL: Swarm intelligence ensemble deep transfer learning for multiple vehicle detection in UAV images., 2022, 34,p. e6726. DOI: https://doi.org/10.1002/cpe.6726.

38. M. Shokouhifar; M.M. Sabbaghi; N. Pilevari Inventory management in blood supply chain considering fuzzy supply/demand uncertainties and lateral transshipment., 2021, 60,p. 103103. DOI: https://doi.org/10.1016/j.transci.2021.103103.

39. R. Naderi; M. Shafiei Nikabadi; A. Alem Tabriz; M.S. Pishvaee Supply chain sustainability improvement using exergy analysis., 2021, 154,p. 107142. DOI: https://doi.org/10.1016/j.cie.2021.107142.

40. M. Sohrabi; M. Zandieh; M. Shokouhifar Sustainable inventory management in blood banks considering health equity using a combined metaheuristic-based robust fuzzy stochastic programming., 2022, DOI: https://doi.org/10.1016/j.seps.2022.101462.. In press

41. K. Sörensen Metaheuristics—The metaphor exposed., 2015, 22,pp. 3-18. DOI: https://doi.org/10.1111/itor.12001.

42. H. Esmaeili; V. Hakami; B.M. Bidgoli; M. Shokouhifar Application-specific clustering in wireless sensor networks using combined fuzzy firefly algorithm and random forest., 2022, 210,p. 118365. DOI: https://doi.org/10.1016/j.eswa.2022.118365.

43. M. Shokouhifar; A. Jalali Optimized sugeno fuzzy clustering algorithm for wireless sensor networks., 2017, 60,pp. 16-25. DOI: https://doi.org/10.1016/j.engappai.2017.01.007.

44. S. Mirjalili; A. Lewis The whale optimization algorithm., 2016, 95,pp. 51-67. DOI: https://doi.org/10.1016/j.advengsoft.2016.01.008.

45. K.M. Wang; K.J. Wang; C.C. Chen Capacitated production planning by parallel genetic algorithm for a multi-echelon and multi-site TFT-LCD panel manufacturing supply chain., 2022, 127,p. 109371. DOI: https://doi.org/10.1016/j.asoc.2022.109371.

46. A.M. Fathollahi-Fard; K. Govindan; M. Hajiaghaei-Keshteli; A. Ahmadi A green home health care supply chain: New modified simulated annealing algorithms., 2019, 240,p. 118200. DOI: https://doi.org/10.1016/j.jclepro.2019.118200.

Figures and Tables

Figure 1: The PFSCM model. [Please download the PDF to view the image]

Figure 2: Overall flowchart of the proposed H-WOA-VNS algorithm. [Please download the PDF to view the image]

Figure 3: Encoding of the direct DVs into a feasible solution. [Please download the PDF to view the image]

Figure 4: Integer local search operators: (a) I-Swap, and (b) I-Exchange. The variables selected for change in SOL[sup.current] and the changed variables in SOL[sup.new] are shown in blue and green, respectively. [Please download the PDF to view the image]

Figure 5: Permutation operators: (a) Exchange, (b) Relocate, (c) OrOpt, (d) TwoOpt, and (e) Reverse. The selected variables for change are shown in colors. [Please download the PDF to view the image]

Figure 6: Best OF versus iteration in WOA and VNS phases: WOA phase (left), VNS phase (right). [Please download the PDF to view the image]

Figure 7: Convergence of the H-WOA-VNS algorithm in terms of the best OF versus NFE. [Please download the PDF to view the image]

Figure 8: The best OF versus NFE obtained by WOA, VNS, and H-WOA-VNS. [Please download the PDF to view the image]

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. [Please download the PDF to view the image]

Figure 10: Comparison of the worst, mean, and best results, obtained by different methods. [Please download the PDF to view the image]

Table 1: List of indices, sets, parameters, and decision variables.

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 3: Raw materials composition (%) for the different P-fertilizers.

Table 4: Capacity of the suppliers (ton/year).

Table 5: Monthly capacity of the manufacturers for the P-fertilizer production (ton).

Table 6: Parameters of the H-WOA-VNS algorithm.

Table 7: Setting of the PFSCM model parameters.

Table 8: Results of the H-WOA-VNS algorithm in 10 successive runs to solve the sustainable PFSCM model for the Case Study.

Table 9: Consumed raw materials for the production of the different fertilizers (ton).

Table 10: Total production of the P-fertilizers by each manufacturer (ton).

Table 11: Consumed raw materials (ton) for the production of the different fertilizers.

Table 12: Total delivered P-fertilizers (ton) and P-uptake (ton) in each distribution center.

Table 13: Distribution of the different P-fertilizers (ton) in the system.

Table 14: PFSCM test examples.

Table 15: Comparison of the average results of H-WOA-VNS over 10 runs against exact search.

Table 16: Comparison of the H-WOA-VNS algorithm against other techniques, over 10 runs.

Table 17: Effect of the weights of the total objective function on the different sub-objectives.

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.

Author Affiliation(s):

[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

Author Note(s):

[*] Correspondence: m_shokouhifar@sbu.ac.ir (M.S.); frank.werner@ovgu.de (F.W.)

DOI: 10.3390/agronomy13020565