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")