Search code examples
pythonoptimizationpyomo

Time-delay response in Pyomo


I am trying to simulate a DAE system with several ODE's, one of them (the controller) showing a huge time lag compared to the simulation time scale. How should I implement this in Pyomo (not considering other packages, already did with Gekko but since apm.exe source is not released as open source, no longer considering the package for my applications).

Currently I have stated the problem as:

odeu = lambda m, t: tau * m.dudt[t] - (m.spu[t] - m.u[t]) == 0
model.odeu = Constraint(model.t, rule=lambda m, t: odeu(m, t))

And what I am trying to produce is something like: Currently I have stated the problem as:

odeu = lambda m, t: tau * m.dudt[t] - (m.spu[t-tde] - m.u[t]) == 0
model.odeu = Constraint(model.t, rule=lambda m, t: odeu(m, t))

The problem is that Pyomo will discretize the equations and use a weird floatting point indexing instead of local evaluations (what is required for optimization, sure, just the floatting point indexing that is weird to me), so the index t-tde does not exist.

I considered searching at each step what would be the closest index to that point, but this will exponentially increase my computation time.

Thanks!


Solution

  • A couple points...

    First, this is invalid code. I'm not sure what the final segment is, but you can't add a positional argument after the named ones and it isn't clear what you are trying to do with t:ode(m, t):

    Constraint(model.t, rule=lambda m, t: odeu(m, t))
    

    Pyomo does evaluate expressions for set indices, and as long as it falls within the set, you're GTG. I'm not sure what you mean by "floating point indexing." I think you are looking for something like this where you have a time index and some offset that you need to use in some constraints:

    # set with lag
    
    from pyomo.environ import * 
    
    mdl = ConcreteModel()    
    mdl.t = Set(initialize=range(5))    # a set index for time    
    mdl.x = Var(mdl.t)                  # a toy variable
    
    tde = 2                             # some offset
    
    # make constraint
    def c1(mdl, t):
        if t - tde < 0:                 # out of bounds
            return Constraint.Skip
        return mdl.x[t - tde] <= 10
    mdl.c1 = Constraint(mdl.t, rule=c1)
    
    mdl.pprint()
    

    yields a model that has the appropriate constraints...

    1 Set Declarations
        t : Dim=0, Dimen=1, Size=5, Domain=None, Ordered=False, Bounds=(0, 4)
            [0, 1, 2, 3, 4]
    
    1 Var Declarations
        x : Size=5, Index=t
            Key : Lower : Value : Upper : Fixed : Stale : Domain
              0 :  None :  None :  None : False :  True :  Reals
              1 :  None :  None :  None : False :  True :  Reals
              2 :  None :  None :  None : False :  True :  Reals
              3 :  None :  None :  None : False :  True :  Reals
              4 :  None :  None :  None : False :  True :  Reals
    
    1 Constraint Declarations
        c1 : Size=3, Index=t, Active=True
            Key : Lower : Body : Upper : Active
              2 :  -Inf : x[0] :  10.0 :   True
              3 :  -Inf : x[1] :  10.0 :   True
              4 :  -Inf : x[2] :  10.0 :   True
    
    3 Declarations: t x c1