The Waste-to-Biomethane Logistic Problem: A Mathematical Optimization Approach

In this paper, we propose a new mathematical optimization approach to make decisions on the optimal design of the complex logistic system required to produce biogas from waste. We provide a novel and flexible decision-aid tool that allows decision makers to optimally determine the locations of different types of plants (pretreatment, anaerobic digestion, and biomethane liquefaction plants) and pipelines involved in the logistic process, according to a given budget, as well as the most efficient distribution of the products (from waste to biomethane) along the supply chain. The method is based on a mathematical optimization model that we further analyze and that, after reducing the number of variables and constraints without affecting the solutions, is able to solve real-size instances in reasonable CPU times. The proposed methodology is designed to be versatile and adaptable to different situations that arise in the transformation of waste to biogas. The results of our computational experiments, both in synthetic and in a case study instance, prove the validity of our proposal in practical applications. Synthetic instances with up to 200 farms and potential locations for pretreatment plants and 100 potential locations for anaerobic digestion and biomethane liquefaction plants were solved, exactly, within <20 min, whereas the larger instances with 500 farms were solved within <2 h. The CPU times required to solve the real-world instance range from 2 min to 6 h, being highly affected by the given budget to install the plants and the percent of biomethane that is required to be injected in the existing gas network.


Introduction
Decarbonization of energy has emerged as a critical global priority in recent years due to the urgent need to combat climate change.Most World Organizations have recognized the significance of transitioning to cleaner and more sustainable energy sources in the next few years to reduce greenhouse gas emissions.Different agreements have been signed to promote this change.Specifically, the United Nations, through its Sustainable Development Goals, has set different targets to ensure affordable, reliable, sustainable, and modern energy for all by 2030 (United Nations, 2015).The International Energy Agency (IEA) has been actively promoting decarbonization efforts by providing policy recommendations, conducting research, and facilitating international cooperation (International Energy Agency, 2021).Additionally, the Intergovernmental Panel on Climate Change (IPCC) (Legg, 2021) has been instrumental in assessing the impacts of climate change and highlighting the importance of decarbonizing the energy sector.In 2019, the European Commission presented the European Green Deal (Fetting, 2020), an initiative that aimed to be a statement of intent on the road-map that Europe should follow for the implementation of the United Nations' Sustainable Development Agendas for 2030 and 2050, designed to mitigate the effects of climate change.Various measures have been taken by these organizations, including advocating for renewable energy investments, promoting energy efficiency, encouraging the use of electric vehicles, and supporting the development of innovative technologies such as carbon capture and storage.These collective efforts by world organizations are crucial in driving the global transition towards a low-carbon future and mitigating the adverse effects of climate change.
One of challenges is to decarbonize the energy system by developing a new power sector based on renewable sources.It is a known fact that one of the main strategies to achieve the proposed goal is the use of biogas as an alternative renewable energy source to carbonbased energies, since it contributes the reduction of greenhouse gases but also to the development of the circular economy through the anaerobic digestion of organic waste from different sources and its transformation into fuel.Since biomethane is the same molecule as natural gas, it can be distributed via the existing gas distribution networks, facilitating the transition from natural gas to biogas energy.Thus, the installation of anaerobic digestion plants for the conversion of organic waste and livestock manure into biomethane has gained prominence in recent years as a sustainable waste-to-energy solution and many countries are making special interest in the assessment of the installation and adaptation of these biogas plants.For instance, the U.S. Farm Bill (National Sustainable Agriculture Coalition, 2023), periodically reauthorized by Congress, includes provisions that support the development and utilization of biogas from agricultural sources.The bill provides funding for biogas projects, such as the installation of anaerobic digesters plants on farms or near them, which capture methane emissions from manure and convert them into usable biogas for electricity generation or transportation fuel.
Analyzing the logistic systems behind the implementation of different modes of energy, particularly biogas, is of paramount importance in ensuring the successful and sustainable integration of renewable energy sources into our global energy landscape.Biogas, derived from organic waste materials through anaerobic digestion, presents a promising alternative to conventional fossil fuels.However, its widespread adoption hinges on addressing intricate logistic challenges.Understanding the logistics involved in the collection, transportation, and processing of organic feedstocks is crucial for optimizing efficiency and minimizing environmental impact.Rigorous analyses on the logistical aspects of biogas production in different countries reveal opportunities for streamlining supply chains, enhancing resource utilization, and reducing overall costs (see, e.g.Muradin and Foltynowicz, 2018;Kwaśny and Balcerzak, 2017).Effective logistic planning also ensures reliable and consistent feedstock availability, a key factor in the stable operation of biogas plants.By delving into the logistics of biogas implementation, we can develop strategies that not only bolster the economic viability of this renewable energy source but also contribute to the broader goal of achieving a more sustainable and resilient energy infrastructure.
Among the decisions to make in the design of an efficient logistic systems, one of the most critical ones is selecting the optimal sites to locate the different types of plants involved in the system (biogas plants, pre-treatment plants, and liquefaction plants).This is a difficult task since it involves different agents, different production and conversion technologies, different types of demand centers, and it has to be coordinated with the different distribution processes (Jafaru et al., 2018).At this point, mathematical optimization plays a very important role since it is the main tool to determine solutions to complex logistics system like the one we are dealing with here.This fact has caused the publication of a large amount of literature related to this topic.For instance, Mathematical optimization has been used in Tampio et al. (2022) for the design of a cost-optimal processing route for a biogas anaerobic digestion plant to produce fertilizer products based on specified regional needs.In Balaman and Selim (2014) a mixed integer linear programming model is developed to determine the appropriate locations for the biogas plants and biomass storages.Scarlat et al. (2018) provides a spatial analysis algorithm that uses data of manure production and collection to assess the spatial distribution of the biogas potential in Europe, in order to decide the location of bioenergy plants.A geographic information system (GIS) based analysis is used in Valenti et al. (2018) to determine the size and location of four biogas plants for the Catania province in Sicily.A multi-objective optimization approach for the optimal location of biogas and biofertilizer plants at the maximum profit and the minimum environmental impact is presented in Díaz-Trujillo and Nápoles-Rivera (2019) and it is applied for a geographical region in Mexico.Mathematical methods have been used to study the location of bioenergy plants in many other regions and countries (Park et al. (2019) in North Dakota, Sultana and Kumar (2012) in the province of Alberta, Egieya et al. (2020) in Slovenia, Amigun and von Blottnitz (2010) in Africa, Silva et al. (2017) in Portugal, Soha et al. (2021) in Hungarian, Igliński et al. (2012) in Poland, Delzeit and Kellner (2013) in Germany).Furthermore, mathematical optimization tools have been already recognized as a crucial tool to make decisions in the design of systems for biogas production.Hernandez and Martin (2017), propose a non linear optimization model to compute the optimal operating conditions of the reactors, the biogas composition and the content of nutrients in the digestate in the production of biodiesel from waste via anaerobic digestion.Hu et al. (2018) give a mathematical optimization framework to analyze the economical and environmental benefits of recovering biogas, caproic acid (C6), and caprylic acid (C8) from different sources of organic waste.Ankathi et al. (2021) consider a mixed integer linear optimization model to determine the location and capacities of biogas plants based only in the location and production of the farms.In terms of locational decisions, some of the classical mathematical approaches have been already adapted to incorporate green goals (see e.g.Martínez and Fransoo, 2017), although more advances are expected within the next few year based on what our society needs.
On the other hand, mathematical optimization is recognized as a fundamental tool for the design and modeling of multiple logistics, transportation, and supply chain problems (see e.g.Hinojosa et al., 2000Hinojosa et al., , 2008;;Blanco et al., 2010Blanco et al., , 2020Blanco et al., , 2022Blanco et al., , 2023;;Blanco and Martínez-Antón, 2024, among many others).Specifically, is particularly useful to construct robust and effective supply chains to integrate bioenergy in any economy (see e.g.De Meyer et al., 2015;Sarker et al., 2018;Tominac et al., 2022;Cambero and Sowlati, 2014;Ghaderi et al., 2016, and the references therein) .In particular, in order to implement an efficient and sustainable biogas distribution system, it is necessary to properly design a robust logistic plan for all the elements involved in the process, such as manure, biomethane, waste, etc.In this phase, it is crucial to decide not only where to locate the anaerobic digestion plants but also, where to locate pre-treatment plants for cleaning and drying the waste or transshipment plants to distribute the products, how to collect the manures from the farms or fields and sending them to the plants, how to link (if needed) the different types of plants, how to distribute the final biomethane to the gas distribution network or to external clients, how to dispatch possible fertilizers back to some of the farms, etc.A few works have already proposed models to make decision in this setting.Jensen et al. (2017) propose minimum cost flow based model for finding the optimal production and investment plan for a biogas supply chain by means of minimizing the transportation cost on an existing supply chain network.Three layers of analysis for designing optimal animal waste supply for anaerobic bio-digestion are detailed in Mayerle and Neiva de Figueiredo ( 2016): (1) Identification of the optimal locations for the anaerobic digestion plants the farms' information; (2) specification of the optimal distribution system; and (3) an operational layer that includes scheduling optimal biomass collection from each farm to minimize biogas loss.Sarker et al. (2019) studies the optimization of the supply chain cost for a biomethane gas production system which is organized in four stages (collecting feedstock to hubs located according to zip code areas, transporting feedstock from hubs to reactor(s), transporting biomethane gas from reactor(s) to condenser(s), and shipping the liquefied biomethane gas from condensers to final demand points).The problem is formulated as a mixed-integer non-linear mathematical optimization problem and, due to the complexity of the non-lineal model, a genetic heuristic algorithm is propose to solve it.
In this paper, we introduce the Waste-to-Biomethane Logistic Problem and provide a general and flexible mathematical optimization model to make decisions of locating biogas, pre-treatment, and liquefaction plants and distributing the different types of elements along a complex network that integrates all the aspects mentioned above.Thus, we consider a supply chain for the biomethane gas production combining six stages: (A) collecting waste and sending it to pre-treatment plants (independently of the zip code); (B) delivering contaminant-free organic waste from pre-treatment plants to anaeorobic digestion plants; (C) constructing pipelines for transporting biomethane from anaerobic digestion plants to injection points of an existing gas pipeline network; (D) constructing pipelines for transporting biomethane from anaerobic digestion plants to biomethane liquefaction plants; (E) dispatching fertilizer from anaerobic digestion plants to waste sources and; (F) shipping liquefied biomethane from biomethane liquefaction plants to customer demand points.We jointly integrate all these complex stages into a mixed integer linear programming (MILP) model seeking to minimize the overall transportation cost of the system by assumming that a limited budget is given to install all the plants and pipelines.Our MILP model can be efficiently solved using off-the-shelf optimization solvers (as Gurobi, CPLEX, or FICO) after proving some theoretical results for reducing the number of variables.This model is applied and tested in a real-word dataset based on the region of upper Yahara Watershed in the state of Wisconsin (see Sampat et al., 2019Sampat et al., , 2021, , for further information about this dataset).
The remainder of this paper is organized as follows.In Section 2 we introduce the Waste-to-Biomethane Logistic Problem (W2BLP) and present our modeling assumptions.Section 3 is devoted to detail the mathematical optimization model that we develop for the problem.We also prove in this section a theoretical result that allows us to significantly reduce the number of variables in the model.In Section 4, we report the results of our computational experiments on realistic synthetic instances.A case study is presented in Section 5 where the model is applied to the real-world dataset based on the region of upper Yahara Watershed.The paper closes in Section 6 with some conclusions and future research.

The Waste-to-Biomethane Logistic Problem
The problem under analysis consists on efficiently use the agricultural waste, as livestock manure, energy crops, municipal waste, etc, to be transformed into biomethane.This transformation requires different phases.First, from a given set of waste sources (WS), as farms or residual storages, each of them producing an amount of waste, it is transported into pre-treatment (PT) plants where the non-organic material is removed, and what remains is dried, pressed and prepared.The location of the PT plants is to be determined among a finite set of potential positions.
The dried contaminant-free organic waste (DOW) obtained after the pre-treatment, is delivered to an anaerobic digestion (AD) plant, where it is transformed into biomethane.The position of the AD plants is also to be decided among a finite set of potential locations.Biomethane is then either directly injected into an existing gas pipeline network (GPN) or processed into liquified natural gas (LNG).We assume that we are given a finite set of injection points in the GPN and that pipelines connecting the AD plants with a selected set of injections point are to be built.A minimum percent of the produced biomethane is to be injected in the GPN, and the remainder biomethane is transformed into LNG and distributed to a given set of extra customers (EC), each of them with a given demand.This minimum percent represents the trade-off between the gas injection to the network with respect to the amount of gas served to the external customers.The liquefaction of biomethane requires the construction of Biomethane Liquefaction (BL) plants as well as pipelines connecting the AD plants with them.Then, the liquified gas is delivered to the customers using tank trucks.Finally, the remaining material in the DOW to biomethane transformation process (digester solids) is then returned back from the AD plants to some of the WS where it can be used as fertilizer.In Figure 1 we illustrate the different elements that appear in this logistic problem.The names are differentiated by color.In blue color we indicate the input data (WS, GPN and EC), whereas in green color we highlight the decisions to be made on: the (PT, AD and BL) plants and the (to GPN and BL) pipelines to be installed.The transportation routes are highlighted with red arrows in the plot.
The following hyphoteses are assumed in our model, based on the technological requirements of the process: • A minimum percentage of the total amount of waste in the WS must be collected and sent to the PT plants to satisfy the gas requirements of the GPN.Since a given percent of the produced biomethane is to be injected in the GPN, this hypothesis implies that a minimum amount of biogas must be produced to be injected to the network.
• Each PT plant sends its whole produced DOW to the AD plants.In the AD plants, the DOW received is transformed into biomethane.The digester solids of this transformation can be delivered to the WS as fertilizers.• We assume that each WS receives a given proportion of the total amount of digester solid produced at the AD plants.In case this proportion is 0, the WS does not receive any digester solid (as in the case of residual storages that may not be interested in the product).We do not force delivering the whole digester solid produced in the AD plants to the WS.
• The whole biomethane production obtained in the AD plants must be sent.It is shared into the GPN and the BL plants but, as already mentioned, a minimum percentage of the total biomethane production is required to be injected into the GPN.
• Each BL plant sends its whole LNG production to the EC.
• The demand of the EC is not assumed to be fully satisfied.Instead, we consider that every EC receives a pre-specified proportion of the total LNG produced at the BL plants.This proportion may depends on its demand or other preference criterion that the agents consider reasonable.
• A percent of the amount (weight and volume) of the different products is assumed to be lost at each phase of the whole process.Specifically, a given percentage of the amount of waste, DOW, and biomethane received in PT, AD and BL, respectively, is lost during the process.This technical hyphotesis is derived from the chemical processes to obtain the different types of products from their raw materials.
The efficient management of this system requires both deciding the location of the different types of plants (PT, AD, and BL) and pipelines (connecting AD with GPN and AD with BL) to be constructed, as well as the distribution of the product along the different phases of this process.
There is a cost for opening each potential location of each different type of plants, and a construction cost for unit length for the pipelines.We assume that a budget is given to install all the plants and the pipelines.
For the distribution of the products we assume that a unit transportation cost is provided for the different phases, namely, delivering from WS to the PT plants (waste), from PT to the AD plants (DOW), from the BL plants to the EC (LNG), and from AD to the WS (digester solid).
The Waste-to-Biomethane Logistic Problem (W2BLP) consists of deciding the location of the different plants and pipelines, and the way the product is distributed, minimizing the overall transportation cost with the given installation budget.

Mathematical Optimization Model
In this section we present the mathematical optimization model that we propose for the W2BLP.First, we define the parameters and variables that we use in our model.− δ: Proportion of waste transformed into DOW at the PT plants (0 < δ < 1).− γ 1 : Proportion of DOW transformed into biomethane at the AD plants (0 < γ 1 < 1).− γ 2 : Proportion of DOW transformed into fertilizers at the AD plants (0 < γ 2 < 1).We assume that γ 1 + γ 2 < 1.
− β: Proportion of biomethane transformed into LNG at the BL plants (0 < β < 1).− D ℓ : Demand of EC ℓ ∈ C. − p: Lower bound percent of the total waste produced at the WS that must be collected and sent to the PT plants (0 < p ≤ 1).
− q: Lower bound percent of the total biomethane produced at the AD plants that must be injected into the GPN (0 < q ≤ 1).− R i : Proportion of the total amount of fertilizer to be deliver to WS i ∈ N ( i∈N R i ≤ 1).
− α ℓ : Proportion of the total amount of LNG at the BL plants that will be received by demand point ℓ ∈ C ( ℓ∈C α ℓ = 1).

Costs. Set-up Costs:
These costs represent the installation costs of the different type of plants and links.They may include not only the construction cost of the different plants and pipelines, but also other costs related with the maintenance and use of them, as labor costs, land costs, etc.
− f 1 j : Set-up cost for opening PT plant j ∈ P 1 .− f 2 j : Set-up cost for opening AD plant j ∈ P 2 .− f 3 j : Set-up cost for opening BL plant j ∈ P 3 .− h jk : Set-up cost for installing a pipeline linking AD plant j ∈ P 2 with BL plant k ∈ P 3 .
− g jℓ : Set-up cost for installing a pipeline linking AD plant j ∈ P 2 with injection point in GPN ℓ ∈ I.
− B: Total budget for installing plants and pipelines.

Transportation Costs:
These costs represent the transportation costs of the different products (waste, DOW, LNG, and digester solid) that are usually transported by trucks from the different geographical locations of the elements in the system.
jℓ ≤ 0: Unit transportation cost (profit) of LNG from BL plant j ∈ P 3 to EC ℓ ∈ C. Note that this non positive cost represent the benefit of delivering to the EC each unit of liquified gas for their particular use.

− c N
ji ≥ 0: Unit transportation cost of digester solid from AD plant j ∈ P 2 to WS i ∈ N .Note that the transportation costs for the biogas through the pipelines are already included in the set-up costs described above.
3.2.Variables.Our model makes decision on the different plants and pipelines that are opened, as well as the amount of product delivered at the different phases of the process.

Binary variables:
The following variables determine whether a plant or a pipeline should be installed.

Continuous Variables:
The following variables of our model decide the amount of product delivered between the different sites in the system.
x 1 ij : amount of waste delivered from WS i to the PT plant j, ∀ i ∈ N , j ∈ P 1 .
x 2 jk : amount of DOW delivered from PT plant j to AD plant k, ∀ j ∈ P 1 , k ∈ P 2 .
x 3 jk : amount of biomethane delivered from AD plant j to BL plant k, ∀ j ∈ P 2 , k ∈ P 3 .
x I jℓ : amount of biomethane delivered from AD plant j to GPN's injection point ℓ, ∀ j ∈ P 2 , ℓ ∈ I.
x N ji : amount of digester solid delivered from AD j to WS i, ∀ j ∈ P 2 , i ∈ N .x C jℓ : amount of LGN delivered from BL plant j to EC ℓ, ∀ j ∈ P 3 , ℓ ∈ C.
3.2.1.Objective Function.The goal of the W2BLP is to minimize the total transportation cost of the system.For that, one must decide the optimal position of the different plants and pipelines with the given budget B. These locations have a direct impact in the transportation cost.
With the above set of variables, the overall transportation cost can be written as follows: (1) 3.3.Constraints.The assumptions of our problem are adequately established by the following constraints: • The budget for installing the different plants and pipelines must be satisfied: • Flow Conservation Constraints: The product must be adequately routed through the intermediate stages of the supply chain and enforced by the following constraints: Constraints (3) ensure that all the waste received at a PT plant is delivered to the AD plants once it is processed (100δ% of the waste is transformed into DOW).100γ 1 % of the amount of DOW received at each AD plant is transformed into biomethane and fully delivered either to the injection points or to the BL plants (Constraint (4)) and 100γ 2 % is converted into digester solid and delivered to the WS (Constraint (5)).Finally, Constraints (6) assure that the amount of biomethane received at the BL plants is converted into 100β% of LNG and fully ddelivered to the EC.
• Demand and Capacity Constraints: The different products obtained in the process are adequately delivered according to demands and capacities.Specifically, the WS must send waste and receive digester solid, the EC must receive LNG and the GPN biomethane.To this end, the following constraints are incorporated to our model: i∈N j∈P 1 x 1 ij ≥ pW (8) Constraints ( 7) assure that the amount delivered from each WS is at most the waste it produces.Constraint (8) forces to send, from the WS's, at least 100p% of the total amount of waste.Constraint ( 9) states that at least 100q% of the produced biomethane must be injected in the GPN.Each EC, by Constraints (10), at most receive its required demand.Note that the transportation cost for the product in this phase is non positive, and then, if possible, this demand will be satisfied.
Note also that, by (3), Constraints ( 8) and ( 9) can be equivalently rewritten as: It may happen that if q is large, the overall demand for the EC that can be satisfied is small.Thus, one can establish a sharing rule for the available LNG, i.e., a vector (α 1 , . . ., α |C| ) ∈ [0, 1] |C| representing the percent of the LNG at the BL plants that will be received by each EC.The following constraints assure the adequate verification of this share: Finally, Constraints (13) force that each WS receives its given percentage of the total production of digester solid: • Distribution through open Plants and Links: The transportation of the product in the different phases forces the plants and/or pipelines to be installed.This is assured by the following constraints: x 3 jk ≤ pδγ 1 W z jk , ∀j ∈ P 2 , ∀k ∈ P 3 , (16) These constraints assure that in case a plant or a link is not open, the product cannot be distributed through the corresponding plant or link.
• Compatibility of open plants and links: The installation of a pipeline linking a plant with other plant or with an injection point in the GPN, is subject to the previous instalation of the plants that the pipelines are connecting.These conditions are imposed in our model with the following linear constraints: Summarizing all the previous descriptions, the W2BLP can be model with the following mathematical optimization problem that we denote by (W 2BLP ) 0 : (2) − ( 7), ( 10) − (20), The above model is a mixed integer linear programming (MILP) problem, which can be difficult to solve in commercial solvers for real-world instances.In the following result we show that the size of the model can be considerably reduced based on the fact that once decided the potential AD plant which is open, it will optimally distribute the biogas to the GPN through its less costly injection point.where ℓ(j) = arg min ℓ∈I g jℓ .
Furthemore, the overall set-up costs of Ŝ is smaller or equal than the ones for S.
Proof.The proof follows straightforward since there are no transportation costs between the AD plants and the injections points but linking them is accounted in the budget constraint.Thus, replacing the link (j, ℓ) by (j, ℓ(j)) does not affect the objective function but reduces the set-up costs.□ The result above allows us to simplify the model by reducing the number of variables modeling the links and the flows between the AD plants and the GPN.Specifically, one can replace the two-index variables t and x I by one-index variables that, abusing of notation, we denote as: 1 if a pipeline linking AD plant j with injection point l(j) is built 0 otherwise ∀ j ∈ P 2 , and x I j : amount of biomethane delivered from AD plant j to GPN's injection point l(j), ∀ j ∈ P 2 .
Thus, if we denote by g j = g jℓ(j) for all j ∈ P 2 , one can replace the two-indices t and x I variables in model (W 2BLP ) 0 by the one-index t and x I variables above.Consequently, Constraints (2), ( 4), ( 11), and ( 17) are replaced by the following inequalities and equations: j∈P 2 x I j ≥ pqδγ 1 W, (23) As one can observe in the next section, the reduced model above is able to solve real instances of the W2BLP in reasonable CPU time.Our approach is then a useful tool for making decisions on the logistic of biogas production system.

Computational Experiments
We have run a series of experiments to analyze the computational performance of our approach.The main goal of the experiments is to determine the computational limitations of our model and its ability to obtain solutions for real-world instances.
We have randomly generated different instances with different sizes and parameters.We generate the coordinates for the WS, the potential location for the different plants (PT, AD,and BL), the injection points in GPN and the EC uniformly in [0, 1000] × [0, 1000].For the sake of simplification, we assume that the number and the location of WS (n) is the same that the number and the location of potential PT plants (m 1 ), and that the number and the location of potential AD plants (m 2 ) and potential BL plants (m 3 ) are also the same.Aditionally, the number of EC (|C|) also coincides with the number Let dist ij denote the Euclidean distance between locations i and j.We considered as unit transportation costs the Euclidean distances between the different points (dist ij ), except for the transportation costs of LNG from BL plants to the EC where the unit transportation cost is considered as a profit and it is defined as: These costs are designed to represent that as closer the extra customer to the open BL plant the larger the profit to satisfy its unit demand.The set-up costs for the PT plants and the BL plants were chosen all equal to one unit (in million of $) and the installation cost for the AD plants wer fixed to 5 units (in million of $).The set-up costs for the pipelines linking the AD plants with the injection points in GPN and the BL plants are 0.1 times the normalized (by the maximum) Euclidean distance between the corresponding points.
For determining adequate budgets for a given instance, we first solve the problem without the budget constraint (2).After solving the problem, we compute its effective set-up cost by adding up the costs of the open plants and pipelines that route a positive flow.We denote by B this set-up cost and consider as budget B = 0.2 B. This budget indicates that one can only use 20% of the cost of makig the best decision with no budget.
The slurry production at each farm (w i for i ∈ N ) has been uniformly generated in [1, 100] ∩ Z + (here Z + stands for the set of nonnegative integer numbers).The percent of the overall produced biomethane that must be injected into the existing GPN, q, ranges in {50%, 70%, 90%}.We have considered parameters p = 1, δ = 0.8, γ 1 = 0.8, γ 2 = 0.15, β = 0.7, α ℓ = D ℓ l∈C D l for all ℓ ∈ C, and The values of the parameters that we consider in our experiments are summarized in Table 1.
The model has been coded in Python 3.7 in an iMac with 3.3GHz with an Intel Core i7 with 4 cores and 16GB 1867 MHz DDR3 RAM.We used Gurobi 9.1.2as optimization solver.A time limit of 6 hours was fixed for all the instances.All the instances with n ≤ 100 were optimally solved.For the larger instances we fixed a MIP Gap limit (the solver stops when reaching such a limit and outputs the best feasible solution).For n = 200 we fix a MIPGap limit of 2%, for n = 500 of 4%.For each combination of parameters n, m, d, q we solved five instances.Thus, we have solved a total of 450 instances.
Figure 2 depicts performance profiles with respect to CPU times and MIP Gaps.In these pictures, we plots the percentage of instances solved up to each value of the CPU or MIP Gap, respectively.We analyze the computational performance of the model for the different values of the parameter q.As one can observe, the model seems to perform slightly better for the the smaller values of parameter q, although the differences are tiny.
Table 1.Summary of the parameters used in our experiments.
(a) CPU times (log scale).For this reason, in what follows we will not differentiate between the different values of the parameter q.
In Table 2 we summarize the results of our computational experiments.The first three columns indicate the size of the instances.Each of the rows summarize the results of a total of 15 instances (five random instances and three values of the parameter q).Column CPUTime (secs) is the average CPU time, in seconds, that the solver required to solve the instances.Column MIPGap indicates the average percent MIP Gap.In case the instance is optimally solved, the MIP Gap is 0%.Otherwise, this number reports either the MIPGap obtained within the time limit (if it is greater than the MIP Gap limit) or the MIP Gap when the MIP Gap limit was reached (that can be slightly smaller than the MIP Gap limit).Finally, column UnSolved indicates the percent instances summarized in the row that were not optimally solved or reaching the MIP Gap limit within the time limit.As mentioned before, we observe that all the instances with up to n = 100 have been optimally solved within the time limit.As expected, the computing time increases with the value of n and m, while it does not seem to depend on the value of d.The average computing time over all these instances is approximately 46.26 seconds, being 1.19 seconds for the instances with n = 25 , 11.11 seconds for n = 50, and 96.44 for n = 100.The average computing time needed to obtain an optimal solution depending on the value of m is 1.5 seconds for m = 5, 19.9 seconds for m = 10, and 104.8 seconds for m = 20.
Regarding the instances with n = 200 we observe that all the instances have been solved, which in this case means that all the instances have reached the MIP gap limit within the time limit, being the average consuming time over all the instances 300 seconds and the average MIP gap 1.57%.One can observe that, for these instances, the computing time is greater as the lower is the value of the parameter d.This may be because the value of parameter d does not affect to the number of binary variables in our improved formulation and it might happen that the smaller the number of EC and injection points, the more difficult it is to decide where to install the different plants to satisfy them at optimal cost.This fact can also be observed for the largest instances with n = 500.For d = 10, 20% of the instances with m = 50 and 6.67% of the instances with m = 100 have not been solved within the time limit, that is, these instances have not reached the MIP gap of 4% in 6 hours.The rest of the instances have been solved within the time limit.The average computing time for all the instances with n = 500 was 3057 seconds and the average MIP gap, 2.55%.
The results shown in Table 2, together with Figure 2, allow us to get a clear empirical evidence that using our improved formulation, the problem can be solved (with the MIP gap limit) within the time limit using a commercial solver.In fact, only 4 of the 450 instances (all of them with n = 500) have reached the time limit with a MIP Gap greater than the MIP Gap limit.This number represents less than 1% of the total instances and only 3.3% of the instances with n = 500.
From our experiments, we conclude that our approach provides a useful tool to make decisions about on the complex logistic problem behind the energy transformation of waste into biogas in realistic instances using reasonable computational resources.Thus, our model can also be used to evaluate different alternatives based on different values for the set-up costs for the plants or different transportation costs.

Case Study
We tested our model in a real-world dataset based on the region of upper Yahara watershed region in the State of Wisconsin (see Figure 3) whose data is available at https: //github.com/zavalab/JuliaBox/tree/master/Graph_S%26C.A detailed description of this dataset can be found in (Ma and Zavala, 2022;Sampat et al., 2019Sampat et al., , 2021;;Tominac et al., 2022).This region consists of 203 dairy farms whose waste production is known and whose locations are predefined.These locations determine our set of WS and the potential positions for the PT plants.The potential positions for the AD plants and the BL were randomly generated in that region.The positions of the injections points on the GPN were also randomly generated in this region and a network connecting this points was constructed.The position of the EC where randomly generated outside that region but in the minimum rectangular box containing it.
In Figure 4 we plot the coordinates of the different agents involved in the W2BLP.Red points in the plot indicate the WS locations (as larger the size of the dot, larger the production of waste in the farm) and the potential locations for PT, blue points are the potential locations for AD plants and for BL plants, green points indicate the EC  The parameters not provided in the above repository were fixed as follows: • The set-up budget was fix to κ × B, for κ ∈ {0.1, 0.2}, where B was computed as described in Section 4: • The rest of parameters have been considered as described in Section 4 (q ranges in {50%, 70%, 90%}, p = 1, δ = 0.8, γ 1 = 0.8, γ 2 = 0.15, β = 0.7, The parameters required to reproduce the obtained results are available at https:// github.com/vblancoOR/w2blp.As in the previous section, we set a time limit of 6 hours for solving the problem.In Table 3 we show the CPU times (in seconds) and MIP Gaps obtained after running our model to the six different configurations of budget and q.In that table, one can observe that except the case B = 0.1 B with q = 0.9, all the instances were optimally solved within the time limit.The case were the optimality is not guaranteed obtained a MIP gap of 0.1% which is insignificant.Furthermore, as expected, a more limited budget has a significative impact in the CPU time required to solve the problem, being the problems with budget 0.1 B more challenging than those with budget 0.2 B.
In Figures 5, 6, and 7 we show the results of our experiment for the Yahara watershed dataset for the different choices of budget, B, and q, and the different decisions that are made by our model.
In Figure 5, we show the results about the optimal locations for the PT plants (thick red points), for the AD plants (blue squares), for the BL plants (blue triangles) and for the links connecting AD plants with BL plants (blue lines) and AD plants with injection points of the GPN (black lines).As can be observed, the budget and the percent of production that must be injected in the GPN has a direct impact in the structure of the obtained solutions.On the one hand, as expected, as smaller the budget, smaller the number of plants and pipelines that are open.In particular, the budget mostly affect to the construction of pipelines to either inject to the GPN or transport to the BL, but also in the number of installed plants.We also observed that the network of the different installed plants and pipelines has more connected components for the larger budget, whereas is almost connected for the small budget.Regarding the value of q, it seems that as larger the value of q closer the new facilities to the GPN.In Figure 6 we show the distribution network between the WS plants, the PT plants, and the AD plants.Specifically, the links for which there is a positive waste flow from WS to PT plants (red lines), a DOW flow from PT plants to AD plants (purple lines) and digester solid flow from AD plants to WS (yellow lines).The main observation that can be drawn from these plots is that in most of the cases there are differentiated clusters of WS plants (farms) that share the same PT plants and AD plants.This observations implies that for larger datasets where the model could not be able to solve the problem, the optimal solution would be adequately approximated by clustering the WS (with an adequate criterion) and solve the problem separately for each cluster.
Finally, in Figure 7, we show the biomethane flow from AD plants to BL plants (blue lines) and to injection points of the GPN (black lines), and LGN flow from BL plants to EC (green lines).Note that some of the AD plants are devoted only to give service to the GPN, others that serve exclusively the BL plants (and then the EC), but one can also  find AD plants that send part of the production to the GPN and what remains, to the BL plants.This behaviour would have never happen if an integrated model as the one we propose would not have been considered.Note also that a single AD plant is allowed to send biomethane to different BL plants.

Conclusions
In this paper, we propose a mathematical optimization approach to design an optimal logistic system to distribute the different products involved in the generation of biogas from waste.Specifically, given a set of waste storage centers, an existing gas pipeline network, and a set of external customers, we provide a decision aid tool to determine the number and optimal locations of pre-treatment plants, anaerobic digestion plants, and biomethane liquefaction plants, as well as the pipelines linking some of these plants.Additionally, we provide a distribution plan to send the different wastes, to serve either the gas pipeline network or the external customers with the produced biogas, and to serve the waste storage centers with the generated fertilizer.All the decisions are made by minimizing the  transportation costs of the different products and restricting the installation of plants and pipelines to a given budget.We develop a new mixed integer linear programming formulation for the problem and prove some results that allow us to reduce the size of the model.With this simplification, we were able to obtain optimal solutions for the problem in real-world intances.
We report the results of an extensive battery of synthetic computational experiments, concluding that our approach is suitable to be applied to different settings and sizes.We also analyze the case study of the Yahara watershed, and study the obtained solutions based on different parameters.
Our future research on this topic includes the incorporation of uncertainty in some parameters of the model.Concretely, the production of waste in the farms or waste storages is known to vary in the different seasons of the year, and this production may has an impact in the solution of the problem.We will design stochastic optimization models that take into account this uncertainty to construct robust solution of the W2BLP.Additionally, we will consider capacity constraints for the different plants and multiple periods for the decisions made in our model.Instead of assuming that the plants and pipelines are installed here-and-now, we will decide in which period of a time horizon each of the installations is constructed and used, assuming that they have a limited capacity and that a budget is available for each period.The integrated model that consider uncertainty, capacity constraints and multiperiod decisions will be closer to reality, but the mathematical programming model for it would be prohibitive even for small instances.Thus, a different solution approach will have to be designed to solve it, that would not be exact, but heuristic.

Figure 1 .
Figure 1.The process modeled in the W2BLP.

Figure 2 .
Figure 2. Performance profile of our computational experience for the different values of parameter q. n

Figure 4 .
Figure 4. Input Data of the Case Study.

Figure 5 .
Figure 5. Optimal location of the different types of plants and pipelines for different values of B and q.

Figure 6 .
Figure 6.Distribution network between WS, PT plants, and AD plants.

Figure 7 .
Figure 7. Distribution network between AD plants, BL plants, injection points of the GPN, and EC.

Table 2 .
Summary of our computational experience.

Table 3 .
D ℓ was randomly generated in [δγ 1 β × min w i , δγ 1 β × max w i ].CPU times (in seconds) and MIP gaps for solving the instances in the case study.