Why does the problem become infeasibly if the bounds of the decision variable are increased?

Question

I am programming an heat exchanger optimization problem as a MINLP problem using Pyomo 5.6.6 and Couenne 0.5.6 as downloaded executable file for MacOS. The model has two decision variables, which are the external supply of heat in form of steam and cooling utility in form of cooling water.

The problem that is occurring is, that the model is only solved if I set boundaries for the use of heat/cooling in a rather small window around the optimal value. Otherwise Couenne says the problem is infeasible.

Does somebody know why the solver has troubles solving with bigger bounds?

The heat exchanger optimization model is based on an expanded-transhipment model.

Right now the heat exchanger network is made up of two hot and two cold streams that are defined by inlet and outlet temperatures as well as their heat capacity flow rates. Additionally there are one free hot stream (hot steam) and one cold stream (cooling water) that should be calculated by the optimization in order to fulfill the heat balance while minimizing the use of the additional hot steam.

I set bounds for the use of heat and cooling utilities as constraints in the form Constraint1: Q <= value1 and Constraint2: Q >= value2.

The problem is that I have set these values quite near to the optimal value, which is impractical in case I don't have a good idea about that value.

The code of the model as well as the example data and solver outputs for working/not-working case can be found on Github here.

I know that the minimal value for additional heat supply has to be 950 in my model. However, this value is only calculated if I set the bounds between 1 and 1100. If I change the bounds to e.g. 0 and 5000 the solver returns infeasible.

In the best case I would like to set the use of additional heating and cooling to Non-Negative and let the rest be calculated.

The code for the model is:

import sys
sys.path.append('Here comes the path of the data.dat file')
from pyomo.environ import *

model = AbstractModel(name="(HEN_MODEL_V1)")

# hot streams
model.I = Set()
# cold streams
model.J = Set()

# k = grid point and kk =temperature interval boundaries
model.K = Set()

# hot/ cold utilities
model.I_HU = Set(within=model.I)
model.J_CU = Set(within=model.J)

# hot process streams
model.I_P = Set(within=model.I)
model.J_P = Set(within=model.J)

# temperature intervals
model.K_I = Set(within=model.K)



#  ----- Parameters ----

# Minimal temperature approach
model.MAT = Param()

# Upper bound for inlet temperature
model.ai = Param(model.I)
model.aj = Param(model.J)

# Upperbound for outlet temperature
model.bi = Param(model.I)
model.bj = Param(model.J)

#Upperbound for heat steam
model.yi = Param(model.I, model.K_I)
model.yj = Param(model.J, model.K_I)

#Upperbounds for R and D
model.yyi = Param(model.I, model.K_I)
model.yyj = Param(model.J, model.K_I)

# Parameter for in and outlet temperature constraints
model.T_H_IN_CON = Param(model.I)
model.T_H_OUT_CON = Param(model.I)

model.T_C_IN_CON = Param(model.J)
model.T_C_OUT_CON = Param(model.J)

model.F_CON = Param(model.I)
model.F2_CON = Param(model.J)


# ----- Variables-----

# Gridtemperature
model.T = Var(model.K, within=NonNegativeReals)

# Residual heat of heat steam i in intervall k
model.R = Var(model.I, model.K_I, within=NonNegativeReals)

# Deficite heat of cold stream j in intervall k
model.D = Var(model.J, model.K_I, within=NonNegativeReals)

# Exchanged Heat from i to j in intervall k
model.Q = Var(model.I, model.J, model.K_I, within=NonNegativeReals)

# Heat in Cold stream j 
model.QC = Var(model.J, model.K_I, within=NonNegativeReals)

#Disaggregates Cold streams j
model.QC1 = Var(model.J, model.K_I, within=NonNegativeReals)
model.QC2 = Var(model.J, model.K_I, within=NonNegativeReals)

# Heat in hot stream i
model.QH = Var(model.I, model.K_I, within=NonNegativeReals)

# Disaggregates hot streams i
model.QH1 = Var(model.I, model.K_I, within=NonNegativeReals)
model.QH2 = Var(model.I, model.K_I, within=NonNegativeReals)

# Heat Capacity flows of i and j
model.FH = Var(model.I_P, within=NonNegativeReals)
model.FC = Var(model.J, within=NonNegativeReals)

# Used additional Steam

model.steam = Var()


# Binary variables X for Inlet-Temp on Gridpoint
# Binary variables Y for Outlet-Temp in Grid k
# Binary variales Z for Stream through intervall k
model.XC = Var(model.J, model.K, within=Binary)
model.XH = Var(model.I, model.K, within=Binary)

model.YC = Var(model.J, model.K_I, within=Binary)
model.YH = Var(model.I, model.K_I, within=Binary)

model.ZC = Var(model.J, model.K_I, within=Binary)
model.ZH = Var(model.I, model.K_I, within=Binary)

# Disaggregated Inlettemperatures from Cold and Hot stream
model.TIN_HD = Var(model.I, model.K, within=NonNegativeReals)
model.TIN_CD = Var(model.J, model.K, within=NonNegativeReals)

# Inlettemperatures from Cold and hot stream
model.TIN_H = Var(model.I, within=NonNegativeReals)
model.TIN_C = Var(model.J, within=NonNegativeReals)

# Outlettemperature from Hot and cold streams 
model.TOUT_HD = Var(model.I, model.K_I, within=NonNegativeReals)
model.TOUT_CD = Var(model.J, model.K_I, within=NonNegativeReals)

# Outlettemperatures from Cold and hot stream
model.TOUT_H = Var(model.I, within=NonNegativeReals)
model.TOUT_C = Var(model.J, within=NonNegativeReals)




# ------ Constraints ------


# DEFINITION OF THE INLET AND OUTLET TEMPERATURES OF THE HOT AND COLD STREAMS
# THAT HAVE TO BE MET


def T1_rule(model, i):
    return model.TIN_H[i] == model.T_H_IN_CON[i]

def T2_rule(model,i):
    return model.TOUT_H[i] == model.T_H_OUT_CON[i]

def T3_rule(model,j):
    return model.TIN_C[j] == model.T_C_IN_CON[j]

def T4_rule(model, j):
    return model.TOUT_C[j] == model.T_C_OUT_CON[j]


model.T1_con = Constraint(model.I, rule=T1_rule)
model.T2_con = Constraint(model.I, rule=T2_rule)
model.T3_con = Constraint(model.J, rule=T3_rule)
model.T4_con = Constraint(model.J, rule=T4_rule)



# DEFINITION OF THE FLOW CAPACITIES OF THE HEAT AND COLD PROCESS STREAMS
# DEFINTION OF THE UPPER AND LOWER BOUNDS OF THE  FLOW CAPACITY OF THE 
# COOLING ULTILITY 




def F1_rule(model, i):
    return model.FH[i] == model.F_CON[i]

def F2_rule_a(model,j):
    return model.FC[j] == model.F2_CON[j]

def F2_rule_b(model,j):
    return model.FC[j] >= 0

def F2_rule_c(model,j):
    return model.FC[j] <= 5000




model.F1_con = Constraint(model.I_P, rule=F1_rule)
model.F2_con_a = Constraint(model.J_P, rule=F2_rule_a)

model.F2_con_b = Constraint(model.J_CU, rule=F2_rule_b)
model.F2_con_c = Constraint(model.J_CU, rule=F2_rule_c)

############ GLG 1 + 2 ################
# DEFINITION OF HEAT CASCADE 


def Rcas_rule(model,i,k):
    if k > 1:
        return model.R[i,k] + sum(model.Q[i,j,k] for j in model.J) == model.R[i,k-1] + model.QH[i,k] 
    else:
        return model.R[i,k] + sum(model.Q[i,j,k] for j in model.J) == model.QH[i,k] 


def Dcas_rule(model,j,k):
    a=len(model.K_I)
    if k < a:
        return model.D[j,k] + sum(model.Q[i,j,k] for i in model.I) == model.D[j,k+1] + model.QC[j,k] 
    else:
        return model.D[j,k] + sum(model.Q[i,j,k] for i in model.I) == model.QC[j,k] 

model.Rcas_con = Constraint(model.I, model.K_I, rule=Rcas_rule)
model.Dcas_con = Constraint(model.J, model.K_I, rule=Dcas_rule)




######### GLG 3 ############
# DEFINITION OF GRIDPOINT TEMPERATURES

def T_k_rule(model,k):
    a = len(model.K)
    if k < a-1:
        return model.T[k] >= model.T[k+1]
    else:
        return model.T[k] <= model.T[k-1]

model.T_k_con = Constraint(model.K, rule=T_k_rule)




######### GLG 4 - 7 ###########
# DISAGGREGATED HOT INLET TEMPERATURES
# ASSIGNMENT TO GRIDPOINT TEMPERATURES

def TIN_Hot(model,i):
    return  sum(model.TIN_HD[i,k] for k in model.K) == model.TIN_H[i]

def THD_1(model,i,k):
    return model.TIN_HD[i,k] >= model.T[k] - model.ai[i] * (1-model.XH[i,k])

def THD_2(model,i,k):
    return model.TIN_HD[i,k] <= model.T[k] 

def THD_3(model,i,k):
    return model.TIN_HD[i,k] <= model.ai[i] * model.XH[i,k]

model.THIN1 = Constraint(model.I, rule=TIN_Hot)
model.THIN2 = Constraint(model.I, model.K, rule=THD_1)
model.THIN3 = Constraint(model.I, model.K, rule=THD_2)
model.THIN4 = Constraint(model.I, model.K, rule=THD_3)
#  




########## GLG 8 - 11 #############
# DISSAGREGATED COLD INLET TEMPERATURES
# ASSIGMENT TO GRIDPOINT TEMPERATUES

def TIN_Cold(model,j):
    return sum(model.TIN_CD[j,k] for k in model.K) == model.TIN_C[j]
#    return model.TIN_C[j] == sum(model.TIN_CD[j,k] for k in model.K)

def TCD_1(model,j,k):
    return model.TIN_CD[j,k] >= model.T[k] - model.MAT - model.aj[j] * (1-model.XC[j,k])

def TCD_2(model,j,k):
    return model.TIN_CD[j,k] <= model.T[k] - model.MAT

def TCD_3(model,j,k):
    return model.TIN_CD[j,k] <= model.aj[j] * model.XC[j,k]

model.TCIN1 = Constraint(model.J, rule=TIN_Cold)
model.TCIN2 = Constraint(model.J, model.K, rule=TCD_1)
model.TCIN3 = Constraint(model.J, model.K, rule=TCD_2)
model.TCIN4 = Constraint(model.J, model.K, rule=TCD_3)
#
######### GLG 12 - 15 ##############
# DISSAGREGATED HOT OUTLET TEMPERATURES
# ASSIGMENT TO TEMPERATURE INTERVALS 





def TOUT_Hot(model,i):
    return sum(model.TOUT_HD[i,k] for k in model.K_I) == model.TOUT_H[i]

def THD_11(model,i,k):
    return model.TOUT_HD[i,k] >= model.T[k] - model.bi[i] * (1-model.YH[i,k])

def THD_22(model,i,k):
    return model.TOUT_HD[i,k] <= model.T[k-1]                                          # -1 Funktuniert sicher nicht

def THD_33(model,i,k):
    return model.TOUT_HD[i,k] <= model.bi[i] * model.YH[i,k]

model.THOUT1 = Constraint(model.I, rule=TOUT_Hot)
model.THOUT2 = Constraint(model.I, model.K_I, rule=THD_11)
model.THOUT3 = Constraint(model.I, model.K_I, rule=THD_22)
model.THOUT4 = Constraint(model.I, model.K_I, rule=THD_33)
#


#
######### GLG 16 - 19 ###########
# DISAGGREGATED COLD OUTLET TEMPERATURES
# ASSIGNMENT TO TEMPERATURE INVERTALS

def TOUT_Cold(model,j):
    return sum(model.TOUT_CD[j,k] for k in model.K_I) == model.TOUT_C[j] 

def TCD_11(model,j,k):
    return model.TOUT_CD[j,k] >= model.T[k] - model.MAT - model.bj[j] * (1-model.YC[j,k])

def TCD_22(model,j,k):
    return model.TOUT_CD[j,k] <= model.T[k-1]                             # Auch hier -1 in der Klammer

def TCD_33(model,j,k):
    return model.TOUT_CD[j,k] <= model.bj[j] * model.YC[j,k]


model.TCOUT1 = Constraint(model.J, rule=TOUT_Cold)
model.TCOUT2 = Constraint(model.J, model.K_I, rule=TCD_11)
model.TCOUT3 = Constraint(model.J, model.K_I, rule=TCD_22)
model.TCOUT4 = Constraint(model.J, model.K_I, rule=TCD_33)




############ GLG 20 - 25 ###############
# BINARY CONSTRAINTS FÜR TEMPERATURE GRID

# Constraints for Binary Variables of Temperature Grid
def XH_rule(model,i):
    return sum(model.XH[i,k] for k in model.K) == 1

def XC_rule(model,j):
    return sum(model.XC[j,k] for k in model.K) == 1

def YH_rule(model,i):
    return sum(model.YH[i,k] for k in model.K_I) == 1

def YC_rule(model,j):
    return sum(model.YC[j,k] for k in model.K_I) == 1

def ZH_rule(model,i,k):
    if k==1:
        return model.ZH[i,k] == model.XH[i,k-1] - model.YH[i,k]    # Auch hier die -1
    else:
        return model.ZH[i,k] == model.ZH[i,k-1] + model.XH[i,k-1] - model.YH[i,k]    # Auch hier die -1


def ZC_rule(model,j,k):
    temp = len(model.K_I)
    if k < temp:
        return model.ZC[j,k] == model.ZC[j,k+1] + model.XC[j,k] - model.YC[j,k]
    else:
        return model.ZC[j,k] == model.XC[j,k] - model.YC[j,k]



model.XH_con = Constraint(model.I, rule=XH_rule)
model.XC_con  = Constraint(model.J, rule=XC_rule)
model.YH_con = Constraint(model.I, rule=YH_rule)
model.YC_con = Constraint(model.J, rule=YC_rule)
model.ZH_con = Constraint(model.I, model.K_I, rule=ZH_rule)
model.ZC_con = Constraint(model.J, model.K_I, rule=ZC_rule)




############### GLG 26 - 32 ################
# DISAGGRATED HEAT LOADS IN Q1 IF HEAT GOES THROUGH k AND Q2 if PARTIAL SPANS k



# Constraints for disaggregated Heat streams 
def QH_rule(model,i,k):
    return model.QH[i,k] == model.QH1[i,k] + model.QH2[i,k]

###### QH1 ######
def QH2_rule(model,i,k):
    return model.QH1[i,k] >= model.FH[i] * (model.T[k-1] - model.T[k]) - model.yi[i,k] * (1-model.ZH[i,k])            # Hier wieder -1  

def QH3_rule(model,i,k):
    return model.QH1[i,k] <= model.FH[i] * (model.T[k-1] - model.T[k])          # hier wieder die -1

def QH4_rule(model,i,k):
    return model.QH1[i,k] <= model.yi[i,k] * model.ZH[i,k] 
###### QH2 ######
def QH5_rule(model,i,k):
    return model.QH2[i,k] >= model.FH[i] * (model.T[k-1] - model.TOUT_HD[i,k]) - model.yi[i,k] * (1-model.YH[i,k])

def QH6_rule(model,i,k):
    return model.QH2[i,k] <= model.FH[i] * (model.T[k-1] - model.TOUT_HD[i,k]) + model.yi[i,k] * (1-model.YH[i,k])

def QH7_rule(model,i,k):
    return model.QH2[i,k] <= model.yi[i,k] * model.YH[i,k]



model.QH_con = Constraint(model.I_P, model.K_I, rule=QH_rule)
model.QH2_con = Constraint(model.I_P, model.K_I, rule=QH2_rule)
model.QH3_con = Constraint(model.I_P, model.K_I, rule=QH3_rule)
model.QH4_con = Constraint(model.I_P, model.K_I, rule=QH4_rule)

model.QH5_con = Constraint(model.I_P, model.K_I, rule=QH5_rule)
model.QH6_con = Constraint(model.I_P, model.K_I, rule=QH6_rule)
model.QH7_con = Constraint(model.I_P, model.K_I, rule=QH7_rule)

############### GLG 33 ###################
# HEAT LOAD FOR HEAT UTILITY

def QH8_rule(model,i,k):
    return model.QH[i,k] <= model.yi[i,k] * model.XH[i,k-1]    

model.QH8_con = Constraint(model.I_HU, model.K_I, rule=QH8_rule)






############## GLG 34 - 40 #####################
# DISAGGREGATED COLD HEAT LOADS; Q1 IF COLD GOES THROUGH k , Q2 IF IT ENDS THERE
# Constraints for disaggregated Cold streams 
def QC_rule(model,j,k):
    return model.QC[j,k] == model.QC1[j,k] + model.QC2[j,k]

####### QC1 #########
def QC2_rule(model,j,k):
    return model.QC1[j,k] >= model.FC[j] * (model.T[k-1] - model.T[k]) - model.yj[j,k] * (1 - model.ZC[j,k])            # Hier wieder -1  

def QC3_rule(model,j,k):
    return model.QC1[j,k] <= model.FC[j] * (model.T[k-1] - model.T[k])          # hier wieder die -1

def QC4_rule(model,j,k):
    return model.QC1[j,k] <= model.yj[j,k] * model.ZC[j,k] 

####### QC2 ###########
def QC5_rule(model,j,k):
    return model.QC2[j,k] >= model.FC[j] * (model.TOUT_CD[j,k] - model.T[k] + model.MAT) - model.yj[j,k] * (1 - model.YC[j,k])

def QC6_rule(model,j,k):
    return model.QC2[j,k] <= model.FC[j] * (model.TOUT_CD[j,k] - model.T[k] + model.MAT) + model.yj[j,k] * (1 - model.YC[j,k])


def QC7_rule(model,j,k):
    return model.QC2[j,k] <= model.yj[j,k] * model.YC[j,k]

model.QC_con = Constraint(model.J, model.K_I, rule=QC_rule)
model.QC2_con = Constraint(model.J, model.K_I, rule=QC2_rule)
model.QC3_con = Constraint(model.J, model.K_I, rule=QC3_rule)
model.QC4_con = Constraint(model.J, model.K_I, rule=QC4_rule)
###
model.QC5_con = Constraint(model.J, model.K_I, rule=QC5_rule)
model.QC6_con = Constraint(model.J, model.K_I, rule=QC6_rule)
model.QC7_con = Constraint(model.J, model.K_I, rule=QC7_rule)





###### GLG 41 + 42 #########
# RESIDUAL HEAT AND DEFIZIT SHOULD BE BELOW A GIVEN UPPER BOUND
#####

def R_rule(model, i, k):
    temp=len(model.K_I)
    if k < temp:
        return model.R[i,k] <= model.yyi[i,k] * sum(model.XH[i,kk] for kk in model.K if kk<= k-1)   
    else:
        return model.R[i,k] == 0

def D_rule(model, j, k):
    if k > 1:
        return model.D[j,k] <= model.yyj[j,k] * (1 - sum(model.XC[j,kk] for kk in model.K if kk <= k-1))  
    else:
        return model.D[j,k] == 0

model.R1_con = Constraint(model.I, model.K_I, rule=R_rule)
model.D1_con = Constraint(model.J, model.K_I, rule=D_rule) 

#####
# DEFINITION OF LOWER AND UPPER BOUND OF STEAM- HEAT THAT CAN BE USED
#####

def Steam_lb(model,):
    return model.QH[3,1] >= 1

def Steam_ub(model,):
    return model.QH[3,1] <= 1500

model.St_lb=Constraint(rule=Steam_lb)
model.St_ub=Constraint(rule=Steam_ub)


##### 
# DEFINITION OF ADDITIONAL VARIABLE TO DISPLAY THE OBJECTIVE VALUE
####

def Steam_use(model):
    return model.steam == sum(model.QH[3,k] for k in model.K_I)

model.Steam_use_con = Constraint(rule=Steam_use)


######
# DEFINITION OF OBJECTIVE RULE: MINIMIZE THE USAGE OF ADDITIONAL STEAM
####


def obj_rule(model):
    return sum(model.QH[3,k] for k in model.K_I) 


model.objective= Objective(rule=obj_rule)


#
#-------------------------Solving---------------------------------------------
#--------------------------------------------------------------------

instance= model.create_instance('data.dat')
solver = SolverFactory('couenne')
results = solver.solve(instance, tee=True)
instance.steam.pprint()

And the text from the data file which is saved as "data.dat" in the path which is specified in the code above:

set K := 0 1 2 3 4 5;
set K_I := 1 2 3 4 5;
set I := 1 2 3;
set J := 1 2 3;
set I_P := 1 2;
set J_P := 1 2;
set I_HU := 3;
set J_CU := 3;


param MAT := 10;

param T_H_IN_CON :=
    1 600
    2 590
    3 680
;

param T_H_OUT_CON :=
    1 370
    2 400
    3 680
;

param T_C_OUT_CON :=
    1 650
    2 490
    3 320
;

param T_C_IN_CON :=
    1 400
    2 350
    3 300
;


param ai :=
    1 1000
    2 1000
    3 1000
;

param aj :=
    1 1000
    2 1000
    3 1000
;

param bi :=
    1 1000
    2 1000
    3 1000
;

param bj :=
    1 1000
    2 1000
    3 1000
;

param yi :=
    1 1 100000
    1 2 100000
    1 3 100000
    1 4 100000
    1 5 100000
    2 1 100000
    2 2 100000
    2 3 100000
    2 4 100000
    2 5 100000
    3 1 100000
    3 2 100000
    3 3 100000
    3 4 100000
    3 5 100000
;

param yj :=
    1 1 100000
    1 2 100000
    1 3 100000
    1 4 100000
    1 5 100000
    2 1 100000
    2 2 100000
    2 3 100000
    2 4 100000
    2 5 100000
    3 1 100000
    3 2 100000
    3 3 100000
    3 4 100000
    3 5 100000
;

param yyi :=
    1 1 100000
    1 2 100000
    1 3 100000
    1 4 100000
    1 5 100000
    2 1 100000
    2 2 100000
    2 3 100000
    2 4 100000
    2 5 100000
    3 1 100000
    3 2 100000
    3 3 100000
    3 4 100000
    3 5 100000
;

param yyj :=
    1 1 100000
    1 2 100000
    1 3 100000
    1 4 100000
    1 5 100000
    2 1 100000
    2 2 100000
    2 3 100000
    2 4 100000
    2 5 100000
    3 1 100000
    3 2 100000
    3 3 100000
    3 4 100000
    3 5 100000
;

param F_CON :=
    1 10
    2 20
;

param F2_CON :=
    1 15
    2 14
;

Answer 1

Couenne is a respectable solver, but as with all solvers, it has weaknesses. I did not look closely at your formulation, but it is possible that near variable values of zero, numerical issues arise, and therefore Couenne terminates unsuccessfully.