increase tolerances from 1e-5 to 2e-5

Author: arght

CHANGELOG.rst

@@ -1,6 +1,10 @@
 Change Log
 =============
 
+[4.18.17RC] - 2026-04-11
+-------------------------
+- [CHANGED] change from 1e-5 to 2e-5 tolerances to avoid extreme values in the optimization problem and to improve the robustness of the model
+
 [4.18.16] - 2026-04-10
 -------------------------
 - [CHANGED] add H2 network utilization in output results

doc/rst/conf.py

@@ -12,7 +12,7 @@
 author = 'Andres Ramos'
 
 # The short X.Y version
-version = 'version 4.18.16'
+version = 'version 4.18.17RC'
 # The full version, including alpha/beta/rc tags
 release = ''
 
@@ -84,13 +84,13 @@
 #
 # html_sidebars = {}
 html_theme = 'alabaster'
-html_title = 'version 4.18.16'
+html_title = 'version 4.18.17RC'
 html_logo  = '../img/openTEPES.png'
 html_last_updated_fmt = ''
 html_show_sphinx = False
 html_theme_options = {
     'analytics_id': 'UA-515200-2',  #  Provided by Google in your dashboard
-    'description': 'version 4.18.16',
+    'description': 'version 4.18.17RC',
     'body_max_width' : 'none',
     'page_width': 'auto',
     'font_family': 'Georgia'

openTEPES/__init__.py

@@ -14,7 +14,7 @@
         >>> import openTEPES as oT
         >>> oT.routine("9n", "C:\\Users\\UserName\\Documents\\GitHub\\openTEPES", "glpk")
 """
-__version__ = "4.18.16"
+__version__ = "4.18.17RC"
 
 from .openTEPES_Main             import main
 from .openTEPES                  import *

openTEPES/openTEPES.py

@@ -39,8 +39,8 @@ def openTEPES_run(DirName, CaseName, SolverName, pIndOutputResults, pIndLogConso
     idxDict['y'  ] = 1
 
     #%% model declaration
-    mTEPES = ConcreteModel('Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.16 - April 10, 2026')
-    print(                 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.16 - April 10, 2026', file=open(f'{_path}/openTEPES_version_{CaseName}.log','w'))
+    mTEPES = ConcreteModel('Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - April 10, 2026')
+    print(                 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - April 10, 2026', file=open(f'{_path}/openTEPES_version_{CaseName}.log','w'))
 
     pIndOutputResults = [j for i,j in idxDict.items() if i == pIndOutputResults][0]
     pIndLogConsole    = [j for i,j in idxDict.items() if i == pIndLogConsole   ][0]

openTEPES/openTEPES_InputData.py

@@ -1,5 +1,5 @@
 """
-Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - March 17, 2026
+Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - April 11, 2026
 """
 
 import time
@@ -931,7 +931,7 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
 
     # variable emission cost [M€/GWh]
     mTEPES.dPar['pVariableEmissionCost'] = mTEPES.dPar['pVariableEmissionCost'].replace(0.0, mTEPES.dPar['pCO2Cost']) #[€/tCO2]
-    mTEPES.dPar['pEmissionVarCost']      = mTEPES.dPar['pEmissionRate'] * 1e-3 * mTEPES.dPar['pVariableEmissionCost']                 #[M€/GWh] = [tCO2/MWh] * 1e-3 * [€/tCO2]
+    mTEPES.dPar['pEmissionVarCost']      = mTEPES.dPar['pEmissionRate'] * 1e-3 * mTEPES.dPar['pVariableEmissionCost']                 # [M€/GWh] = [tCO2/MWh] * 1e-3 * [€/tCO2]
     mTEPES.dPar['pEmissionVarCost']      = mTEPES.dPar['pEmissionVarCost'].reindex(sorted(mTEPES.dPar['pEmissionVarCost'].columns), axis=1)
 
     # minimum up- and downtime and maximum shift time converted to an integer number of time steps
@@ -1122,9 +1122,9 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
     mTEPES.dPar['pDemandElecPeak']          = pd.Series([0.0 for p,ar in mTEPES.par], index=mTEPES.par)
     mTEPES.dPar['pDemandHeatPeak']          = pd.Series([0.0 for p,ar in mTEPES.par], index=mTEPES.par)
     for p,ar in mTEPES.par:
-        # values < 1e-5 times the maximum demand for each area (an area is related to operating reserves procurement, i.e., country) are converted to 0
+        # values < 2e-5 times the maximum demand for each area (an area is related to operating reserves procurement, i.e., country) are converted to 0
         mTEPES.dPar['pDemandElecPeak'][p,ar] = mTEPES.dPar['pDemandElec'].loc[p,:,:][[nd for nd in d2a[ar]]].sum(axis=1).max()
-        mTEPES.dPar['pEpsilonElec']          = mTEPES.dPar['pDemandElecPeak'][p,ar]*1e-5
+        mTEPES.dPar['pEpsilonElec']          = mTEPES.dPar['pDemandElecPeak'][p,ar]*2e-5
 
         # these parameters are in GW
         mTEPES.dPar['pDemandElecPos']     [mTEPES.dPar['pDemandElecPos'] [[nd for nd in d2a[ar]]] <  mTEPES.dPar['pEpsilonElec']] = 0.0
@@ -1186,7 +1186,7 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
 
         if mTEPES.dPar['pIndHeat']:
             mTEPES.dPar['pDemandHeatPeak'][p,ar] = mTEPES.dPar['pDemandHeat'].loc[p,:,:][[nd for nd in d2a[ar]]].sum(axis=1).max()
-            mTEPES.dPar['pEpsilonHeat']          = mTEPES.dPar['pDemandHeatPeak'][p,ar]*1e-5
+            mTEPES.dPar['pEpsilonHeat']          = mTEPES.dPar['pDemandHeatPeak'][p,ar]*2e-5
             mTEPES.dPar['pDemandHeat']             [mTEPES.dPar['pDemandHeat']    [[nd for nd in   d2a[ar]]] <  mTEPES.dPar['pEpsilonHeat']] = 0.0
             mTEPES.dPar['pHeatPipeNTCFrw'][mTEPES.dPar['pHeatPipeNTCFrw'] < mTEPES.dPar['pEpsilonElec']] = 0
             mTEPES.dPar['pHeatPipeNTCBck'][mTEPES.dPar['pHeatPipeNTCBck'] < mTEPES.dPar['pEpsilonElec']] = 0
@@ -1211,7 +1211,7 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
     mTEPES.dPar['pEmissionVarCost']   = mTEPES.dPar['pEmissionVarCost'].loc  [:,mTEPES.g ]
 
     # replace < 0.0 by 0.0
-    mTEPES.dPar['pEpsilon'] = 1e-6
+    mTEPES.dPar['pEpsilon'] = 1e-6 # measure of how close two values are
     mTEPES.dPar['pMaxPower2ndBlock']  = mTEPES.dPar['pMaxPower2ndBlock'].where (mTEPES.dPar['pMaxPower2ndBlock']  > mTEPES.dPar['pEpsilon'], 0.0)
     mTEPES.dPar['pMaxCharge2ndBlock'] = mTEPES.dPar['pMaxCharge2ndBlock'].where(mTEPES.dPar['pMaxCharge2ndBlock'] > mTEPES.dPar['pEpsilon'], 0.0)
 
@@ -1316,7 +1316,7 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
         mTEPES.dPar['pHeatPipeUpInvest']           = mTEPES.dPar['pHeatPipeUpInvest'].loc          [mTEPES.hc]
 
     # replace very small costs by 0
-    mTEPES.dPar['pEpsilon'] = 1e-5           # this value in EUR/GWh is lower than the O&M variable cost of any technology, independent of the area
+    mTEPES.dPar['pEpsilon'] = 2e-5           # this value in EUR/GWh is lower than the O&M variable cost of any technology, independent of the area
 
     mTEPES.dPar['pLinearVarCost']  [mTEPES.dPar['pLinearVarCost']  [[g for g in mTEPES.g]] < mTEPES.dPar['pEpsilon']] = 0.0
     mTEPES.dPar['pConstantVarCost'][mTEPES.dPar['pConstantVarCost'][[g for g in mTEPES.g]] < mTEPES.dPar['pEpsilon']] = 0.0

openTEPES/openTEPES_Main.py

@@ -660,7 +660,7 @@
 # For more information on this, and how to apply and follow the GNU AGPL, see
 # <https://www.gnu.org/licenses/>.
 
-# Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - April 10, 2026
+# Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - April 11, 2026
 # simplicity and transparency in power systems planning
 
 # Developed by
@@ -693,7 +693,7 @@
 BLUE   = "\033[34m"
 RESET  = "\033[0m"
 
-print(GREEN + 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.16 - April 10, 2026' + RESET)
+print(GREEN + 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - April 11, 2026' + RESET)
 print(BLUE  + '#### Academic research license - for non-commercial use only ####' + RESET + '\n')
 
 parser = argparse.ArgumentParser(description='Introducing main parameters...')

openTEPES/openTEPES_ModelFormulation.py

@@ -1,5 +1,5 @@
 """
-Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - April 05, 2026
+Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - April 11, 2026
 """
 
 import time
@@ -169,8 +169,8 @@ def eTotalGCost(OptModel,n):
                                                 mTEPES.pLoadLevelDuration[p,sc,n]() * sum(mTEPES.pLinearOMCost   [       re] * OptModel.vTotalOutput    [p,sc,n,re] for re in mTEPES.re if (p,re) in mTEPES.pre) )
     setattr(OptModel, f'eTotalGCost_{p}_{sc}_{st}', Constraint(mTEPES.n, rule=eTotalGCost, doc='system variable generation operation cost [MEUR]'))
 
-    # the small tolerance pEpsilon=1e-5 is added to avoid pumping/charging with curtailment/spillage
-    pEpsilon = 1e-5
+    # the small tolerance pEpsilon=2e-5 is added to avoid pumping/charging with curtailment/spillage
+    pEpsilon = 2e-5
     def eTotalCCost(OptModel,n):
         return OptModel.vTotalCCost    [p,sc,n] == mTEPES.pLoadLevelDuration[p,sc,n]() * sum((mTEPES.pLinearVarCost[p,sc,n,eh]+pEpsilon) * OptModel.vESSTotalCharge[p,sc,n,eh] for eh in mTEPES.eh if (p,eh) in mTEPES.peh and eh not in mTEPES.el)
     setattr(OptModel, f'eTotalCCost_{p}_{sc}_{st}', Constraint(mTEPES.n, rule=eTotalCCost, doc='system variable consumption operation cost [MEUR]'))
@@ -184,8 +184,8 @@ def eTotalECost(OptModel,n):
     def eTotalEmissionArea(OptModel,n,ar):
         if mTEPES.pEmission[p,ar] == math.inf or sum(mTEPES.pEmissionRate[g] for g in mTEPES.g if g in g2a[ar] and (p,g) in mTEPES.pg) == 0.0:
             return Constraint.Skip
-        return OptModel.vTotalEmissionArea[p,sc,n,ar] == (mTEPES.pLoadLevelDuration[p,sc,n]() * 1e-3 * (sum(mTEPES.pEmissionRate[nr] * OptModel.vTotalOutput    [p,sc,n,nr] for nr in mTEPES.nr if nr in g2a[ar] and (p,nr) in mTEPES.pnr)    #1e-3 to change from tCO2/MWh to MtCO2/GWh
-                                                                                                     +  sum(mTEPES.pEmissionRate[bo] * OptModel.vTotalOutputHeat[p,sc,n,bo] for bo in mTEPES.bo if bo in g2a[ar] and (p,bo) in mTEPES.pbo)))  #1e-3 to change from tCO2/MWh to MtCO2/GWh
+        return OptModel.vTotalEmissionArea[p,sc,n,ar] == (mTEPES.pLoadLevelDuration[p,sc,n]() * 1e-3 * (sum(mTEPES.pEmissionRate[nr] * OptModel.vTotalOutput    [p,sc,n,nr] for nr in mTEPES.nr if nr in g2a[ar] and (p,nr) in mTEPES.pnr)    # 1e-3 to change from tCO2/MWh to MtCO2/GWh
+                                                                                                     +  sum(mTEPES.pEmissionRate[bo] * OptModel.vTotalOutputHeat[p,sc,n,bo] for bo in mTEPES.bo if bo in g2a[ar] and (p,bo) in mTEPES.pbo)))  # 1e-3 to change from tCO2/MWh to MtCO2/GWh
     setattr(OptModel, f'eTotalEmissionArea_{p}_{sc}_{st}', Constraint(mTEPES.n*mTEPES.ar, rule=eTotalEmissionArea, doc='area total emission [MtCO2 eq]'))
 
     def eTotalECostArea(OptModel,n,ar):