[Python] How to use the cold-start options with GEKKO correctly?

I was testing some options with Gekko for an optimization with MINLP , but I still couldn't figure out how the option COLDSTART = 2 work. Using the variable values from a prior solution as initial values to a new solve still somehow shows that infeasibilities exist:

Could not perform presolve. Failed at block 0 with infeasible variables and equations reported in presolve.txt

Do i need to manually put extra equations so that the number of variables is equal to the number of equations? And why does sometimes using the prior solution as initial values to a new solve lead to worse solution with COLDSTART = 1? Shouldn't the first integer solution found be at least as good as the initial values (or solution from the previous solve)

Update: here's how the problem looks like. Basically I just wanted to find stream matches for heat recovery

from gekko import GEKKO
import numpy as np
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
Ts_h = np.array([100,105,110,115,120,125,130,135,140])
Tt_h = np.array([80,85,90,95,100,105,110,115,120])
htc_h = np.array([0.5,0.5,0.5,1.5,1.5,1.5,1.5,1.5,0.5])
C_h = np.array([10,10,10,10,10,10,10,10,10])
Ts_c = np.array([100,110,120,130])
Tt_c = np.array([110,120,130,140])
htc_c = np.array([1.5,0.5,1.5,1.5])
C_c = np.array([40,40,40,40])
I = len(C_h)
J = len(C_c)
eps = 1e-5
s = 2
ALMTD = 10
CAPEX_F = 4
Cost_he_fix =20000*CAPEX_F
Cost_he_area =200*CAPEX_F
Cost_cu = 50
Cost_hu = np.array([500,1000,1000,1000])
discount_rate = 0.1
lifetime = 20
CRF = discount_rate*(1+discount_rate)**lifetime/((1+discount_rate)**lifetime-1)
U=0.5
Tmin = 5

#bounds
omega = np.ones((I, J))
gamma = np.ones((I, J))
for i in range(I):
for j in range(J):
omega[i,j]=np.min([C_h*(Ts_h-Tt_h),C_c[j]*(Tt_c[j]-Ts_c[j])])
gamma[i,j]=np.max([0,Tt_c[j]-Tt_h,Tt_c[j]-Ts_h,Ts_c[j]-Tt_h,Ts_c[j]-Ts_h])

#Variables for temperatures of stream
T_h = m.Array(m.Var, (I, s+1))
T_c = m.Array(m.Var, (J, s+1))

for i in range(I):
m.Equation(T_h[i,0] == Ts_h)
for k in range(1,s+1):
T_h[i,k].lower = Tt_h
T_h[i,k].upper = Ts_h
m.Equation(T_h[i,k] <= T_h[i,k-1])
for j in range(J):
m.Equation(T_c[j,s] == Ts_c[j])
for k in range(0,s):
T_c[j,k].lower = Ts_c[j]
T_c[j,k].upper = Tt_c[j]
m.Equation(T_c[j,k+1] <= T_c[j,k])


#Variables for HEx existence, heat load, and approach temps
z_he = m.Array(m.Var, (I, J, s), integer=True, lb=0, ub=1)
Q_he = m.Array(m.Var, (I, J, s), lb=0)
Q_cu = m.Array(m.Var, (I), lb=0, value=0)
Q_hu = m.Array(m.Var, (J), lb=0, value=0)
dT1_he = m.Array(m.Var, (I, J, s), lb=Tmin)
dT2_he = m.Array(m.Var, (I, J, s), lb=Tmin)
#Explicit calculation of Areas
LMTD = np.empty((I, J, s), dtype='object')
A_he = np.empty((I, J, s), dtype='object')
for i in range(I):
for j in range(J):
for k in range(s):
m.Equation(dT1_he[i,j,k]==z_he[i,j,k]*(T_h[i,k]-T_c[j,k])+(1-z_he[i,j,k])*(gamma[i,j]*0+Tmin))
m.Equation(dT2_he[i,j,k]==z_he[i,j,k]*(T_h[i,k+1]-T_c[j,k+1])+(1-z_he[i,j,k])*(gamma[i,j]*0+Tmin))
m.Equation(Q_he[i,j,k] < z_he[i,j,k]*omega[i,j]) LMTD[i,j,k]=m.Intermediate((0.5*dT1_he[i,j,k]*dT2_he[i,j,k]*(dT1_he[i,j,k]+dT2_he[i,j,k]))**(1/3))
A_he[i,j,k]=m.Intermediate(Q_he[i,j,k]*(1/htc_h+1/htc_c[j])/LMTD[i,j,k])

#energy balance for stream
for i in range(I):
m.Equation(Q_cu == C_h*(T_h[i,s]-Tt_h))
m.Equation(sum([Q_he[i,j,k] for j in range(J) for k in range(s)])+Q_cu == C_h*(Ts_h-Tt_h))

for j in range(J):
m.Equation(Q_hu[j] == C_c[j]*(Tt_c[j]-T_c[j,0]))
m.Equation(sum([Q_he[i,j,k] for i in range(I) for k in range(s)])+Q_hu[j] == C_c[j]*(Tt_c[j]-Ts_c[j]))


#Energy balance for stage
for k in range(s):
for i in range(I):
m.Equation(sum([Q_he[i,j,k] for j in range(J)]) == C_h*(T_h[i,k]-T_h[i,k+1]))
for k in range(s):
for j in range(J):
m.Equation(sum([Q_he[i,j,k] for i in range(I)]) == C_c[j]*(T_c[j,k]-T_c[j,k+1]))


sum_Cost_cu = m.Intermediate(sum(Q_cu)*Cost_cu)
sum_Cost_hu = m.Intermediate(sum([Q_hu[j]*Cost_hu[j] for j in range(J)]))
sum_Cost_he_area = m.Intermediate(sum([A_he[i,j,k] for i in range(I) for j in range(J) for k in range(s)])*Cost_he_area)
sum_Cost_he_fix = m.Intermediate(sum([z_he[i,j,k] for i in range(I) for j in range(J) for k in range(s)])*Cost_he_fix)
TAC = m.Intermediate(sum_Cost_cu+sum_Cost_hu+CRF*(sum_Cost_he_fix+sum_Cost_he_area))

m.solver_options = ['minlp_maximum_iterations 10000', \
'minlp_max_iter_with_int_sol 2000', \
'minlp_gap_tol 1e-2',
'objective_convergence_tolerance 1.0e-6',
#1=depth first (find integer solution faster), 2=breadth first,
#3=lowest objective leaf, 4=highest objective leaf
'minlp_branch_method 2']
m.options.SOLVER= 1
m.options.COLDSTART = 1
m.Minimize(TAC)
m.solve(disp=True)
print(TAC[0])


For example, I can then save the solution and use it again for the next solve as initial values:

#Fill var arrays for intial values
z_he_init = np.empty((I,J,s))
Q = np.empty((I,J,s))
Q_cu_init = np.empty((I))
Q_hu_init = np.empty((J))
T_h_init = np.empty((I,s+1))
T_c_init = np.empty((J,s+1))
dT1_he_init = np.empty((I,J,s))
dT2_he_init = np.empty((I,J,s))
for i in range(I):
for j in range(J):
for k in range(s):
z_he_init[i,j,k] = z_he[i,j,k][0]
Q[i,j,k] = Q_he[i,j,k][0]
Q_cu_init = Q_cu[i][0]
Q_hu_init[j] = Q_hu[j][0]
T_h_init[i,k] = T_h[i,k][0]
T_c_init[j,k] = T_c[j,k][0]
dT1_he_init[i,j,k] = dT1_he[i,j,k][0]
dT2_he_init[i,j,k] = dT2_he[i,j,k][0]
T_h_init[i,s] = T_h[i,s][0]
T_c_init[j,s] = T_c[j,s][0]

# Initialization all var
for i in range(len(Ts_h)):
for k in range(1,s+1):
T_h[i,k].value = T_h_init[i,k] #initialise
for j in range(len(Ts_c)):
for k in range(0,s):
T_c[j,k].value = T_c_init[j,k] #initialise
for j in range(len(Ts_c)):
Q_hu[j].value = Q_hu_init[j]
for i in range(len(Ts_h)):
Q_cu[i].value = Q_cu_init[i]
for j in range(len(Ts_c)):
for k in range(s):
Q_he[i,j,k].value = Q[i,j,k]
z_he[i,j,k].value = z_he_init[i,j,k]
dT1_he[i,j,k].value = dT1_he_init[i,j,k]
dT2_he[i,j,k].value = dT2_he_init[i,j,k]

[url="https://stackoverflow.com/questions/79523417/how-to-use-the-cold-start-options-with-gekko-correctly"]Continue reading...[/url][/i][/i][/i]