Pyomo: How to write a constraint for each (i,j) pair

Question

Edit: I have also added a constraint 1.5 to illustrate maybe a different way of approaching the constraint.

I am trying to write the following constraints in Pyomo for each (i,j) pair on an MxN grid:

The code that I have thus far is as follows, and I am just hoping I can get some feedback on whether or not the constraint definition is written properly to meet the intention.The idea is that each (i,j) cell on the 6x6 grid will have the following two constraints.

model = AbstractModel()

#Define the index sets for the grid, time horizions, and age classes:
model.Iset = RangeSet(6)
model.Jset = RangeSet(6)
model.Tset = RangeSet(7)
model.Kset = RangeSet(50)

#Define model parameters:
model.s = Param(within=NonNegativeIntegers)

#Define model variables:
model.juvenille = Var(model.Iset, model.Jset, model.Tset, model.Kset,
              within=NonNegativeReals, initialize = "some expression"

#Constraints: 

# Constraint #1
def juv_advance(model, i, j, t, k):
    return model.juvenille[i,j,t+1,k+1] == model.juvenille[i,j,t,k]*model.juvsurv

# Constraint #1.5
def juv_advance(model, t, k):
    return model.juvenille[t+1,k+1] == model.juvenille[t,k]*model.s \\
           for i in model.Iset for j in model.Jset

# Constraint #2
def juv_total(model, i, j, t, k):
    return sum(model.juvenille[k] for k in range(1,50))

Additionally, if anybody feels like answering this... how could you save the calculated j_t+1 value to use as the initial value in the next time period.

Your mathematically written constraints don't seem to be well-defined. In the first equation you list out the values that k can take but don't list the values i, j, or t can take. Please update your constraint descriptions so that I can give you accurate guidance on your implementation. — Bethany Nicholson, Aug 14 '19 at 14:51
I updated the values in the main image. Is that more clear? essentially on each cell of the grid will have those two constraints laid out and then ideally the model would solve at time t = 1, update the variable, solve at t = 2, update the variable, ... , and repeat that process until the time index is complete. — GrayLiterature, Aug 14 '19 at 14:58
I am also wondering if it would be possible to define an indexed set of coordinate pairs? So you could have something like... model.set = Set( (1,1), (1,2), (1,3) ... (6,6) ) — GrayLiterature, Aug 14 '19 at 16:52

PhilippK · Accepted Answer · 2019-08-14T20:44:00.253

I would try something like this:

model = AbstractModel()

#Define the index sets for the grid, time horizions, and age classes:
model.Iset = RangeSet(6)
model.Jset = RangeSet(6)
model.Tset = RangeSet(7)
model.Kset = RangeSet(50)

#Define model parameters:
model.s = Param(within=NonNegativeIntegers)

#Define model variables:
model.juvenille = Var(model.Iset, model.Jset, model.Tset, model.Kset,
              within=NonNegativeReals, initialize="some expression")

# As far as I see your problem in you second constraint the big J is a new variable ? 
 If that is the case than you have to create it:

model.J_big =Var(model.Iset, model.Jset, model.Tset, within=NonNegativeReals)

#Constraints: 

# Constraint #1
def juv_advance(model, i, j, t, k):
k_len = len(model.Kset)
t_len = len(model.Tset)
    if k == 1 and t == 1:
         return "some expression"
    elif t < t_len and k < k_len:
         return model.juvenille[i,j,t+1,k+1] == model.juvenille[i,j,t,k]*model.s
    else:
         return "Here has to come a statement what should happen with the last index (because if you are iterating to k=50 there is no k=51) " 


model.ConstraintNumber1 = Constraint(model.Iset, model.Jset, model.Tset, model.Kset, rule=juv_advance)


# Constraint #2
def juv_total(model, i, j, t, k):
    return model.J_big[i,k,j] == sum(model.juvenille[i,j,t,k] for k in model.Kset)

model.ConstraintNumber2 = Constraint(model.Iset, model.Jset, model.Tset, rule=juv_total)

It is important that you not only define the rule of the constraint but also the constraint itself. Also you have to have in mind that you K and T sets are ending somewhere and that an expression of k+1 does not work if there is no k+1. Another point that could be mentioned is that if you start with k+1 == something the first k-value that is taken into account is k = 2.

I hope this helps, maybe someone also knows something smarter, I am also quite new to pyomo.

In regards to your comment about the "Big J", it could itself be a variable that is constructed from the small j's. It just wasn't clear to me whether one would call that a constraint or a variable. — GrayLiterature, Aug 14 '19 at 20:28
As far as I understood it you have to define every variable that exists in your "real mathematical model" (on paper) has to be constructed as a variable in pyomo first. Than in the second step you define via constraints how those variables are connected to each other. So basically you are setting up boundaries and also the "equations" from you mathematical model as constraints. So indeed J is a variable which is first constructed and then constrainted to take the value of the sum of j(i,j,k,t) over k. — PhilippK, Aug 14 '19 at 20:40
Cool yeah that is very helpful to think about. Thanks a lot! — GrayLiterature, Aug 14 '19 at 20:41

Pyomo: How to write a constraint for each (i,j) pair

1 Answers1