•  Gurobi Staff

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ł

•    from gurobipy import *import numpy as npfrom math import *# indicesmyList = list(range(0, 4))vctr = np.array(myList)I = vctr  # existing generatormyList = list(range(0, 13))vctr = np.array(myList)NI = vctr  # new generatormyList = list(range(0, 3))vctr = np.array(myList)J = vctr  # loadmyList = list(range(0, 7))vctr = np.array(myList)EL = vctr  # existing linemyList = list(range(0, 7))vctr = np.array(myList)NL = vctr  # new linemyList = list(range(0, 6))vctr = np.array(myList)N = vctr  # busmyList = list(range(0, 5))vctr = np.array(myList)T = vctr  # period# parameterscpe = np.array([10, 5, 5, 10])  # capacity of existing generatorscpn = np.array([10, 7, 5, 3, 3, 3, 2, 5, 3, 10, 8, 5, 2])  # capacity of new generatorsce = np.array([25, 35, 37, 25])  # cost of existing generatorcn = np.array(    [22, 30, 35, 40, 40, 40, 55, 35, 40, 22, 29, 33, 55])  # cost of new generatorccg = 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,    })  # reactanceLF = np.array([0.5, 0.65, 0.8, 0.9, 1])  # load factorweight = np.array([0.4, 0.3, 0.3])  # load weightdu = np.array([1510, 2800, 2720, 1120, 610])  # duration of each periodfmaxel = np.array([10, 7, 7, 7, 7, 7, 7])  # maximum flow of existing linesfmaxnl = np.array([10, 7, 7, 7, 7, 7, 7])  # maximum flow of new lineswind = 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 = 0m = Model()# variablesge = m.addVars(    I, T, vtype=GRB.CONTINUOUS, lb=ss, name="ge")  # output power of existing generatorsgge = m.addVars(    I, N, T, vtype=GRB.CONTINUOUS, lb=ss, name="gge")  # output power of existing generatorsgn = m.addVars(    NI, T, vtype=GRB.CONTINUOUS, lb=ss, name="gn")  # output power of new generatorsggn = m.addVars(    NI, N, T, vtype=GRB.CONTINUOUS, lb=ss, name="ggn")  # output power of new generatorscost_inv = m.addVar(vtype=GRB.CONTINUOUS, lb=ss, name="cost_inv")  # investment costcost_gen = m.addVar(vtype=GRB.CONTINUOUS, lb=ss, name="cost_gen")  # generation costtheta = m.addVars(N, T, vtype=GRB.CONTINUOUS, name="theta")  # angle of voltagefel = m.addVars(EL, T, vtype=GRB.CONTINUOUS, name="fel")  # powe flow of existing linesfnl = m.addVars(NL, T, vtype=GRB.CONTINUOUS, name="fnl")  # powe flow of new linesfell = m.addVars(    N, N, T, vtype=GRB.CONTINUOUS, name="fell")  # powe flow of existing linesfnll = m.addVars(N, N, T, vtype=GRB.CONTINUOUS, name="fnll")  # powe flow of new linesPl = m.addVars(J, N, T, vtype=GRB.CONTINUOUS, name="Pl")  # Demandunl = m.addVars(NL, vtype=GRB.BINARY, name="unl")  # binary of new linesung = m.addVars(NI, vtype=GRB.BINARY, name="ung")  # binary of new generatorsobjf = m.addVar(vtype=GRB.CONTINUOUS, name="objf")  # objective# constraintsc1 = 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.

•  Gurobi Staff

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] = 22c15[1,3,4]: Pl[1,3,4] = 16.5c15[2,4,4]: Pl[2,4,4] = 16.5

which 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ł

•    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#indicesmyList = list(range(0,4))vctr=np.array(myList)I=vctr# generatormyList = list(range(0,3))vctr=np.array(myList)J=vctr#loadmyList = list(range(0,7))vctr=np.array(myList)EL=vctr# linemyList = list(range(0,6))vctr=np.array(myList)N=vctr#busmyList = list(range(0,5))vctr=np.array(myList)T=vctr#period#parameterscpe=np.array([10,5,5,10])# capacity of generatorsce=np.array([25,35,37,25])# cost of  generatorline,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})#reactanceLF=np.array([0.5,0.65,0.8,0.9,1])#load factorweight=np.array([0.4,0.3,0.3])#load weightdu=np.array([1510,2800,2720,1120,610])#duration of each periodfmaxel=np.array([10,7,7,7,7,7,7])#maximum flow of  linespeak=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=0m=Model()m.setParam('FeasibilityTol', 10e-9)m.setParam('IntFeasTol', 10e-9)#variablesge=m.addVars(I,T,vtype=GRB. CONTINUOUS,lb=ss,name='ge') #output power of  generatorsgge=m.addVars(I,N,T,vtype=GRB. CONTINUOUS,lb=ss,name='gge') #output power of  generatorscost_gen=m.addVar(vtype=GRB. CONTINUOUS,lb=ss,name='cost_gen') #generation costtheta=m.addVars(N,T,vtype=GRB. CONTINUOUS,name='theta') #angle of voltagefel=m.addVars(EL,T,vtype=GRB. CONTINUOUS,name='fel') #power flow of existing linesfell=m.addVars(N,N,T,vtype=GRB. CONTINUOUS,name='fell') #power flow of existing linesPl=m.addVars(J,N,T,vtype=GRB. CONTINUOUS,name='Pl') #Demandobjf=m.addVar(vtype=GRB. CONTINUOUS,name='objf') #objective#constraintsc2=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)
•  Gurobi Staff

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] = 0

Now, 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] = 6c15[1,3,4]: Pl[1,3,4] = 4.5c15[2,4,4]: Pl[2,4,4] = 4.5

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

•    Thank you, I figured out the issue. I simplified code and removed all fixed variables.

from gurobipy import *import numpy as np#indicesmyList = list(range(0,4))vctr=np.array(myList)I=vctr# generatormyList = list(range(0,3))vctr=np.array(myList)J=vctr#loadmyList = list(range(0,7))vctr=np.array(myList)EL=vctr# linemyList = list(range(0,6))vctr=np.array(myList)N=vctr#busmyList = list(range(0,5))vctr=np.array(myList)T=vctr#period#parameterscpe=np.array([10,5,5,10])# capacity of generatorsce=np.array([25,35,37,25])# cost of  generatorline,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})#reactanceLF=np.array([0.5,0.65,0.8,0.9,1])#load factorweight=np.array([0.4,0.3,0.3])#load weightdu=np.array([1510,2800,2720,1120,610])#duration of each periodfmaxel=np.array([10,7,7,7,7,7,7])#maximum flow of  linespeak=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=0m=Model()#variablesge=m.addVars(I,T,vtype=GRB. CONTINUOUS,lb=ss,name='ge') #output power of  generatorscost_gen=m.addVar(vtype=GRB. CONTINUOUS,lb=ss,name='cost_gen') #generation costtheta=m.addVars(N,T,vtype=GRB. CONTINUOUS,name='theta') #angle of voltagefel=m.addVars(EL,T,vtype=GRB. CONTINUOUS,name='fel') #power flow of existing linesobjf=m.addVar(vtype=GRB. CONTINUOUS,name='objf') #objective#constraintsc1=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.

•  Gurobi Staff

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ł

•    Thanks a lot for your help.