Merge branch 'main' into water_network

LICENSE

@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2020, Ignacio Grossmann Research Group
+Copyright (c) 2024, SECQUOIA Research Group
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

benchmark.py

@@ -79,7 +79,8 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"):
         # "hda",
         "jobshop",
         # "kaibel",
-        # "logical",
+        # "positioning",
+        # "spectralog",
         # "med_term_purchasing",
         # "methanol",
         # "mod_hens",

gdplib/__init__.py

@@ -1,7 +1,8 @@
 import gdplib.mod_hens
 import gdplib.modprodnet
 import gdplib.biofuel
-import gdplib.logical  # Requires logical expression system
+import gdplib.positioning
+import gdplib.spectralog
 import gdplib.stranded_gas  # Requires logical expression system
 import gdplib.gdp_col
 import gdplib.hda
@@ -13,3 +14,6 @@
 import gdplib.med_term_purchasing
 import gdplib.syngas
 import gdplib.water_network
+import gdplib.ex1_linan_2023
+import gdplib.small_batch
+import gdplib.cstr

gdplib/batch_processing/batch_processing.py

@@ -69,9 +69,9 @@ def build_model():
 
     # Sets
 
-    model.PRODUCTS = Set(doc='Set of Products')
-    model.STAGES = Set(doc='Set of Stages', ordered=True)
-    model.PARALLELUNITS = Set(doc='Set of Parallel Units', ordered=True)
+    model.PRODUCTS = Set(doc="Set of Products")
+    model.STAGES = Set(doc="Set of Stages", ordered=True)
+    model.PARALLELUNITS = Set(doc="Set of Parallel Units", ordered=True)
 
     # TODO: this seems like an over-complicated way to accomplish this task...
     def filter_out_last(model, j):
@@ -101,27 +101,27 @@ def filter_out_last(model, j):
 
     # Parameters
 
-    model.HorizonTime = Param(doc='Horizon Time')
-    model.Alpha1 = Param(doc='Cost Parameter of the units')
-    model.Alpha2 = Param(doc='Cost Parameter of the intermediate storage tanks')
-    model.Beta1 = Param(doc='Exponent Parameter of the units')
-    model.Beta2 = Param(doc='Exponent Parameter of the intermediate storage tanks')
+    model.HorizonTime = Param(doc="Horizon Time")
+    model.Alpha1 = Param(doc="Cost Parameter of the units")
+    model.Alpha2 = Param(doc="Cost Parameter of the intermediate storage tanks")
+    model.Beta1 = Param(doc="Exponent Parameter of the units")
+    model.Beta2 = Param(doc="Exponent Parameter of the intermediate storage tanks")
 
-    model.ProductionAmount = Param(model.PRODUCTS, doc='Production Amount')
+    model.ProductionAmount = Param(model.PRODUCTS, doc="Production Amount")
     model.ProductSizeFactor = Param(
-        model.PRODUCTS, model.STAGES, doc='Product Size Factor'
+        model.PRODUCTS, model.STAGES, doc="Product Size Factor"
     )
-    model.ProcessingTime = Param(model.PRODUCTS, model.STAGES, doc='Processing Time')
+    model.ProcessingTime = Param(model.PRODUCTS, model.STAGES, doc="Processing Time")
 
     # These are hard-coded in the GAMS file, hence the defaults
     model.StorageTankSizeFactor = Param(
-        model.STAGES, default=StorageTankSizeFactor, doc='Storage Tank Size Factor'
+        model.STAGES, default=StorageTankSizeFactor, doc="Storage Tank Size Factor"
     )
     model.StorageTankSizeFactorByProd = Param(
         model.PRODUCTS,
         model.STAGES,
         default=StorageTankSizeFactorByProd,
-        doc='Storage Tank Size Factor by Product',
+        doc="Storage Tank Size Factor by Product",
     )
 
     # TODO: bonmin wasn't happy and I think it might have something to do with this?
@@ -148,27 +148,27 @@ def get_log_coeffs(model, k):
         return log(model.PARALLELUNITS.ord(k))
 
     model.LogCoeffs = Param(
-        model.PARALLELUNITS, initialize=get_log_coeffs, doc='Logarithmic Coefficients'
+        model.PARALLELUNITS, initialize=get_log_coeffs, doc="Logarithmic Coefficients"
     )
 
     # bounds
     model.volumeLB = Param(
-        model.STAGES, default=VolumeLB, doc='Lower Bound of Volume of the Units'
+        model.STAGES, default=VolumeLB, doc="Lower Bound of Volume of the Units"
     )
     model.volumeUB = Param(
-        model.STAGES, default=VolumeUB, doc='Upper Bound of Volume of the Units'
+        model.STAGES, default=VolumeUB, doc="Upper Bound of Volume of the Units"
     )
     model.storageTankSizeLB = Param(
-        model.STAGES, default=StorageTankSizeLB, doc='Lower Bound of Storage Tank Size'
+        model.STAGES, default=StorageTankSizeLB, doc="Lower Bound of Storage Tank Size"
     )
     model.storageTankSizeUB = Param(
-        model.STAGES, default=StorageTankSizeUB, doc='Upper Bound of Storage Tank Size'
+        model.STAGES, default=StorageTankSizeUB, doc="Upper Bound of Storage Tank Size"
     )
     model.unitsInPhaseUB = Param(
-        model.STAGES, default=UnitsInPhaseUB, doc='Upper Bound of Units in Phase'
+        model.STAGES, default=UnitsInPhaseUB, doc="Upper Bound of Units in Phase"
     )
     model.unitsOutOfPhaseUB = Param(
-        model.STAGES, default=UnitsOutOfPhaseUB, doc='Upper Bound of Units Out of Phase'
+        model.STAGES, default=UnitsOutOfPhaseUB, doc="Upper Bound of Units Out of Phase"
     )
 
     # Variables
@@ -214,13 +214,16 @@ def get_volume_bounds(model, j):
         return (model.volumeLB[j], model.volumeUB[j])
 
     model.volume_log = Var(
-        model.STAGES, bounds=get_volume_bounds, doc='Logarithmic Volume of the Units'
+        model.STAGES, bounds=get_volume_bounds, doc="Logarithmic Volume of the Units"
     )
     model.batchSize_log = Var(
-        model.PRODUCTS, model.STAGES, doc='Logarithmic Batch Size of the Products'
+        model.PRODUCTS,
+        model.STAGES,
+        bounds=(0, 10),
+        doc="Logarithmic Batch Size of the Products",
     )
     model.cycleTime_log = Var(
-        model.PRODUCTS, doc='Logarithmic Cycle Time of the Products'
+        model.PRODUCTS, doc="Logarithmic Cycle Time of the Products"
     )
 
     def get_unitsOutOfPhase_bounds(model, j):
@@ -244,7 +247,7 @@ def get_unitsOutOfPhase_bounds(model, j):
     model.unitsOutOfPhase_log = Var(
         model.STAGES,
         bounds=get_unitsOutOfPhase_bounds,
-        doc='Logarithmic Units Out of Phase',
+        doc="Logarithmic Units Out of Phase",
     )
 
     def get_unitsInPhase_bounds(model, j):
@@ -266,7 +269,7 @@ def get_unitsInPhase_bounds(model, j):
         return (0, model.unitsInPhaseUB[j])
 
     model.unitsInPhase_log = Var(
-        model.STAGES, bounds=get_unitsInPhase_bounds, doc='Logarithmic Units In Phase'
+        model.STAGES, bounds=get_unitsInPhase_bounds, doc="Logarithmic Units In Phase"
     )
 
     def get_storageTankSize_bounds(model, j):
@@ -291,15 +294,15 @@ def get_storageTankSize_bounds(model, j):
     model.storageTankSize_log = Var(
         model.STAGES,
         bounds=get_storageTankSize_bounds,
-        doc='Logarithmic Storage Tank Size',
+        doc="Logarithmic Storage Tank Size",
     )
 
     # binary variables for deciding number of parallel units in and out of phase
     model.outOfPhase = Var(
-        model.STAGES, model.PARALLELUNITS, within=Binary, doc='Out of Phase Units'
+        model.STAGES, model.PARALLELUNITS, within=Binary, doc="Out of Phase Units"
     )
     model.inPhase = Var(
-        model.STAGES, model.PARALLELUNITS, within=Binary, doc='In Phase Units'
+        model.STAGES, model.PARALLELUNITS, within=Binary, doc="In Phase Units"
     )
 
     # Objective
@@ -336,7 +339,7 @@ def get_cost_rule(model):
         )
 
     model.min_cost = Objective(
-        rule=get_cost_rule, doc='Minimize the Total Cost of the Plant Design'
+        rule=get_cost_rule, doc="Minimize the Total Cost of the Plant Design"
     )
 
     # Constraints
@@ -371,7 +374,7 @@ def processing_capacity_rule(model, j, i):
         model.STAGES,
         model.PRODUCTS,
         rule=processing_capacity_rule,
-        doc='Processing Capacity',
+        doc="Processing Capacity",
     )
 
     def processing_time_rule(model, j, i):
@@ -402,7 +405,7 @@ def processing_time_rule(model, j, i):
         )
 
     model.processing_time = Constraint(
-        model.STAGES, model.PRODUCTS, rule=processing_time_rule, doc='Processing Time'
+        model.STAGES, model.PRODUCTS, rule=processing_time_rule, doc="Processing Time"
     )
 
     def finish_in_time_rule(model):
@@ -424,7 +427,7 @@ def finish_in_time_rule(model):
             for i in model.PRODUCTS
         )
 
-    model.finish_in_time = Constraint(rule=finish_in_time_rule, doc='Finish in Time')
+    model.finish_in_time = Constraint(rule=finish_in_time_rule, doc="Finish in Time")
 
     # Disjunctions
 
@@ -554,7 +557,7 @@ def no_batch_rule(disjunct, i):
         [0, 1],
         model.STAGESExceptLast,
         rule=storage_tank_selection_disjunct_rule,
-        doc='Storage Tank Selection Disjunct',
+        doc="Storage Tank Selection Disjunct",
     )
 
     def select_storage_tanks_rule(model, j):
@@ -581,7 +584,7 @@ def select_storage_tanks_rule(model, j):
     model.select_storage_tanks = Disjunction(
         model.STAGESExceptLast,
         rule=select_storage_tanks_rule,
-        doc='Select Storage Tanks',
+        doc="Select Storage Tanks",
     )
 
     # though this is a disjunction in the GAMs model, it is more efficiently formulated this way:
@@ -612,7 +615,7 @@ def units_out_of_phase_rule(model, j):
         )
 
     model.units_out_of_phase = Constraint(
-        model.STAGES, rule=units_out_of_phase_rule, doc='Units Out of Phase'
+        model.STAGES, rule=units_out_of_phase_rule, doc="Units Out of Phase"
     )
 
     def units_in_phase_rule(model, j):
@@ -641,7 +644,7 @@ def units_in_phase_rule(model, j):
         )
 
     model.units_in_phase = Constraint(
-        model.STAGES, rule=units_in_phase_rule, doc='Units In Phase'
+        model.STAGES, rule=units_in_phase_rule, doc="Units In Phase"
     )
 
     def units_out_of_phase_xor_rule(model, j):
@@ -665,7 +668,7 @@ def units_out_of_phase_xor_rule(model, j):
     model.units_out_of_phase_xor = Constraint(
         model.STAGES,
         rule=units_out_of_phase_xor_rule,
-        doc='Exclusive OR for Units Out of Phase',
+        doc="Exclusive OR for Units Out of Phase",
     )
 
     def units_in_phase_xor_rule(model, j):
@@ -689,16 +692,16 @@ def units_in_phase_xor_rule(model, j):
     model.units_in_phase_xor = Constraint(
         model.STAGES,
         rule=units_in_phase_xor_rule,
-        doc='Exclusive OR for Units In Phase',
+        doc="Exclusive OR for Units In Phase",
     )
 
-    return model.create_instance(join(this_file_dir(), 'batch_processing.dat'))
+    return model.create_instance(join(this_file_dir(), "batch_processing.dat"))
 
 
 if __name__ == "__main__":
     m = build_model()
-    TransformationFactory('gdp.bigm').apply_to(m)
-    SolverFactory('gams').solve(
-        m, solver='baron', tee=True, add_options=['option optcr=1e-6;']
+    TransformationFactory("gdp.bigm").apply_to(m)
+    SolverFactory("gams").solve(
+        m, solver="baron", tee=True, add_options=["option optcr=1e-6;"]
     )
     m.min_cost.display()

gdplib/biofuel/model.py

@@ -52,7 +52,8 @@ def build_model():
     [2] Chen, Q., & Grossmann, I. E. (2019). Effective generalized disjunctive programming models for modular process synthesis. Industrial & Engineering Chemistry Research, 58(15), 5873-5886. https://doi.org/10.1021/acs.iecr.8b04600
     """
     m = ConcreteModel('Biofuel processing network')
-    m.bigM = Suffix(direction=Suffix.LOCAL, initialize=7000)
+    m.bigM = Suffix(direction=Suffix.LOCAL)
+    m.bigM[None] = 7000
     m.time = RangeSet(0, 120, doc="months in 10 years")
     m.suppliers = RangeSet(10)  # 10 suppliers
     m.markets = RangeSet(10)  # 10 markets

gdplib/cstr/__init__.py

@@ -0,0 +1,3 @@
+from .gdp_reactor import build_model
+
+__all__ = ['build_model']

gdplib/cstr/gdp_reactor.py

@@ -17,7 +17,7 @@
 from pyomo.opt.base.solvers import SolverFactory
 
 
-def build_cstrs(NT: int = 5) -> pyo.ConcreteModel():
+def build_model(NT: int = 5) -> pyo.ConcreteModel():
     """
     Build the CSTR superstructure model of size NT.
     NT is the number of reactors in series.
@@ -915,7 +915,7 @@ def obj_rule(m):
 
 
 if __name__ == "__main__":
-    m = build_cstrs()
+    m = build_model()
     pyo.TransformationFactory("core.logical_to_linear").apply_to(m)
     pyo.TransformationFactory("gdp.bigm").apply_to(m)
     pyo.SolverFactory("gams").solve(

gdplib/ex1_linan_2023/README.md

@@ -0,0 +1,18 @@
+## Example 1 Problem of Liñán (2023)
+
+The Example 1 Problem of Liñán (2023) is a simple optimization problem that involves two Boolean variables, two continuous variables, and a nonlinear objective function. The problem is formulated as a Generalized Disjunctive Programming (GDP) model.
+
+The Boolean variables are associated with disjuncts that define the feasible regions of the continuous variables. The problem also includes logical constraints that ensure that only one Boolean variable is true at a time. Additionally, there are two disjunctions, one for each Boolean variable, where only one disjunct in each disjunction must be true. A specific logical constraint also enforces that `Y1[3]` must be false, making this particular disjunct infeasible.
+
+The objective function is `-0.9995999999999999` when the continuous variables are alpha = 0 (`Y1[2]=True `) and beta=-0.7 (`Y2[3]=True`).
+
+The objective function originates from Problem No. 6 of Gomez's paper, and Liñán introduced logical propositions, logical disjunctions, and the following equations as constraints.
+
+### References
+
+> [1] Liñán, D. A., & Ricardez-Sandoval, L. A. (2023). A Benders decomposition framework for the optimization of disjunctive superstructures with ordered discrete decisions. AIChE Journal, 69(5), e18008. https://doi.org/10.1002/aic.18008
+>
+> [2] Gomez, S., & Levy, A. V. (1982). The tunnelling method for solving the constrained global optimization problem with several non-connected feasible regions. In Numerical Analysis: Proceedings of the Third IIMAS Workshop Held at Cocoyoc, Mexico, January 1981 (pp. 34-47). Springer Berlin Heidelberg. https://doi.org/10.1007/BFb0092958
+
+
+---

gdplib/ex1_linan_2023/__init__.py

@@ -0,0 +1,3 @@
+from .ex1_linan_2023 import build_model
+
+__all__ = ['build_model']