Pyomo / pyomo

An object-oriented algebraic modeling language in Python for structured optimization problems.
https://www.pyomo.org
Other
1.92k stars 498 forks source link

Support for Integrals on hierarchical models #764

Open blnicho opened 5 years ago

blnicho commented 5 years ago

A bug report submitted to the forum found that Integral components cannot be declared on Blocks:

Here is a minimal example :

from pyomo.dae.integral import Integral
from pyomo.environ import *
from pyomo.dae.contset import ContinuousSet

m = ConcreteModel(name='m')
m.time = ContinuousSet(bounds=(0,1), name='time')

m.block = Block()
m.block.x = Var(m.time)

def _instant_obj(m, t):
    return m.x[t]
m.block.int_cost = Integral(m.time, wrt=m.time, rule=_instant_obj)   # this line raise an error

Here is the Traceback of the Error :

Traceback (most recent call last):
  File "/home/admin/Documents/02-Recherche/02-Python/lms2/lms2/template/bug_Integrals_blocks.py", line 20, in <module>
    m.block.int_cost = Integral(m.time, wrt=m.time, rule=_instant_obj)
  File "/home/admin/Documents/02-Recherche/02-Python/lms2/venv/lib/python3.6/site-packages/pyomo/core/base/block.py", line 540, in __setattr__
    self.add_component(name, val)
  File "/home/admin/Documents/02-Recherche/02-Python/lms2/venv/lib/python3.6/site-packages/pyomo/core/base/block.py", line 980, in add_component
    val.construct(data)
  File "/home/admin/Documents/02-Recherche/02-Python/lms2/venv/lib/python3.6/site-packages/pyomo/core/base/expression.py", line 410, in construct
    self.add(None, _init_rule(self._parent()))
  File "/home/admin/Documents/02-Recherche/02-Python/lms2/venv/lib/python3.6/site-packages/pyomo/dae/integral.py", line 135, in _trap_rule
    ds = sorted(m.find_component(wrt.local_name))
TypeError: 'NoneType' object is not iterable

It seems that m.find_component(wrt.local_name) is empty. So an error is raised when trying to sort it. It is empty cause he is trying to find the component 'time' in the block 'b', whereas it is a component of the model 'm'.

ReinboldV commented 5 years ago

Hi there, I am adding here some additional information about this issue. And describing new strange behaviour, I think this might be linked but may be worth a new issue (?)

Case where time set is within the root model (issue#764 small example)

from pyomo.environ import ConcreteModel, Var, Block, Expression, TransformationFactory
from pyomo.dae import ContinuousSet, Integral

# an expression for test
def _instant_obj(m, t):
    return m.x[t]

m = ConcreteModel(name='m')
m.b = Block()
m.t = ContinuousSet(bounds=(0, 1))
m.b.x = Var(m.t)
m.b.exp = Expression(m.t, rule=_instant_obj)

# m.b.int_cost = Integral(m.t, wrt=m.t, rule=lambda m, t: m.exp[t]) # Not working --> issue#764
m.int_cost = Integral(m.t, wrt=m.t, rule=lambda m, t: m.b.exp[t]) # Working cause the integral is declared within the block that holds the "time"
TransformationFactory('dae.finite_difference').apply_to(m, nfe=4)
m.pprint()
del m

Case where the time is hold by the block :

m = ConcreteModel(name='m')
m.b = Block()
m.b.t = ContinuousSet(bounds=(0, 1))
m.b.x = Var(m.b.t)
m.b.exp = Expression(m.b.t, rule=_instant_obj)

m.b.int_cost = Integral(m.b.t, wrt=m.b.t, rule=lambda m, t: m.exp[t]) # Seems to work cause the block 'b' holds the time set.
# m.int_cost = Integral(m.b.t, wrt=m.b.t, rule=lambda m, t: m.b.exp[t]) # Not working --> same behaviour as in issue#764

TransformationFactory('dae.finite_difference').apply_to(m, nfe=4)
m.pprint()

Output :

int_cost : Size=1, Index=None
    Key  : Expression
    None : 0.5*(b.x[1] + b.x[0])

In this case, the integral is defined within the block as expected but the Transformation does not behave the same as the first case, and depends on whether the integral is declared within m or m.b.
I am not sure we are facing one or multiple issues, I hope it will help...

Greetings