Different platforms and different results
Answered
Hi,
I wrote a problem to be solved by Gurobi (Python), before increasing a parameter, I got correct results by Gurobi and another platform, when I increase a parameter I will get infeasible problem in Gurobi, while another platform gave me optimal solution and it should be feasible. I have academic license for Gurobi, too. Could someone help me what is wrong with Gurobi solver?
-
Hello Arash,
It is possible that your model is on the verge of (in)feasibility. This is often caused by shaky numerics. I would recommend having a look at our Guidelines for Numerical Issues. You should also try experimenting with the NumericFocus parameter and try using a lower FeasibilityTol and IntFeasTol. You could try running the model with different Seed parameter values on the two machines to see whether it behaves differently. You could also share the model. Note that uploading files in the Community Forum is not possible but we discuss an alternative in Posting to the Community Forum.
The Knowledge Base article Why does Gurobi perform differently on different machines? might also be of interest to you.
Best regards,
Jaromił0 -
Thanks for your reply,
from gurobipy import *
import numpy as np
from math import *
# indices
myList = list(range(0, 4))
vctr = np.array(myList)
I = vctr # existing generator
myList = list(range(0, 13))
vctr = np.array(myList)
NI = vctr # new generator
myList = list(range(0, 3))
vctr = np.array(myList)
J = vctr # load
myList = list(range(0, 7))
vctr = np.array(myList)
EL = vctr # existing line
myList = list(range(0, 7))
vctr = np.array(myList)
NL = vctr # new line
myList = list(range(0, 6))
vctr = np.array(myList)
N = vctr # bus
myList = list(range(0, 5))
vctr = np.array(myList)
T = vctr # period
# parameters
cpe = np.array([10, 5, 5, 10]) # capacity of existing generators
cpn = np.array([10, 7, 5, 3, 3, 3, 2, 5, 3, 10, 8, 5, 2]) # capacity of new generators
ce = np.array([25, 35, 37, 25]) # cost of existing generator
cn = np.array(
[22, 30, 35, 40, 40, 40, 55, 35, 40, 22, 29, 33, 55]
) # cost of new generator
ccg = np.array(
[100e3, 80e3, 60e3, 30e3, 40e3, 45e3, 20e3, 70e3, 35e3, 110e3, 85e3, 50e3, 15e3]
) # capital cost for new generators($/kw)
ccl = np.array(
[5e3, 8e3, 12e3, 10e3, 7e3, 5e3, 6e3]
) # capital cost for new lines($/kw)
line, X = multidict(
{
(0, 1): 0.229,
(1, 0): 0.229,
(1, 2): 0.229,
(2, 1): 0.229,
(0, 3): 0.229,
(3, 0): 0.229,
(1, 3): 0.229,
(3, 1): 0.229,
(3, 4): 0.229,
(4, 3): 0.229,
(4, 5): 0.229,
(5, 4): 0.229,
(2, 5): 0.229,
(5, 2): 0.229,
}
) # reactance
LF = np.array([0.5, 0.65, 0.8, 0.9, 1]) # load factor
weight = np.array([0.4, 0.3, 0.3]) # load weight
du = np.array([1510, 2800, 2720, 1120, 610]) # duration of each period
fmaxel = np.array([10, 7, 7, 7, 7, 7, 7]) # maximum flow of existing lines
fmaxnl = np.array([10, 7, 7, 7, 7, 7, 7]) # maximum flow of new lines
wind = np.array([0, 0, 2, 2, 0, 0])
PV = np.array([0, 0, 2, 2, 0, 0])
EV = np.array([0, 0, 1, 1, 0, 0])
f = np.array([0, 0, 1, 1, 1, 0])
peak = 55
# map(i,n)
I1 = [0, 1, 2, 3]
N1 = [1, 2, 5, 0]
# map1(el,n,m)
EL1 = [0, 1, 2, 3, 4, 5, 6]
N2 = [0, 1, 0, 1, 3, 4, 2]
M2 = [1, 2, 3, 3, 4, 5, 5]
# map2(in,n)
NI1 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
M3 = [0, 0, 1, 1, 3, 2, 2, 4, 4, 5, 5, 5, 5]
# loc-d(j,n)
J1 = [0, 1, 2]
N3 = [2, 3, 4]
# map3(nl,n,m)
NL1 = [0, 1, 2, 3, 4, 5, 6]
M4 = [0, 1, 0, 1, 3, 4, 2]
M5 = [1, 2, 3, 3, 4, 5, 5]
ss = 0
m = Model()
# variables
ge = m.addVars(
I, T, vtype=GRB.CONTINUOUS, lb=ss, name="ge"
) # output power of existing generators
gge = m.addVars(
I, N, T, vtype=GRB.CONTINUOUS, lb=ss, name="gge"
) # output power of existing generators
gn = m.addVars(
NI, T, vtype=GRB.CONTINUOUS, lb=ss, name="gn"
) # output power of new generators
ggn = m.addVars(
NI, N, T, vtype=GRB.CONTINUOUS, lb=ss, name="ggn"
) # output power of new generators
cost_inv = m.addVar(vtype=GRB.CONTINUOUS, lb=ss, name="cost_inv") # investment cost
cost_gen = m.addVar(vtype=GRB.CONTINUOUS, lb=ss, name="cost_gen") # generation cost
theta = m.addVars(N, T, vtype=GRB.CONTINUOUS, name="theta") # angle of voltage
fel = m.addVars(EL, T, vtype=GRB.CONTINUOUS, name="fel") # powe flow of existing lines
fnl = m.addVars(NL, T, vtype=GRB.CONTINUOUS, name="fnl") # powe flow of new lines
fell = m.addVars(
N, N, T, vtype=GRB.CONTINUOUS, name="fell"
) # powe flow of existing lines
fnll = m.addVars(N, N, T, vtype=GRB.CONTINUOUS, name="fnll") # powe flow of new lines
Pl = m.addVars(J, N, T, vtype=GRB.CONTINUOUS, name="Pl") # Demand
unl = m.addVars(NL, vtype=GRB.BINARY, name="unl") # binary of new lines
ung = m.addVars(NI, vtype=GRB.BINARY, name="ung") # binary of new generators
objf = m.addVar(vtype=GRB.CONTINUOUS, name="objf") # objective
# constraints
c1 = m.addConstr(
cost_inv
== (
sum(unl[nl] * ccl[nl] * fmaxnl[nl] for nl in NL)
+ sum(ung[ni] * ccg[ni] * cpn[ni] for ni in NI)
)
)
c2 = m.addConstr(
cost_gen
== (
sum(ge[i, t] * ce[i] * du[t] for i in I for t in T)
+ sum(gn[ni, t] * cn[ni] * du[t] for ni in NI for t in T)
)
)
c3 = m.addConstr(objf == cost_inv + cost_gen)
c4 = m.addConstrs(ge[i, t] <= cpe[i] for i in I for t in T)
c5 = m.addConstrs(gn[ni, t] <= ung[ni] * cpn[ni] for ni in NI for t in T)
c7 = m.addConstrs(
fel[l, t] == (1 / X[i, j]) * (theta[i, t] - theta[j, t]) * 100
for t in T
for l, i, j in zip(EL1, N2, M2)
)
c9 = m.addConstrs(fel[l, t] == fell[i, j, t] for t in T for l, i, j in zip(EL1, N2, M2))
c99 = m.addConstrs(
fell[i, j, t] == 0 for i in N for j in N for t in T if (i, j) not in zip(N2, M2)
)
c10 = m.addConstrs(
fnl[nl, t] == fnll[i, j, t] for t in T for nl, i, j in zip(NL1, M4, M5)
)
c100 = m.addConstrs(
fnll[i, j, t] == 0 for i in N for j in N for t in T if (i, j) not in zip(M4, M5)
)
c11 = m.addConstrs(ge[i, t] == gge[i, n, t] for i, n in zip(I1, N1) for t in T)
c12 = m.addConstrs(
gge[i, n, t] == 0 for i in I for n in N for t in T if (i, n) not in zip(I1, N1)
)
c13 = m.addConstrs(gn[ni, t] == ggn[ni, n, t] for ni, n in zip(NI1, M3) for t in T)
c14 = m.addConstrs(
ggn[ni, n, t] == 0 for ni in NI for n in N for t in T if (ni, n) not in zip(NI1, M3)
)
c15 = m.addConstrs(
Pl[j, n, t] == peak * LF[t] * weight[j] for j, n in zip(J1, N3) for t in T
)
c16 = m.addConstrs(
Pl[j, n, t] == 0 for j in J for n in N for t in T if (j, n) not in zip(J1, N3)
)
c17 = m.addConstrs(
gge.sum("*", n, t)
+ ggn.sum("*", n, t)
- fell.sum(n, "*", t)
+ fell.sum("*", n, t)
- fnll.sum(n, "*", t)
+ fnll.sum("*", n, t)
>= (Pl.sum("*", n, t) + EV[n] - PV[n] - wind[n])
- 0.05 * (Pl.sum("*", n, t) + EV[n] - PV[n] - wind[n])
+ 0.05 * 6 * (Pl.sum("*", n, t) + EV[n] - PV[n] - wind[n])
for n in N
for t in T
)
c18 = m.addConstrs(fel[l, t] <= fmaxel[l] for l in EL for t in T)
c19 = m.addConstrs(fel[l, t] >= -fmaxel[l] for l in EL for t in T)
c20 = m.addConstrs(fnl[nl, t] <= fmaxnl[nl] * unl[nl] for nl in NL for t in T)
c21 = m.addConstrs(fnl[nl, t] >= -fmaxnl[nl] * unl[nl] for nl in NL for t in T)
c22 = m.addConstrs(
fnl[nl, t] * X[i, j] - (theta[i, t] - theta[j, t]) * 100 >= -(1 - unl[nl]) * 1e6
for t in T
for nl, i, j in zip(NL1, M4, M5)
)
c23 = m.addConstrs(
fnl[nl, t] * X[i, j] - (theta[i, t] - theta[j, t]) * 100 <= (1 - unl[nl]) * 1e6
for t in T
for nl, i, j in zip(NL1, M4, M5)
)
m.setObjective(objf, GRB.MINIMIZE)
# m.params.NonConvex = 2
# opt.options['NonConvex'] = 2
# orignumvars = m.NumVars
# m.feasRelaxS(0, True, False, True)
m.optimize()
# m.optimize()
# m.computeIIS()
m.printAttr("x")
print("runtime is", m.Runtime)This is my model and it should be feasible, I have checked FeasibilityTol ,but I do not know how to use it. when I increase parameter ''peak'' in my model, it will be infeasible. I appreciate your help.
0 -
First, you should give all your constraints unique names to make it easier to distinguish them. For example
c15 = m.addConstrs(
(Pl[j, n, t] == peak * LF[t] * weight[j] for j, n in zip(J1, N3) for t in T),
name ="c15")You can then compute an IIS as described in How do I determine why my model is infeasible?
m.optimize()
m.computeIIS()
m.write("iis.ilp")You can open the file \(\texttt{iis.ilp}\) in any standard text editor and see which subset of your constraints makes the model infeasible. For \(\texttt{peak=55}\), you can see that constraints \(\texttt{c15}\) fix the values of your \(\texttt{Pl}\) variables to
c15[0,2,4]: Pl[0,2,4] = 22
c15[1,3,4]: Pl[1,3,4] = 16.5
c15[2,4,4]: Pl[2,4,4] = 16.5which is somehow responsible for the infeasibility of your model. With the IIS, you can try to find out why your model is infeasible.
Best regards,
Jaromił0 -
Hi,
I have following code , but I do not know when I increase ''peak'' model is infeasible, while it should be feasible. For peak=10, I got optimal solution, but for more than 10, it will be infeasible that is wrong. why Gurobi does not work here, while for peaks more than 10 , I got optimal solution in another platform which is correct. I decreased FeasibilityTol and IntFeasTol, but still infeasible, and checking "iis.ilp" says that c17 constraint has problem, while it does not make sense.
from gurobipy import *
import numpy as np
#indices
myList = list(range(0,4))
vctr=np.array(myList)
I=vctr# generator
myList = list(range(0,3))
vctr=np.array(myList)
J=vctr#load
myList = list(range(0,7))
vctr=np.array(myList)
EL=vctr# line
myList = list(range(0,6))
vctr=np.array(myList)
N=vctr#bus
myList = list(range(0,5))
vctr=np.array(myList)
T=vctr#period
#parameters
cpe=np.array([10,5,5,10])# capacity of generators
ce=np.array([25,35,37,25])# cost of generator
line,X=multidict({(0,1):0.229,(1,0):0.229,(1,2):0.229,(2,1):0.229,(0,3):0.229,(3,0):0.229,(1,3):0.229,(3,1):0.229,(3,4):0.229,(4,3):0.229,(4,5):0.229,(5,4):0.229,(2,5):0.229,(5,2):0.229})#reactance
LF=np.array([0.5,0.65,0.8,0.9,1])#load factor
weight=np.array([0.4,0.3,0.3])#load weight
du=np.array([1510,2800,2720,1120,610])#duration of each period
fmaxel=np.array([10,7,7,7,7,7,7])#maximum flow of lines
peak=15
#map(i,n)
I1=[0,1,2,3]
N1=[1,2,5,0]
#map1(el,n,m)
EL1=[0,1,2,3,4,5,6]
N2=[0,1,0,1,3,4,2]
M2=[1,2,3,3,4,5,5]
#loc-d(j,n)
J1=[0,1,2]
N3=[2,3,4]
ss=0
m=Model()
m.setParam('FeasibilityTol', 10e-9)
m.setParam('IntFeasTol', 10e-9)
#variables
ge=m.addVars(I,T,vtype=GRB. CONTINUOUS,lb=ss,name='ge') #output power of generators
gge=m.addVars(I,N,T,vtype=GRB. CONTINUOUS,lb=ss,name='gge') #output power of generators
cost_gen=m.addVar(vtype=GRB. CONTINUOUS,lb=ss,name='cost_gen') #generation cost
theta=m.addVars(N,T,vtype=GRB. CONTINUOUS,name='theta') #angle of voltage
fel=m.addVars(EL,T,vtype=GRB. CONTINUOUS,name='fel') #power flow of existing lines
fell=m.addVars(N,N,T,vtype=GRB. CONTINUOUS,name='fell') #power flow of existing lines
Pl=m.addVars(J,N,T,vtype=GRB. CONTINUOUS,name='Pl') #Demand
objf=m.addVar(vtype=GRB. CONTINUOUS,name='objf') #objective
#constraints
c2=m.addConstr(cost_gen==(sum(ge[i,t]*ce[i]*du[t] for i in I for t in T)),name="c2")
c3=m.addConstr(objf==cost_gen,name="c3")
c4=m.addConstrs((ge[i,t]<= cpe[i] for i in I for t in T),name="c4")
c7=m.addConstrs((fel[l,t]==(1/X[i,j])*(theta[i,t]-theta[j,t])*100 for t in T for l,i,j in zip(EL1,N2,M2)),name="c7" )
c9=m.addConstrs((fel[l,t]==fell[i,j,t] for t in T for l,i,j in zip(EL1,N2,M2)),name="c9" )
c99=m.addConstrs((fell[i,j,t]==0 for i in N for j in N for t in T if (i,j) not in zip(N2,M2)),name="c99" )
c11=m.addConstrs((ge[i,t]==gge[i,n,t] for i,n in zip(I1,N1) for t in T),name="c11")
c12=m.addConstrs((gge[i,n,t]==0 for i in I for n in N for t in T if (i,n) not in zip(I1,N1)),name="c12" )
c15=m.addConstrs((Pl[j,n,t] ==peak*LF[t]*weight[j]for j,n in zip(J1,N3) for t in T),name="c15")
c16=m.addConstrs((Pl[j,n,t]==0 for j in J for n in N for t in T if (j,n) not in zip(J1,N3)),name="c16" )
c17=m.addConstrs(((gge.sum('*',n,t)-fell.sum(n,'*',t)+fell.sum('*',n,t)==(Pl.sum('*',n,t) ))for n in N for t in T),name="c17")
c18=m.addConstrs((fel[l,t]<=fmaxel[l] for l in EL for t in T),name="c18")
c19=m.addConstrs((fel[l,t]>=-fmaxel[l] for l in EL for t in T),name="c19")
m.setObjective(objf, GRB. MINIMIZE)
#orignumvars = m.NumVars
#m.feasRelaxS(0, True, False, True)
#m.computeIIS()
#m.write("iis.ilp")
m.optimize()
#m.optimize()
m.printAttr('x')
print ('runtime is',m.Runtime)0 -
When you look at the IIS generated by Gurobi, you will see that most variables are fixed to 0 or some other value. If you substitute these values and use all constraints to substitute left over variables, you will end up with
c4[3,4]: ge[3,4] <= 10
c17[0,4]: ge[3,4] - 873.3624454148472 theta[0,4] + 436.6812227074236 theta[1,4] + 436.6812227074236 theta[3,4] = 0
c17[2,4]: ge[1,4] + 436.6812227074236 theta[1,4] - 873.3624454148472 theta[2,4] + 436.6812227074236 theta[5,4] - 6 = 0
c17[3,4]: 436.6812227074236 theta[0,4] - 1310.043668122271 theta[3,4] + 436.6812227074236 theta[1,4] + 436.6812227074236 theta[4,4] - 4.5 = 0
c17[4,4]: 436.6812227074236 theta[3,4] - 873.3624454148472 theta[4,4] + 436.6812227074236 theta[5,4] - 4.5 = 0
c17[5,4]: ge[2,4] + 436.6812227074236 theta[2,4] - 873.3624454148472 theta[5,4] + 436.6812227074236 theta[4,4] = 0Now, you have to figure out, why this set of constraints is infeasible. Note that the values 6 and 4.5 come from the fixings
c15[0,2,4]: Pl[0,2,4] = 6
c15[1,3,4]: Pl[1,3,4] = 4.5
c15[2,4,4]: Pl[2,4,4] = 4.5You should probably try to substitute the theta variables, i.e., solve the equality constraint for, e.g., theta[0,4], and substitute it in other constraints. Then proceed until you cannot substitute anymore.
0 -
Thank you, I figured out the issue. I simplified code and removed all fixed variables.
from gurobipy import *
import numpy as np
#indices
myList = list(range(0,4))
vctr=np.array(myList)
I=vctr# generator
myList = list(range(0,3))
vctr=np.array(myList)
J=vctr#load
myList = list(range(0,7))
vctr=np.array(myList)
EL=vctr# line
myList = list(range(0,6))
vctr=np.array(myList)
N=vctr#bus
myList = list(range(0,5))
vctr=np.array(myList)
T=vctr#period
#parameters
cpe=np.array([10,5,5,10])# capacity of generators
ce=np.array([25,35,37,25])# cost of generator
line,X=multidict({(0,1):0.229,(1,0):0.229,(1,2):0.229,(2,1):0.229,(0,3):0.229,(3,0):0.229,(1,3):0.229,(3,1):0.229,(3,4):0.229,(4,3):0.229,(4,5):0.229,(5,4):0.229,(2,5):0.229,(5,2):0.229})#reactance
LF=np.array([0.5,0.65,0.8,0.9,1])#load factor
weight=np.array([0.4,0.3,0.3])#load weight
du=np.array([1510,2800,2720,1120,610])#duration of each period
fmaxel=np.array([10,7,7,7,7,7,7])#maximum flow of lines
peak=15
#map(i,n)
I1=[0,1,2,3]
N1=[1,2,5,0]
#map1(el,n,m)
EL1=[0,1,2,3,4,5,6]
N2=[0,1,0,1,3,4,2]
M2=[1,2,3,3,4,5,5]
#loc-d(j,n)
J1=[0,1,2]
N3=[2,3,4]
ss=0
m=Model()
#variables
ge=m.addVars(I,T,vtype=GRB. CONTINUOUS,lb=ss,name='ge') #output power of generators
cost_gen=m.addVar(vtype=GRB. CONTINUOUS,lb=ss,name='cost_gen') #generation cost
theta=m.addVars(N,T,vtype=GRB. CONTINUOUS,name='theta') #angle of voltage
fel=m.addVars(EL,T,vtype=GRB. CONTINUOUS,name='fel') #power flow of existing lines
objf=m.addVar(vtype=GRB. CONTINUOUS,name='objf') #objective
#constraints
c1=m.addConstr(cost_gen==(sum(ge[i,t]*ce[i]*du[t] for i in I for t in T)),name="c1")
c2=m.addConstr(objf==cost_gen,name="c2")
c3=m.addConstrs((ge[i,t]<= cpe[i] for i in I for t in T),name="c3")
c4=m.addConstrs((fel[l,t]==(1/X[n,m])*(theta[n,t]-theta[m,t])*100 for t in T for (l,n,m) in zip(EL1,N2,M2)),name="c4" )
c5=m.addConstrs((sum(ge[i,t] for i in I if (i,n) in zip(I1,N1))-sum(fel[l,t] for l in EL if (l,n) in zip(EL1,N2))+sum(fel[l,t] for l in EL if (l,n) in zip(EL1,M2))==(sum(weight[j]*peak*LF[t] for j in J if (j,n) in zip(J1,N3)))for n in N for t in T),name="c5")
c6=m.addConstrs((fel[l,t]<=fmaxel[l] for l in EL for t in T),name="c6")
c7=m.addConstrs((fel[l,t]>=-fmaxel[l] for l in EL for t in T),name="c7")
m.setObjective(objf, GRB. MINIMIZE)
m.optimize()
m.printAttr('x')Issue is related to 'c4' and 'c5'. I want to define that :
c4=m.addConstrs((fel[l,t]==(1/X[n,m])*(theta[n,t]-theta[m,t])*100 for t in T for (l,n,m) in zip(EL1,N2,M2)),name="c4" )---> from n to m
also it ca be:
c55=m.addConstrs((fel[l,t]==(1/X[m,n])*(theta[m,t]-theta[n,t])*100 for t in T for (l,m,n) in zip(EL1,M2,N2)),name="c55" )--->from m to n
and based on 'c6' and 'c7', fel[l,t] can be positive or negative. In 'c5', i am going to define for each n and t:
summation of ge[i,t] for all i related to n (zip(I1,N1))+summation of fel[l,t] for all l related to n(zip(EL1,N2) or zip(EL1,M2) n and m define same thing )(if l from n to m, take fel[l,t] negative. if l from m to n, take fel[l,t] positive) ==summation (peak*LF[t]*weight[j] for all j related to n(zip(J1,N3))
How can i define these constraints?
Based on my constraints for 'c4' and 'c5', I only get fel[l,t]>=0 in my results, while it can be negative too, but it considers only ''from n to m'' for fel[l,t], so we do not have ''from m to n'' for fel[l,t] and that is why it will be infeasible.
0 -
Hi Arash,
Good job on identifying the source of infeasibility.
I think that your constraints are correct. You state that \(\texttt{fel}\) can and should possibly take negative values. Please note that the default lower bound for variables in Gurobi is \(0\). Thus, to allow for negative \(\texttt{fel}\) values, you have to explicitly state so when creating the variables.
fel=m.addVars(EL,T,lb=-GRB.INFINITY,vtype=GRB. CONTINUOUS,name='fel') #power flow of existing lines
The above change makes your model feasible.
Best regards,
Jaromił0 -
Thanks a lot for your help.
0
Please sign in to leave a comment.
Comments
8 comments