Search code examples
optimizationpyomooperations-research

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

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:

enter image description here

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.

Solution

  • 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.