Add Logic-Based Discrete-Steepest Descent Algorithm in GDPOpt

[ZedongPeng]

Summary/Motivation:

This PR introduces the implementation of the Logic-Based Discrete Steepest Descent algorithm in GDPOpt.

The Logic-based Discrete-Steepest Descent Algorithm (LD-SDA) is a solution method for GDP problems involving ordered Boolean variables. The LD-SDA reformulates these ordered Boolean variables into integer decisions called external variables. The LD-SDA solves the reformulated GDP problem using a two-level decomposition approach where the upper-level subproblem determines external variable configurations. Subsequently, the remaining continuous and discrete variables are solved as a subproblem only involving those constraints relevant to the given external variable arrangement, effectively taking advantage of the structure of the GDP problem.

More details in the paper https://arxiv.org/abs/2405.05358 .

@emma58 @bernalde

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the BSD license.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[codecov]

Codecov Report

Attention: Patch coverage is 26.90355% with 144 lines in your changes missing coverage. Please review.

Project coverage is 88.37%. Comparing base (1ada528) to head (1144dc1).
Report is 746 commits behind head on main.

Files	Patch %	Lines
pyomo/contrib/gdpopt/ldsda.py	23.40%	144 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3331      +/-   ##
==========================================
- Coverage   88.50%   88.37%   -0.13%     
==========================================
  Files         868      869       +1     
  Lines       98418    98614     +196     
==========================================
+ Hits        87100    87154      +54     
- Misses      11318    11460     +142     

Flag	Coverage Δ
linux	85.90% <26.90%> (-0.13%)	⬇️
osx	75.19% <26.90%> (-0.10%)	⬇️
other	86.40% <26.90%> (-0.12%)	⬇️
win	83.71% <26.90%> (-0.12%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[ZedongPeng]

I will add tests to increase the code coverage.

[ZedongPeng]

Hi @dovallev and @David-Linan,

This PR includes a general implementation of LDSDA. When you have some time, could you please review it? Thanks in advance!

[blnicho]

@ZedongPeng is this ready for review now?

[emma58]

@ZedongPeng could you please run black on this so we can see if other tests are passing?

[ZedongPeng]

Hi @blnicho and @jsiirola . Do you know how to resolve the following 'gams' not found issue?

==================================== ERRORS ====================================
__________ ERROR collecting pyomo/contrib/gdpopt/tests/test_ldsda.py ___________
pyomo/contrib/gdpopt/tests/test_ldsda.py:7: in <module>
    class TestGDPoptLDSDA(unittest.TestCase):
pyomo/contrib/gdpopt/tests/test_ldsda.py:11: in TestGDPoptLDSDA
    SolverFactory('gams').available() and SolverFactory('gams').license_is_valid(),
pyomo/solvers/plugins/solvers/GAMS.py:660: in available
    raise NameError(
E   NameError: No 'gams' command found on system PATH - GAMS shell solver functionality is not available.

[jsiirola]

Hi @blnicho and @jsiirola . Do you know how to resolve the following 'gams' not found issue?

Yup: use SolverFactory('gams').available(False) (That is passing exception_flag=False, which suppresses the exception. In the future, we would like to change the default behavior to NOT raise exceptions, but that is a complicated change / deprecation path.

[emma58]

Hi @ZedongPeng, it also looks like the Jenkins failures are real. Here's the stack trace:

self = <pyomo.contrib.gdpopt.tests.test_ldsda.TestGDPoptLDSDA testMethod=test_solve_four_stage_dynamic_model>

    @unittest.skipUnless(
        SolverFactory('gams').available() and SolverFactory('gams').license_is_valid(),
        "gams solver not available",
    )
    def test_solve_four_stage_dynamic_model(self):
    
        model = build_model(mode_transfer=True)
    
        # Discretize the model using dae.collocation
        discretizer = TransformationFactory('dae.collocation')
        discretizer.apply_to(model, nfe=10, ncp=3, scheme='LAGRANGE-RADAU')
        # We need to reconstruct the constraints in disjuncts after discretization.
        # This is a bug in Pyomo.dae. https://github.com/Pyomo/pyomo/issues/3101
        for disjunct in model.component_data_objects(ctype=Disjunct):
            for constraint in disjunct.component_objects(ctype=Constraint):
                constraint._constructed = False
                constraint.construct()
    
        for dxdt in model.component_data_objects(ctype=Var, descend_into=True):
            if 'dxdt' in dxdt.name:
                dxdt.setlb(-300)
                dxdt.setub(300)
    
        for direction_norm in ['L2', 'Linf']:
            result = SolverFactory('gdpopt.ldsda').solve(
                model,
                direction_norm=direction_norm,
                minlp_solver='gams',
                minlp_solver_args=dict(solver='knitro'),
                starting_point=[1, 2],
                logical_constraint_list=[
                    model.mode_transfer_lc1.name,
                    model.mode_transfer_lc2.name,
                ],
                time_limit=100,
            )
>           self.assertAlmostEqual(value(model.obj), -23.305325, places=4)

pyomo/pyomo/contrib/gdpopt/tests/test_ldsda.py:46: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
pyomo/pyomo/common/numeric_types.py:382: in value
    tmp = obj(exception=True)
pyomo/pyomo/core/base/objective.py:430: in __call__
    return super().__call__(exception)
pyomo/pyomo/core/base/expression.py:59: in __call__
    return arg(exception=exception)
pyomo/pyomo/core/expr/base.py:118: in __call__
    return visitor.evaluate_expression(self, exception)
pyomo/pyomo/core/expr/visitor.py:1301: in evaluate_expression
    ans = visitor.dfs_postorder_stack(exp)
pyomo/pyomo/core/expr/visitor.py:919: in dfs_postorder_stack
    flag, value = self.visiting_potential_leaf(_sub)
pyomo/pyomo/core/expr/visitor.py:1202: in visiting_potential_leaf
    return True, value(node, exception=self.exception)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

obj = <pyomo.core.base.var.VarData object at 0x7f6f349403c0>, exception = True

    def value(obj, exception=True):
        """
        A utility function that returns the value of a Pyomo object or
        expression.
    
        Args:
            obj: The argument to evaluate. If it is None, a
                string, or any other primitive numeric type,
                then this function simply returns the argument.
                Otherwise, if the argument is a NumericValue
                then the __call__ method is executed.
            exception (bool): If :const:`True`, then an exception should
                be raised when instances of NumericValue fail to
                evaluate due to one or more objects not being
                initialized to a numeric value (e.g, one or more
                variables in an algebraic expression having the
                value None). If :const:`False`, then the function
                returns :const:`None` when an exception occurs.
                Default is True.
    
        Returns: A numeric value or None.
        """
        if obj.__class__ in native_types:
            return obj
        #
        # Test if we have a duck typed Pyomo expression
        #
        if not hasattr(obj, 'is_numeric_type'):
            #
            # TODO: Historically we checked for new *numeric* types and
            # raised exceptions for anything else.  That is inconsistent
            # with allowing native_types like None/str/bool to be returned
            # from value().  We should revisit if that is worthwhile to do
            # here.
            #
            if check_if_numeric_type(obj):
                return obj
            else:
                if not exception:
                    return None
                raise TypeError(
                    "Cannot evaluate object with unknown type: %s" % obj.__class__.__name__
                )
        #
        # Evaluate the expression object
        #
        if exception:
            #
            # Here, we try to catch the exception
            #
            try:
                tmp = obj(exception=True)
                if tmp is None:
>                   raise ValueError(
                        "No value for uninitialized NumericValue object %s" % (obj.name,)
                    )
E                   ValueError: No value for uninitialized NumericValue object x1[0.015505]

pyomo/pyomo/common/numeric_types.py:384: ValueError

[AlbertLee125]

The error in Jenkins is related to a lack of initialization of the variables after the discretization. There is a related bug report in issue #3101
Should we skip this test until that issue is resolved?