DAE transformation and Arc expansion are order dependent

[andrewlee94]

I have been working with a a model using both a ContinuousSet and Arcs, and have discovered that the order in which the transformations are applied is critical in whether the model is constructed properly.

If I apply the DAE transformation first, and then expand the Arcs, I get the behaviour I would expect, where the Arc transformation writes equality constraints at all point in the ContinuousSet. However, if I apply the Arc expansion first, the equality constraints are only written at the points in the ContunuousSet that already exist, and no additional Constraints are written during the DAE expansion. This then results in unexpected degrees of freedom.

It seems clear that the DAE transformation is not interacting fully with the expanded Arcs for some reason, and I think this either needs to be fixed, or something done to make it clear that this will not work (an Exception or warning message?). Ultimately, I don't see a reason why these two should be order dependent, so I would express a preference for getting the trasnformations to work together properly.

[jsiirola]

@andrewlee94: OK - i think this is indeed a bug. I can see why DAE then Expand works (the DAE has updated all the relevant ContinuousSet objects). I thought that the problem you are seeing was because pyomo.network was not preserving the corresponding index variables, however, I just spent an hour digging through the code, and @gseastream was very careful to preserve indexing everywhere.

At this point, a simple example that we can dig through would be most helpful.

[andrewlee94]

@jsiirola I have an example on a branch of my IDAES fork you can use (https://github.com/andrewlee94/idaes-dev/tree/issue_33). If you look at ../examples/tutorials/Tutorial_3_Dynamic_Flowsheet.py, you can reverse the order of the transformations and see the problem (you might want to reduce the number of finite elements in the DAE transformation to 4 for testing).

[andrewlee94]

I have tried to reproduce this issue with a simpler model, and so far had no success. I tried the following, which I thought would roughly reproduce the operations in my model:

from pyomo.environ import *

from pyomo.dae import ContinuousSet

from pyomo.network import Arc, Port

m = ConcreteModel()

m.s = ContinuousSet(initialize=[0,10])

def block_rule(b, i):
    b.v = Var(initialize=1)
    return

m.b1 = Block(m.s, rule=block_rule)
m.b2 = Block(m.s, rule=block_rule)

m.p1 = Port()
m.p2 = Port()

m.r1 = Reference(m.b1[:].v)
m.r2 = Reference(m.b2[:].v)

m.p1.add(m.r1, "var")
m.p2.add(m.r2, "var")

m.arc = Arc(source=m.p1, destination=m.p2)

TransformationFactory("network.expand_arcs").apply_to(m)

m.discretizer = TransformationFactory('dae.finite_difference')
m.discretizer.apply_to(m,
                       nfe=4,
                       wrt=m.s,
                       scheme="BACKWARD")

m.arc_expanded.display()

However, so far everything has worked as would be expected with the Arc being properly expanded by the DAE transformation.

[andrewlee94]

@jsiirola One thought that just crossed my mind is that this might be an order of operations thing. If the DAE transformation was applied to the Arc before the References for some reason, it might be possible for this to happen.

[blnicho]

It looks like the problem arises when the Reference has multiple indexing sets. This example seems to replicate the behavior:

from pyomo.environ import *
from pyomo.dae import ContinuousSet
from pyomo.network import Arc, Port

m = ConcreteModel()

m.s = ContinuousSet(initialize=[0,10])

def block_rule(b, i):
    b.v = Var(['a','b'], initialize=1)
    return

m.b1 = Block(m.s, rule=block_rule)
m.b2 = Block(m.s, rule=block_rule)

m.p1 = Port()
m.p2 = Port()

m.r1 = Reference(m.b1[:].v[:])
m.r2 = Reference(m.b2[:].v[:])

m.p1.add(m.r1, "var")
m.p2.add(m.r2, "var")

m.arc = Arc(source=m.p1, destination=m.p2)

TransformationFactory("network.expand_arcs").apply_to(m)

m.discretizer = TransformationFactory('dae.finite_difference')
m.discretizer.apply_to(m,
                       nfe=4,
                       wrt=m.s,
                       scheme="BACKWARD")

m.arc_expanded.display()

I think this is because the discretization transformation looks for components that are explicitly indexed by the ContinuousSet and the Reference component is obscuring the individual indexing sets.

In particular, the discretization transformation relies on checking the dim of a component and then uses that to determine if the component's index or implicit_subsets should be checked for ContinuousSets. It looks like the dim attribute is not set correctly on the Reference and m.r1._implicit_subsets is None. The equality constraint added in the Arc expansion transformation gets its indexing set from the Reference component which is why the discretization transformation is not finding it.