Optimization problem

Hi! I have to maximize an expression that is the sum of many sub-expressions (in the code they are called espressione ) that have some constraints like even their variables contained in the vector var[].

When i try to run the code i have this error:

Solution count 0

Model is infeasible or unbounded

But if i print GRB.OPTIMAL it is 2!. And if i try to do:

fob.computeIIS()
fob.write("model.ilp")

after optimize() my model, i recieve:

Computing Irreducible Inconsistent Subsystem (IIS)...

without having a solution

Here the code:

d=0
var=[]
k=1
indice=0
s=0
vincoli=[]
n=[]
vettore_finale=[]

fob=Model()
fob.params.NonConvex=2
for h in range (0,int(matrice1[i][j])-1):
if h==0:
vettore_somma=[]
for x in range (0,int(matrice1[i][j])-1):

   var.append(fob.addVar(vtype='C',name=str(matrice3[i][j][17][x])))
   vincoli.append(fob.addConstr(var[x+indice]<=vmax)
   vincoli.append(fob.addConstr(var[x+indice]>=vmin))
   var.append(fob.addVar(vtype='C',name=str(matrice3[i][j][18][x])))
   vincoli.append(fob.addConstr(var[x+1+indice]<=sigmamax))
   vincoli.append(fob.addConstr(var[x+1+indice]-sigma<=(sigmamax)))
   vincoli.append(fob.addConstr(-var[x+1+indice]+sigma<=(sigmamax)))
   var.append(fob.addVar(vtype='C',name=str(matrice3[i][j][19][x])))
   vincoli.append(fob.addConstr(var[x+2+indice]==var[x+indice]*var[x+indice]))
   vincoli.append(fob.addConstr(((var[x+2+indice]+e-2*var[x+indice]*e)/matrice6[i][j][0][4][x+1][0])<=imax)

   indice=indice+2

fob.update()
if h==0:
l=len(var)
o=0
u=0
if int(matrice1[i][j])>=3 and h==0:
for x in range (1,l-3,3):

for y in range (x,l-3,3):
vincoli.append(fob.addConstr((var[x]-var[y+3])<=sigmamax))
vincoli.append(fob.addConstr((-var[x]+var[y+3])<=sigmamax))
var.append(fob.addVar(vtype='C',name=str(matrice3[i][j][20][o])))
o=o+1
dim=len(var)
fob.update()

u=0
for y in range (x+1,l-3,3): #vincoli sulle vij^2
vincoli.append(fob.addConstr(var[l+u]==var[x-1]*var[y+1]) )
vincoli.append(fob.addConstr(((var[x+1]+var[y+3]-2*var[l+u])/matrice6[i][j][0][4][x][u+2])<=imax))
u=u+1
print(x)
fob.update()


sin=fob.addVar(vtype='C',name='sin'+str(var[(d+1)])+str(sigma))#sini0
if h!=0 and matrice1[i][j]>=3:
for g in range(0,len(n)):
a=str((var[d+1]))
    b=str(n[g])
    if ((a[13])) in b:
vettore_somma.append(n[g])
q=fob.addVar(vtype='C',name=str(var[(d+1)])+str(sigma))
vincoli.append(fob.addConstr(q==var[d+1]-sigma))
cos=fob.addVar(vtype='C',name='cos'+str(var[(d+1)])+str(sigma))
fob.update()
vincoli.append(fob.addGenConstrSin(q,sin,name='sin'+str(var[d+1])+str(sigma)))
vincoli.append(fob.addGenConstrCos(q,cos,name='cos'+str(var[d+1])+str(sigma)))
fob.update()

u=(var[d+2])*matrice6[i][j][0][8][k][k]
vettore_somma.append(u)
u=(var[d])*e*(matrice6[i][j][0][8][k][0]*cos+matrice6[i][j][0][9][k][0]*sin)
vettore_somma.append(u)

p=0

f=int(matrice1[i][j])-1-k

for y in range (s+l,l+f+s):
if p!=0:
n1=3*p
else:
n1=4*p


sin1=fob.addVar(vtype='C',name='sin'+str(var[(d+1)])+str(var[(d+4+n1)])) q1=fob.addVar(vtype='C',name=str(var[(d+1)])+str(var[(d+4+n1)]))
vincoli.append(fob.addConstr(q1==var[d+1]-var[d+4+n1]))
cos1=fob.addVar(vtype='C',name='cos'+str(var[(d+1)])+str(var[(d+4+n1)]))
fob.update()
                                   vincoli.append(fob.addGenConstrSin(q1,sin1,name='sin'+str(var[d+1])+str(var[d+4+n1])))
vincoli.append(fob.addGenConstrCos(q1,cos1,name='cos'+str(var[d+1])+str(var[d+4+n1])))
fob.update()
u=(var[y])*(matrice6[i][j][0][8][k][p+k+1]*(cos1)+matrice6[i][j][0][9][k][p+k+1]*(sin1))
vettore_somma.append(u)


sin2=fob.addVar(vtype='C',name='sin'+str(var[d+4+n1])+str(var[(d+1)]))
q2=fob.addVar(vtype='C',name=str(var[d+4+n1])+str(var[d+1]))
vincoli.append(fob.addConstr(q2==-var[d+1]+var[d+4+n1]))
cos2=fob.addVar(vtype='C',name='cos'+str(var[(d+4+n1)])+str(var[(d+1)]))
fob.update()
vincoli.append(fob.addGenConstrSin(q2,sin2,name='sin'+str(var[d+4+n1])+str(var[d+1])))
vincoli.append(fob.addGenConstrCos(q2,cos2,name='cos'+str(var[d+4+n1])+str(var[d+1])))
fob.update()
n.append((var[y])*(matrice6[i][j][0][8][k][p+k+1]*cos2+matrice6[i][j][0][9][k][p+k+1]*sin2))
p=p+1

espressione=0
for z in range (0,len(vettore_somma)):
espressione=espressione+vettore_somma[z]
vettore_finale.append(espressione)
print(espressione)
if ass_er==1:
fob.addConstr(espressione<=matrice2[i][j][h+1][17]/Prif)

fob.addConstr(espressione>=(matrice2[i][j][h+1][4]*matrice2[i][j][h+1][9]/(Prif*100)))
d=d+3
k=k+1
s=p
vettore_somma=[]

fob.setObjective(quicksum(vettore_finale[i] for i in range (0,len(vettore_finale))))
fob.update()
fob.optimize()
print(var)
print(vettore_finale)
print (GRB.OPTIMAL)
fob.computeIIS()
fob.write("model.ilp")