pyomo
pyomo copied to clipboard
Pyomo
Pyomo DAE flatten ignores whether blocks are active
Summary
In July, I ran into unexpected behavior when the PETSc interface in IDAES uses pyomo.dae.flatten to create the block with variables/equations to be solved. It appears that the flattening process ignores whether the block it's aggregating variables from is active or not, and so variables and constraints from deactivated blocks end up being used in the problem because they are not individually deactivated.
@jsiirola suggested the following fix, which has worked for my applications:
diff --git a/pyomo/dae/flatten.py b/pyomo/dae/flatten.py
index 922ffa9..13f56ee 100644
--- a/pyomo/dae/flatten.py
+++ b/pyomo/dae/flatten.py
@@ -305,7 +305,10 @@ def generate_sliced_components(b, index_stack, slice_, sets, ctype, index_map):
if type(descend_data) is IndexedComponent_slice:
try:
# Attempt to find a data object matching this slice
- descend_data = next(iter(descend_data))
+ _data = next(iter(descend_data))
+ while not _data.active:
+ _data = next(iter(descend_data))
+ descend_data = _data
except StopIteration:
# For this particular idx (and given indices), no
# block data object exists to descend into.
@@ -313,6 +316,8 @@ def generate_sliced_components(b, index_stack, slice_, sets, ctype, index_map):
continue
else:
descend_data = sub
+ if not descend_data.active:
+ continue
# Recursively generate sliced components from this data object
for st, v in generate_sliced_components(
I didn't open an issue then, because I assumed our conversation meant it would be fixed without it, but it probably got forgotten about as soon as a more important issue came to the fore.
Edit: Also, it would be good to make including inactive components an option. I can imagine contexts (like scaling methods) where including them would be the desired behavior.
@dallan-keylogic Thanks for opening this issue. This patch looks good, although I would also like an active=True optional argument that, by default, only descends into active blocks. I think the original reason I excluded this argument was that "active" can be ill-defined for a slice, i.e. do we check the active flag of all or any of the blocks/data objects included in the slice. That said, a flag to not descend into inactive blocks should be well-defined.
@dallan-keylogic Can you open a PR with this patch and a test indicating the desired (currently failing) behavior.
@Robbybp I can. It will take a few days to get around to because I don't have a Pyomo dev environment set up yet.
Returning to this issue finally and am running into an issue when trying to implement the suggested fix. The following code (adapted from the current tests) triggers an infinite loop with the updated function.
from pyomo.environ import ConcreteModel, Var, Set
from pyomo.common.collections import ComponentSet
from pyomo.dae.flatten import flatten_components_along_sets
m = ConcreteModel()
m.time = Set(initialize=[1,2,3])
m.space = Set(initialize=[0.0, 0.5, 1.0])
m.comp = Set(initialize=['a','b'])
m.v0 = Var()
m.v1 = Var(m.time)
m.v2 = Var(m.time, m.space)
m.v3 = Var(m.time, m.space, m.comp)
m.v_tt = Var(m.time, m.time)
m.v_tst = Var(m.time, m.space, m.time)
@m.Block()
def b(b):
@b.Block(m.time)
def b1(b1):
b1.v0 = Var()
b1.v1 = Var(m.space)
b1.v2 = Var(m.space, m.comp)
@b1.Block(m.space)
def b_s(b_s):
b_s.v0 = Var()
b_s.v1 = Var(m.space)
b_s.v2 = Var(m.space, m.comp)
@b.Block(m.time, m.space)
def b2(b2):
b2.v0 = Var()
b2.v1 = Var(m.comp)
b2.v2 = Var(m.time, m.comp)
# Deactivating b1 should get rid of both variables directly on it
# as well as those on the subblock b_s
m.b.b1.deactivate()
sets = ComponentSet((m.time,))
sets_list, comps_list = flatten_components_along_sets(m, sets, Var)
@dallan-keylogic thanks for trying this out. See my first attempt at this functionality in #2643. Can you test this out and see if it behaves as you would expect? It has the following behavior with a slight modification of your example script:
from pyomo.environ import ConcreteModel, Var, Set, ComponentUID
from pyomo.common.collections import ComponentSet
from pyomo.dae.flatten import flatten_components_along_sets
from pyomo.core.base.set import UnindexedComponent_set
m = ConcreteModel()
m.time = Set(initialize=[1,2,3])
m.space = Set(initialize=[0.0, 0.5, 1.0])
m.comp = Set(initialize=['a','b'])
m.v0 = Var()
m.v1 = Var(m.time)
m.v2 = Var(m.time, m.space)
m.v3 = Var(m.time, m.space, m.comp)
m.v_tt = Var(m.time, m.time)
m.v_tst = Var(m.time, m.space, m.time)
@m.Block()
def b(b):
@b.Block(m.time)
def b1(b1):
b1.v0 = Var()
b1.v1 = Var(m.space)
b1.v2 = Var(m.space, m.comp)
@b1.Block(m.space)
def b_s(b_s):
b_s.v0 = Var()
b_s.v1 = Var(m.space)
b_s.v2 = Var(m.space, m.comp)
@b.Block(m.time, m.space)
def b2(b2):
b2.v0 = Var()
b2.v1 = Var(m.comp)
b2.v2 = Var(m.time, m.comp)
# Deactivating b1 should get rid of both variables directly on it
# as well as those on the subblock b_s
m.b.b1.deactivate()
sets = ComponentSet((m.time,))
sets_list, comps_list = flatten_components_along_sets(
m, sets, Var, active=True
)
for sets, comps in zip(sets_list, comps_list):
if len(sets) == 1 and sets[0] is UnindexedComponent_set:
print([s.name for s in sets])
for c in comps:
print(" %s" % c.name)
else:
print([s.name for s in sets])
for c in comps:
print(" %s" % ComponentUID(c.referent))
['UnindexedComponent_set']
v0
['time']
v1[*]
v2[*,0.0]
v2[*,0.5]
v2[*,1.0]
v3[*,0.0,a]
v3[*,0.0,b]
v3[*,0.5,a]
v3[*,0.5,b]
v3[*,1.0,a]
v3[*,1.0,b]
b.b2[*,0.0].v0
b.b2[*,0.0].v1[a]
b.b2[*,0.0].v1[b]
b.b2[*,0.5].v0
b.b2[*,0.5].v1[a]
b.b2[*,0.5].v1[b]
b.b2[*,1.0].v0
b.b2[*,1.0].v1[a]
b.b2[*,1.0].v1[b]
['time', 'time']
v_tt[*,*]
v_tst[*,0.0,*]
v_tst[*,0.5,*]
v_tst[*,1.0,*]
b.b2[*,0.0].v2[*,a]
b.b2[*,0.0].v2[*,b]
b.b2[*,0.5].v2[*,a]
b.b2[*,0.5].v2[*,b]
b.b2[*,1.0].v2[*,a]
b.b2[*,1.0].v2[*,b]
That snippet of code works as expected. Made some comments on that PR.