diff --git a/doc/OnlineDocs/explanation/solvers/gdpopt.rst b/doc/OnlineDocs/explanation/solvers/gdpopt.rst index a6bee3c60fd..f11a6e3650b 100644 --- a/doc/OnlineDocs/explanation/solvers/gdpopt.rst +++ b/doc/OnlineDocs/explanation/solvers/gdpopt.rst @@ -12,11 +12,12 @@ The main advantage of these techniques is their ability to solve subproblems in a reduced space, including nonlinear constraints only for ``True`` logical blocks. As a result, GDPopt is most effective for nonlinear GDP models. -Three algorithms are available in GDPopt: +Four algorithms are available in GDPopt: 1. Logic-based outer approximation (LOA) [`Turkay & Grossmann, 1996`_] 2. Global logic-based outer approximation (GLOA) [`Lee & Grossmann, 2001`_] 3. Logic-based branch-and-bound (LBB) [`Lee & Grossmann, 2001`_] +4. Logic-based discrete steepest descent algorithm (LD-SDA) [`Ovalle et al., 2025`_] Usage and implementation details for GDPopt can be found in the PSE 2018 paper (`Chen et al., 2018`_), or via its @@ -28,6 +29,7 @@ Credit for prototyping and development can be found in the ``GDPopt`` class docu .. _Lee & Grossmann, 2001: https://doi.org/10.1016/S0098-1354(01)00732-3 .. _Lee & Grossmann, 2000: https://doi.org/10.1016/S0098-1354(00)00581-0 .. _Chen et al., 2018: https://doi.org/10.1016/B978-0-444-64241-7.50143-9 +.. _Ovalle et al., 2025: https://doi.org/10.1016/j.compchemeng.2024.108993 GDPopt can be used to solve a Pyomo.GDP concrete model in two ways. The simplest is to instantiate the generic GDPopt solver and specify the desired algorithm as an argument to the ``solve`` method: @@ -63,27 +65,27 @@ An example that includes the modeling approach may be found below. :skipif: not glpk_available Required imports - >>> from pyomo.environ import * - >>> from pyomo.gdp import * + >>> import pyomo.environ as pyo + >>> from pyomo.gdp import Disjunct, Disjunction Create a simple model - >>> model = ConcreteModel(name='LOA example') + >>> model = pyo.ConcreteModel(name='LOA example') - >>> model.x = Var(bounds=(-1.2, 2)) - >>> model.y = Var(bounds=(-10,10)) - >>> model.c = Constraint(expr= model.x + model.y == 1) + >>> model.x = pyo.Var(bounds=(-1.2, 2)) + >>> model.y = pyo.Var(bounds=(-10,10)) + >>> model.c = pyo.Constraint(expr= model.x + model.y == 1) >>> model.fix_x = Disjunct() - >>> model.fix_x.c = Constraint(expr=model.x == 0) + >>> model.fix_x.c = pyo.Constraint(expr=model.x == 0) >>> model.fix_y = Disjunct() - >>> model.fix_y.c = Constraint(expr=model.y == 0) + >>> model.fix_y.c = pyo.Constraint(expr=model.y == 0) >>> model.d = Disjunction(expr=[model.fix_x, model.fix_y]) - >>> model.objective = Objective(expr=model.x + 0.1*model.y, sense=minimize) + >>> model.objective = pyo.Objective(expr=model.x + 0.1*model.y, sense=pyo.minimize) Solve the model using GDPopt - >>> results = SolverFactory('gdpopt.loa').solve( + >>> results = pyo.SolverFactory('gdpopt.loa').solve( ... model, mip_solver='glpk') # doctest: +IGNORE_RESULT Display the final solution @@ -158,25 +160,25 @@ To use the GDPopt-LBB solver, define your Pyomo GDP model as usual: :skipif: not baron_available Required imports - >>> from pyomo.environ import * + >>> import pyomo.environ as pyo >>> from pyomo.gdp import Disjunct, Disjunction Create a simple model - >>> m = ConcreteModel() - >>> m.x1 = Var(bounds = (0,8)) - >>> m.x2 = Var(bounds = (0,8)) - >>> m.obj = Objective(expr=m.x1 + m.x2, sense=minimize) + >>> m = pyo.ConcreteModel() + >>> m.x1 = pyo.Var(bounds = (0,8)) + >>> m.x2 = pyo.Var(bounds = (0,8)) + >>> m.obj = pyo.Objective(expr=m.x1 + m.x2, sense=pyo.minimize) >>> m.y1 = Disjunct() >>> m.y2 = Disjunct() - >>> m.y1.c1 = Constraint(expr=m.x1 >= 2) - >>> m.y1.c2 = Constraint(expr=m.x2 >= 2) - >>> m.y2.c1 = Constraint(expr=m.x1 >= 3) - >>> m.y2.c2 = Constraint(expr=m.x2 >= 3) + >>> m.y1.c1 = pyo.Constraint(expr=m.x1 >= 2) + >>> m.y1.c2 = pyo.Constraint(expr=m.x2 >= 2) + >>> m.y2.c1 = pyo.Constraint(expr=m.x1 >= 3) + >>> m.y2.c2 = pyo.Constraint(expr=m.x2 >= 3) >>> m.djn = Disjunction(expr=[m.y1, m.y2]) Invoke the GDPopt-LBB solver - >>> results = SolverFactory('gdpopt.lbb').solve(m) + >>> results = pyo.SolverFactory('gdpopt.lbb').solve(m) WARNING: 09/06/22: The GDPopt LBB algorithm currently has known issues. Please use the results with caution and report any bugs! @@ -186,9 +188,70 @@ To use the GDPopt-LBB solver, define your Pyomo GDP model as usual: >>> print(results.solver.termination_condition) optimal - >>> print([value(m.y1.indicator_var), value(m.y2.indicator_var)]) + >>> print([pyo.value(m.y1.indicator_var), pyo.value(m.y2.indicator_var)]) [True, False] +Logic-based Discrete-Steepest Descent Algorithm (LD-SDA) +-------------------------------------------------------- + +The GDPopt-LDSDA solver exploits the ordered Boolean variables in the disjunctions to solve the GDP model. +It requires an **exclusive OR (XOR) logical constraint** to ensure that exactly one disjunct is active in each disjunction. +The solver also requires a **starting point** for the discrete variables and allows users to choose between two **direction norms**, `'L2'` and `'Linf'`, to guide the search process. + +.. note:: + + The current implementation of the GDPopt-LDSDA requires an explicit LogicalConstraint to enforce the exclusive OR condition for each disjunction. + +To use the GDPopt-LDSDA solver, define your Pyomo GDP model as usual: + +.. doctest:: + :skipif: not baron_available + + Required imports + >>> import pyomo.environ as pyo + >>> from pyomo.gdp import Disjunct, Disjunction + + Create a simple model + >>> m = pyo.ConcreteModel() + + Define sets + >>> I = [1, 2, 3, 4, 5] + >>> J = [1, 2, 3, 4, 5] + + Define variables + >>> m.a = pyo.Var(bounds=(-0.3, 0.2)) + >>> m.b = pyo.Var(bounds=(-0.9, -0.5)) + + Define disjuncts for Y1 + >>> m.Y1_disjuncts = Disjunct(I) + >>> for i in I: + ... m.Y1_disjuncts[i].y1_constraint = pyo.Constraint(expr=m.a == -0.3 + 0.1 * (i - 1)) + + Define disjuncts for Y2 + >>> m.Y2_disjuncts = Disjunct(J) + >>> for j in J: + ... m.Y2_disjuncts[j].y2_constraint = pyo.Constraint(expr=m.b == -0.9 + 0.1 * (j - 1)) + + Define disjunctions + >>> m.y1_disjunction = Disjunction(expr=[m.Y1_disjuncts[i] for i in I]) + >>> m.y2_disjunction = Disjunction(expr=[m.Y2_disjuncts[j] for j in J]) + + Logical constraints to enforce exactly one selection + >>> m.Y1_limit = pyo.LogicalConstraint(expr=pyo.exactly(1, [m.Y1_disjuncts[i].indicator_var for i in I])) + >>> m.Y2_limit = pyo.LogicalConstraint(expr=pyo.exactly(1, [m.Y2_disjuncts[j].indicator_var for j in J])) + + Define objective function + >>> m.obj = pyo.Objective( + ... expr=4 * m.a**2 - 2.1 * m.a**4 + (1 / 3) * m.a**6 + m.a * m.b - 4 * m.b**2 + 4 * m.b**4, + ... sense=pyo.minimize + ... ) + + Invoke the GDPopt-LDSDA solver + >>> results = pyo.SolverFactory('gdpopt.ldsda').solve(m, + ... starting_point=[1,1], + ... logical_constraint_list=[m.Y1_limit, m.Y2_limit], + ... direction_norm='Linf', + ... ) GDPopt implementation and optional arguments -------------------------------------------- @@ -204,4 +267,5 @@ GDPopt implementation and optional arguments ~pyomo.contrib.gdpopt.gloa.GDP_GLOA_Solver ~pyomo.contrib.gdpopt.ric.GDP_RIC_Solver ~pyomo.contrib.gdpopt.branch_and_bound.GDP_LBB_Solver + ~pyomo.contrib.gdpopt.ldsda.GDP_LDSDA_Solver