Using variable to index size of array, element of array and range of for loop
Answered
I am using Gurobi for an optimization model.
Within the optimization, I want to use 'Np' and 'interval', which depend on variable Ns, to define the length of an array, set the range of a for loop and select certain elements of an array. Below the relevant part of my code:
Ns =m.addVar(vtype=gp.GRB.INTEGER, name='Ns')
X = m.addVar(vtype=gp.GRB.CONTINUOUS, name='X') #1/Ns
Np = Nt*Ns
Npdiv = X/Nt #1/Np
interval = len(excesscumulative)*Npdiv
pickups = np.empty(Np, dtype=object)
for i in range(Np):
pickups[i] = excesscumulative[i*interval]
peakstorage = np.max(pickups)+ capacity
Some notes:
- Nt is a (constant) float
- X is equal to 1/Ns (using the constraint 'm.addConstr(Nships*X == 1)')
- 'excesscumulative' is an array of variables (elements depend on variables of the optimization)
- 'peakstorage' is later used as a minimum for the actual storage with constraint 'm.addConstr(Xst >= peakstorage)'
I want to use the variables Np and interval here in this way because this calculation must be part of my optimization; Np influences the interval and therefore the peak storage and the needed storage capacity. Optimizing Ns (which determines Np) against Xst (the needed storage capacity) is the main goal of this part of my optimization.
I am getting errors when trying to use Np and interval in the way I do here. How can I solve these errors while still making sure this calculation remains part of my optimization?
Thank you in advance for your help.
Kind regards,
Tycho
-
A small addition just to be sure. This is the error I am getting: TypeError: expected a sequence of integers or a single integer, got '<gurobi.LinExpr: 34.73041761101102 <gurobi.Var *Awaiting Model Update*>>'
0 -
Hi Tycho,
It is not possible to use Gurobi Var objects as indices because a Gurobi Var object does not have any value associated with it until the model is solved.
To model the relationship as the one you described, the common approach is to define auxiliary binary variables.
I did not fully understand your use-case. Could you explain what \(N_s\), \(N_p\), and \(\mbox{pickups}\) represent in English? This would help to see how you should define the auxiliary variables.
Best regards,
Maliheh
0 -
Hi Malileh,
Thank you for your response!
What I'm trying to model is a fuel production facility with storage capacity Xst and 'Ns' ships transporting the produced fuel to the usage location. As I mentioned, finding the (economic) optimum between Ns and Xst is the main goal here. If more ships are used, pickups are more frequent and less storage is needed.
Below an explanation of the used variables:
- Ns = number of ships
- X = 1/Ns (with constraint m.addConstr(Nships*X == 1))(needed because you cannot divide by a variable)
- Nt = number of trips a ship can do in a year (calculated earlier in the model primarily based on the distance and ship speed)
- Np = the total number of pickups done in a year (number of ships x number of trips per ship)
- Npdiv = 1/Np
- interval = the amount of elements between two pickups in the array giving the cumulative 'overproduction' over the whole year (so for example if interval=6, there will be a pickup at the instance corresponding to every 6th element in the array giving the cumulative 'overproduction' over the whole year ('excesscumulative'))('overproduction' in this case means more production in a period between two pickups then what could be transported in one ship)
- pickups = the cumulative 'overproduction' at the pickup moments
- peakstorage = the highest storage needed in the year (the biggest 'overproduction' at a pickup moment in the year + the capacity of a ship)
I hope this makes it more clear!
Kind regards,
Tycho
0 -
Hi Tycho,
- Let us define upper bounds \(U_s\) and \(U_p\) for the number of ships and total number of pickups. You would need to use your problem knowledge to come up with a simple heuristic to find an upper bound for the total number of ships required and then \(U_p = N_tU_s\).
- You can then define binary decision variables as below:
\[N_s = \begin{cases} 1, \mbox{if the number of ships is}~ s \\ 0, \mbox{otherwise} \end{cases}\]
\[N_p = \begin{cases} 1, \mbox{if the number of pickups is}~ p \\ 0, \mbox{otherwise} \end{cases}\]
\[\mbox{pickups}_p = \mbox{The cumulative overproduction if the number of pickups is}~ p\]
- Below are some of the constraints you need:
- \(\sum_{s=1}^{U_s} N_s = 1\): This constraint ensures a unique value for the number of ships.
- \(\sum_{p=1}^{U_p} N_p= 1\): This constraint ensures a unique value for the number of pickups.
- \(\sum_{p=1}^{U_p} p N_p = \sum_{s=1}^{U_s} N_t \times s \times N_s\): This constraint ensures the number of pickups equals the multiplication of number of ships and number of trips.
I think the binary variables \(N_s\) and \(N_p\) should allow modeling the rest of your constraints.
Best regards,
Maliheh
0 -
Hi Maliheh,
Thank you again for your response!
I think I understand how to code this now. Just be sure, I have two more questions:
1. First (just to be sure that was clear), the vector 'excesscumulative' has the overproduction on an hourly basis over time. 'Pickups' should also be a vector with values over time (but then in less time instances, i.e. only at the pickup moments), which would mean 'pickups' in the model you propose would be a matrix with a vector for every value of 'p'. Is that indeed what you meant?
2. For the third constraint you formulated, shouldn't the left side of the equation be 'SUM(Np*p)' with the sum from p=1 to Up?
0 -
Hi Maliheh,
One more question that came up regarding actually coding the constraints in Python (my apologies if this is a trivial question; it's my first time working with Gurobi):
3. If there is a total of 'Us' values for variable Ns (of which only one is equal to 1), how do I define this variable? Can I just define it as a binary variable (m.addVar(vtype=gp.GRB.BINARY, name='Nships')) or should it be an array of variables of some sort?
0 -
Hi Maliheh,
I tried to code what you explained today, but ran into some trouble when trying to formulate the final constraint stating that the storage must be bigger than the peak storage needed (Xst >= peakstorage)
Below I copied the code for the constraints that you formulated (including my suggested adjustment for the third constraint):
"""
m.addConstr(gp.quicksum(Ns[s] for s in range(1, Us + 1)) == 1, name='Constraint1')
m.addConstr(gp.quicksum(Np] for p in range(1, Up + 1)) == 1, name='Constraint2')
m.addConstr(gp.quicksum(Np[p] * p for p in range(1, Up+1)) == Nt * gp.quicksum(Ns[s] * s for s in range(1, Us + 1)), name=f'Constraint3_{p}')
"""
Below I have also copied the relevant piece of code doing the 'pre-calculations'. As you can see I try to create a matrix 'pickups' where every row is a vector pickups for a different number of ships.
"""
Nt =int(365*24/triptime) #the number of trips one ship can do in a year# Create binary decision variables Ns and Np
Us = 40 #max number of shipsif Us > (len(produced)/Nt):
print('Too many ships, maximum is',len(produced)/Nt)
sys.exit()Up = int(Nt*Us) #max number of pickups
Ns = np.empty(Us+1,dtype=object)
Np = np.empty(Up+1,dtype=object)
for s in range(1, Us + 1):
Ns[s] = m.addVar(vtype=GRB.BINARY, name=f'Ns_{s}')
for p in range(1, Up + 1):
Np[p] = m.addVar(vtype=GRB.BINARY, name=f'Np_{p}')excess = np.empty((Us+1,len(produced)), dtype=object)
excesscumulative = np.empty((Us+1,len(produced)), dtype=object)for i in range(1,Us+1):
excess[i,:] = produced - i*shipcapacity*Nt/(365*24) #produced - picked up hourly = excess
excesscumulative[i,:] = np.cumsum(excess[i,:]) #cumulative surplus
pickups = np.empty((Us+1,len(produced)), dtype=object)for i in range(1,Us):
interval = int(math.floor(len(excesscumulative)/(i*Nt)))
for k in range(i*Nt):
pickups[i,k] = excesscumulative[i,k*interval]"""
Then I want to have a constraint that picks the right vector from my matrix 'pickups' that corresponds to the chosen amount of ships. This is where I get stuck. Below some of the main constraints that I tried with the errors they gave me:
"""
1) m.addConstrs(Xst >= np.max(np.dot(Ns,pickups)+shipcapacity))
Error: TypeError: unsupported operand type(s) for *: 'NoneType' and 'NoneType'
2) m.addConstr(Xst-shipcapacity >= np.max(Ns[s] * pickups[s, :] for s in range(1, Us + 1)))Error: TypeError: unsupported operand type(s) for -: 'gurobipy.LinExpr' and 'generator'
3) m.addConstrs(Xst >= np.max(pickups[gp.quicksum(Ns[s]*s for s in range(1, Us + 1)),:])+shipcapacity)Error: IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices
"""
As everything that I tried is not working, I think I should approach this constraint in a different way but I do not really see how. Could you give me some advice on how to code this?
It would be greatly appreciated and thank you in advance!
Kind regards,
Tycho
0 -
Hi Maliheh,
Sorry for the many messages in a row, but I continued working on this and I believe I solved the issues with the last constraint and have something now that should be working, which I copied below. 'Produced' is the array with the hourly production volumes (which are 'gurobi.LinExpr').
"""
for j in range(0,len(produced+1)):
for s in range(Us+1):
m.addConstr(a[j] == Ns[s]*pickups[s,j])
alist = [a[j] for j in range(len(produced+1))]
m.addGenConstrMax(Xst, alist, name='maxconstr')"""
The other constraints stayed the same. The 'pre-calculations' did change slightly:
"""
Nt =int(365*24/triptime) #the number of trips one ship can do in a year
# Create binary decision variables Ns and Np
Us = 40 #max number of shipsif Us > (len(produced)/Nt):
print('Too many ships, maximum is',len(produced)/Nt)
sys.exit()Up = int(Nt*Us) #max number of pickups
Ns = np.empty(Us+1,dtype=object)
Np = np.empty(Up+1,dtype=object)
for s in range(1, Us + 1):
Ns[s] = m.addVar(vtype=GRB.BINARY, name=f'Ns_{s}')
for p in range(1, Up + 1):
Np[p] = m.addVar(vtype=GRB.BINARY, name=f'Np_{p}')excess = np.empty((Us+1,len(produced)), dtype=object)
excesscumulative = np.empty((Us+1,len(produced)), dtype=object)for i in range(1,Us+1):
excess[i,:] = produced - i*shipcapacity*Nt/(365*24) #produced - picked up hourly = excess
excesscumulative[i,:] = np.cumsum(excess[i,:]) #cumulative surplus
pickups = np.empty((Us+1,len(produced)), dtype=object)for i in range(1,Us):
interval = int(math.floor(len(produced)/(i*Nt)))
for k in range(i*Nt):
pickups[i,k] = excesscumulative[i,k*interval]a = m.addVars(len(produced),vtype=GRB.CONTINUOUS, name='a', lb=-20000)
"""
The code does not given any more errors now. However, when I look at the results, something weird still seems to be happening. The constraint I added last somehow influences the values of 'excesscumulative', 'excess' and 'produced'. If I see correctly, the solver seems to centralize production completely in the instances where the fuel is picked up by a ship. I have examined all the results, but it is not clear to me why the solver does this. If I take out that last constraint, I have 'normal' results again.
Do you maybe have an idea why this could be happening? If it is easier for you to view my whole code, I would be happy to share it in a private message.
I hope you are able to help with this issue!
Kind regards,
Tycho
0 -
Sorry, I missed your previous messages.
For the third constraint you formulated, shouldn't the left side of the equation be 'SUM(Np*p)' with the sum from p=1 to Up?
Yes, you are right. I edited my previous message.
I do not have a clear understanding of what you are trying to model because your model has so many pieces and I do not know all the input and variable definitions in detail.
I can give you a general approach to understand why Gurobi does not return your expected solution \(y\) as optimal and returns solution \(x\) instead.
Fix the decision variables to the values of solution \(y\) and optimize the model. Two scenarios can happen:
- The model is feasible but the objective value of solution \(y\) is worse than the objective value of solution \(x\). So, it is clear why Gurobi returned the solution \(x\) as optimal.
- The model is infeasible meaning that your expected solution \(y\) does not satisfy some of the constraints. You can then follow the steps in the article "How do I determine why my model is infeasible?" to figure out which constraints are responsible for the infeasibility of solution \(y\).
For me to be able to help further, I need to know what you are trying to model. Is there any reference paper that you are trying to formulate your model based on? If so, could you please refer me to that paper?
Best regards,
Maliheh
0 -
Hi Maliheh,
Thank you for your message. I checked the model status and the model is feasible; the problem is not that the solver does not find a solution, but it is impossible that the solution it is finding is realistic/the optimal one. This means something with that last constraint I added must be causing this, but I cannot seem to figure out what.
To help in understanding, I plotted some figures. Below you can see the plots of 'excess', 'excesscumulative' and 'produced' the way they should look like (and when that last constraint is not in the model they look like this). Ofcourse they would differ based on the optimization inputs/results, but they should like somewhat like this (except for the fact that excesscumulative should end at zero but I still have to implement a constraint for that).
Below you can see the same figures, but now with the last constraint implemented in the code. As you can see, for some reason the solver now decides to only produce when there is a pickup and then do the entire needed production at once. This is causing the production capacity to be unrealistically high (about 40x times higher than without the constraint in this case)(the production 'machines' are also variables in the optimization which determine the amount of produced fuel per hour (aka the elements of array 'produced') together with some other input parameters).
As all production is done in the pickup moments in this optimization, 'Xst' (the needed storage for excess production) and therefore also the values of 'a' remain zero. 'Xst and 'a' are defined as follows (I added lb and ub just to make sure there are no default boundaries):
m.addVar(vtype=gp.GRB.CONTINUOUS, name='Xstorage',lb=-20000,ub=20000)
a = m.addVars(len(produced),vtype=GRB.CONTINUOUS, name='a', lb=-20000,ub=20000)
The fact that the solver is keeping these at zero must be the issue, as this causes the array 'pickups' to consist of zeros and because of that also the elements of 'excesscumulative' must be zero. It is however not clear to me why the solver does this. What I then also do not understand is why that leads to the 'excess' graph as shown above, because I have defined the relationship between 'excess' and 'excesscumulative' as follows:
excesscumulative[i,:] = np.cumsum(excess[i,:])
As far as I see, this would mean excess should also consist of zeros for excesscumulative to consist of zeros, instead of the spikes that can be seen now. That would then lead to a constant production, as 'excess' is defined as follows:
excess[i,:] = produced - s*shipcapacity*Nt/(365*24)
(where s is the number of ships so we substract the average hourly pickup from the hourly production to find the excess)
The only reason I could think of for the solver keeping Xst and therefore 'a' at zero is that the storage would be incredibly expensive (as the objective of the model is to minimize the total costs), but the same happens when I put the costs of storage at zero.
Lastly, the spikes at the end in the second group of graphs are also surprising.
Long story short, I believe the issues are caused by the solver keeping either 'Xst' or 'a' at zero, which is causing some peculiar results. I think the main question to figure out is why the solver is doing that when that last constraint is implemented.
I unfortunately do not have a reference paper that I am directly basing my model upon. My apologies for the complicated questions, I would love to have a call to explain/show you the model if that is easier, but I'm not sure if that is an option.
Kind regards,
Tycho
0 -
Thanks for the further clarification.
Sorry, we do not offer one-on-one support to our academic users beyond licensing and installation issues. Even in-depth modeling support is not part of our typical support offering to our commercial users simply because modeling requires a complete and in-depth understanding of the application domain.
Here are my final comments:
- If the \(\texttt{excesscumulative}_t\) decision variable is defined to be the sum of the \(\texttt{excess}_{t^{\prime}}\) variables over all \(t^{\prime} \leq t\), it is indeed strange that the \(\texttt{excesscumulative}\) figure remains at zero despite having spikes at the \(\texttt{excess}\) figure. I would suggest checking the implementation of this constraint or maybe there is an error in your plotting implementation.
- You can fix the decision variables \(Xst\) and \(a\) at values that you expect to see (the values achieved before adding the last constraint) and optimize the model again. This can hopefully give you some insights.
- To fix the variable \(x\) at a value \(a\), you can add a new constraint in the form \(x = a\) or fix the lower and upper bounds of the variable \(x\) at \(a\).
Best regards,
Maliheh
0
Please sign in to leave a comment.
Comments
11 comments