A Flexible Robust Possibilistic Programming Approach for Sustainable Second-Generation Biogas Supply Chain Design under Multiple Uncertainties

Abstract

The goal of this research is to develop a novel second-generation-based biogas supply chain network design (BG-SCND) model that takes into account the triple bottom line approach. Biogas is a promising renewable energy source that can be obtained from a variety of easily accessible second-generation wastes, including animal manure, municipal waste, and agricultural leftovers. Integrated optimization of the biogas generation system is essential for a speedy and environmentally friendly transition to sustainable biodiesel production. The dynamic environment of the energy market significantly impairs the decisions of the BG-SCND model; therefore, a hybrid solution approach using flexible programming and possibilistic programming is suggested. To verify the suggested model and approach for solving the problem, a thorough computational analysis of a case study is conducted. The case study findings demonstrate that considerable investment is necessary to attain social and environmental well-being goals and safeguard decisions against epistemic uncertainty. Policymakers involved in the planning of biogas production and distribution projects may find the proposed approach useful.

1. Introduction

The world is currently witnessing the worst energy shortage in history, which is threatening both the sustainability of the energy supply and the advancement of the economy. Coal, petroleum products, and other conventional resources such as methane are all nonrenewable. Moreover, the need to investigate alternative sources has been reinforced by declining environmental health, rising import costs, and exponential growth in energy prices. To address the global energy issue, alternative energy supplies that are abundant, eco-friendly, and socially sustainable are required. Biofuels have been suggested as a substitute by experts and policymakers due to their wide commercial applications. Biogas, one of the available biofuels, is particularly significant, since it has a very high calorific value and may be used directly as fuel or indirectly to produce power [1,2]. Biogas is the combination of gases made through the breakdown of organic materials with the nonappearance of oxygen (anaerobically) [3]. Anaerobic digestion has been rated as one of the most energy-efficient and ecologically friendly methods for producing biogas and has many benefits [4]. By using locally accessible biomass, it may significantly cut carbon emissions when compared to fossil fuels, which have spurred significant public and private interest globally [5]. Biogas can be produced from different types of waste, such as animal waste, municipal solid waste, and various types of crop residuals. Hence, using these wastes as biomass for the production of biogas can be both an excellent waste control strategy and a cheap source of energy generation [6].

In addition to employing inexpensive feedstock and processing technologies, an effective system of the biogas supply chain (SC) is crucial in enabling large-scale biodiesel production to be financially, environmentally, and socially viable. In this context, the increasing development of the business has prompted policymakers to employ SC network design optimization models by adopting a sustainable approach [7]. All supply chain network echelons must be arranged in a way that ensures the effective flow of resources and information to effectively turn biomass into biogas. Waste collection, pretreatment processing, storage, and transferring biomass, as well as converting biomass to biogas, and its distribution in the market are different stages of the biogas SC. An effective SC network design necessitates centralized strategic, tactical, and operational-level planning to maximize the economic feasibility of the biogas supply chain. In the past, awful impacts of commercialization on the environment have led all industries toward sustainable growth, while concurrently considering socio-economic aspects [8,9].

In this context, it is vital to establish a decision support system (DSS) that can aid policymakers in achieving both commercialization and sustainable growth goals. Additionally, due to the highly complex market, biogas SCs are more vulnerable to parameter uncertainty than commercial SCs, since biomass supplies are wastes, which are seasonal and dependent on the use of primary items, while the price of biogas and biomass, as well as transportation costs, are determined by global fluctuations in crude oil prices. Looking at the research related to biogas SC design, most of them have ignored the crucial aspect of operational uncertainty. Although some of them have taken a stochastic approach to address this issue, stochastic techniques require previous data to build probability distributions for imprecise parameters, which are not always accessible for all parameters. In light of this, experts have suggested employing a fuzzy technique that does not require previous data and is based on the possibility distribution. There are only a couple of studies that have employed a fuzzy solution approach to designing a biogas SC. Furthermore, there is no study in the field of biogas supply chain network design (SCND) that has provided a DSS for managers and policymakers to achieve an efficient trade-off among the three sustainability aspects of economics, environment, and society by using a fuzzy optimization technique. To bridge the research gap, this research proposes a multi-objective mixed-integer linear programming (MILP)-based biogas supply chain network design (BG-SCND) model, considering operational environment uncertainty. The BG-SCND model’s goal is to maximize social benefits and reduce carbon emissions from biogas SC-related operations, while also reducing the system’s overall SC cost. Against this background, the following research questions are addressed in this study.

• How can efficient compromise among sustainability dimensions in a biogas SC network be attained based on the preferences of policymakers?
• How can efficient and robust results for a multi-objective BG-SCND model be obtained in a highly dynamic environment?
• How can guidelines for the social life cycle assessment of products (GSLCAP) theory be employed to quantify the social impact of biogas SC operations?

The remainder of the research is structured as follows. In Section 2, a thorough overview of the related literature is provided. The problem definition and the BG-SCND model are presented in Section 3. The FRPP approach is explained in Section 4. A case study to validate the BG-SCND model and FRPP solution approach is detailed in Section 5. Finally, Section 6 concludes this research.

2. Literature Review

This section begins by summarizing the most recent SCND models and solution techniques proposed for various types of biofuels, such as methanol, biodiesel, and bioethanol. Following this, a thorough analysis of the SCND models and solution approaches for biogas research is provided. Mathematical modeling-based approaches have been used in several studies. Among them, MILP has been widely used to design methanol, biodiesel, and bioethanol distribution networks.

In the category of bioethanol SCND, Akgul et al. [10] provided a static modeling-based optimization model for the development of a bioethanol SC network, analyzing factors such as the influence of carbon duty on economic and environmental goals, as well as the trade-off of economic and environmental dimensions. An optimization model was presented by Ghaderi et al. [11] to examine the economic, social, and environmental impact of a bioethanol SC. To deal with the SC’s unpredictable environment, the authors developed a possibilistic programming (PP)-based method that maximizes the mean value of the SC’s performance. Gonela et al. [12] suggested a hybrid bioethanol SC network design model that accounts for sustainability dimensions using a stochastic uncertainty technique. Marvin et al. [13] also worked on the maximization of the net present value of a bioethanol SC by formulating the network design problem as a MILP model in the USA region. The provided model assists in choosing the best site for biorefineries, their ideal capacity, and their feedstock cultivation.

Many researchers have worked on developing optimization models for biodiesel SCs. In this regard, Babazadeh [14] created a multi-product biodiesel SCND model using a deterministic solution approach to minimize the total cost of the system. Ghelichi et al. [15] suggested a scenario-based flexible stochastic approach to reduce the present value cost for a Jatropha-based biodiesel SCND model. Mirhashemi et al. [16] developed a two-stage optimization approach to enhance the commercialization of oleifera-based biodiesel. In the first step, the proposed study uses data envelopment analysis (DEA) to rank appropriate oleifera growing sites, followed by a MILP model to derive strategic and tactical-level decisions of investigated SCs. Habib et al. [17] assessed the feasibility of waste animal fat to produce biodiesel for the region of Pakistan. For this purpose, the author developed an SC model intended to minimize overall costs. An RPP-based solution technique was provided to address the operational uncertainty in the SC model. This study was further extended by Habib et al. [18] by taking into account disruption risk linked with the previously investigated biodiesel production and distribution model. Mohtashami et al. [19] proposed a two-stage optimization model for a Jatropha-based biodiesel SC. In the first step, a common-weight DEA approach was used to identify potential Jatropha-growing locations, and then a multi-objective model based on the sustainability paradigm was used to obtain decisions. An ε-constraint approach was applied to obtain an efficient trade-off solution among sustainability dimensions.

For methanol SCNDs, a few authors have proposed decision-making models. In this regard, Villicaña-García et al. [20] developed a mathematical model to minimize water usage and carbon emissions, and maximize profit for the strategic planning of the methanol SC, while accounting for incentives to use local resources. Nugroho et al. [21] proposed an agent-based simulation model to investigate the techno-economic and life cycle evaluation of the methanol supply chain. The proposed methanol SC model’s primary considerations include the biomass and methanol ordering strategy, capital budgeting of the methanol SC, and CO₂ emissions/ton of methanol. Leduc et al. [22] used spatial modeling to analyze the cost of the methanol supply chain in the German region. For this reason, a variety of costs related to the methanol SC, including biomass collection, handling, shipping, and methanol production, were considered.

Several authors have proposed mathematical models for biogas SCNDs. In this aspect, Sarker et al. [23] used a non-linear mixed-integer model to optimize a biomethane gas SC network design with varied feedstock crops, grass, wood wastes, and animals. The suggested solution does not account for the uncertainty associated with the model input parameters and was solved using a genetic algorithm. Jensen et al. [24] proposed a deterministic MILP model for biogas SC network design based on animal manure and crop waste, and calculated operational and capital investments from the place of production to the market zone. Egieya et al. [25] also used a deterministic MILP technique for biogas SC design and used agricultural waste as a feedstock. Three objectives—reducing overall costs, reducing greenhouse gas emissions, and maximizing profit—were considered. Garbs and Geldermann [26] provided a deterministic model and used a linear cost optimization technique to investigate the environmental and economic effects of manure transportation to a biogas plant. Díaz-Trujillo and Nápoles-Rivera [27] suggested a deterministic MILP model for biogas SC design, intending to maximize profit and minimize environmental impact. A similar deterministic approach for biogas SCNDs was adopted by Park et al. [28], Balaman and Selim [29], Cheraghalipour and Roghanian [30], and Galvez et al. [31]. All of the biogas SCND research mentioned above uses a deterministic modeling technique to provide optimized solutions; nevertheless, the robustness of their solutions cannot be guaranteed.

A few publications, such as Zirngast et al. [32] and Khishtandar [33], used a stochastic strategy to cope with operational uncertainty. However, stochastic techniques need a large quantity of past data in order to accurately predict the distribution. Therefore, fuzzy-based approaches have been found very effective when used to deal with uncertain circumstances. In this aspect, Franco et al. [34] adopted a fuzzy-based weighted overlap dominance approach to obtain location decisions of biogas-processing facilities. Yilmaz Balaman and Selim [35] presented a MILP model to build an SC network for the anaerobic digestion of biomass for biogas generation. The model efficiently determined the quantity and location of biogas production facilities, as well as biomass supply and distribution decisions. Fuzzy goal programming was employed as a solution approach. This study was further extended by Yılmaz Balaman and Selim [36] by accounting for thermal energy shortages in terms of service-level objectives for biogas markets. The authors also considered a fuzzy environment for the proposed biogas SC model and used fuzzy goal programming to find a solution. A comparative analysis of the existing biogas research is provided in Table 1.

Examining research related to biogas production indicates that most of them employed a deterministic MILP modeling strategy, while only a small number of them used a fuzzy approach to take into account uncertainty in the biogas SC parameters. The failure to account for these uncertainties in the biogas SC model may result in unrealistic or suboptimal results. To address the challenge, fuzzy PP-based uncertainty is used to represent the uncertain behavior of BG-SCND model parameters because there is a lack of credible historical data regarding these parameters. According to Pishvaee et al. [37] and Torabi. et al. [38], the most effective method for addressing parameter uncertainty when there is a lack of historical data is the PP technique, which incorporates both subjective and objective data. Additionally, by employing the flexible programming (FP) technique, the idea of soft constraints is also used, and supply, demand, and capacity constraints are all modeled as soft constraints. Furthermore, the final results’ robustness is critical, especially for strategic planning (i.e., those about the network and facility location). Therefore, a flexible robust possibilistic programming (FRPP) technique is developed that significantly benefits from both FP and robust possibilistic programming (RPP) advantages in order to obtain robust solutions while addressing various kinds of uncertainties.

To the best of the authors’ knowledge, this is the first study in the domain of biogas SC network design to use a flexible RPP technique to acquire decisions that concurrently take into account economic, environmental, and social factors. Further, this study is the first to measure the social effects of biogas SCND operations using the most recent social life cycle assessment approach, known as GSLCAP. Another drawback of current studies is that existing biogas SC models only take into account the biomass supply to the biorefinery and ignore the biorefinery to digestate markets phase, which might be crucial for the economic viability of biogas. Based on the identified gap, the following are contributions of this research:

• Developing a multi-period biogas SC optimization model to obtain integrated strategic, tactical, and operational-level decisions, while simultaneously considering the economic, environmental, and social dimensions of sustainability.
• Proposing an interactive flexible robust possibilistic programming technique that significantly benefits from both flexible programming and robust programming and provides robust strategic and tactical-level decisions under multiple uncertainties.
• Using the GSLCAP technique to evaluate the social dimension of sustainability, which takes a systemic approach by considering job creation potential, unemployment rate, economic investment decisions, and development level to uplift the underprivileged regions of the biogas SC.

Table 1. A comparative analysis of the existing biogas SCND research.

Author	Type of Feedstock	Decision Levels	Method/Analysis	Uncertainty Type		Sustainability Dimension Considered			Major Supply Chain Decisions Considered
				Stochastic	Fuzzy	Environmental	Economic	Social
Sarker et al. [23]	Crops, grass, wood wastes, and animal manure	Strategic, tactical	Genetic algorithm.	🗴	🗴		✓		Location–allocation, reactor demand, and natural loss of feedstock
Jensen et al. [24]	Crop waste	Strategic, tactical	Deterministic MILP	🗴	🗴		✓		Location–allocation, capacity, and material flow
Egieya et al. [25]	Crop residue, manure, grass	Operational, strategic, tactical	Deterministic MILP	🗴	🗴	✓	✓	✓	Location, material flow, transport modes, and conversion technologies
Garbs and Geldermann [26]	Poultry	Strategic, tactical	Microsoft Excel-based modeling	🗴	🗴	✓	✓		Transport model, location, and material flow
Díaz-Trujillo and Nápoles-Rivera [27]	Manure, organic waste	Strategic, tactical	Deterministic MILP	🗴	🗴	✓	✓		Location, material flow, and technology,
Park et al. [28]	Manure	Strategic	Deterministic MILP	🗴	🗴	✓	✓		Location–allocation and capacity
Balaman and Selim [29]	Agricultural, animal, and organic wastes	Operational, strategic, tactical	Deterministic MILP	🗴	🗴		✓		Location–allocation, capacities, biomass storage
Cheraghalipour and Roghanian [30]	Agricultural waste	Strategic and operational	Genetic algorithm and stochastic fractal search algorithm	🗴	🗴		✓		Location–allocation and capacity
Galvez et al. [31]	Domestic waste, agriculture, and livestock waste	Strategic and operational	Deterministic MILP and analytical hierarchy process	🗴	🗴		✓		Location–allocation and material flow
Zirngast et al. [32]	Animal manures and crop residual	Strategic and tactical	Monte Carlo simulation	✓		✓	✓		Location–allocation and material flow, capacity, technology, and transportation mode
Khishtandar [33]	Agricultural waste, domestic waste, and animal waste	Strategic and operational	Monte Carlo simulation and genetic algorithm	✓			✓		Location–allocation, material flow
Yilmaz Balaman and Selim [35]	Animal wastes and energy crops	Strategic, tactical	Fuzzy goal programming		✓	✓	✓		Location–allocation, storage, biomass cultivation, capacity
Yılmaz Balaman and Selim [36]	Animal manure and energy crops	Operational, strategic, tactical	Fuzzy goal programming		✓		✓		Location–allocation, capacity, material flow
This study	Sawdust, wheat straw, bagasse, rice husk	Operational, strategic, tactical	Fuzzy flexible robust possibilistic programming-based approach		✓	✓	✓	✓	Location–allocation, storage, capacity, material flow, digestate, GSLCAP-based social aspect estimation

Table 1. A comparative analysis of the existing biogas SCND research.

Author	Type of Feedstock	Decision Levels	Method/Analysis	Uncertainty Type		Sustainability Dimension Considered			Major Supply Chain Decisions Considered
				Stochastic	Fuzzy	Environmental	Economic	Social
Sarker et al. [23]	Crops, grass, wood wastes, and animal manure	Strategic, tactical	Genetic algorithm.	🗴	🗴		✓		Location–allocation, reactor demand, and natural loss of feedstock
Jensen et al. [24]	Crop waste	Strategic, tactical	Deterministic MILP	🗴	🗴		✓		Location–allocation, capacity, and material flow
Egieya et al. [25]	Crop residue, manure, grass	Operational, strategic, tactical	Deterministic MILP	🗴	🗴	✓	✓	✓	Location, material flow, transport modes, and conversion technologies
Garbs and Geldermann [26]	Poultry	Strategic, tactical	Microsoft Excel-based modeling	🗴	🗴	✓	✓		Transport model, location, and material flow
Díaz-Trujillo and Nápoles-Rivera [27]	Manure, organic waste	Strategic, tactical	Deterministic MILP	🗴	🗴	✓	✓		Location, material flow, and technology,
Park et al. [28]	Manure	Strategic	Deterministic MILP	🗴	🗴	✓	✓		Location–allocation and capacity
Balaman and Selim [29]	Agricultural, animal, and organic wastes	Operational, strategic, tactical	Deterministic MILP	🗴	🗴		✓		Location–allocation, capacities, biomass storage
Cheraghalipour and Roghanian [30]	Agricultural waste	Strategic and operational	Genetic algorithm and stochastic fractal search algorithm	🗴	🗴		✓		Location–allocation and capacity
Galvez et al. [31]	Domestic waste, agriculture, and livestock waste	Strategic and operational	Deterministic MILP and analytical hierarchy process	🗴	🗴		✓		Location–allocation and material flow
Zirngast et al. [32]	Animal manures and crop residual	Strategic and tactical	Monte Carlo simulation	✓		✓	✓		Location–allocation and material flow, capacity, technology, and transportation mode
Khishtandar [33]	Agricultural waste, domestic waste, and animal waste	Strategic and operational	Monte Carlo simulation and genetic algorithm	✓			✓		Location–allocation, material flow
Yilmaz Balaman and Selim [35]	Animal wastes and energy crops	Strategic, tactical	Fuzzy goal programming		✓	✓	✓		Location–allocation, storage, biomass cultivation, capacity
Yılmaz Balaman and Selim [36]	Animal manure and energy crops	Operational, strategic, tactical	Fuzzy goal programming		✓		✓		Location–allocation, capacity, material flow
This study	Sawdust, wheat straw, bagasse, rice husk	Operational, strategic, tactical	Fuzzy flexible robust possibilistic programming-based approach		✓	✓	✓	✓	Location–allocation, storage, capacity, material flow, digestate, GSLCAP-based social aspect estimation

3. Biogas Supply Chain Model

This paper presents a supply chain planning problem for second-generation biogas. For this purpose, a BG-SCND model is developed that takes into account economic, social, and environmental aspects, while dealing with different operational uncertainties. The considered SC network comprises five echelons: biomass supply points, biomass collection centers, biogas production plants, biogas distribution and storage centers, and digestate and biogas market zones. The framework of the proposed BG-SCND model is provided in Figure 1.

Biomass from feedstock supply locations is collected and sent to biomass collection centers, where it is stored and then delivered to biogas production plants based on demand. At a biogas plant, biomass is processed, biogas is produced through anaerobic digestion, and digestate is obtained as a byproduct. Biogas is transported from biogas production plants to biogas distribution and storage facilities, whereas digestate is transported straight to digestate market zones. Finally, biogas is delivered from distribution and storage centers to market zones (customers) to meet market demand. Figure 2 outlines the comprehensive working of the considered SC and the indices used for decision making in the BG-SCND model.

3.1. Notations

Indices

• f  biomass type
• p biomass supply centers
• c biomass collection centers
• b biogas production plants
• m biogas market zones
• n digestate market zones
• s biogas distribution and storage centers
• t  time period
• a production capacity of biogas production plants
• u capacity level of biogas distribution and storage centers

Parameters

• $I\tilde{C}{B}_{ba}$ installation cost of biogas production plant b with capacity level a ($)
• $I\tilde{C}{C}_c$  installation cost of feedstock collection center c ($)
• $I\tilde{C}{S}_{su}$  installation cost of biogas distribution and storage center s with capacity u ($)
• $C\tilde{F}{T}_{pft}$  cost of feedstock type f obtained at feedstock supply center p in time period t ($/ton)
• $H\tilde{C}{C}_{cft}$  handling and storage cost of biomass f at collection center c in period t ($/ton)
• $C\tilde{E}{H}_c$  carbon emissions in inventory handling at feedstock collection center c (kg of CO₂/ton)
• $P\tilde{C}{B}_b$  production cost per ton of biogas at production plant b ($/ton)
• $TC{F}_{pc}$  transportation cost of feedstock from feedstock supply center p to collection center c ($/ton.km)
• $TC{F}_{cb}$  transportation cost of feedstock from feedstock collection center c to biogas production plant b ($/ton.km)
• $TC{B}_{bs}$  transportation cost of biogas from biogas production plant b to biogas storage center s ($/ton.km)
• $TC{B}_{sm}$  transportation cost of biogas from biogas storage center s to market m ($/ton.km)
• $TC{D}_{bn}$  transportation cost of digestate from biogas production plant b to market zone n ($/ton.km)
• $C_{tax}$  carbon emissions tax ($/kg of CO₂)
• $C\tilde{E}{P}_b$  quantity of carbon emissions per kg of biomass pretreatment at biogas production plant b (kg of CO₂/kg of biomass)
• $CE{T}_{pc}$  carbon emissions during feedstock transportation from supply center p to feedstock collection center c (kg of CO₂/ton)
• $CE{T}_{cb}$  carbon emissions during feedstock transportation from feedstock collection center c to biogas production plant b (kg of CO₂/ton)
• $CE{T}_{bs}$  carbon emissions during biogas transportation from biogas production plant b to biogas distribution and storage center s (kg of CO₂/ton)
• $CE{T}_{sm}$  carbon emissions during biogas transportation from biogas distribution and storage centers to market zone m (kg of CO₂/ton)
• $CE{T}_{bn}$  carbon emissions during digestate transportation from biogas production plant b to market zone n (kg of CO₂/ton)
• $D_{mt}$  demand of market zone m for biogas in time period t (ton/period)
• $D_{nt}$  demand of market zone n for digestate in time period t (ton/period)
• $Q{F}_{pft}$  quantity available for feedstock type f at feedstock supply center p in time t
• ${\eta }_f$  yield factor of feedstock type f
• $C{P}_{cft}$   capacity of collection center c for feedstock f in period t
• $C{P}_{ba}$  capacity of biogas plant b with level a
• $C{P}_{su}$  capacity of biogas distribution center s with level u
• $In{v}_c$  investment in region c to establish a feedstock collection center ($)
• $Dev{l}_c$  development level of region c
• $In{v}_{ba}$  investment in region b to establish a biogas plant with capacity a ($)
• $Dev{l}_b$  development level of region b
• $In{v}_{su}$  investment in region s to establish a distribution center with capacity u ($)
• $Dev{l}_s$  development level of region s
• $Un{e}_c$  rate of unemployment in region c
• $Jo{b}_c$  job creation potential of feedstock collection center c
• $Un{e}_b$  rate of unemployment in region b
• $Jo{b}_{ba}$  job creation potential of a biogas plant in region b with capacity a
• $Un{e}_s$  rate of unemployment in region s
• $Jo{b}_{su}$  job creation potential of distribution center in region s with capacity u
• ${\mathrm{\Phi }}^{EGI}$  priority weight of economic growth index
• ${\mathrm{\Phi }}^{JPI}$  priority weight of job creation potential index

Decision  variables

• $F{T}_{pcft}$  amount of feedstock type f transported from feedstock supply center p to biomass collection center c in time period t (ton/period)
• $F{T}_{cbft}$  amount of biomass type f supplied from collection center c to biogas production plant b in time period t (ton/period)
• $B{T}_{bst}$  amount of biogas transported from biogas production plant b to biogas distribution and storage center s in time period t (ton/period)
• $B{T}_{smt}$  amount of biogas transported from biogas distribution and storage center s to market zone m in time period t (ton/period)
• $D{T}_{bnt}$  amount of digestate transported from biogas production plant b to digestate market zone n in time period t (ton/period)
• $AB{P}_{bt}$  amount of biogas produced at biogas production plant b in time period t (ton/period)
• $AD{P}_{bqt}$  amount of digestate produced at biogas production plant b in time period t (ton/period)
• $w_c$  1 if feedstock collection center c is selected; otherwise, 0
• $y_{ba}$  1 if biogas production plant b with capacity level a is selected; otherwise, 0
• $z_{su}$  1 if biogas storage center s with storage capacity level u is selected; otherwise, 0

3.2. Assumptions

Following are the assumptions for the proposed BG-SCND model:

• The transportation distances between each echelon are known.
• Biogas and digestate demand and transportation costs are considered uncertain.
• Carbon emission tax is enacted following regional government policy for all stakeholders.
• Each terminal in the BG-SCND model has a homogenous fleet of vehicles.

3.3. Economic Objective Function

The economic objective function seeks to reduce the overall costs of the suggested supply chain model. The first three elements of the objective function represent the installation costs of biomass collection centers, biogas production facilities, and biogas distribution and storage centers, respectively. The inclusion of binary variables ensures that the echelons are only developed if they are economically viable. The fourth term indicates the overall purchase cost at the individual feedstock supply centers. The fifth term represents the feedstock handling cost at the feedstock collection center and its carbon emissions cost. The sixth term reflects the total cost of biogas generation, as well as the expenses of carbon emissions. The seventh term represents the cost of transporting feedstock from feedstock supply centers to feedstock collection sites, as well as the cost of carbon emissions. The eighth term refers to the cost of transporting feedstock from feedstock collection locations to biogas production plants, coupled with the cost of carbon emissions. The ninth term depicts the cost of transporting biogas from production plants to distribution and storage centers, as well as the cost of carbon emissions. The tenth term indicates the cost of transporting biogas from distribution and storage centers to biogas demand zones, including the cost of carbon emissions. The eleventh term displays the cost of transporting digestate from biogas production facilities to digestate demand zones, as well as the costs of carbon emissions.

$$\begin{array}{l}Minimize E(f_{.TC})={\displaystyle \sum _c^{C}I\tilde{C}{C}_c}\times w_c+{\displaystyle \sum _b^B{\displaystyle \sum _a^{A}I\tilde{C}{B}_{ba}\times y_{ba}}}+{\displaystyle \sum _s^S{\displaystyle \sum _u^{U}I\tilde{C}{S}_{su}\times z_{su}}}+\\ {\displaystyle \sum _p^P{\displaystyle \sum _c^C{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[C\tilde{F}{T}_{pft}\right]}}}}F{T}_{pcft}+{\displaystyle \sum _p^P{\displaystyle \sum _c^C{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[H\tilde{C}{C}_{cft}+\left(C\tilde{E}{H}_c\times C_{tax}\right)\right]}}}}F{T}_{pcft}+\\ {\displaystyle \sum _b^B{\displaystyle \sum _t^T\left[P\tilde{C}{B}_b+\left(C\tilde{E}{P}_b\times C_{tax}\right)\right]}}AB{P}_{bt}+{\displaystyle \sum _p^P{\displaystyle \sum _c^C{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[T\tilde{C}{F}_{pc}+\left(C\tilde{E}{T}_{pc}\times C_{tax}\right)\right]}}}}F{T}_{pcft}+\\ {\displaystyle \sum _c^C{\displaystyle \sum _b^B{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[T\tilde{C}{F}_{cb}+\left(C\tilde{E}{T}_{cb}\times C_{tax}\right)\right]}}}}F{T}_{cbft}+{\displaystyle \sum _b^B{\displaystyle \sum _s^S{\displaystyle \sum _t^T\left[T\tilde{C}{B}_{bs}+\left(C\tilde{E}{T}_{bs}\times C_{tax}\right)\right]}}}B{T}_{bst}+\\ {\displaystyle \sum _s^S{\displaystyle \sum _m^M{\displaystyle \sum _t^T\left[T\tilde{C}{B}_{sm}+\left(C\tilde{E}{T}_{sm}\times C_{tax}\right)\right]}}}B{T}_{smt}+{\displaystyle \sum _b^B{\displaystyle \sum _n^N{\displaystyle \sum _t^T\left[T\tilde{C}{D}_{bn}+\left(C\tilde{E}{T}_{bn}\times C_{tax}\right)\right]}}}D{T}_{bnt}\end{array}$$

(1)

3.4. Social Objective Function

The provided social objective of the BG-SCND consists of two primary components: an economic growth index and a job creation potential index. The economic growth index (EGI) has two key components. The first is the overall investment in prospective SC facilities located in specific regions, and the second is the level of development in the region where investment is being made. The EGI strives to construct SC facilities with the maximum possible investment in underdeveloped regions. The EGI is expressed in Equation (2):

$$EGI={\displaystyle \sum _c^{C}I\tilde{n}{v}_c\times w_c(1-Dev{l}_c)}+{\displaystyle \sum _b^B{\displaystyle \sum _a^{A}I\tilde{n}{v}_{ba}\times y_{ba}(1-Dev{l}_b)+}}{\displaystyle \sum _s^S{\displaystyle \sum _u^{U}I\tilde{n}{v}_{su}\times z_{su}(1-Dev{l}_s)}}$$

(2)

The second component of the social objective is the job potential index (JPI). The JPI considers two important parameters: the first is the potential for job creation by an SC facility in a prospective region, and the second is the rate of unemployment in the region where SC facilities are established. The JPI seeks to enhance job creation potential by focusing on locations with low unemployment rates. The following equation represents the JPI formulation:

$$JPI={\displaystyle \sum _c^{C}U\tilde{n}{e}_c\times J\tilde{o}{b}_c\times w_c}+{\displaystyle \sum _b^B{\displaystyle \sum _a^{A}U\tilde{n}{e}_b\times J\tilde{o}{b}_{ba}\times y_{ba}+}}{\displaystyle \sum _s^S{\displaystyle \sum _u^{U}U\tilde{n}{e}_s\times J\tilde{o}{b}_{su}\times z_{su}}}$$

(3)

Normalization is necessary for both components, the EGI and PGI, to obtain the overall value of the social objective function, which is represented by Equations (4) and (5).

$$EG{I}^{nrm}=\left[\frac{EGI-EG{I}^{\min }}{EG{I}^{\max }-EG{I}^{\min }}\right]$$

(4)

$$JP{I}^{nrm}=\left[\frac{JPI-JP{I}^{\min }}{JP{I}^{\max }-JP{I}^{\min }}\right]$$

(5)

The normalized values of the EGI and PGI are organized as follows to obtain the total social impact objective value of the BG-SCND model activities.

$$Minimize E(f_{soc})={\mathrm{\Phi }}^{EGI}EG{I}^{nrm}+{\mathrm{\Phi }}^{JPI}JP{I}^{nrm}$$

(6)

3.5. Constraints of the BG-SCND Model

Equation (7) fulfills the biogas demand for each market zone in period t.

$${\displaystyle \sum _s^{S}B{T}_{smt}}\ge {\tilde{D}}_{mt} \forall m,t$$

(7)

Equation (8) fulfills the digestate demand for all markets in each period.

$${\displaystyle \sum _b^{B}D{T}_{bnt}}\ge {\tilde{D}}_{nt} \forall n,t$$

(8)

Equation (9) assures that the amount of feedstock type f transported to the feedstock collection center should be less than the quantity existing at the feedstock supply center p.

$$\tilde{Q}{F}_{pft}\ge {\displaystyle \sum _c^{C}F{T}_{pcft} \forall p,f,t}$$

(9)

Equation (10) assures that the amount of biomass f sent from collection center c to biogas production plan b must be less than the quantity of feedstock supplied from supply center p to collection center c during each time period t.

$${\displaystyle \sum _p^{P}F{T}_{pcft} \ge {\displaystyle \sum _b^{B}F{T}_{cbft}} \forall c,f,t}$$

(10)

Equations (11) and (12) represent biomass conversion to biogas and digestate.

$${\displaystyle \sum _c^C{\displaystyle \sum _f^{F}F{T}_{cbft}\times {\eta }_f}}=AB{P}_{bt} \forall b,t$$

(11)

$${\displaystyle \sum _c^C{\displaystyle \sum _f^{F}F{T}_{cbft}\times (1-{\eta }_f})}=AD{P}_{bt} \forall b,t$$

(12)

Equations (13) and (14) require that the amount of biogas and digestate supplied from the biogas plant cannot exceed the total amount produced.

$$AB{P}_{bt}\ge {\displaystyle \sum _s^{S}B{T}_{bst}} \forall b,t$$

(13)

$$AD{P}_{bt}\ge {\displaystyle \sum _n^{N}D{T}_{bnt}} \forall b,t$$

(14)

Equation (15) ensures that the total amount of biogas carried from the biogas distribution center to the market is always less than the total amount of biogas delivered from the biogas plant to the biogas storage and distribution center.

$${\displaystyle \sum _b^{B}B{T}_{bst}-{\displaystyle \sum _m^{M}B{T}_{smt}}\ge 0 \forall s,t}$$

(15)

Equations (16)–(18) limit the processing capacities of functional biomass collection centers, biogas plants, and biogas distribution centers, respectively.

$${\displaystyle \sum _p^{P}F{T}_{pcft}\le C{P}_{cft}}\times w_c \forall c,f,t$$

(16)

$${\displaystyle \sum _c^C{\displaystyle \sum _f^{F}F{T}_{cbft}}}\times {\eta }_f\le {\displaystyle \sum _a^{A}C{P}_{ba}{y}_{ba}} \forall b,t$$

(17)

$${\displaystyle \sum _b^{B}B{T}_{bst}\le {\displaystyle \sum _u^{U}C{P}_{su}{z}_{su}} \forall s,t}$$

(18)

Equation (19) allows the operating biogas plant to select one capacity level and Equation (20) restricts the selected storage and distribution facility to a single capacity level.

$${\displaystyle \sum _a^A{y}_{ba}\le 1 \forall b}$$

(19)

$${\displaystyle \sum _u^U{z}_{su}\le 1 \forall s}$$

(20)

Constraints (21)–(23) assure that all kinds of considered decision variables are non-negative.

$$F{T}_{pcft}, F{T}_{cbft}, B{T}_{bst}, B{T}_{smt}, D{T}_{bnt}, AB{P}_{bt}, AD{P}_{bt}\ge 0$$

(21)

Binary and non-negativity.

$$w_c\in \left[0,1\right] \forall c$$

(22)

$$y_{ba},z_{su}\in \left[0,1\right] \forall b,a,s,u$$

(23)

4. Proposed Solution Methodology

Due to the prevalence of multiple types of uncertainties in the SC network design, a combination of several methodologies to uncertainty programming has been increasingly employed in recent years. Against this background, a hybrid solution approach, FRPP, which is developed by merging flexible programming and robust probabilistic programming, is used to deal with the diverse nature of uncertainties involved in the BG-SCND model. The generic formulation of the FRPP technique is systematically described here.

4.1. Possibilistic Programming (PP) Model

A generic uncertain SC mathematical model is provided in Equation (24).

$$\begin{array}{l}\mathrm{Minimize} q=\tilde{a}m+\tilde{b}n\\ \mathrm{Such} \mathrm{that} Cn\le \tilde{r}\\ En=0\\ Fn\le \tilde{G}m\\ Hn\le 1\\ m\in \{0,1\}, n\ge 0\end{array}$$

(24)

where a, b, and r represent uncertain installation and shipping expenses, and demands, respectively, which follow the triangular fuzzy number. The constraints’ coefficients are represented by the matrices C, E, G, F, and H. Furthermore, vector m represents binary parameters, and vector n represents continuous variables. According to Habib [39], if D is an uncertain number with a triangular possibility distribution $\tilde{D} = (D^{opt},D^{most},D^{pes})$, its optimistic $(D^{opt})$, most likely $(D^{most})$, and pessimistic $(D^{pes})$ values for TFN can be calculated using Equations (25) and (26) [40,41].

$$D^{opt}=(1-{\partial }_1)D^{most}$$

(25)

$$D^{pes}=(1+{\partial }_1)D^{most}$$

(26)

A graphical representation of triangular possibility distribution for the uncertain parameter $\tilde{D} = (D^{opt},D^{most},D^{pes})$ is provided in Figure 3.

The membership function for the uncertain parameter D following TFN can be developed as follows:

$${\mu }_{\tilde{D}}(v)=\left\{\begin{array}{l}\frac{v-D^{pes}}{D^{mos}-D^{pes}}, if D^{pes}\le v\le D^{mos}\\ 1, if v=D^{mos}\\ \frac{D^{opt}-v}{D^{opt}-D^{mos}}, if D^{mos}\le v\le D^{opt}\\ 0, if v\le D^{pes} or v\ge D^{opt}\end{array}\right.$$

(27)

(a)

Expected value (ExpV) operator

To transform uncertain parameters of the objective function, the ExpV operator proposed by Jiménez et al. [42] is employed. The EV for an uncertain parameter can be obtained as follows:

$$ExpV(\tilde{D})=\frac{D^{pes}+2D^{mos}+D^{opt}}{4}$$

(28)

(b)

Me-measure

To obtain the corresponding certain form of uncertain constraints of the BG-SCND model, the Me-measure provided by Xu and Zhou [43] is employed. The Me-measure efficiently allows a manager to pick any point on the spectrum of the Nec-measure and Pos-measure. If $\tilde{\partial }$ is an uncertain parameter of constraint, then its equivalent Me-measure form can be obtained as follows:

$$Me\left\{\tilde{\partial }\ge v\right\}=Nec\left\{\tilde{\partial }\ge v\right\}+\pi \left[Pos\left\{\tilde{\partial }\ge v\right\}-Nec\left\{\tilde{\partial }\ge v\right\}\right]$$

(29)

where $\pi$ controls the behavior of the SC manager on the spectrum of Nec-measure and Pos-measure. Me-measure transformations for the following cases are given below:

$$Me\left\{\tilde{\partial }\ge v\right\}= \left\{\begin{array}{l} 1 v\le {\partial }_1\\ \pi +(1-\pi )\frac{{\partial }_2-v}{{\partial }_2-{\partial }_1} {\partial }_1\le v\le {\partial }_2\\ \pi \frac{{\partial }_3-v}{{\partial }_3-{\partial }_2} {\partial }_2\le v\le {\partial }_3\\ 0 v\ge {\partial }_3 \end{array}\right\}$$

(30)

$$Me\left\{\tilde{\partial }\le v\right\}= \left\{\begin{array}{l} 0 v\le {\partial }_1\\ \pi \frac{v-{\partial }_1}{{\partial }_2-{\partial }_1} {\partial }_1\le v\le {\partial }_2\\ \pi +(1-\pi )\frac{v-{\partial }_2}{{\partial }_3-{\partial }_2} {\partial }_2\le v\le {\partial }_3\\ 1 v\ge {\partial }_3 \end{array}\right\}$$

(31)

Based on Equations (30) and (31), the Me-measure for the cases of $Me\left\{\tilde{\partial }\le v\right\}$ and $Me\left\{\tilde{\partial }\ge v\right\}$ can be obtained as below:

$$Me\left\{\tilde{\partial }\le v\right\}\ge \beta \iff \pi +(1-\pi )\frac{v-{\partial }_2}{{\partial }_3-{\partial }_2}\ge \beta \iff v\ge \frac{(\beta -\pi ){\partial }_3+(1-\pi ){\partial }_2}{1-\pi }$$

(32)

$$Me\left\{\tilde{\partial }\ge v\right\}\ge \beta \iff \pi +(1-\pi )\frac{{\partial }_2-v}{{\partial }_2-{\partial }_1}\ge \beta \iff v\le \frac{(\beta -\pi ){\partial }_1+(1-\pi ){\partial }_2}{1-\pi }$$

(33)

The mathematical model provided in Equation (24) is converted into an equivalent form by using the ExpV operator provided in Equation (28) and the Me-measure provided in Equations (32) and (33) as below:

$$\begin{array}{l}\mathrm{Minimize} ExpV\left[q\right]=\left(\frac{a^{pes}+2a^{most}+a^{opt}}{4}\right)m+\left(\frac{b^{pes}+2b^{most}+b^{opt}}{4}\right)n\\ \mathrm{Such} \mathrm{that} Cn\le \frac{({\beta }_1-\pi )r^{pes}+(1-\pi )r^{most}}{1-\pi }\\ En=0\\ Fn\le \frac{({\beta }_2-\pi )G^{pes}+(1-\pi )G^{most}}{1-\pi }m\\ Hn\le 1\\ m\in \{0,1\}, n\ge 0, 0.5\le {\beta }_1,{\beta }_2\le 1\end{array}$$

(34)

where β₁ and β₂ are the minimum levels of confidence for uncertain constraints.

4.2. Flexible Possibilistic Programming (FPP)

Herein, the PP formulation obtained in Equation (34) is further modified by integrating flexible programming. The modified form after FP integration is as follows:

$$\begin{gathered} \mathrm{Minimize} ExpV\left[q\right]=\left(\frac{a^{pes}+2a^{most}+a^{opt}}{4}\right)m+\left(\frac{b^{pes}+2b^{most}+b^{opt}}{4}\right)n \\ \begin{array}{l}\mathrm{Such} \mathrm{that} Cn\tilde{\le }\frac{({\beta }_1-\pi )r^{pes}+(1-\pi )r^{most}}{1-\pi }\\ En=0\\ Fn\tilde{\le }\frac{({\beta }_2-\pi )G^{pes}+(1-\pi )G^{most}}{1-\pi }m\\ Hn\le 1\\ m\in \{0,1\}, n\ge 0, 0.5\le {\beta }_1,{\beta }_2\le 1\end{array} \end{gathered}$$

(35)

In Equation (35), the symbol of $\tilde{\le }$ is added for constraints containing uncertain parameters.

$$\begin{array}{l}\mathrm{Minimize} ExpV\left[q\right]=\left(\frac{a^{pes}+2a^{most}+a^{opt}}{4}\right)m+\left(\frac{b^{pes}+2b^{most}+b^{opt}}{4}\right)n\\ \mathrm{Such} \mathrm{that} Cn\le \left[\frac{\left({\beta }_1-\pi \right)r^{pes}+\left(1-\pi \right)r^{most}}{1-\pi }\right]+\left[\frac{{\xi }^{pes}+2{\xi }^{most}+{\xi }^{opt}}{4}\left(1-{\epsilon }_1\right)\right]\\ En=0\\ Fn\le \left[\frac{\left({\beta }_2-\pi \right)G^{pes}+\left(1-\pi \right)G^{most}}{1-\pi }\right]m+\left[\frac{{\theta }^{pes}+2{\theta }^{most}+{\theta }^{pes}}{4}\left(1-{\epsilon }_2\right)\right]m\\ Hn\le 1\\ m\in \left\{0,1\right\},n\ge 0,0.5\le {\beta }_1,{\beta }_2\le 1,0\le {\epsilon }_1,{\epsilon }_2\le 1\end{array}$$

(36)

In Equation (36) $\left[\frac{{\xi }^{pes}+2{\xi }^{most}+{\xi }^{opt}}{4}(1-{\epsilon }_1)\right]$ and $\left[\frac{{\theta }^{pes}+2{\theta }^{most}+{\theta }^{pes}}{4}(1-{\epsilon }_2)\right]m$ are added in constraints containing $\tilde{\le }$ symbols, which makes the constraints more flexible. These terms define the extent of a possible violation of each constraint based on the preference of the SC manager. Here, ${\epsilon }_1 \mathrm{and} {\epsilon }_2$ are the confidence level (CL) of the manager for the possible violation of the constraint.

4.3. Flexible Robust Possibilistic Programming (FRPP)

The FPP formulation provided in Equation (36) can efficiently tackle the operational uncertainty which arises due to lack of information and also incorporates flexibility in the uncertain constraints. However, FPP cannot control the deviation over and above the ExpV of the objective function. This drawback can be effectively overcome by the incorporation of RPP. In this research, the RPP-II form proposed by Pishvaee et al. [37] is adopted. The equivalent FRPP formulation of Equation (36) is as below:

$$\begin{gathered} \begin{array}{l}\mathrm{Minimize} ExpV\left[q\right]+\mathrm{\Psi }(q_{...\max }-ExpV\left[q\right])+{\kappa }_1\left[\left(\frac{({\beta }_1-\pi )r^{pes}+(1-\pi )r^{most}}{1-\pi }\right)-r^{pes}\right] +\\ {\kappa }_2\left[\left(\frac{({\beta }_2-\pi )G^{pes}+(1-\pi )G^{most}}{1-\pi }\right)-G^{pes}\right] +{\wp }_1\left[\frac{{\xi }^{pes}+2{\xi }^{most}+{\xi }^{opt}}{4}(1-{\epsilon }_1)\right] + {\wp }_2\left[\frac{{\theta }^{pes}+2{\theta }^{most}+{\theta }^{pes}}{4}(1-{\epsilon }_2)\right]m\\ \mathrm{Such} \mathrm{that} ExpV\left[q\right]=\left(\frac{a^{pes}+2a^{most}+a^{opt}}{4}\right)m+\left(\frac{b^{pes}+2b^{most}+b^{opt}}{4}\right)n \\ q_{\max }=a^{most}m+b^{most}n \end{array} \\ \begin{array}{l} Cn\le \left[\frac{({\beta }_1-\pi )r^{pes}+(1-\pi )r^{most}}{1-\pi }\right]+\left[\frac{{\xi }^{pes}+2{\xi }^{most}+{\xi }^{opt}}{4}(1-{\epsilon }_1)\right]\\ En=0\\ Fn\le \left[\frac{({\beta }_2-\pi )G^{pes}+(1-\pi )G^{most}}{1-\pi }\right]m+\left[\frac{{\theta }^{pes}+2{\theta }^{most}+{\theta }^{pes}}{4}(1-{\epsilon }_2)\right]m\\ Hn\le 1\\ m\in \{0,1\}, n\ge 0, 0.5\le {\beta }_1,{\beta }_2\le 1, 0\le {\epsilon }_1,{\epsilon }_2\le 1 \end{array} \end{gathered}$$

(37)

where $\mathrm{\Psi }$ is the scaling factor of the second term that provides optimality robustness (OR), ${\kappa }_1$ and ${\kappa }_2$ are the penalties for violating the uncertain constraints due to operational uncertainty that enhance feasibility robustness (FR), and ${\wp }_1$ and ${\wp }_2$ are penalties for integrating the flexibility in constraint goals.

4.4. Transformation of the BG-SCND Model into an Equivalent Crisp Form Using FRPP Methodology

In this section, the FRPP methodology is implemented on the proposed BG-SCND optimization model through the following three phases:

4.4.1. Stage I: Conversation of the BG-SCND Optimization Model into a PP Formulation

In the first stage, the uncertain BG-SCND model is transformed into a PP formulation as follows:

$$\begin{array}{l}Minimize ExpV(f_{TC})={\displaystyle \sum _c^C\left(\frac{IC{C}_c^{pes}+2IC{C}_c^{most}+IC{C}_c^{opt}}{4}\right)}\times w_c+{\displaystyle \sum _b^B{\displaystyle \sum _a^A\left(\frac{IC{B}_{ba}^{pes}+2IC{B}_{ba}^{most}+IC{B}_{ba}^{opt}}{4}\right)\times y_{ba}}}+\\ {\displaystyle \sum _s^S{\displaystyle \sum _u^U\left(\frac{IC{S}_{su}^{pes}+2IC{S}_{su}^{most}+IC{S}_{su}^{opt}}{4}\right)\times z_{su}}}+{\displaystyle \sum _p^P{\displaystyle \sum _c^C{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[\left(\frac{CF{T}_{pft}^{pes}+2CF{T}_{pft}^{most}+CF{T}_{pft}^{opt}}{4}\right)\right]}}}}F{T}_{pcft}+\\ {\displaystyle \sum _p^P{\displaystyle \sum _c^C{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[\left(\frac{HC{C}_{cft}^{pes}+2HC{C}_{cft}^{most}+HC{C}_{cft}^{opt}}{4}\right)+\left(\left(\frac{CE{H}_c^{pes}+2CE{H}_c^{most}+CE{H}_c^{opt}}{4}\right)\times C_{tax}\right)\right]}}}}F{T}_{pcft}+\\ {\displaystyle \sum _b^B{\displaystyle \sum _t^T\left[\left(\frac{PC{B}_b^{pes}+2PC{B}_b^{most}+PC{B}_b^{opt}}{4}\right)+\left(\left(\frac{CE{P}_b^{pes}+2CE{P}_b^{most}+CE{P}_b^{opt}}{4}\right)\times C_{tax}\right)\right]}}AB{P}_{bt}+\\ {\displaystyle \sum _p^P{\displaystyle \sum _c^C{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[\left(\frac{TC{F}_{pc}^{pes}+2TC{F}_{pc}^{most}+TC{F}_{pc}^{opt}}{4}\right)+\left(\left(\frac{CE{T}_{pc}^{pes}+2CE{T}_{pc}^{most}+CE{T}_{pc}^{opt}}{4}\right)\times C_{tax}\right)\right]}}}}F{T}_{pcft}+\\ {\displaystyle \sum _c^C{\displaystyle \sum _b^B{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[\left(\frac{TC{F}_{cb}^{pes}+2TC{F}_{cb}^{most}+TC{F}_{cb}^{opt}}{4}\right)+\left(\left(\frac{CE{T}_{cb}^{pes}+2CE{T}_{cb}^{most}+CE{T}_{cb}^{opt}}{4}\right)\times C_{tax}\right)\right]}}}}F{T}_{cbft}+\\ {\displaystyle \sum _b^B{\displaystyle \sum _s^S{\displaystyle \sum _t^T\left[\left(\frac{TC{B}_{bs}^{pes}+2TC{B}_{bs}^{most}+TC{B}_{bs}^{opt}}{4}\right)+\left(\left(\frac{CE{T}_{bs}^{pes}+2CE{T}_{bs}^{most}+CE{T}_{bs}^{opt}}{4}\right)\times C_{tax}\right)\right]}}}B{T}_{bst}+\\ {\displaystyle \sum _s^S{\displaystyle \sum _m^M{\displaystyle \sum _t^T\left[\left(\frac{TC{B}_{sm}^{pes}+2TC{B}_{sm}^{most}+TC{B}_{sm}^{opt}}{4}\right)+\left(\left(\frac{CE{T}_{sm}^{pes}+2CE{T}_{sm}^{most}+CE{T}_{sm}^{opt}}{4}\right)\times C_{tax}\right)\right]}}}B{T}_{smt}+\\ {\displaystyle \sum _b^B{\displaystyle \sum _n^N{\displaystyle \sum _t^T\left[\left(\frac{TC{D}_{bn}^{pes}+2TC{D}_{bn}^{most}+TC{D}_{bn}^{opt}}{4}\right)+\left(\left(\frac{CE{T}_{bn}^{pes}+2CE{T}_{bn}^{most}+CE{T}_{bn}^{opt}}{4}\right)\times C_{tax}\right)\right]}}}D{T}_{bnt}\end{array}$$

(38)

$$Maximize ExpV(f_{soc})={\mathrm{\Phi }}^{EGI}EG{I}^{nrm}+{\mathrm{\Phi }}^{JPI}JP{I}^{nrm}$$

(39)

Economic growth index (EGI)

$$\begin{array}{l}EGI={\displaystyle \sum _c^C\left[\frac{In{v}_c^{pes}+2In{v}_c^{most}+In{v}_c^{opt}}{}\right]\times w_c(1-Dev{l}_c)}+{\displaystyle \sum _b^B{\displaystyle \sum _a^A\left[\frac{In{v}_{ba}^{pes}+2In{v}_{ba}^{most}+In{v}_{ba}^{opt}}{}\right]\times y_{ba}(1-Dev{l}_b)+}}\\ {\displaystyle \sum _s^S{\displaystyle \sum _u^U\left[\frac{In{v}_{su}^{pes}+2In{v}_{su}^{most}+In{v}_{su}^{opt}}{}\right]\times z_{su}(1-Dev{l}_s)}}\end{array}$$

(40)

$$\begin{array}{l}JPI={\displaystyle {\displaystyle \sum }_c^C}\left[\left\{\frac{Un{e}_c^{pes}+2Un{e}_c^{most}+Un{e}_c^{opt}}{4}\right\}\times \left\{\frac{Jo{b}_c^{pes}+2Jo{b}_c^{most}+Jo{b}_c^{opt}}{4}\right\}\right]\times w_c+\\ {\displaystyle {\displaystyle \sum }_b^B}{\displaystyle {\displaystyle \sum }_a^A}\left[\left\{\frac{Un{e}_b^{pes}+2Un{e}_b^{most}+Un{e}_b^{opt}}{4}\right\}\times \left\{\frac{Jo{b}_{ba}^{pes}+2Jo{b}_{ba}^{most}+Jo{b}_{ba}^{opt}}{4}\right\}\right]\times y_{ba}+\\ {\displaystyle {\displaystyle \sum }_s^S}{\displaystyle {\displaystyle \sum }_u^U}\left[\left\{\frac{Un{e}_s^{pes}+2Un{e}_s^{most}+Un{e}_s^{opt}}{4}\right\}\times \left\{\frac{Jo{b}_{su}^{pes}+2Jo{b}_{su}^{most}+Jo{b}_{su}^{opt}}{4}\right\}\right]\times z_{su}\end{array}$$

(41)

$${\displaystyle \sum _s^{S}B{T}_{smt}}\tilde{\ge }\frac{({\delta }_1-\phi )D_{mt}^{opt}+(1-\phi )D_{mt}^{most}}{(1-\phi )} \forall m,t|0.5\le {\delta }_1\le 1$$

(42)

$${\displaystyle \sum _b^{B}D{T}_{bnt}}\tilde{\ge } \frac{({\delta }_2-\phi )D_{nt}^{opt}+(1-\phi )D_{nt}^{most}}{(1-\phi )} \forall n,t|0.5\le {\delta }_2\le 1$$

(43)

$$\frac{({\delta }_3-\phi )Q{F}_{pft}^{pes}+(1-\phi )Q{F}_{pft}^{most}}{1-\phi }\tilde{\ge }{\displaystyle \sum _c^{C}F{T}_{pcft} \forall p,f,t}|0.5\le {\delta }_3\le 1$$

(44)

and Equations (10)–(23).

4.4.2. Stage II: Conversation of the BG-SCND Optimization Model into an FPP Formulation

The second phase of the suggested solution approach involves converting the PP formulation of the BG-SCND model into its equivalent FPP form as follows:

$${\displaystyle \sum _s^{S}B{T}_{smt}}\ge \left[\frac{({\delta }_1-\phi )D_{mt}^{opt}+(1-\phi )D_{mt}^{most}}{(1-\phi )}\right]-\left[\left(\frac{{\varpi }^{opt}+2{\varpi }^{most}+{\varpi }^{pes}}{4}\right)(1-{\vartheta }^{Bdm})\right] \forall m,t|0.5\le {\delta }_1\le 1$$

(45)

$${\displaystyle \sum _b^{B}D{T}_{bnt}}\ge \left[\frac{({\delta }_2-\phi )D_{nt}^{opt}+(1-\phi )D_{nt}^{most}}{(1-\phi )}\right]-\left[\left(\frac{{\mu }^{pes}+{\mu }^{most}+{\mu }^{opt}}{4}\right)(1-{\vartheta }^{Ddm})\right] \forall n,t|0.5\le {\delta }_2\le 1$$

(46)

$$\left[\frac{({\delta }_3-\phi )Q{F}_{pft}^{pes}+(1-\phi )Q{F}_{pft}^{most}}{1-\phi }\right]+\left[\left(\frac{{\sigma }^{opt}+2{\sigma }^{most}+{\sigma }^{pes}}{4}\right)\left(1-{\vartheta }^{b\sup }\right)\right]\ge {\displaystyle \sum _c^{C}F{T}_{pcft} \forall p,f,t} 0.5\le {\delta }_3\le 1$$

(47)

and Equations (10)–(23), (38) and (39).

4.4.3. Stage III: Conversation of the BG-SCND Optimization Model into an FRPP Formulation

In the third stage, the BG-SCND model’s FPP form is converted into FRPP form as follows:

$$\begin{array}{l}{F}_{TC}= ExpV(f_{TC})+ \mathrm{\Psi }\left\{f_{TC}{}_{\max }-ExpV(f_{TC})\right\}+{\kappa }_{TC}^{{}^{Bdm}}{\displaystyle \sum _m^M{\displaystyle \sum _t^T\left[\left\{D_{mt}^{opt}-\frac{({\delta }_1-\phi )D_{mt}^{opt}+(1-\phi )D_{mt}^{most}}{(1-\phi )}\right\}\right]}} +\\ {\kappa }_{TC}^{{}^{Ddm}}{\displaystyle \sum _n^N{\displaystyle \sum _t^T\left[\left\{D_{nt}^{opt}-\frac{({\delta }_2-\phi )D_{nt}^{opt}+(1-\phi )D_{nt}^{most}}{(1-\phi )}\right\}\right]}} +{\kappa }_{TC}^{b\sup }{\displaystyle \sum _p^P{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[\left\{\frac{({\delta }_3-\phi )Q{F}_{pft}^{pes}+(1-\phi )Q{F}_{pft}^{most}}{1-\phi }\right\}-Q{F}_{pft}^{pes}\right]}}}+\\ {\wp }_{TC}^{{}^{Bdm}}\left[\left(\frac{{\varpi }^{opt}+2{\varpi }^{most}+{\varpi }^{pes}}{4}\right)(1-{\vartheta }^{Bdm})\right]+{\wp }_{TC}^{{}^{Ddm}}\left[\left(\frac{{\mu }^{pes}+{\mu }^{most}+{\mu }^{opt}}{4}\right)(1-{\vartheta }^{Ddm})\right]+{\wp }_{TC}^{b\sup }\left[\left(\frac{{\sigma }^{opt}+2{\sigma }^{most}+{\sigma }^{pes}}{4}\right)\left(1-{\vartheta }^{b\sup }\right)\right]\end{array}$$

(48)

$$\begin{array}{l}{F}_{SOC}= ExpV(f_{SOC})+ \mathrm{\Psi }\left\{f_{SOC}{}_{\max }-ExpV(f_{SOC})\right\}+{\kappa }_{SOC}^{Bdm}{\displaystyle \sum _m^M{\displaystyle \sum _t^T\left[\left\{D_{mt}^{opt}-\frac{({\delta }_1-\phi )D_{mt}^{opt}+(1-\phi )D_{mt}^{most}}{(1-\phi )}\right\}\right]}} +\\ {\kappa }_{SOC}^{Ddm}{\displaystyle \sum _n^N{\displaystyle \sum _t^T\left[\left\{D_{nt}^{opt}-\frac{({\delta }_2-\phi )D_{nt}^{opt}+(1-\phi )D_{nt}^{most}}{(1-\phi )}\right\}\right]}} +{\kappa }_{SOC}^{b\sup }{\displaystyle \sum _p^P{\displaystyle \sum _f^F{\displaystyle \sum _t^T\left[\left\{\frac{({\delta }_3-\phi )Q{F}_{pft}^{pes}+(1-\phi )Q{F}_{pft}^{most}}{1-\phi }\right\}-Q{F}_{pft}^{pes}\right]}}}+\\ {\wp }_{SOC}^{Bdm}\left[\left(\frac{{\varpi }^{opt}+2{\varpi }^{most}+{\varpi }^{pes}}{4}\right)(1-{\vartheta }^{Bdm})\right]+{\wp }_{SOC}^{Ddm}\left[\left(\frac{{\mu }^{pes}+{\mu }^{most}+{\mu }^{opt}}{4}\right)(1-{\vartheta }^{Ddm})\right]+{\wp }_{SOC}^{b\sup }\left[\left(\frac{{\sigma }^{opt}+2{\sigma }^{most}+{\sigma }^{pes}}{4}\right)\left(1-{\vartheta }^{b\sup }\right)\right]\end{array}$$

(49)

and Equations (10)–(23) and (45)–(49).

In Equation (48), ${\kappa }_{TC}^{Bdm},{\kappa }_{TC}^{Ddm},{\kappa }_{TC}^{b\sup }$, and in Equation (49), ${\kappa }_{SOC}^{Bdm},{\kappa }_{SOC}^{Ddm},{\kappa }_{SOC}^{b\sup }$, are penalties for the violation of uncertain constraint targets for economic and social objectives, respectively. While ${\wp }^{Bdm},{\wp }^{Ddm},\mathrm{and}{\wp }^{b\sup }$ are the soft constraint violation penalties.

4.4.4. Stage IV: Obtaining Extreme Solutions for Each Objective and Conversion of the Multi-Objective Model into a Single Objective

After obtaining the FRPP form for each objective, the next step is to acquire extreme solutions: positive extreme solutions (PES) and negative extreme solutions (NES) for each objective of the optimization model. By using these extreme solutions, membership functions for both objectives are developed as follows:

$${\mu }_{F_{TC}}=\left\{\begin{array}{l}1,\\ \frac{F_{TC}^{NES}-F_{TC}}{F_{TC}^{NES}-F_{TC}^{PES}}\\ 0,\end{array}\right. ,\begin{array}{c} \mathrm{If}{F}_{TC}\le F_{TC}^{PES}\\ \mathrm{If}{F}_{TC}^{PES}<F_{TC}<F_{TC}^{NES}\\ \mathrm{If}{F}_{TC}\ge F_{TC}^{NES}\end{array}$$

(50)

$${\mu }_{F_{SOC}}=\left\{\begin{array}{l}0,\\ \frac{F_{SOC}-F_{SOC}^{NES}}{F_{SOC}^{PBS}-F_{SOC}^{NBS}}\\ 1,\end{array}\right. , \begin{array}{c} If{F}_{SOC}\le F_{SOC}^{NES}\\ \mathrm{If}{F}_{SOC}^{NES}<F_{SOC}<F_{SOC}^{PBS}\\ \mathrm{If}{F}_{SOC}\ge F_{PBS}^{soc}\end{array}$$

(51)

Using the membership functions provided in Equations (50) and (51), the multi-objective BG-SCND model is merged into a single objective using the Torabi and Hassini (TH) approach (for details of the TH method see [44,45]).

$$\begin{array}{l}\mathrm{Maximize} \varphi {\wp }^{comp}+\{(1-\varphi )({\nu }_{TC}\times {\mu }_{F_{TC}}+{\nu }_{SOC}\times {\mu }_{F_{SOC}})\}\\ {\mu }_{F_{TC}}\ge {\wp }^{comp}\\ {\mu }_{F_{SOC}}\ge {\wp }^{comp}\\ {\wp }^{comp}\in \left[0,1\right]\end{array}$$

(52)

Equations (10)–(23) and (45)–(49).

In Equation (52), $\varphi$ is the compensation coefficient, ${\wp }^{comp}$ is the mandatory lowest satisfaction, and ${\mu }_{F_{TC}}$ and ${\mu }_{F_{SOC}}$ are satisfaction levels for BG-SCND model objectives.

Figure 4 provides the systematic procedure of the employed solution methodology.

5. Case Study

To validate the proposed feedstock-based supply chain model, the largest province by population in Pakistan, Punjab, is considered. Biomass is collected from multiple districts of Punjab and sent to a supply center, where it is transported to biomass collection centers established in major cities. After that, biomass is supplied to production plants to obtain biogas. Biogas from these plants is delivered to consumer market zones via storage and distribution channels, while leftover digestate is delivered directly to digestate demand zones.

Facility site decisions are integrated as a binary decision variable in the BG-DSCND model, since transportation costs are a significant portion of the overall biogas SC costs and are, therefore, necessary to optimize. The provided model chooses the operational biogas processing facility’s ideal capacity, besides deciding the location of the facility. The shortlisted locations in the BG-SCND model are selected using a multi-criteria evaluation approach that satisfies all requirements. In this case, there are 26 feedstock supply points, eight locations of preprocessing plants, five locations of production plant facilities, and four market zone hubs of biogas. Due to easy accessibility, dung, bagasse, rice husk, and organic municipal waste are used as feedstock. The potential biogas-processing facilities are depicted in Figure 5 on a district-based geographical map.

5.1. Data Acquisition for the BG-SCND Model

Herein, the required data for the BG-SCND model are provided. The demand for biogas in each demand zone was calculated based on methane usage during the previous three years. Statistics for the last three years were obtained from the ministry of energy (Power division) (http://www.mowp.gov.pk/ (accessed on 5 January 2021)). Based on the data for the last three years, 0.2% of methane demand was considered as demand for biogas for each market zone. Data for municipal waste for different districts of Punjab were obtained from the local government and community development department (https://lgcd.punjab.gov.pk/ (accessed on 10 January 2021)). Input parameters for the social objective were obtained from the Pakistan bureau of the statistics department. (https://www.pbs.gov.pk/publication/pakistan-social-standards-measurement-survey-2019-20-provincial (accessed on 20 January 2021)). Data for bagasse and rice husk were obtained from the agriculture marketing information service department (http://www.amis.pk/ (accessed on 20 January 2021)), while data for animal dung were obtained from the livestock and dairy department (https://lddb.org.pk/ (accessed on 1 December 2020)). The statistics on biogas plant installation, biomass procurement, biomass material handling, and production expenses were obtained from biogas production firms, published studies, and government agency reports. CO₂ emissions during biogas production and transportation among processing facilities were adopted from published articles by Gonela et al. [12] and Gürü et al. [46]. An emission tax of $10/ton of CO₂ was assumed based on the report published by the World Bank [47]. Since the adopted values of various parameters are tainted with epistemic uncertainty, these parameters are considered uncertain in the BG-SCND model. The main strategic and tactical-level decisions made in the BG-SCND model to attain a trade-off among the sustainability dimensions are as follows:

• Number of biomass collection centers to be made operational against specified sustainability dimension preferences.
• Number and capacity level of the biogas production plant and storage and distribution center to be made operational.
• The quantity of biomass to be allocated from operational biomass supply centers to biomass collection centers and from operational biomass collection centers to biogas production plants.
• The quantity of biogas to be supplied from the production plant to the storage and distribution center and from biogas distribution centers to biogas market zones.
• The quantity of digestate to be allocated from the biogas production plant to digestate demand zones.

Appendix A contains the most probable practical estimates for the BG-SCND model parameters.

5.2. Validation of the BG-SCND Optimization Model Results and Computational Analysis for a Small-Scale Case

In the next stage, collected data and the equivalent FRPP form of the BG-SCND model provided in Equations (10)–(23), (45)–(49), and (52) were coded in LINGO software and the results were obtained. Since the provided FRPP approach is interactive, the managers must first pick the values of the following parameters before obtaining results: ${\delta }_1,{\delta }_2,{\delta }_3$ (CL of manager), $\mathrm{\Psi }$(scaling parameter), ${\kappa }_{TC}^{Bdm},{\kappa }_{TC}^{Ddm},{\kappa }_{TC}^{b\sup },{\kappa }_{SOC}^{Bdm},{\kappa }_{SOC}^{Ddm},{\kappa }_{SOC}^{b\sup }$ (uncertain constraint violation penalty), and ${\wp }^{Bdm},{\wp }^{Ddm},{\wp }^{b\sup }$ (flexible constraint violation penalty). The value of the scaling parameter $(\mathrm{\Psi })$ provides robustness in terms of optimality and reduces the maximum deviation of the ExpV of the BG-SCND model objective. Further, the values of uncertain constraint violation and soft constraint violation penalty values $({\kappa }_{TC}^{Bdm},{\kappa }_{TC}^{Ddm},{\kappa }_{TC}^{b\sup },{\kappa }_{SOC}^{Bdm},{\kappa }_{SOC}^{Ddm},{\kappa }_{SOC}^{b\sup },{\wp }^{Bdm},{\wp }^{Ddm}, \mathrm{and} {\wp }^{b\sup })$ are decided based on the CL of the manager. For a low CL value, higher penalty values are decided, while for a high CL value, lower penalty levels are decided. It is also found that the values ${\kappa }_{TC}^{Bdm},{\kappa }_{TC}^{Ddm},{\kappa }_{TC}^{b\sup },{\kappa }_{SOC}^{Bdm},{\kappa }_{SOC}^{Ddm},{\kappa }_{SOC}^{b\sup }$ are highly sensitive. A little increase in these values has a large influence on the entire cost of the considered biogas SC. The RPP-II formulation proposed by Pishvaee et al. [37] is used in this work; however, if the CL value is very low, the RPP-II formulation is not appropriate due to its sensitivity. Instead, the HWRPP formulation is suggested (for details see Pishvaee et al. [37]).

A small-scale scenario is initially presented to thoroughly investigate the decisions made by the BG-SCND model. For this illustration, eight biomass supply facilities, four potential locations for collection and preprocessing facilities, three potential locations for biogas plants, two potential locations for storage facilities, three locations for biogas markets, and two locations for digestate sites are considered. Considering a CL of 80% for the biomass supply and biogas demand constraints $({\delta }_1,{\delta }_2,{\delta }_3)$, the PES and NES are obtained, which are illustrated in

Table 2. The PES and NES for the BG-SCND model objectives.

BG-SCND Model Objectives	Economic Aspect (USD)	Social Aspect (%)
Minimize FTC	(PES) 33,518,210	(NES) 0.050
Maximize FSOC	(NES) 38,710,140	(PES) 1.000

δ₁ = 0.8, δ₂ = 0.8, δ₃ = 0.8.

. The results show that both objectives conflict with one another. For example, at the lowest value of the economic objective (33,518,210), which is the PES for the economic objective, the NES of the social objective (0.050) is attained, and vice versa.

Following that, using Equations (50) and (51), the membership functions for economic and social objectives are established as follows:

$${\mu }_{F_{TC}}=\left\{\begin{array}{l}1, if F_{TC}\le 33,518,210\\ \frac{38,710,140-F_{TC}}{38,710,140-33,518,210} , if 33,518,210 < F_{TC}<38,710,140\\ 0, if F_{TC}\ge 38,710,140\end{array}\right.$$

(53)

$${\mu }_{F_{SOC}}=\left\{\begin{array}{l}0, if F_{SOC}\le 0.050\\ \frac{F^{soc}-0.050}{1.00-0.050} , if 0.050 < F_{SOC}<1.00\\ 1, if F_{SOC}\ge 1.00\end{array}\right.$$

(54)

Using the membership function provided in Equations (53) and (54) with objective priority weights of ${\nu }_{TC}=0.8 \mathrm{and} {\nu }_{SOC}=0.2$, the values of ${\mu }_{F_{TC}}=0.98$(33,535,960) and ${\mu }_{F_{SOC}}=0.20$(0.2487) for the economic and social objectives are obtained, respectively. Figure 6 and Figure 7 show the BG-SCND model decisions derived at these parameters during the first and second planning periods, correspondingly. According to the results, it is observed that all of the four collection and preprocessing centers are made operational. In the following tier of biogas plants (ba), two of three potential biogas plant sites are made operational: Nankana and Rahimyar Khan. The optimization model selects a lower capacity level (B1-A1) for the Nankana biogas plant and a higher capacity level (B2-A2) for the Rahimyar Khan biogas plant. Both prospective sites are made operational at higher capacity levels in the next tier of biogas distribution centers (su). Overall, it is noticed that, because economic objectives are given more importance over social objectives, the BG-SCND model successfully utilizes economies of scale.

5.3. Results Analysis for the Large-Scale BG-SCND Optimization Model Case

This section discusses the outcomes of the BG-SCND optimization model’s large-scale scenario. As discussed in Section 4, this research utilizes the TH approach to convert the multi-objective BG-SCND model into a single objective. The TH approach is interactive and allows managers to obtain BG-SCND model decisions based on their preferences by setting the priority weights of the objectives $({\nu }_{TC}, {\nu }_{SOC})$. Herein, the results of the BG-SCND model are analyzed for two extreme cases of ${\nu }_{TC}=1.0, {\nu }_{SOC}=0.0$ and ${\nu }_{TC}=0.0, {\nu }_{SOC}=1.0$. For all uncertain constraints, a confidence level of 1 and δ₁ = 0.8, δ₂ = 0.8, and δ₃ = 0.8 are assumed in order to make the findings of both scenarios comparable.

Results of the first case only consider the economic aspect and depict that, among 26 biomass supply points, 18 are supplying biomass to preprocessing and collection centers. In the following echelon, all eight prospective sites for a biomass collection facility have been made operational. According to detailed results, the Gujrat facility covers the eastern sections of Punjab, while the Sheikhupura and Faisalabad biomass collection facilities serve the middle areas of Punjab. Multan and Okara, on the other hand, cover the central and southern portions of Punjab, while Mianwali and DG Khan biomass collecting facilities serve the western regions. Since no binary decision is associated with the biomass supply center tier, all of the supply centers are made operational and feedstock is supplied to the seven selected collection centers. Further, it has been noted that biomass collection facilities are strategically positioned close to each location’s greatest potential biomass supply to reduce total SC costs and carbon emissions during transportation. In the next echelon of the BG-SCND model, three out of five biogas production plants are made operational in Rahimyar Khan, Jhang, and Sargodha. Among the operational biogas production sites, the Rahimyar Khan site is made operational with a capacity of 20,000 m³/period, while Jhang and Jehlum biogas plants are made operational with a biomass processing capacity of 150,000 m³/period. The biomass for the Rahimyar Khan biogas plant primarily comes from the DG khan and Multan biomass collection centers, the biomass for the Jhang biogas plant arrives from the Faisalabad and Mianwali collection centers, and most of the biomass for the Jehlum biogas plant originates from the Gujrat and Sheikhupura facilities. In the next tier, biogas produced at the Rahimyar Khan plant is largely supplied to the Sahiwal distribution center, the Jhang plant supplies biogas to the Jhang and Khushab distribution centers, while the Jehlum plant supplies its biogas to Gujranwala and Khushab distribution centers.

Finally, in the last tier of the biogas SC, the Gujranwala distribution center meets the demand of the Lahore zone, while Rawalpindi receives biogas produced from the Khushab plant. The Sahiwal facility largely feeds biogas to the Faisalabad region, and the Jhang distribution center meets the demand for biogas in the Multan zone. Strategic decisions for the case ${\nu }_{TC}=1.0, {\nu }_{SOC}=0.0$ are provided in Figure 8.

The second scenario favors the social component over the economic side; hence, priority weights ${\nu }_{TC}=0.0, {\nu }_{SOC}=1.0$ are specified in the TH method. The obtained findings show that the results of this case are completely different from the first case in terms of location and capacity decisions. The results show that all of the preprocessing and collection facilities, biogas plants, and distribution sites are now operating. Furthermore, greater capacity levels are chosen for all of these location decisions, resulting in a significant rise in the total cost of the BG-SCND model. This is due to the fact that the social objective aims at maximizing the number of job opportunities and economic investments in regions with a low level of development that have higher unemployment rates. Detailed decisions for this case are provided in Figure 9.

5.4. Impact of the BG-SCND Model’s Objective Priority Weights on Strategic-Level Decisions

This section provides a comprehensive computational investigation of the effect of varying the priority weights of the BG-SCND model objectives on the model’s strategic-level decisions. For this purpose, various combinations of priority weights between 0 and 1 are developed for the economic and social objectives. As discussed earlier, the absolute priority of the economic objective over the social objective provides decisions based on the total cost aspect. In this scenario, the aim of minimizing costs is achieved by balancing economies of scale and overall transportation costs. Due to this, the capacity decisions made for this scenario choose a combination of higher and lower-level capacities for the biogas processing facilities. On the other hand, when the priority increasingly shifts to the social dimension, the degree of satisfaction with the economic objective begins to deteriorate. This is because, when the contribution of the social aspect is incorporated into the model, it attains the goal of social objective maximization by selecting numerous biogas facilities in regions that have low development levels and high unemployment rates. BG-SCND model strategic-level decisions involving location and capacity decisions against each combination of objective priority weight are provided in

Table 3. Impact of the BG-SCND model’s objective priority weights on strategic-level decisions.

Objective Weight		Objective Satisfaction Level		Operational Supply Terminals (p)	Operational Collection Center (c)	Operational Biogas Plant (ba)		Operational Distribution Center (su)
Economic	Social	Economic	Social			Location (b)	Capacity (a) (m3/period)	Location (s)	Capacity (u) (m3/period)
0	1	0.00%	100%	P2-P5, P7-P17, P20, P22, P24, P25	A1–A8	Nankana	150,000	Gujranwala	170,000
						Rahimyar Khan	150,000	Khushab	170,000
						Jhang	150,000	Sahiwal	170,000
						Sargodha	150,000	Muzaffargarh	170,000
						Jehlum	150,000	Jhang	170,000
0.2	0.8	3.99%	100%	P2, P4, P5, P7-P17, P20, P22	A1–A8	Rahimyar Khan	150,000	Gujranwala	170,000
						Jhang	150,000	Khushab	170,000
						Sargodha	150,000	Sahiwal	170,000
						Jehlum	150,000	Muzaffargarh	170,000
						Nankana	150,000	Jhang	170,000
0.4	0.6	75.38%	56.25%	P2-P5, P7-P18, P26	A1–A8	Rahimyar Khan	150,000	Gujranwala	170,000
						Jhang	150,000	Khushab	170,000
								Sahiwal	170,000
						Jehlum	150,000	Jhang	170,000
0.5	0.5	93.84%	42.12%	P2-P5, P7-P18, P26	A1–A8	Rahimyar Khan	150,000	Gujranwala	170,000
						Jhang	100,000	Khushab	170,000
								Sahiwal	170,000
						Jehlum	100,000	Jhang	170,000
0.6	0.4	93.84%	42.12%	P2-P5, P7-P18, P26	A1–A8	Rahimyar Khan	150,000	Gujranwala	170,000
						Jhang	100,000	Khushab	170,000
								Sahiwal	170,000
						Jehlum	100,000	Jhang	170,000
0.8	0.2	97.43%	35.90%	P2-P5, P7-P18, P26	A1–A8	Rahimyar Khan	150,000	Gujranwala	170,000
						Jhang	100,000	Khushab	170,000
								Sahiwal	170,000
						Jehlum	100,000	Jhang	120,000
1	0	100%	20.16%	P2-P5, P7-P18, P25-P26	A1–A8	Rahimyar Khan	150,000	Gujranwala	120,000
						Jhang	100,000	Khushab	170,000
								Sahiwal	170,000
						Jehlum	100,000	Jhang	120,000

.

5.5. The Effect of the Compensation Coefficient on the Objectives of the BG-SCND Model

To demonstrate the effect of the compensation coefficient (${\wp }^{comp}$) on the BG-SCND model’s objectives’ satisfaction levels ${\mu }_{F_{TC}}$ and ${\mu }_{F_{SOC}}$, a detailed analysis is provided by altering the ${\wp }^{comp}, {\nu }_{TC}, \mathrm{and} {\nu }_{SOC}$ values between 0 to 1. As discussed above, ${\wp }^{comp}$ is an interactive parameter and decided by the manager based on the real-time environment; hence, multiple efficient solutions can be obtained. The impact of variation on the compensation coefficient and priority weights on the BG-SCND model’s objectives’ attainment levels is graphically shown in Figure 10a–e. In the mentioned figure, the BG-SCND model’s objective satisfaction levels are drawn along the vertical axis, while objective preference weights are drawn along the horizontal axis. From the analysis of the results, it is observed that a balanced solution for the model’s objectives develops as the value of the compensation coefficient grows. This is because the model’s purpose for greater levels of ${\wp }^{comp}$ in the TH method is to augment the minimal level of satisfaction for all model objectives. As a result, in Figure 10e, satisfaction levels have the least amount of dispersion, whereas, in Figure 10a, goal satisfaction levels have the most dispersion.

6. Conclusions, Limitations, and Future Research Directions

Rising energy needs, as well as the rapid depletion of natural resources, have increased the desire for sustainable energy supplies. Against this background, this research suggests a BG-SCND model for biogas that uses multiple feedstocks (bagasse, organic MSW, dung, and rice husk). The following echelons are taken into account in the proposed model: feedstock supply centers, feedstock collection centers, biogas production facilities, biogas market zones, and digestate market zones. In addition to maximizing the social impact of biogas SC decisions, the presented model is concerned with minimizing the overall cost and environmental impact. To deal with the various types of epistemic uncertainties involved in the BG-SCND model, the FRPP approach is proposed, which not only efficiently deals with epistemic uncertainty, but also incorporates robustness in the obtained results and allows managers to make decisions based on their preferences. The major findings of the research are the following:

• Results of the proposed model demonstrate that, with a minor increase in the total cost, the strategic-level decisions can be made secure.
• It is also found that a fair inclusion of the social aspect in the BG-SCND model’s decisions may be accomplished at a marginal cost, but the absolute achievement of the social objective exponentially raises the overall cost of the proposed biogas SC.
• The sensitivity analysis was performed for various interactive parameters, such as ${\wp }^{comp}, {\nu }_{TC}, \mathrm{and} {\nu }_{SOC}$, depicting that the strategic-level decisions provided by the FRPP approach are robust.
• Greater economic preferences result in fewer facilities being operational, leading to the development of a centralized network, whereas greater social objective preference result in a large number of facilities being operational, creating a decentralized network configuration.

This research can be further extended in several directions. For example, the suggested BG-SCND model’s scope can be expanded by including the notion of SC coordination and contract, in which supply centers transport biomass to biogas production facilities and acquire biogas at a reasonable cost. This does not only encourage biogas supply chain partners to be more collaborative, but also helps to commercialize biogas. This study exclusively uses anaerobic digestion as a method of producing biogas; biogas production technology is not considered a decision. The research will be more useful if production technology is taken into account as a decision variable. However, increasing the scope of the optimization model will increase the overall number of decision variables. In such cases, the suggested solution approaches will take longer to provide global optimal decisions. To cope with this, intelligent and innovative algorithms and heuristics will be really valuable. Furthermore, this study employs an RPP-based approach that can effectively tackle operational uncertainty; however, incorporating disruption risk in the planning phase will further enhance the effectiveness of this approach.