Using OMLT block (Linear Tree) to estimate derivative w.r.t. time leads to error when applying pyomo.environ.TransformationFactory

[tomasvanoyen]

Hi,

first of all let me thank you for this nice work.

I am trying to use the library to provide a data-driven surrogate model of a time derivative of quantity (say speed) and use this optimize the total cost of a route, considering also that the cost of the accelerating depends on time.

As such, I am defining the following:

`t_sim = np.array(range(0, 100)) model1 = pyo.ConcreteModel() model1.t = pyo_dae.ContinuousSet(initialize = t_sim) model1.speed = pyo.Var(model1.t, bounds=(0, 40)) model1.modelled_acceleration = pyo.Var(model1.t, bounds=(-10, 10)) model1.forc = pyo.Var(model1.t, bounds=(0, 10)) model1.cost_petrol = pyo.Var(model1.t) model1.speed_dot = pyo_dae.DerivativeVar(model1.speed, wrt=model1.t) model1.lt = OmltBlock() formulation1_lt = LinearTreeGDPFormulation(ltmodel, transformation="hull") # ltmodel being a LinearTreeDefinition of a regression model model1.lt.build_formulation(formulation1_lt)

model1.connection_constrain_1 = pyo.Constraint( model1.t, rule=lambda m, t: m.forc[t] == model1.lt.inputs[0] ) model1.connection_constrain_2 = pyo.Constraint( model1.t, rule=lambda m, t: m.modelled_acceleration[t] == model1.lt.outputs[0] ) model1.speed_dot_ode = pyo.Constraint( model1.t, rule=lambda m, t: m.speed_dot[t] == m.modelled_acceleration[t] - 1 ) model1.speed_min = pyo.Constraint( range(0, len(t_sim)), rule=lambda m, k: m.speed[t_sim[k]] >= 10 ) model1.control = sum([model1.forc[t_sim[k]] for k in range(0, len(t_sim))]) model1.obj = pyo.Objective( expr = model1.control, sense=pyo.minimize )`

However, applying then the TransformationFactory pyo.TransformationFactory( "dae.collocation" ).apply_to(model1, nfe=len(t_sim))

provides the error below.

Could you provide any guidance on how to resolve this error?

Thanks in advance!

Tomas

---

TypeError Traceback (most recent call last) Cell In[37], line 3 1 pyo.TransformationFactory( 2 "dae.collocation" ----> 3 ).apply_to(model1, nfe=len(t_sim))

File ~/mambaforge/envs/fs_ems_ops/lib/python3.11/site-packages/pyomo/core/base/transformation.py:77, in Transformation.apply_to(self, model, kwds) 75 if not hasattr(model, '_transformation_data'): 76 model._transformation_data = TransformationData() ---> 77 reverse_token = self._apply_to(model, kwds) 78 timer.report() 80 return reverse_token

File ~/mambaforge/envs/fs_ems_ops/lib/python3.11/site-packages/pyomo/dae/plugins/colloc.py:464, in Collocation_Discretization_Transformation._apply_to(self, instance, **kwds) 461 elif self._scheme_name == 'LAGRANGE-LEGENDRE': 462 self._get_legendre_constants(currentds) --> 464 self._transformBlock(instance, currentds)

File ~/mambaforge/envs/fs_ems_ops/lib/python3.11/site-packages/pyomo/dae/plugins/colloc.py:501, in Collocation_Discretization_Transformation._transformBlock(self, block, currentds) 498 disc_info['afinal'] = self._afinal[currentds] 499 disc_info['scheme'] = self._scheme_name --> 501 expand_components(block) 503 for d in block.component_objects(DerivativeVar, descend_into=True): 504 dsets = d.get_continuousset_list()

File ~/mambaforge/envs/fs_ems_ops/lib/python3.11/site-packages/pyomo/dae/misc.py:124, in expand_components(block) 118 # Record the missing BlockData before expanding components. This is for 119 # the case where a ContinuousSet indexed Block is used in a Constraint. 120 # If the Constraint is expanded before the Block then the missing 121 # BlockData will be added to the indexed Block but will not be 122 # constructed correctly. 123 for blk in block.component_objects(Block, descend_into=True): --> 124 missing_idx = set(blk.index_set()) - set(blk._data.keys()) 125 if missing_idx: 126 blk._dae_missing_idx = missing_idx

File ~/mambaforge/envs/fs_ems_ops/lib/python3.11/site-packages/pyomo/core/base/set.py:572, in SetData.iter(self) 564 def iter(self) -> Iterator[typingAny]: 565 """Iterate over the set members 566 567 Raises AttributeError for non-finite sets. This must be (...) 570 underlying indexing set). 571 """ --> 572 raise TypeError( 573 "'%s' object is not iterable (non-finite Set '%s' " 574 "is not iterable)" % (self.class.name, self.name) 575 )

TypeError: 'GlobalSet' object is not iterable (non-finite Set 'NonNegativeIntegers' is not iterable)`

[rmisener]

@bammari -- Possible to take a look at this one?

[tomasvanoyen]

Hi @rmisener @bammari ,

Thanks for this nice library.

Yet I am wondering how to solve my problem. Any take on this? Guidance on how to implement / circumvent this problem ? Or is this just not possible and we should switched to a different method?

Thanks in any case!

Regards,

Tomas

[bammari]

@tomasvanoyen Thank you for raising this issue! Another student and I are looking into this and we will respond shortly. Thank you.

[bammari]

Hi @emma58! I'm hoping I can get your input here. I believe that because we're using Pyomo.GDP for the linear tree formulations, we get this error when a NonNegativeInteger set is introduced during transformation. I haven't had the opportunity to look into this further but is there a way around this?

@tomasvanoyen You can get around this error by following Emma's Solution Below!

Please let me know if you have any additional questions in the meantime.

Bashar

[emma58]

@tomasvanoyen, @bammari is correct that the DAE transformation is getting tripped up by what the gdp.hull transformation is doing: You can fix this by calling dae.collocation first: When you construct LinearTreeGDPFormulation, set transformation="custom". This will return you a GDP formulation of the linear model tree rather than a MILP. Then, after you've built the model, call:

pyo.TransformationFactory( "dae.collocation" ).apply_to(model1, nfe=len(t_sim))
pyo.TransformationFactory("gdp.hull").apply_to(model1)

[tomasvanoyen]

Hi @emma58, @bammari thank you for your response.

I am stuck in second gear, and therefore hadn't have the chance yet to check in on your solution.

I hope to get back to you by Friday.

Thanks again, Tomas