Revised vapor composition constraints fix ZeroDivisionError when gdp.hull reformulation is applied

benchmark.py

@@ -1,7 +1,7 @@
 import os
 import json
 import time
-import sys
+from contextlib import redirect_stdout
 from datetime import datetime
 from importlib import import_module
 from pyomo.environ import *
@@ -28,14 +28,13 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"):
     None
     """
     model = model.clone()
-    stdout = sys.stdout
+    log_path = os.path.join(result_dir, f"{strategy}_{subsolver}.log")
+    json_path = os.path.join(result_dir, f"{strategy}_{subsolver}.json")
     if strategy in ["gdp.bigm", "gdp.hull"]:
         transformation_start_time = time.time()
         TransformationFactory(strategy).apply_to(model)
         transformation_end_time = time.time()
-        with open(
-            result_dir + "/" + strategy + "_" + subsolver + ".log", "w"
-        ) as sys.stdout:
+        with open(log_path, "w") as log_file, redirect_stdout(log_file):
             results = SolverFactory(subsolver).solve(
                 model, tee=True, timelimit=timelimit
             )
@@ -50,9 +49,7 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"):
         "gdpopt.lbb",
         "gdpopt.ric",
     ]:
-        with open(
-            result_dir + "/" + strategy + "_" + subsolver + ".log", "w"
-        ) as sys.stdout:
+        with open(log_path, "w") as log_file, redirect_stdout(log_file):
             results = SolverFactory(strategy).solve(
                 model,
                 tee=True,
@@ -63,9 +60,10 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"):
                 time_limit=timelimit,
             )
             print(results)
+    else:
+        raise ValueError(f"Unknown benchmark strategy: {strategy}")
 
-    sys.stdout = stdout
-    with open(result_dir + "/" + strategy + "_" + subsolver + ".json", "w") as f:
+    with open(json_path, "w") as f:
         json.dump(results.json_repn(), f)
     return None
 

gdplib/biofuel/model.py

@@ -4,7 +4,7 @@
 
 The model enforces constraints to ensure that raw material supplies do not exceed available amounts, product shipments meet market demands exactly, and production at each site matches outgoing shipments and available resources.
 It also optimizes transportation costs by managing both variable and fixed costs associated with active transportation routes.
-The disjunctions in the model define the operational modes for facility sites (modular, conventional, or inactive) and the activity status of supply and product routes (active or inactive). 
+The disjunctions in the model define the operational modes for facility sites (modular, conventional, or inactive) and the activity status of supply and product routes (active or inactive).
 These elements allow the model to simulate different operational scenarios and strategic decisions, optimizing the network's layout and logistics based on economic and market conditions.
 The objective of the model is to optimize the network layout and production allocation to minimize total costs, which include setup and teardown of facilities, production costs, and transportation costs.
 

gdplib/gdp_col/initialize.py

@@ -72,9 +72,9 @@ def set_value_if_not_fixed(var, val):
             for i in range(1, num_active_trays)
         ]
         for tray in range(2, num_active_trays):
-            indx = new_indices[tray - 1]
-            lower = floor(indx)
-            frac_above = indx - lower
+            index = new_indices[tray - 1]
+            lower = floor(index)
+            frac_above = index - lower
             # Take linear combination of values
             tray_indexed_data.loc[tray] = (
                 tray_indexed_data.loc[lower] * (1 - frac_above)

gdplib/hda/HDA_GDP_gdpopt.py

@@ -20,7 +20,7 @@ def HDA_model():
 
     m.alpha = Param(initialize=0.3665, doc="compressor coefficient")
     m.compeff = Param(initialize=0.750, doc="compressor efficiency")
-    m.gam = Param(initialize=1.300, doc="ratio of cp to cv")
+    m.heat_capacity_ratio = Param(initialize=1.300, doc="ratio of cp to cv")
     m.abseff = Param(initialize=0.333, doc="absorber tray efficiency")
     m.disteff = Param(initialize=0.5000, doc="column tray efficiency")
     m.uflow = Param(initialize=50, doc="upper bound - flow logicals")
@@ -1312,7 +1312,7 @@ def Compelec(_m, comp_):
                     * m.f[stream]
                     / 60.0
                     * (1.0 / m.compeff)
-                    * (m.gam / (m.gam - 1.0))
+                    * (m.heat_capacity_ratio / (m.heat_capacity_ratio - 1.0))
                     for (comp1, stream) in m.icomp
                     if comp_ == comp1
                 )
@@ -1324,9 +1324,13 @@ def Compelec(_m, comp_):
 
         def Ratio(_m, comp_):
             if comp == comp_:
-                return m.presrat[comp_] ** (m.gam / (m.gam - 1.0)) == sum(
+                return m.presrat[comp_] ** (
+                    m.heat_capacity_ratio / (m.heat_capacity_ratio - 1.0)
+                ) == sum(
                     m.p[stream] for (comp1, stream) in m.ocomp if comp_ == comp1
-                ) / sum(m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_)
+                ) / sum(
+                    m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_
+                )
             return Constraint.Skip
 
         b.ratio = Constraint([comp], rule=Ratio, doc='pressure ratio (out to in)')
@@ -1416,7 +1420,7 @@ def Relvol(_m, dist_):
                 return m.avevlt[dist] == sqrt(divided1 * divided2)
             return Constraint.Skip
 
-        b.relvol = Constraint([dist], rule=Relvol, doc='average relative volatilty')
+        b.relvol = Constraint([dist], rule=Relvol, doc='average relative volatility')
 
         def Undwood(_m, dist_):
             if dist_ == dist:
@@ -2248,11 +2252,19 @@ def Valcmb(_m, valve, compon):
 
         def Valt(_m, valve):
             return sum(
-                m.t[stream] / (m.p[stream] ** ((m.gam - 1.0) / m.gam))
+                m.t[stream]
+                / (
+                    m.p[stream]
+                    ** ((m.heat_capacity_ratio - 1.0) / m.heat_capacity_ratio)
+                )
                 for (valv, stream) in m.oval
                 if valv == valve
             ) == sum(
-                m.t[stream] / (m.p[stream] ** ((m.gam - 1.0) / m.gam))
+                m.t[stream]
+                / (
+                    m.p[stream]
+                    ** ((m.heat_capacity_ratio - 1.0) / m.heat_capacity_ratio)
+                )
                 for (valv, stream) in m.ival
                 if valv == valve
             )

gdplib/kaibel/kaibel_init.py

@@ -1,28 +1,28 @@
 """
-        Calculation of the theoretical minimum number of trays and initial
-                             temperature values.
-             (written by E. Soraya Rawlings, esoraya@rwlngs.net)
-
- The separation of four components require a sequence of at least three distillation
- columns. Here, we calculate the minimum number of theoretical trays for the three
- columns. The sequence is shown in Figure 2.
-
-          COLUMN 1              COLUMN 2              COLUMN 3
-              -----                 ----                 -----
-             |     |               |    |               |     |
-           -----   |   A        -----   |             -----   |
-          |     |<---> B --    |     |<----> A --    |     |<---> A
-          |     |      C   |   |     |       B   |   |     |
-    A     |     |          |   |     |           |   |     |
-    B     |     |          |   |     |           |   |     |
-    C --->|     |           -->|     |            -->|     |
-    D     |     |              |     |               |     |
-          |     |              |     |               |     |
-          |     |<-            |     |<-             |     |<-
-           -----   |            -----   |             -----   |
-             |     |              |     |               |     |
-              -------> D           -------> C            -------> B
-   Figure 2. Sequence of columns for the separation of a quaternary mixture
+       Calculation of the theoretical minimum number of trays and initial
+                            temperature values.
+            (written by E. Soraya Rawlings, esoraya@rwlngs.net)
+
+The separation of four components require a sequence of at least three distillation
+columns. Here, we calculate the minimum number of theoretical trays for the three
+columns. The sequence is shown in Figure 2.
+
+         COLUMN 1              COLUMN 2              COLUMN 3
+             -----                 ----                 -----
+            |     |               |    |               |     |
+          -----   |   A        -----   |             -----   |
+         |     |<---> B --    |     |<----> A --    |     |<---> A
+         |     |      C   |   |     |       B   |   |     |
+   A     |     |          |   |     |           |   |     |
+   B     |     |          |   |     |           |   |     |
+   C --->|     |           -->|     |            -->|     |
+   D     |     |              |     |               |     |
+         |     |              |     |               |     |
+         |     |<-            |     |<-             |     |<-
+          -----   |            -----   |             -----   |
+            |     |              |     |               |     |
+             -------> D           -------> C            -------> B
+  Figure 2. Sequence of columns for the separation of a quaternary mixture
 """
 
 from __future__ import division

gdplib/kaibel/kaibel_prop.py

@@ -1,4 +1,4 @@
-""" Properties of the system """
+"""Properties of the system"""
 
 from pyomo.environ import ConcreteModel
 

gdplib/kaibel/kaibel_side_flash.py

@@ -1,4 +1,4 @@
-""" Side feed flash """
+"""Side feed flash"""
 
 from pyomo.environ import (
     ConcreteModel,