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

PyomoNLP scaling factors on sub-blocks #3295

Merged
merged 19 commits into from
Aug 16, 2024
Merged

Conversation

bknueven
Copy link
Contributor

@bknueven bknueven commented Jun 19, 2024

Fixes # N/A

Summary/Motivation:

Scaling factors in the PyomoNLP interface will not be applied if the scaling_factor Suffix is on a sub-block.

Changes proposed in this PR:

  • Change PyomoNLP scaling determination methods to use SuffixFinder instead (23d316e, a7791e1)
  • Change PyomoGreyBoxNLP and PyomoNLPWithGreyBoxes to use SuffixFinder (8215ef8, c7b1bca)
  • Maintain backwards compatibility (9ca3928)
  • Update the tests so they fail with nested scaling_factors (0b4d16b)

As-is this PR would change the behavior of these methods in that they would always return scaling factors. I cannot seem to exactly replicate the prior behavior which would allow all existing tests to pass without modification.

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.

@bknueven bknueven changed the title [WIP] PyomoNLP scaling PyomoNLP scaling factors on sub-blocks Jun 21, 2024
@blnicho blnicho requested a review from Robbybp June 25, 2024 18:59
Copy link
Contributor

@Robbybp Robbybp left a comment

Choose a reason for hiding this comment

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

Thanks for the PR, this is a good change. I have a minor concern about backward incompatibility, see below. Also, can you add a small unit test to TestPyomoNLP in test_nlp.py that exercises this functionality?

Comment on lines 231 to 236
scaling_suffix_finder = SuffixFinder('scaling_factor')
for i, v in enumerate(self._pyomo_model_var_datas):
v_scaling = scaling_suffix_finder.find(v)
if v_scaling is not None:
need_scaling = True
self._primals_scaling[i] = v_scaling
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 a minor edge case, but should we enforce that we never use a suffix that is "outside the scope" (above) the pyomo_model that was used to construct the PyomoNLP? Maybe SuffixFinder should support a context argument.

Copy link
Contributor

@Robbybp Robbybp Jul 30, 2024

Choose a reason for hiding this comment

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

@blnicho @jsiirola @michaelbynum Can I get a second opinion here? The question is: what should happen when a PyomoNLP is constructed from a subblock that has relevant scaling factors on the parent block:

import pyomo.environ as pyo
from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP

m = pyo.ConcreteModel()
m.b = pyo.Block()
m.b.x = pyo.Var([1, 2], initialize={1: 100, 2: 20})

# Components so we don't have an empty NLP
m.b.eq = pyo.Constraint(expr=m.b.x[1]*m.b.x[2] == 2000)
m.b.obj = pyo.Objective(expr=m.b.x[1]**2 + m.b.x[2]**2)

m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT)
m.scaling_factor[m.b.x[1]] = 1e-2
m.scaling_factor[m.b.x[2]] = 1e-1

nlp = PyomoNLP(m.b)
scaling = nlp.get_primals_scaling()
print(scaling)

Current behavior:

None

Proposed by this PR:

[0.01, 0.1 ]

Copy link
Member

Choose a reason for hiding this comment

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

I think that I would advocate the old behavior (although I can see value in both approaches!).

My rationale is the following: if we define a "model" as all "active components reachable through active blocks contained within the reference block," then we should exclude Suffixes outside of the subtree, as Suffix is an active component.

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 inclined to agree.

Comment on lines 307 to 311
scaling_suffix_finder = SuffixFinder('scaling_factor', 1.0)
primals_scaling = np.ones(self.n_primals())
for i, v in enumerate(self.get_pyomo_variables()):
primals_scaling[i] = scaling_suffix_finder.find(v)
return primals_scaling
Copy link
Contributor

Choose a reason for hiding this comment

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

Returning ones instead of None here when no scaling factors have been defined is the backward incompatibility I'm concerned about. Maybe something like this would be better?

scaling_suffix_finder = SuffixFinder("scaling_factor")
scaling_factors = [scaling_suffix_finder.find(v) for v in self.get_pyomo_variables()]
has_scaling_factors = any(sf is not None for sf in scaling_factors)
if has_scaling_factors:
    return [sf if sf is not None else 1.0 for sf in scaling_factors]
else:
    return None

It's not exactly the same as having any non-empty scaling_factor suffix used to be enough to avoid returning None. Now we need to actually contain one of the primal variables. Maybe we could make the SuffixFinder._get_suffix_list method public, and return None if all these suffixes are empty?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I implemented something equivalent originally (see 23d316e) but this will break the existing "hacks" in the implicit function solver:

# I need scaling_factor so Pyomo NLPs I create from these blocks
# don't break when ProjectedNLP calls get_primals_scaling
block.scaling_factor = Suffix(direction=Suffix.EXPORT)
# HACK: scaling_factor just needs to be nonempty
block.scaling_factor[block._obj] = 1.0

To support both cases this function could an optional argument always_return_numeric=False or something similar to determine which behavior is desired. Maybe you have a better idea about how to de-conflict?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, if we do this, I'll have to explicitly set scaling factors to 1 in the implicit function solver, which wouldn't be the worst thing ever. I think explicitly checking SuffixFinder._get_suffix_list for at least one non-empty suffix would get around this as well.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

In 9ca3928 I have opted just to replicate the original check. This seems the safest and doesn't involve accessing private methods or modifying other parts of Pyomo.

@bknueven
Copy link
Contributor Author

bknueven commented Jun 27, 2024

Thanks for the PR, this is a good change. I have a minor concern about backward incompatibility, see below. Also, can you add a small unit test to TestPyomoNLP in test_nlp.py that exercises this functionality?

Thanks for the review. I believe I have addressed all comments except for potential modifications to SuffixFinder to support a context argument. I believe the same potential issue would also affect the NL writer.

Copy link
Contributor

@Robbybp Robbybp left a comment

Choose a reason for hiding this comment

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

Thanks for making this backwards compatible. That part looks good. Just want to get a second opinion about suffixes on parent components.

Comment on lines 231 to 236
scaling_suffix_finder = SuffixFinder('scaling_factor')
for i, v in enumerate(self._pyomo_model_var_datas):
v_scaling = scaling_suffix_finder.find(v)
if v_scaling is not None:
need_scaling = True
self._primals_scaling[i] = v_scaling
Copy link
Contributor

@Robbybp Robbybp Jul 30, 2024

Choose a reason for hiding this comment

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

@blnicho @jsiirola @michaelbynum Can I get a second opinion here? The question is: what should happen when a PyomoNLP is constructed from a subblock that has relevant scaling factors on the parent block:

import pyomo.environ as pyo
from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP

m = pyo.ConcreteModel()
m.b = pyo.Block()
m.b.x = pyo.Var([1, 2], initialize={1: 100, 2: 20})

# Components so we don't have an empty NLP
m.b.eq = pyo.Constraint(expr=m.b.x[1]*m.b.x[2] == 2000)
m.b.obj = pyo.Objective(expr=m.b.x[1]**2 + m.b.x[2]**2)

m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT)
m.scaling_factor[m.b.x[1]] = 1e-2
m.scaling_factor[m.b.x[2]] = 1e-1

nlp = PyomoNLP(m.b)
scaling = nlp.get_primals_scaling()
print(scaling)

Current behavior:

None

Proposed by this PR:

[0.01, 0.1 ]

Copy link

codecov bot commented Jul 31, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 88.52%. Comparing base (f46da7f) to head (3a083c3).

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #3295      +/-   ##
==========================================
- Coverage   88.52%   88.52%   -0.01%     
==========================================
  Files         868      868              
  Lines       98398    98413      +15     
==========================================
+ Hits        87109    87116       +7     
- Misses      11289    11297       +8     
Flag Coverage Δ
linux 86.03% <100.00%> (+<0.01%) ⬆️
osx 75.31% <100.00%> (+<0.01%) ⬆️
other 86.53% <100.00%> (-0.02%) ⬇️
win 83.84% <100.00%> (-0.02%) ⬇️

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.

@blnicho
Copy link
Member

blnicho commented Aug 15, 2024

@bknueven the PR adding the context option to SuffixFinder has been merged. I think you were planning to use that in this PR based on the discussion during the dev call. Leave a comment when this is ready to re-review.

@bknueven
Copy link
Contributor Author

@Robbybp @blnicho this is ready for re-review after integrating the new context option for the SuffixFinder.

Copy link
Member

@jsiirola jsiirola left a comment

Choose a reason for hiding this comment

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

This is fine, but there is a cleaner / more efficient way to make use of the SuffixFinder

Comment on lines 302 to 308
val = SuffixFinder('scaling_factor', context=self._pyomo_model).find(obj)
# maintain backwards compatibility
scaling_suffix = self._pyomo_model.component('scaling_factor')
if scaling_suffix and scaling_suffix.ctype is pyo.Suffix:
if obj in scaling_suffix:
return scaling_suffix[obj]
return 1.0
return None
return 1.0 if val is None else val
else:
return val
Copy link
Member

Choose a reason for hiding this comment

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

This is fine, but could be significantly simpler:

scaling_finder = SuffixFinder('scaling_factor', default=1.0, context=self._pyomo_model)
val = suffix_finder.find(self.get_pyomo_objective())
if not scaling_finder.all_suffixes:
    return None
return val

Comment on lines 311 to 327
def get_primals_scaling(self):
scaling_suffix_finder = SuffixFinder(
'scaling_factor', context=self._pyomo_model
)
primals_scaling = np.ones(self.n_primals())
ret = None
for i, v in enumerate(self.get_pyomo_variables()):
val = scaling_suffix_finder.find(v)
if val is not None:
primals_scaling[i] = val
ret = primals_scaling
# maintain backwards compatibility
scaling_suffix = self._pyomo_model.component('scaling_factor')
if scaling_suffix and scaling_suffix.ctype is pyo.Suffix:
primals_scaling = np.ones(self.n_primals())
for i, v in enumerate(self.get_pyomo_variables()):
if v in scaling_suffix:
primals_scaling[i] = scaling_suffix[v]
return primals_scaling
return None
else:
return ret
Copy link
Member

Choose a reason for hiding this comment

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

This is fine, but could be significantly simpler:

scaling_finder = SuffixFinder('scaling_factor', default=1.0, context=self._pyomo_model)
primals_scaling = numpy.fromiter(
    (scaling_finder.find(v) for v in self.get_pyomo_variables()), 
    count=self.n_primals(),
)
if not scaling_finder.all_suffixes:
    return None
return primals_scaling

Comment on lines 330 to 346
def get_constraints_scaling(self):
scaling_suffix_finder = SuffixFinder(
'scaling_factor', context=self._pyomo_model
)
constraints_scaling = np.ones(self.n_constraints())
ret = None
for i, c in enumerate(self.get_pyomo_constraints()):
val = scaling_suffix_finder.find(c)
if val is not None:
constraints_scaling[i] = val
ret = constraints_scaling
# maintain backwards compatibility
scaling_suffix = self._pyomo_model.component('scaling_factor')
if scaling_suffix and scaling_suffix.ctype is pyo.Suffix:
constraints_scaling = np.ones(self.n_constraints())
for i, c in enumerate(self.get_pyomo_constraints()):
if c in scaling_suffix:
constraints_scaling[i] = scaling_suffix[c]
return constraints_scaling
return None
else:
return ret
Copy link
Member

Choose a reason for hiding this comment

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

This is fine, but could be significantly simpler:

scaling_finder = SuffixFinder('scaling_factor', default=1.0, context=self._pyomo_model)
constraints_scaling = numpy.fromiter(
    (scaling_finder.find(v) for v in self.get_pyomo_constraints()), 
    count=self.n_constriaints(),
)
if not scaling_finder.all_suffixes:
    return None
return constraints_scaling

Comment on lines 627 to 615
self._primals_scaling = np.ones(self.n_primals())
scaling_suffix = self._pyomo_nlp._pyomo_model.component('scaling_factor')
scaling_suffix_finder = SuffixFinder(
'scaling_factor', context=self._pyomo_model
)
for i, v in enumerate(self.get_pyomo_variables()):
v_scaling = scaling_suffix_finder.find(v)
if v_scaling is not None:
need_scaling = True
self._primals_scaling[i] = v_scaling
# maintain backwards compatibility
scaling_suffix = self._pyomo_model.component('scaling_factor')
if scaling_suffix and scaling_suffix.ctype is pyo.Suffix:
need_scaling = True
for i, v in enumerate(self.get_pyomo_variables()):
if v in scaling_suffix:
self._primals_scaling[i] = scaling_suffix[v]
Copy link
Member

Choose a reason for hiding this comment

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

This is fine, but could be significantly simpler:

scaling_finder = SuffixFinder('scaling_factor', default=1.0, context=self._pyomo_model)
self._primals_scaling = numpy.fromiter(
    (scaling_finder.find(v) for v in self.get_pyomo_variables()), 
    count=self.n_primals(),
)
need_scaling = bool(scaling_finder.all_suffixes)

@bknueven
Copy link
Contributor Author

This is fine, but there is a cleaner / more efficient way to make use of the SuffixFinder

Thanks for showing me that; I implemented your suggestion.

Copy link
Contributor

@Robbybp Robbybp left a comment

Choose a reason for hiding this comment

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

Can you add a test that makes sure we don't pull a scaling factor from a parent block?

@bknueven
Copy link
Contributor Author

Can you add a test that makes sure we don't pull a scaling factor from a parent block?

Sure, I just copy+pasted your example above and made it a test in 7e79e19.

Copy link
Contributor

@Robbybp Robbybp left a comment

Choose a reason for hiding this comment

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

Thanks, looks good.

@jsiirola jsiirola merged commit 29f8ede into Pyomo:main Aug 16, 2024
25 of 29 checks passed
lbianchi-lbl added a commit to lbianchi-lbl/watertap that referenced this pull request Aug 20, 2024
lbianchi-lbl added a commit to lbianchi-lbl/watertap that referenced this pull request Aug 21, 2024
lbianchi-lbl added a commit to watertap-org/watertap that referenced this pull request Aug 21, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
No open projects
Development

Successfully merging this pull request may close these issues.

5 participants