Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GDP rewrite #354

Merged
merged 273 commits into from
Feb 28, 2018
Merged

GDP rewrite #354

merged 273 commits into from
Feb 28, 2018

Conversation

jsiirola
Copy link
Member

Fixes #318, #241

Summary/Motivation:

This is a rewrite of the pyomo.gdp package.

Changes proposed in this PR:

  • Rewrite of the Disjunct and Disjunction components to bring them in line with current Component conventions.
  • Adding more examples
  • Rewrites of the bigm and chull transformations
  • Adding a cuttingplane reformulation
  • Adding GDPopt solver (to pyomo.contrib)

Legal Acknowledgement

By contributing to this software project, I 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.

emma58 and others added 30 commits February 7, 2017 17:22
…d troubles), fixes variables in stickies to validate model, changes add slacks transformation to modify the body of the constraint directly
…esults from first layout on minlp.org (and it's rightgit add stickies.pygit add stickies.py)
… begins to add bounds in the batch processing model
…ing model, and notes Francisco's suggestion for second term of disjunction for whether or not there is a storage tank (but no conclusion yet)
…essions by digging into the _SumExpression objects
…essions by digging into the _SumExpression objects
…ession tests pass, adding unit test file for bigm
…ansformation block and xor constraints added to parent block and using original indicator variables
…tions, indexed constraints in the disjunction)
GDPopt solver initial implementation
Substituting he disaggregated variables into a constraint should
preserve the "equality nature" of the constraint.  In addition, catch
the special case of "x == 0" in a disjunct so we can fix the
disaggregated variable to 0.
This supports declaring "empty" disjunctions; fixes #241.
Refactor the gdp.cuttingplane transformation
  - avoid the use of CUIDs
  - improved initialization of separation problems
  - graceful exit when subproblem solves fail
Most of the changes are due to:
 - equality constraints are no longer split into inequalities
 - BigM constraints on disaggregated vars that are redundant to bounds
   are omitted.
@jsiirola jsiirola self-assigned this Feb 16, 2018
return suffix_list

def _apply_to(self, instance, targets=None, **kwds):
config = self.CONFIG().set_value(kwds.pop('options', {}))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm very confused how this is passing tests, because on my system, set_value() has no return value (so it's always None), which causes config = None, which causes issues when you try to do config.set_value(kwds) below.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need the current pyutilib master.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, yes. That seems to have done the trick.

CONFIG = ConfigBlock("gdp.bigm")
CONFIG.declare('bigM', ConfigValue(
default=None,
domain=_to_dict,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does domain mean with respect to a ConfigBlock?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that the domain is passed to the ConfigValue constructor. domain is a callable object (could be a type, function, or functor). When a value is set, the incoming value is passed as the single argument to __call__ and the return value is what is stored in the ConfigValue. At it's simplest, it can be a type like int so that all incoming values are automatically cast to int. In this case it is a little more complicated, as it accepts dicts or "values" (int/float/tuple). The former is preserved and the latter is cast to {None: val}.

To answer your question, ConfigBlocks can not have domains, although they accept an implicit_domain, which is used when the ConfigBlock accepts undeclared (implicit) keys - in which case, any implicit values are passed through the implicit_domain callable.

Copy link
Contributor

@qtothec qtothec left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few questions.
Also, I should note that bigM doesn't currently support big-E Expression objects attached to disjuncts. I assume that chull does not as well.

I got tired of looking through all the examples.


# bounds for box coordinates
# TODO: What are they trying to accomplish here? Why is it always index
# 2 for the first two??
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an error in the original GAMS file. The intent is to find the min and max x and y coordinates that the boxes can have, based on the circle locations and radii. The code as written works because circle 2 has the higher x and y coordinate values.

def no_overlap_rule(model, i, j):
return [model.no_overlap_disjunct[i, j, direction] for direction in \
model.BoxRelations]
model.no_overlap_disjunction = Disjunction(model.BoxPairs, rule=no_overlap_rule)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are we going to support the shortcut notation?

@model.Disjunction(model.BoxPairs)
def no_overlap_rule(model, i, j):
    ...

Similarly with Disjunct.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. That notation should (already) be supported.


model = AbstractModel()

# TODO: it looks like they set a bigM for each j. Which I need to look up how to do...
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

j's are stages?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, j's are stages. I think the gams file uses 3 different Big-M values, one for each disjunction. They are indexed by stage, but I don't think they actually depend on it. I never got as far as making my Big-M values match those ones.

UnitsInPhaseUB = log(6)
UnitsOutOfPhaseUB = log(6)
# TODO: YOU ARE HERE. YOU HAVEN'T ACTUALLY MADE THESE THE BOUNDS YET, NOR HAVE YOU FIGURED OUT WHOSE
# BOUNDS THEY ARE. AND THERE ARE MORE IN GAMS.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems unfinished?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that is likely... I ran into some mistakes in the gams file for this problem, and so I know I wasn't ever able to match their results for this one.

model.STAGES = Set(ordered=True)
model.PARALLELUNITS = Set(ordered=True)

# TODO: this seems like an over-complicated way to accomplish this task...
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given what's written, the filter seems reasonable.

###############

def storage_tank_selection_disjunct_rule(disjunct, selectStorageTank, j):
model = disjunct.model()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this line is not strictly necessary? :0)

return [model.storage_tank_selection_disjunct[selectTank, j] for selectTank in [0,1]]
model.select_storage_tanks = Disjunction(model.STAGESExceptLast, rule=select_storage_tanks_rule)

# though this is a disjunction in the GAMs model, it is more efficiently formulated this way:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jsiirola's thought was that in both these disjunctions you are just choosing a value for unitsInPhase and unitsOutOfPhase, so you could just formulate it that way rather than having the disjunctions. But I guess that might make it less useful as a GDP test case?

"""

alias('core.add_slack_variables', \
doc="Create a model where we had slack variables to every constraint "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo. "where we add"

elif hasattr(arg, '__getitem__'):
return (_Initializer.dict_like, arg)
else:
# Hopefully this thing is castable to teh type that is desired
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo: "the"


class GDP_Error(Exception):
"""Exception raised while processing GDP Models"""


# THe following should eventually be promoted so that all
# IndexedComponents can use it
class _Initializer(object):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some short description?

@@ -1,16 +1,16 @@
[ 0.00] Setting up Pyomo environment
[ 0.00] Applying Pyomo preprocessing actions
[ 0.00] Creating model
[ 0.01] Creating model
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This broad comment applies to all the book examples. Since changes were made to the book examples, we need to make sure that something is added to the Pyomo website. I recommend we add a link and some description below the Pyomo book on the Documentation page. I will add an issue, and we can discuss it in the pyomo meeting.

@@ -9,15 +9,16 @@
# ___________________________________________________________________________

import pyomo.core.plugins.transform.relax_integrality
#import pyomo.core.plugins.transform.eliminate_fixed_vars
#import pyomo.core.plugins.transform.standard_form
# import pyomo.core.plugins.transform.eliminate_fixed_vars
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the reason for these commented out imports?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think at some point some of the transformations were deprecated and thus commented out. This PR doesn't change that fact, other than a cosmetic addition of a space after the #.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I see. In my view it looked like they were added and then commented out.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that @whart222 "deprecated" them by removing their registration and seeing if anyone complained.

# packages are optional and/or under development.
#
_optional_packages = set([
'pyomo.contrib.example',
'pyomo.contrib.preprocessing'])
'pyomo.contrib.preprocessing',
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What are these _optional_packages used for?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The _optional_packages import the necessary contrib modules, such as the transformation plugins that I have under pyomo.contrib.preprocessing. The comment lines above describe its functionality.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, the comment is confusing since it says they "need to be loaded", but also says that they "are optional". No biggie. I am fine with it.

@@ -20,7 +20,7 @@
from pyomo.gdp import *
from pyomo.repn import generate_canonical_repn

logger = logging.getLogger('pyomo.core')
logger = logging.getLogger('pyomo.gdp')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great.

@@ -39,7 +39,7 @@ def _apply_solver(self):
xfrm.apply_to(self._instance)

xfrm = TransformationFactory('gdp.bigm')
xfrm.apply_to(self._instance, default_bigM=self.options.get('bigM',10**6))
xfrm.apply_to(self._instance, bigM=self.options.get('bigM',10**6))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Totally minor, but why 10**6 vs 100000 used in an earlier file as the default? (1) I think it was 1e6 vs 1e5, but I am not sure. Also, what is wrong with 1e5? or 1e6?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that @whart222 wrote this code and that was his style. I agree that 1e6 would be better.

Key : Lower : Body : Upper : Active
None : 1.0 : cc.expr1.indicator_var + cc.expr2.indicator_var + cc.expr3.indicator_var : 1.0 : True
Key : Disjuncts : Active : XOR
None : ['cc.expr1', 'cc.expr2', 'cc.expr3'] : True : True
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This move to list from a sum is expected?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a move from a Constraint-like to a list, with a flag for whether it's an XOR disjunction vs. an OR disjunction. It introduces support for OR disjunctions, whereas the old approach only supported XOR disjunctions.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes: this is expected. The original GDP implementation hijacked standard Constraints to implement the Disjunction (because at the time deriving Components was even more convoluted than it is now). One of the big points of the rewrite was to make GDP use "standard" components. The move to a list of disjuncts instead of the MIP equivalent is more intuitive / informative from the modeling perspective.

@@ -159,7 +159,7 @@ def _warm_start(self, instance):
mst_file.write("<header/>\n")
mst_file.write("<quality/>\n")
mst_file.write("<variables>\n")
for var in instance.component_data_objects(Var, active=True):
for var in instance.component_data_objects(Var):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we need to loop through variables that are not active?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It turns out that Vars aren't active (they are not ActiveComponents and cannot be deactivated). The active flag only controlled if the generator descended into all blocks or only active blocks. This change allows the generator to descend into all blocks - not just the active ones.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right. (We just had this conversation :)

@@ -179,7 +179,7 @@ def _warm_start(self, instance):
smap = instance.solutions.symbol_map[self._smap_id]
byObject = smap.byObject
with open(self._warm_start_file_name, 'w') as mst_file:
for vdata in instance.component_data_objects(Var, active=True):
for vdata in instance.component_data_objects(Var):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See above.

@@ -139,7 +139,7 @@ def _populate_glpk_instance ( self, model ):
if objective.is_minimizing(): sense = GLP_MIN

constraint_list = model.component_map(Constraint, active=True)
variable_list = model.component_map(Var, active=True)
variable_list = model.component_map(Var)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above.

@codecov-io
Copy link

Codecov Report

Merging #354 into master will increase coverage by <.01%.
The diff coverage is 72.07%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #354      +/-   ##
==========================================
+ Coverage   65.71%   65.72%   +<.01%     
==========================================
  Files         467      474       +7     
  Lines       60729    62289    +1560     
==========================================
+ Hits        39906    40937    +1031     
- Misses      20823    21352     +529
Impacted Files Coverage Δ
pyomo/repn/plugins/baron_writer.py 5.49% <ø> (ø) ⬆️
pyomo/solvers/plugins/solvers/GUROBI.py 29.11% <0%> (ø) ⬆️
pyomo/solvers/plugins/solvers/CPLEX.py 19.27% <0%> (ø) ⬆️
pyomo/solvers/plugins/solvers/glpk_direct.py 7.6% <0%> (ø) ⬆️
pyomo/core/base/symbolic.py 99.09% <100%> (ø) ⬆️
pyomo/gdp/plugins/bilinear.py 80.41% <100%> (ø) ⬆️
pyomo/bilevel/plugins/solver2.py 95% <100%> (ø) ⬆️
pyomo/mpec/plugins/solver2.py 100% <100%> (ø) ⬆️
pyomo/util/modeling.py 100% <100%> (ø)
pyomo/repn/plugins/cpxlp.py 93.79% <100%> (-0.32%) ⬇️
... and 37 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f9d4c48...393e327. Read the comment docs.

@jsiirola jsiirola merged commit 983a3a5 into master Feb 28, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants