Add the DAE transformation support for Disjunct

ZedongPeng commented 9 months ago

Summary

Currently, the DAE transformation ignores the derivative equations in Disjunct. It would be highly beneficial if we can support it.

Rationale

If we try the following code, the derivative equation inside Disjunct will not be transformed.

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2022
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

# Sample Problem 1 (Ex 1 from Dynopt Guide)
#
#   min X2(tf)
#   s.t.    X1_dot = u          X1(0) = 1
#       X2_dot = X1^2 + u^2     X2(0) = 0
#       tf = 1

from pyomo.environ import *
from pyomo.dae import *
from pyomo.gdp import Disjunct

m = ConcreteModel()

m.t = ContinuousSet(bounds=(0, 1))

m.x1 = Var(m.t, bounds=(0, 1))
m.x2 = Var(m.t, bounds=(0, 1))
m.u = Var(m.t, initialize=0)

m.x1dot = DerivativeVar(m.x1)
m.x2dot = DerivativeVar(m.x2)

m.obj = Objective(expr=m.x2[1])

def _x1dot(M, i):
    if i == 0:
        return Constraint.Skip
    return M.x1dot[i] == M.u[i]

m.x1dotcon = Constraint(m.t, rule=_x1dot)

m.bb = Disjunct()

# original code
# def _x2dot(M, i):
#     if i == 0:
#         return Constraint.Skip
#     return M.x2dot[i] == M.x1[i] + M.u[i] ** 2

# m.x2dotcon = Constraint(m.t, rule=_x2dot)

# replace by the following
def _x2dot(M, i):
    model = M.model()
    if i == 0:
        return Constraint.Skip
    return model.x2dot[i] == model.x1[i] + model.u[i] ** 2

m.bb.x2dotcon = Constraint(m.t, rule=_x2dot)

def _init(M):
    yield M.x1[0] == 1
    yield M.x2[0] == 0
    yield ConstraintList.End

m.init_conditions = ConstraintList(rule=_init)

discretizer = TransformationFactory('dae.collocation')
discretizer.apply_to(m, nfe=5, ncp=6, scheme='LAGRANGE-RADAU')

m.bb.x2dotcon.pprint()

Output

x2dotcon : Size=1, Index=t, Active=True
    Key : Lower : Body                         : Upper : Active
      1 :   0.0 : x2dot[1] - (x1[1] + u[1]**2) :   0.0 :   True

Description

Maybe we just need to change some component_objects(Block, descend_into=True) into component_objects([Block, Disjunct], descend_into=True). Any ideas? @blnicho

ZedongPeng commented 9 months ago

If we reconstruct the constraint, the result will be correct as follows.

x2dotcon : Size=30, Index=t, Active=True
    Key      : Lower : Body                                              : Upper : Active
    0.007962 :   0.0 : x2dot[0.007962] - (x1[0.007962] + u[0.007962]**2) :   0.0 :   True
    0.039603 :   0.0 : x2dot[0.039603] - (x1[0.039603] + u[0.039603]**2) :   0.0 :   True
    0.087595 :   0.0 : x2dot[0.087595] - (x1[0.087595] + u[0.087595]**2) :   0.0 :   True
    0.139093 :   0.0 : x2dot[0.139093] - (x1[0.139093] + u[0.139093]**2) :   0.0 :   True
    0.180293 :   0.0 : x2dot[0.180293] - (x1[0.180293] + u[0.180293]**2) :   0.0 :   True
         0.2 :   0.0 :                x2dot[0.2] - (x1[0.2] + u[0.2]**2) :   0.0 :   True
    0.207962 :   0.0 : x2dot[0.207962] - (x1[0.207962] + u[0.207962]**2) :   0.0 :   True
    0.239603 :   0.0 : x2dot[0.239603] - (x1[0.239603] + u[0.239603]**2) :   0.0 :   True
    0.287595 :   0.0 : x2dot[0.287595] - (x1[0.287595] + u[0.287595]**2) :   0.0 :   True
    0.339093 :   0.0 : x2dot[0.339093] - (x1[0.339093] + u[0.339093]**2) :   0.0 :   True
    0.380293 :   0.0 : x2dot[0.380293] - (x1[0.380293] + u[0.380293]**2) :   0.0 :   True
         0.4 :   0.0 :                x2dot[0.4] - (x1[0.4] + u[0.4]**2) :   0.0 :   True
    0.407962 :   0.0 : x2dot[0.407962] - (x1[0.407962] + u[0.407962]**2) :   0.0 :   True
    0.439603 :   0.0 : x2dot[0.439603] - (x1[0.439603] + u[0.439603]**2) :   0.0 :   True
    0.487595 :   0.0 : x2dot[0.487595] - (x1[0.487595] + u[0.487595]**2) :   0.0 :   True
    0.539093 :   0.0 : x2dot[0.539093] - (x1[0.539093] + u[0.539093]**2) :   0.0 :   True
    0.580293 :   0.0 : x2dot[0.580293] - (x1[0.580293] + u[0.580293]**2) :   0.0 :   True
         0.6 :   0.0 :                x2dot[0.6] - (x1[0.6] + u[0.6]**2) :   0.0 :   True
    0.607962 :   0.0 : x2dot[0.607962] - (x1[0.607962] + u[0.607962]**2) :   0.0 :   True
    0.639603 :   0.0 : x2dot[0.639603] - (x1[0.639603] + u[0.639603]**2) :   0.0 :   True
    0.687595 :   0.0 : x2dot[0.687595] - (x1[0.687595] + u[0.687595]**2) :   0.0 :   True
    0.739093 :   0.0 : x2dot[0.739093] - (x1[0.739093] + u[0.739093]**2) :   0.0 :   True
    0.780293 :   0.0 : x2dot[0.780293] - (x1[0.780293] + u[0.780293]**2) :   0.0 :   True
         0.8 :   0.0 :                x2dot[0.8] - (x1[0.8] + u[0.8]**2) :   0.0 :   True
    0.807962 :   0.0 : x2dot[0.807962] - (x1[0.807962] + u[0.807962]**2) :   0.0 :   True
    0.839603 :   0.0 : x2dot[0.839603] - (x1[0.839603] + u[0.839603]**2) :   0.0 :   True
    0.887595 :   0.0 : x2dot[0.887595] - (x1[0.887595] + u[0.887595]**2) :   0.0 :   True
    0.939093 :   0.0 : x2dot[0.939093] - (x1[0.939093] + u[0.939093]**2) :   0.0 :   True
    0.980293 :   0.0 : x2dot[0.980293] - (x1[0.980293] + u[0.980293]**2) :   0.0 :   True
           1 :   0.0 :                      x2dot[1] - (x1[1] + u[1]**2) :   0.0 :   True

ZedongPeng commented 9 months ago

Hi @blnicho . Any idea about this?

Pyomo / pyomo