Pyomo: Building Up A Constraint Over An Index

Question

I am trying to create a constraint for model.z that builds up as a summation during a new time index. I tried to use the += but this does not work with Constraints it seems. I am also trying to use an Expression, because I think that this might be the way to go here, but I am not having much success just yet. I included a M.W.E below and separated the output. The output that I am trying to accomplish is also included.

from __future__ import division
from pyomo.environ import *
from pyomo import environ as pym
model = ConcreteModel()

Imax = 1
Jmax = 1
Tmax = 3
model.Iset = RangeSet(1, Imax)
model.Jset = RangeSet(1, Jmax)
model.Tset = RangeSet(0, Tmax)
model.Tset2 = RangeSet(1, Tmax)

model.x = Var(model.Iset, model.Jset, model.Tset, initialize=5)
model.y = Var(model.Iset, model.Jset, model.Tset)
model.y[1,1,0].fix(1)
model.z = Var(model.Iset, model.Jset, model.Tset)

def first_rule(model,i,j,t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.y[i,j,t] == model.x[i,j,t-1]
model.rule1 = Constraint(model.Iset,model.Jset,model.Tset, rule=first_rule)

def second_rule(model,i,j,t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.z[i,j,t] == sum(model.y[i,j,t-1] for t in model.Tset2)
model.rule2 = Constraint(model.Iset,model.Jset,model.Tset, rule=second_rule)

model.x.pprint()
model.rule1.pprint()
model.rule2.pprint()

Current Output:

x : Size=4, Index=x_index
    Key       : Lower : Value : Upper : Fixed : Stale : Domain
    (1, 1, 0) :  None :     5 :  None : False : False :  Reals
    (1, 1, 1) :  None :     5 :  None : False : False :  Reals
    (1, 1, 2) :  None :     5 :  None : False : False :  Reals
    (1, 1, 3) :  None :     5 :  None : False : False :  Reals
rule1 : Size=3, Index=rule1_index, Active=True
    Key       : Lower : Body                : Upper : Active
    (1, 1, 1) :   0.0 : y[1,1,1] - x[1,1,0] :   0.0 :   True
    (1, 1, 2) :   0.0 : y[1,1,2] - x[1,1,1] :   0.0 :   True
    (1, 1, 3) :   0.0 : y[1,1,3] - x[1,1,2] :   0.0 :   True
rule2 : Size=3, Index=rule2_index, Active=True
    Key       : Lower : Body                                        : Upper : Active
    (1, 1, 1) :   0.0 : z[1,1,1] - (y[1,1,0] + y[1,1,1] + y[1,1,2]) :   0.0 :   True
    (1, 1, 2) :   0.0 : z[1,1,2] - (y[1,1,0] + y[1,1,1] + y[1,1,2]) :   0.0 :   True
    (1, 1, 3) :   0.0 : z[1,1,3] - (y[1,1,0] + y[1,1,1] + y[1,1,2]) :   0.0 :   True

Desired Output:

x : Size=4, Index=x_index
    Key       : Lower : Value : Upper : Fixed : Stale : Domain
    (1, 1, 0) :  None :     5 :  None : False : False :  Reals
    (1, 1, 1) :  None :     5 :  None : False : False :  Reals
    (1, 1, 2) :  None :     5 :  None : False : False :  Reals
    (1, 1, 3) :  None :     5 :  None : False : False :  Reals
rule1 : Size=3, Index=rule1_index, Active=True
    Key       : Lower : Body                : Upper : Active
    (1, 1, 1) :   0.0 : y[1,1,1] - x[1,1,0] :   0.0 :   True
    (1, 1, 2) :   0.0 : y[1,1,2] - x[1,1,1] :   0.0 :   True
    (1, 1, 3) :   0.0 : y[1,1,3] - x[1,1,2] :   0.0 :   True
rule2 : Size=3, Index=rule2_index, Active=True
    Key       : Lower : Body                                        : Upper : Active
    (1, 1, 1) :   0.0 : z[1,1,1] - (y[1,1,0])                       :   0.0 :   True
    (1, 1, 2) :   0.0 : z[1,1,2] - (y[1,1,0] + y[1,1,1])            :   0.0 :   True
    (1, 1, 3) :   0.0 : z[1,1,3] - (y[1,1,0] + y[1,1,1] + y[1,1,2]) :   0.0 :   True

Answer 1

You just need to sum over another index variable than $t$. Here is the correct code:

from __future__ import division
from pyomo.environ import *
from pyomo import environ as pym
model = ConcreteModel()

Imax = 1
Jmax = 1
Tmax = 3
model.Iset = RangeSet(1, Imax)
model.Jset = RangeSet(1, Jmax)
model.Tset = RangeSet(0, Tmax)
model.Tset2 = RangeSet(1, Tmax)

model.x = Var(model.Iset, model.Jset, model.Tset, initialize=5)
model.y = Var(model.Iset, model.Jset, model.Tset)
model.y[1,1,0].fix(1)
model.z = Var(model.Iset, model.Jset, model.Tset)

def first_rule(model,i,j,t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.y[i,j,t] == model.x[i,j,t-1]
model.rule1 = Constraint(model.Iset,model.Jset,model.Tset, rule=first_rule)

def second_rule(model,i,j,t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.z[i,j,t] == sum(model.y[i,j,k] for k in model.Tset2 if k<= t)
model.rule2 = Constraint(model.Iset,model.Jset,model.Tset, rule=second_rule)

model.x.pprint()
model.rule1.pprint()
model.rule2.pprint()

which produces

x : Size=4, Index=x_index
    Key       : Lower : Value : Upper : Fixed : Stale : Domain
    (1, 1, 0) :  None :     5 :  None : False : False :  Reals
    (1, 1, 1) :  None :     5 :  None : False : False :  Reals
    (1, 1, 2) :  None :     5 :  None : False : False :  Reals
    (1, 1, 3) :  None :     5 :  None : False : False :  Reals
rule1 : Size=3, Index=rule1_index, Active=True
    Key       : Lower : Body                : Upper : Active
    (1, 1, 1) :   0.0 : y[1,1,1] - x[1,1,0] :   0.0 :   True
    (1, 1, 2) :   0.0 : y[1,1,2] - x[1,1,1] :   0.0 :   True
    (1, 1, 3) :   0.0 : y[1,1,3] - x[1,1,2] :   0.0 :   True
rule2 : Size=3, Index=rule2_index, Active=True
    Key       : Lower : Body                                        : Upper : Active
    (1, 1, 1) :   0.0 :                         z[1,1,1] - y[1,1,1] :   0.0 :   True
    (1, 1, 2) :   0.0 :            z[1,1,2] - (y[1,1,1] + y[1,1,2]) :   0.0 :   True
    (1, 1, 3) :   0.0 : z[1,1,3] - (y[1,1,1] + y[1,1,2] + y[1,1,3]) :   0.0 :   True
```

Answer 2

The best code was answered above, however, I think that this is another way to implement what I intended as well. Although, it is more complicated but it does provide another way of thinking about the problem for anyone who sees it down the line.

from __future__ import division
from pyomo.environ import *
from pyomo import environ as pym
model = ConcreteModel()

Imax = 1
Jmax = 1
Tmax = 5
model.Iset = RangeSet(1, Imax)
model.Jset = RangeSet(1, Jmax)
model.Tset = RangeSet(0, Tmax)
model.Tset2 = RangeSet(1, Tmax)

model.x = Var(model.Iset, model.Jset, model.Tset, initialize=5)
model.y = Var(model.Iset, model.Jset, model.Tset)
model.y[1,1,0].fix(1)
model.z = Var(model.Iset, model.Jset, model.Tset)

def first_rule(model,i,j,t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.y[i,j,t] == model.x[i,j,t-1]
model.rule1 = Constraint(model.Iset,model.Jset,model.Tset, rule=first_rule)

def second_rule(model,i,j,t):
    if t == 0:
        return Constraint.Skip
    if t == 1:
        return model.z[i,j,t] == model.y[i,j,t-1]
    else:
        return model.z[i,j,t] == model.z[i,j,t-1] + model.y[i,j,t-1]
model.rule2 = Constraint(model.Iset,model.Jset,model.Tset, rule=second_rule)

model.pprint()

Edited Output:

2 Constraint Declarations
    rule1 : Size=5, Index=rule1_index, Active=True
        Key       : Lower : Body                : Upper : Active
        (1, 1, 1) :   0.0 : y[1,1,1] - x[1,1,0] :   0.0 :   True
        (1, 1, 2) :   0.0 : y[1,1,2] - x[1,1,1] :   0.0 :   True
        (1, 1, 3) :   0.0 : y[1,1,3] - x[1,1,2] :   0.0 :   True
        (1, 1, 4) :   0.0 : y[1,1,4] - x[1,1,3] :   0.0 :   True
        (1, 1, 5) :   0.0 : y[1,1,5] - x[1,1,4] :   0.0 :   True
    rule2 : Size=5, Index=rule2_index, Active=True
        Key       : Lower : Body                             : Upper : Active
        (1, 1, 1) :   0.0 :              z[1,1,1] - y[1,1,0] :   0.0 :   True
        (1, 1, 2) :   0.0 : z[1,1,2] - (z[1,1,1] + y[1,1,1]) :   0.0 :   True
        (1, 1, 3) :   0.0 : z[1,1,3] - (z[1,1,2] + y[1,1,2]) :   0.0 :   True
        (1, 1, 4) :   0.0 : z[1,1,4] - (z[1,1,3] + y[1,1,3]) :   0.0 :   True
        (1, 1, 5) :   0.0 : z[1,1,5] - (z[1,1,4] + y[1,1,4]) :   0.0 :   True