DOI: 10.3303/CET25117176 Paper Received: 15 December 2024; Revised: 8 April 2025; Accepted: 30 April 2025 Please cite this article as: Villamizar-Villamizar A.E., Rosso-Ceron A.M., Mayo-Florez D.J., 2025, Mathematical Programming for Sustainable Heat Alternatives Selection in a High-altitude Region: a Case Study of Pamplona-north of Santander, Chemical Engineering Transactions, 117, 1051- 1056 DOI:10.3303/CET25117176 CHEMICAL ENGINEERING TRANSACTIONS VOL. 117, 2025 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Fabrizio Bezzo, Flavio Manenti, Gabriele Pannocchia, Almerinda di Benedetto Copyright © 2025, AIDIC Servizi S.r.l. ISBN 979-12-81206-17-5; ISSN 2283-9216 Mathematical Programming for Sustainable Heat Alternatives Selection in a High-Altitude Region: a Case Study of Pamplona-North of Santander Alvaro E. Villamizar-Villamizar, Ana M. Rosso-Cerón*, David J. Mayo-Florez Chemical Engineering Program, Universidad de Pamplona, km 1 vía Bucaramanga; Pamplona-Colombia ana.rosso@unipamplona.edu.co This study introduces an innovative mathematical programming methodology to determine optimal sustainable heating alternatives for high-altitude regions. Employing integer linear programming, a multi-objective optimization model was developed to simultaneously minimize the total present value and CO₂ emissions. The model integrates a comprehensive framework of design, operational, and financial constraints to address the heating demands over a 19-year period starting in 2021. The analysis is focussed on Pamplona, North of Santander, a region characterized by challenging climatic conditions and specific energy needs. The model evaluates diverse heat generation technologies, including gasification, anaerobic digestion, solar collectors for heating water, and conventional fuel-based systems and natural gas grid connection. These considerations offer valuable insights for policymakers and stakeholders aiming to advance sustainable energy solutions in comparable settings. 1. Introduction Pamplona, located in the Noth of Santander, Colombia, within the Páramo de Santurbán, lies at an altitude of 2,350 metres above sea level. The region is characterised by a cold and humid climate, with annual average temperatures ranging from 5°C to 16°C. This páramo ecosystem plays a vital role in capturing and regulating water resources, which are essential for the downstream agricultural and urban areas. Given the cold climate, the demand for heating is significant, particularly during colder months, thus necessitating efficient and sustainable energy solutions for both residential and institutional buildings. Furthermore, the energy requirements are increased due to the population of over 56,894 people, whose daily activities, education, and commerce depend on reliable heating systems (Sistema Único de Información, 2022). Unlike many countries with established district heating networks, Colombia lacks this infrastructure. Instead, the supply of heat energy relies on a combination of decentralised and centralised systems. In some urban areas, centralised natural gas networks provide a reliable source of heat for both heating and cooking (UPME, 2022). In contrast, other sources, such as hot water systems, solar radiation, liquefied petroleum gas (LPG), and biomass for processes like anaerobic digestion and pyrolysis, behave in a decentralised manner. While these solutions cater to local needs, they lack the integration and efficiency of district heating networks, underscoring the need for scalable, innovative approaches to meet heating demands, especially in high-altitude regions like Pamplona. In the region, heating distribution systems are generally organised into three types: the first is decentralised, where individual users install their own heating units, typically in isolated homes. The second is centralised, where a heating plant connected to the natural gas network distributes heat to small communities or residential complexes. The third approach involves expanding the natural gas grid to reach wider urban and rural areas, providing heat to high-altitude zones. In response to these challenges, this study proposes the adoption of kWh as the standard unit for measuring heat supply. Standardising heat measurement in kWh offers several advantages: it allows for integration with electrical and renewable energy systems, simplifies billing by applying charges like those used for electricity, 1051 and facilitates comparisons across various energy sources. This approach provides a solid foundation for the development of efficient, scalable, and user-friendly heating solutions tailored to local conditions. The need for sustainable heating solutions becomes even more urgent when considering the recent changes in energy demand observed in Pamplona. These shifts have been significantly influenced by the COVID-19 pandemic, which altered residential and commercial energy consumption patterns. Villamizar-Villamizar et al. (2022) employed LEAP software to project long-term energy demand in the region, noting an increase in residential consumption as people spent more time indoors during lockdowns. This trend underscores the growing necessity for energy-efficient and sustainable heating technologies to meet both current and future demands. In response to these emerging needs, researchers have investigated strategies to optimise energy systems. Stamatellou et al. (2021) examined various building ventilation techniques to improve indoor air quality in a post- pandemic context. Their findings highlighted that while enhanced ventilation increases heating costs, it is essential for maintaining healthy indoor environments, particularly in cold climates like Pamplona. These insights have contributed to further research into cost-effective and energy-efficient heating solutions for regions with similar climatic challenges. Subsequent studies have explored innovative approaches to heat generation in areas with high demand. Jahangir et al. (2022) studied hybrid power and heat generation systems during the COVID-19 crisis in Italy, demonstrating the benefits of integrating renewable sources such as solar heating and biomass to reduce reliance on conventional fuels. These hybrid systems effectively addressed fluctuating energy demands and reduced greenhouse gas emissions, providing promising solutions for regions like Pamplona. Karasu (2023) further advanced this field by optimising thermal requirements in smart district heating systems. His research demonstrated the potential of advanced control strategies to improve energy efficiency and reduce costs through real-time demand adaptation. This approach is particularly relevant to Pamplona, where flexible and efficient heating systems are critical due to the region's climate conditions. Moreover, optimised heating systems not only address immediate energy demands but also contribute to broader sustainability goals outlined in Colombia's national energy strategies, such as the PROURE programme (Ministerio de Minas y Energía, 2023). From this perspective, this study proposes a novel multi-objective mixed integer linear programming (MOMILP) approach to identify optimal sustainable heating alternatives for Pamplona. The MOMILP model aims to minimise net present value (NPV) while reducing CO₂ emissions, by incorporating various heat generation technologies such as biomass gasification, pyrolysis, LPG, and solar collectors. These technologies are designed to meet thermal demand in both centralised and decentralised areas through the extension of natural gas networks. The model accounts for specific design, operational, and budgetary constraints. This research provides valuable insights for policymakers and stakeholders, contributing to the advancement of energy efficiency and sustainability in Colombia's high-altitude regions. 2. Mathematical Formulation This section presents the MOMILP model, including the objective functions and the design, operational, and budgetary constraints for the heating system. 2.1 Objective Functions The objective functions are designed to minimize the NPV and CO₂ emissions. The NPV is calculated as the sum of the total capital costs of each heating plant utilized, along with fixed and variable operation and maintenance costs, fossil fuel expenses, and the costs associated with extending gas grids (Eq. (1)). CO₂ emissions include those produced by fuel-based plants, such as biomass gasifiers and diesel plants (Eq. (2)). NPV = ∑ 1 (1+𝑅)𝑡𝑡 [(∑ 𝐼𝑛𝐶𝑜𝑝,𝑡 ∑ 𝐶𝐴𝑧,𝑝,𝑡𝑧𝑝 ) + (∑ 𝑉𝑎𝑟𝐶𝑜𝑡,𝑝 ∑ 𝐸𝐺𝑧,𝑝,𝑡𝑧𝑝 ) + (∑ 𝐹𝑖𝑥𝐶𝑜𝑝,𝑡 ∑ 𝐶𝑎𝑇𝑧,𝑝,𝑡𝑧𝑝 ) + (∑ 𝐹𝑢𝐶𝑜𝑟𝑒.𝑡 ∑ 𝐹𝑧,𝑟𝑒,𝑡𝑧𝑟𝑒|𝑟𝑒∈𝑟𝑒𝑐 ) + (∑ 𝐸𝑇𝑟𝑡,𝑧 ∗ 𝑇𝑟𝑎𝐶𝑜𝑡,𝑧𝑧|𝑧∈𝑐 + ∑ 𝐼𝑛𝑇𝑟𝐶𝑜𝑡,𝑧𝑧|𝑧∈𝑐 ∗ 𝐿𝑧 ∗ 𝑇𝑟𝑧,𝑡 )] (1) 𝐶𝑂2 emissions = ∑ [ ∑ (𝐸𝑧 ∗ 𝐸𝑇𝑟𝑡,𝑧 + ∑ 𝐹𝐸𝑧,𝑝 ∗ 𝐸𝐺𝑧,𝑝,𝑡𝑝 )𝑧 ]𝑡 (2) 2.2 Constraints The constraints of the model are grouped into various categories, such as operational, design, and financial limitations for heating plants. 1052 2.2.1 Operational Constraints The operational constraints ensure the fulfilment of heating requirements for both water and household applications. Energy Balance: The heating produced by each plant p, during period t must satisfy the corresponding heating demand, for both household applications (Eq. (3)) and hot water (Eq. (4)). ∑ (𝐸𝐺𝑧,𝑝,𝑡) + 𝐸𝑇𝑟𝑧,𝑡 =𝑝|𝑝∈𝑝𝑐 𝐷𝐸𝑧,𝑡 ∗ (1 − 𝑅𝐴𝑧,𝑡); ∀ 𝑧, 𝑡 (3) ∑ 𝐸𝐺𝑧,𝑝,𝑡 = 𝑅𝐴𝑧,𝑡𝑝|𝑝∈𝑝𝑎 𝐷𝐸𝑧,𝑡; ∀ 𝑧, 𝑡 (4) Thermal Generation Limit Condition: Reflects that the energy generated by the plant is limited by the product of operating hours and the amount of heat generated, ensuring operational capacities are not exceeded (Eq. (5)). 𝐸𝐺𝑧,𝑝,𝑡 ≤ 𝐵𝐷𝑇𝑧,𝑝 ∗ Q𝑧,𝑝,𝑡; ∀ 𝑧, 𝑝, 𝑡 (5) Installed Capacity Efficiency: Reflect that the heat flow is limited by the product of efficiency and installed thermal capacity, ensuring that the heat production does not exceed the plant's available capacity (Eq. (6)). Q𝑧,𝑝,𝑡 ≤ 𝑈𝑧,𝑝 ∗ 𝐶𝑎𝑇𝑧,𝑝,𝑡; ∀ 𝑧, 𝑝, 𝑡 (6) Extension lines of natural gas: The heat transmitted through the gas lines cannot exceed the amount of heat supplied if the line is operational (Eq. (7)). Additionally, the pipelines capacity cannot be activated until the decision to implement the gas networks has been made (Eq. (8)). 𝐸𝑇𝑟𝑧,𝑡 ≤ 𝑇𝑟𝑜𝑧,𝑡 ∗ 𝑇𝐷𝑧 ∗ 𝐿𝑧 ∗ 𝐶𝑇𝑧,𝑡 ; ∀ 𝑧|𝑧 ∈ 𝑧𝑐, 𝑡 (7) 𝑇𝑟𝑜𝑧,𝑡 ≤ ∑ 𝑇𝑟𝑧,𝑣𝑣≤𝑡−𝑇𝑡 ; ∀ 𝑧|𝑧 ∈ 𝑧𝑐, 𝑡 (8) 2.2.2 Design Constraints These constraints incorporate installed capacities as decision variables, considering resource availability, annual limits on the number of plants that can be installed, and additional capacity restrictions. Available Installed Capacity: The previously installed capacity plus the newly added capacity (Eq. (9)). 𝐶𝑎𝑇𝑧,𝑝,𝑡 = CAO𝑧,𝑝 ∗ 𝐼𝐶𝐴𝑧,𝑝 + ∑ 𝐶𝐴𝑧,𝑝,𝑡 ; ∀ 𝑧, 𝑝, 𝑡𝑣≤𝑡−𝐵𝑇 (9) Added Capacity Limits: Limits are set for the added capacity (Eq. (10)), ensuring that these do not exceed the maximum capacity throughout the entire planning horizon (Eq. (11)). 𝑃𝑎𝑚𝑖𝑛𝑧,𝑝 ∗ 𝐵𝑖𝑛𝑃𝑧,𝑝,𝑡 ≤ 𝑃𝐴𝑧,𝑝,𝑡 ≤ 𝑃𝑎𝑚𝑎𝑥𝑧,𝑝 ∗ 𝐵𝑖𝑛𝑃𝑧,𝑝,𝑡 ; ∀ 𝑧, 𝑝, 𝑡 (10) ∑ ∑ 𝑃𝐴𝑧,𝑝,𝑡 𝑧𝑡 ≤ 𝑃𝑎𝑚𝑎𝑥𝑧,𝑝 ; ∀ 𝑧, 𝑝 (11) Conversion of Primary Energy into Heat: The primary sources used for heating generation must not exceed the available flow of resources (Eq. (12)). 𝐸𝐺𝑧,𝑝,𝑡 ≤ 𝑁𝑟𝑒,𝑝 ∗ 𝑃𝑟𝑖𝐹 𝑧,𝑟𝑒,𝑝,𝑡 ∗ 𝐶𝑝 𝑟𝑒,𝑝 ; ∀ 𝑧, 𝑟𝑒, 𝑝, 𝑡 (12) Specific Consumption of Thermal Energy Sources: The total amount of available thermal energy sources to be consumed during a period can be expressed as the sum of the available fuel in each plant (Eq. (13)). Besides, the total amount of available fuel to be consumed locally, is restricted by the local energy resource (Eq. (14)). 𝐹𝑧,𝑟𝑒,𝑡 = ∑ 𝑃𝑟𝑖𝐹𝑧,𝑟𝑒,𝑝,𝑡 ; ∀ 𝑧, 𝑟𝑒, 𝑡 𝑝 (13) 𝐹𝑧,𝑟𝑒,𝑡 ≤ 𝑅𝐸𝑧,𝑟𝑒,𝑡 ; ∀ 𝑧, 𝑟𝑒, 𝑡 (14) 1053 2.2.3 Budgetary Constraint The total cost of the hybrid system to be installed is limited by the financial support allocated for the new capacities and the expansion of natural gas network coverage (Eq. (15)). ∑ 1 1+𝑅𝑡𝑡 [(∑ 𝐶𝐼𝑝,𝑡 ∗ 𝑃𝐴𝑧,𝑝,𝑡𝑝 ) + 𝐼𝑛𝑇𝑟𝐶𝑜𝑧,𝑡 ∗ 𝐿𝑧 ∗ 𝑇𝑟𝑧,𝑡] ≤ 𝐵𝑈𝐺𝑧; ∀ 𝑧 (15) 2.3 Case Study The MOMILP model was applied to Pamplona, which spans 29,840 hectares, 24.2% of which are designated as páramo zones, crucial for biodiversity and the water supply to 2.5 million people. The municipality, which is connected to the natural gas grid, has a heating demand of 648,100 kWh for approximately 56,894 residents (Villamizar et al., 2022). This demand is influenced by climate, population density, and socio-economic factors. Additionally, the model inputs included cost data from the International Energy Agency (2024), energy capacities from the Sistema Único de Información (2022), energy resources from NASA POWER (2024), and emission factors from the Ministry of Mines and Energy, Colombia (2023). 3. Results and Analysis The MOMILP model, developed in GAMS using the MIP SIMPLEX algorithm and solved through the ε-constraint method, provides valuable insights into the trade-offs between maximizing profits and minimizing CO₂ emissions. A critical finding is the significant impact of adjusting the NPV tolerance gap, which directly influences emission outcomes. Near-optimal cost solutions yield rapid CO₂ savings with minimal profit loss; however, this efficiency diminishes as the tolerance gap widens. The analysis of the Pareto front (Figure 1) highlights that maximizing capital recovery inevitably leads to higher emissions, as observed in alternative scenario (A1), which achieves the highest profit (3,353.72 thousand USD) but also the greatest emissions (11,008.05 thousand kg CO₂). Conversely, alternative scenario (A10), with the lowest profit (2,939.89 thousand USD), corresponds to the lowest emissions (7,670.89 thousand kg CO₂). These variations are attributed to the capacity and technological choices of the plants included in each alternative. The Pareto front follows a convex trend, where initial emission reductions are cost-effective but become progressively more expensive as deeper reductions are pursued. Early reductions rely on economical strategies, such as integrating low-emission technologies like LPG, digesters, and solar thermal accumulators, which provide substantial gains at relatively low costs. However, as emission targets become more ambitious, advanced measures such as solar thermal collectors or biomass systems must be deployed, significantly raising marginal costs. This trend, driven by the transition to cleaner technologies, aligns with the findings of Riedmüller et al. (2024) and Gbadeyan et al. (2024), who observed that costlier measures dominate as low-cost opportunities are exhausted. Furthermore, the phenomenon of diminishing returns steepens the marginal cost curve, a hallmark of a convex Pareto frontier (Vannoni et al., 2021). Additionally, CO₂ emissions are not solely determined by the level of investment but also by operational factors, highlighting the variability of outcomes for identical portfolios. This underscores the importance of precise planning and optimization in designing sustainable energy systems. The MOMILP model effectively illustrates how achieving emission reductions involves balancing economical strategies with investments in advanced technologies. As such, the results emphasize the need for a nuanced approach to decision-making, particularly in energy transitions, where both economic feasibility and environmental sustainability are critical. Figure 1. Pareto front of CO₂ emissions and net present value for ten alternative scenarios The energy matrix participation analysis for heat generation from 2021 to 2040 (Figure 2.), with a total demand of 14,288,000 kWh, shows a consistent role for solar collection and biodigestion. Solar collection remains fixed at 0.7%, indicating its stable but modest contribution to meeting the heat demand. Similarly, biodigestion 1054 contributes 1.6% throughout the period, highlighting its reliable role in sustainable energy generation. However, the natural gas grid shows a steady increase in its share, rising from 24.6% in the early years to 67.1% by 2040, reflecting the growing reliance on conventional energy sources as the primary solution for heat generation. In contrast, LPG combustion experiences a gradual decline, from 73.1% in 2021 to 30.6% by 2040. This decrease suggests a shift away from LPG as a primary heat source, driven by the increased reliance on natural gas and the fixed contributions from renewable sources. This trend indicates a transition to more sustainable and cleaner energy options while meeting the rising heat demand over the planning horizon. Figure 2. Percentage Contribution of Heat Generated Over the Planning Horizon for 10 Alternative Scenarios 4. Conclusions The application of the MOMILP model to the municipality of Pamplona, Colombia, reveals a significant trade-off between NPV and CO₂ emissions in the energy generation system. The Pareto front analysis illustrates a clear pattern where initial emission reductions can be achieved through cost-effective strategies, such as the integration of solar thermal collectors (fixed at 0.7% contribution) and biodigesters (1.6% contribution). These renewable technologies provide substantial emissions reductions with relatively low costs. However, as the demand for deeper emissions reductions grows, the marginal cost of achieving further reductions increases due to the need for more advanced technologies. For instance, while the natural gas grid steadily increases its share from 24.6% in 2021 to 67.1% by 2040, the contribution from LPG combustion declines from 73.1% to 30.6% over the same period. This shift reflects a transition towards cleaner energy options. The convex nature of the Pareto front aligns with findings in the literature, such as Riedmüller et al. (2024) and Gbadeyan et al. (2024), where initial cost-effective emission reduction opportunities are exhausted, requiring higher investments for deeper reductions. Furthermore, the results underscore the importance of balancing economic feasibility with environmental sustainability, as both elements play critical roles in energy system planning. As a recommendation for future research, it is suggested that the model be expanded to include the integration of energy storage systems. This would allow for a more comprehensive analysis of the region's potential for decarbonizing the energy matrix while meeting growing heat demands. Additionally, exploring policy incentives and regulatory frameworks would provide further insights into how government intervention could support the adoption of clean technologies and facilitate the transition to a low-carbon energy system in regions with similar environmental and socio-economic conditions. Nomenclature Sets Model parameter (Units) re - Resources R - Discount rate p - Plants InCop,t - Unit investment cost, USD/kW t - Period VarCot,p - Unit variable cost, USD/kWh z - Zones FixCop,t - Unit fuel cost, USD/kW Subsets: FuCore.t - Unit fuel cost, USD/kg D ⊆ z − Decentralized zones TraCot,z - Unit operation cost of natura gas lines, USD/kWh C ⊆ z - Centralized zones InTrCot,z - Unit investment cost of new lines of natural gas, USD/km pa ⊆ p - Hot water plants Lz - Length of natural gas of pipelines, km pb ⊆ p – Fuel thermal plants Ez - Emissions factor of natural gas, kg CO2/kWh 1055 Continuos Variables FEz,p - Emission factor of local resources, kg CO2/kWh CaTz,p,t - Total installed capacity, kW DEz,t - Heat demand, kWh Qz,p,t - Operational capacity, kW RAz,t - Fraction of demand for hot water CAz,p,t - Additional capacity to install, kW BDTz,p - Annual operating time of heat plants, h EGz,p,t - Generated heat, kWh Uz,p - Heat generation capacity efficiency PriFz,re,p,t - Total fuel flow consumes, kg/y CTz,t - Heat generated per pipeline length, kW/km Fz,re,t - Total fuel flow available, kg/y TDz - Time factor, h ETrz,t - Heat Transfer, kWh Pamaxz,p - Maximum capacity to be added, kW Binary Variables Paminz,p - Minimum capacity to be added, kW BinPz,p,t - If a new capacity is installed, otherwise Cpre,p - Heating value of primary source BinCz,p - If an old capacity is operated, otherwise Nre,p - Primary thermal resource to heat conversion efficiency Troz,t - If new natural gas pipelines operate, otherwise REz,re,t - Primary resource available, Unit of source Trz,t - If new natural gas pipelines are installed, otherwise) BUGz - Total budget for new heating plants and networks over the entire period, USD Acknowledgments The authors would like to acknowledge the financial support from the Ministry of Science, Technology, and Innovation of Colombia for the project "Development of a computational tool and renewable energy technologies for energy transition in high mountain areas post-pandemic," CD 82605 CT ICETEX 2022-0644. References Gbadeyan O.J., Muthivhi J., Linganiso L.Z., Deenadayalu N., 2024, Decoupling economic growth from carbon emissions: A transition toward low-carbon energy systems—A critical review, Clean Technologies, Vol 6(3), 1076–1113, doi:10.3390/cleantechnol6030054. International Energy Agency, 2024, "Cost Data for Energy Systems," International Energy Agency, https://www.iea.org, accessed 28.07.2024. Jahangir M.H., Zavvari S., Parvin M., Yanıkoğlu İ., Aykut A., Ay D., Özpoyraz İ., 2022, Optimizing the Size of the Hybrid Power and Heat Generation System during COVID-19 Crisis (Case Study: Italy), Energy Conversion and Management, Vol 259, 115569, doi:10.1016/j.enconman.2022.115569. Karasu M.B., 2023, Thermal Request Optimization of a Smart District Heating System, Energy and Buildings, Vol 275, 112405, doi:10.1016/j.enbuild.2023.112405. Ministerio de Minas y Energía, Colombia, 2023, PROURE: Promoting Rational and Efficient Use of Energy in Colombia’s Highlands, Energy Policy and Planning Journal. NASA POWER, 2024, "Energy Resources Data," NASA POWER Data Viewer, https://power.larc.nasa.gov/data-access-viewer/, accessed 25.08.2024. Riedmüller S., Rivetta F., Zittel J., 2024, Long-Term Multi-Objective Optimization for Integrated Unit Commitment and Investment Planning for District Heating Networks, arXiv Preprint, doi:10.48550/arXiv.2410.06673. Stamatellou A.-M., Zogou O., Stamatelos A., 2021, Energy Cost Assessment and Optimization of Post-COVID- 19 Building Ventilation Strategies, Applied Energy, Vol 289, 116658, doi:10.1016/j.apenergy.2021.116658. SUI, 2022, Energy Consumption Trends in High Altitude Regions, Journal of Sustainable Energy Systems. UPME, 2022, National Energy Plan 2020-2050, UPME Report Vannoni A., Sorce A., Traverso A., Massardo A.F., 2021, Techno-Economic Analysis of Power-to-Heat Systems, E3S Web of Conferences, Vol 238, 03003, doi:10.1051/e3sconf/202123803003. Villamizar-Villamizar Á.E., Rosso-Cerón A.M., Becerra-Rodríguez J., 2024, Post-pandemic long-term energy demand forecasting using LEAP software in the Páramo de Santurbán: The case of Pamplona North of Santander, Chemical Engineering Transactions, Vol 109, 121–130, doi:10.3303/CET241109. 1056 https://power.larc.nasa.gov/data-access-viewer/ CET-vol117-b.pdf 180villamizarvillamizar.pdf Mathematical Programming for Sustainable Heat Alternatives Selection in a High-Altitude Region: a Case Study of Pamplona-North of Santander