How do I correctly solve this MINLP problem and ensure that a feasible solution is found using GEKKO?

Question

I am trying to solve an MINLP problem to generate a heat exchanger (HE) network with heat pumps (HP) using a stage-wise superstructure with 9 stages which has about 7.7k variables. I find it very confusing because changing the initial values of some variables like the heat loads of exchangers/heat pumps by trial and error somehow determines if a solution is found or not after relatively long solution times. Is there a proper approach to tackle this problem? And which possible methods can I use to improve my model? So far, I've tried to minimize the hot utility needed Q_hu which will yield a solution that can be used to minimize TAC. This has worked sometimes, but only by pure coincidence when I somehow plug in specific inital values.

I've heard that solving it using IPOPT first could work, but using IPOPT also didn't work in my case. However, using the 2 hot and 2 cold streams almost always yielded at least a suboptimal solution which changes according to the initial values that I plugged in.

As for the heat pump calculations, I have reformulated the COP formula so that it avoids singularity when the temperature difference is zero. I know that COP is actually dependent on source temperature, but I decided it might be better to make it dependent solely on temperature lift.

from gekko import GEKKO
import numpy as np
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
# Input data with 2 Hot 2 Cold
# Ts_h = np.array([90,100])
# Tt_h = np.array([89,99])
# Ts_c = np.array([120,130])
# Tt_c = np.array([121,131])
# C_h = np.array([500,250])
# C_c = np.array([500,500])
# alt input data with more streams
Ts_h = np.array([100,100,100.0,108.0,91.0,130.0,135.0,150.0,120.0])
Tt_h = np.array([99,99,99,107.0,90.0,30.0,30,65.0,30.0])
Ts_c = np.array([30,120.0,100.0,127,159,155,189,182,172,187,158])
Tt_c = np.array([120.0,121.0,101.0,128,160,156,190,183,173,188,159])
C_h = np.array([650,1300,200.0,250.0,900,5.2,6.6,3.6,2.2])
C_c = np.array([13,380,240,1000,800,800,1000,230,350,230,170])
I = len(Ts_h) #hot stream
J = len(Ts_c) #cold stream

omega = 4000.0
gamma = 1000.0
eps = 1e-5
s = 9
Tmin = eps
ALMTD = 10
Cost_he_fix =29000.0
Cost_he_area =270.0
Cost_hp = 800
Cost_cu = 0.005*8000
Cost_hu = 0.085*8000
Cost_el = 0.1*8000
discount_rate = 0.07
lifetime = 20
CRF = discount_rate*(1+discount_rate)**lifetime/((1+discount_rate)**lifetime-1)
U=0.5

#Variables for temp intervals and constraints
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[i])
    for k in range(1,s+1):
        T_h[i,k].lower = Tt_h[i]
        T_h[i,k].upper = Ts_h[i]
        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])
    

#Binary variables for heatex existence, heat flow, and approach temps
z_he = m.Array(m.Var, (I, J, s), integer=True, lb=0, ub=1)
z_hp = m.Array(m.Var, (I, J, s), integer=True, lb=0, ub=1)
Q_he = m.Array(m.Var, (I, J, s), lb=0, value=1000)
Q_hp = m.Array(m.Var, (I, J, s), lb=0, value=1)
Q_cu = m.Array(m.Var, (I), lb=0, value=eps)
Q_hu = m.Array(m.Var, (J), lb=0, value=eps)
dT1_he = m.Array(m.Var, (I, J, s), lb=Tmin, ub=gamma)
dT2_he = m.Array(m.Var, (I, J, s), lb=Tmin, ub=gamma)
dT1_hp = m.Array(m.Var, (I, J, s), lb=0, ub=gamma)
dT2_hp = m.Array(m.Var, (I, J, s), lb=0, ub=gamma)
#Explicit calculation of Areas
COP = np.empty((I, J, s), dtype='object')
P_hp = np.empty((I, J, s), dtype='object')
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])*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])*Tmin)
            m.Equation(Q_he[i,j,k]*(z_he[i,j,k]-1) == 0) #logical constraint Hall et al
            #LMTD approx by CHEN
            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]/U/LMTD[i,j,k])
            # A_he[i,j,k]=m.Intermediate(Q_he[i,j,k]/U/ALMTD)
            #constraints HP
            m.Equation(dT1_hp[i,j,k]==-z_hp[i,j,k]*(T_h[i,k]-T_c[j,k]))
            m.Equation(dT2_hp[i,j,k]==-z_hp[i,j,k]*(T_h[i,k+1]-T_c[j,k+1]))
            m.Equation(Q_hp[i,j,k]*(z_hp[i,j,k]-1) == 0)
            COP[i,j,k]=m.Intermediate(0.5*(1-293/(303+0.5*(dT1_hp[i,j,k]+dT2_hp[i,j,k])))**-1)
            # COP[i,j,k]=m.Intermediate(3) #take a constant COP value
            P_hp[i,j,k]=m.Intermediate(Q_hp[i,j,k]/(COP[i,j,k]-1))


#no stream split for every stream
for i in range(I):
    for k in range(s):
        m.Equation(sum([z_he[i,j,k]+z_hp[i,j,k] for j in range(J)]) <= 1)
for j in range(J):
    for k in range(s):
        m.Equation(sum([z_he[i,j,k]+z_hp[i,j,k] for i in range(I)]) <= 1)


#energy balance for each stream
for i in range(I):
    m.Equation(sum([Q_he[i,j,k] for j in range(J) for k in range(s)])
               +sum([Q_hp[i,j,k] for j in range(J) for k in range(s)])
               +Q_cu[i] == C_h[i]*(Ts_h[i]-Tt_h[i]))
    
for j in range(J):
    m.Equation(sum([Q_he[i,j,k] for i in range(I) for k in range(s)])
               +sum([Q_hp[i,j,k] for i in range(I) for k in range(s)])
               +sum([P_hp[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 each stage
for k in range(s):
    for i in range(I):
        m.Equation(sum([Q_he[i,j,k] for j in range(J)])
                   +sum([Q_hp[i,j,k] for j in range(J)])== C_h[i]*(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)])
                   +sum([Q_hp[i,j,k] for i in range(I)])
                   +sum([P_hp[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[i]*Cost_hu for i in range(J)]))
sum_Cost_el = m.Intermediate(sum([P_hp[i,j,k] for i in range(I) for j in range(J) for k in range(s)])*Cost_el)
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)
sum_Cost_hp1 = m.Intermediate(sum([Q_hp[i,j,k] for i in range(I) for j in range(J) for k in range(s)])*Cost_hp)
sum_Cost_hp2 = m.Intermediate(sum([P_hp[i,j,k] for i in range(I) for j in range(J) for k in range(s)])*Cost_hp)
TAC = m.Intermediate(sum_Cost_cu+sum_Cost_hu+sum_Cost_el+CRF*(sum_Cost_he_fix+sum_Cost_he_area+sum_Cost_hp1+sum_Cost_hp2))
Obj1 = m.Intermediate(sum_Cost_cu+sum_Cost_hu+1*(sum_Cost_he_fix+sum_Cost_he_area+sum_Cost_hp1+sum_Cost_hp2))
Obj2 = m.Intermediate(sum(Q_hu))   

m.options.SOLVER= 1
m.options.COLDSTART = 0
m.solver_options = ['minlp_maximum_iterations 500', \
            # minlp iterations with integer solution
            'minlp_max_iter_with_int_sol 500', \
                    'minlp_gap_tol 0.1',
                    '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 1']
m.Minimize(Obj2)
m.solve(disp=True)
print('heat recovery / kW:', sum(Q_he[i,j,k][0] for i in range(I) for j in range(J) for k in range(s)))
print('Area / m2:', sum([A_he[i,j,k][0] for i in range(I) for j in range(J) for k in range(s)]))
print(f'Total annual operating costs /t EUR:', (sum_Cost_hu[0]+sum_Cost_cu[0]+sum_Cost_el[0])/1000)
print(f'Total investment costs /t EUR:', (sum_Cost_he_fix[0]+sum_Cost_he_area[0]+sum_Cost_hp1[0]+sum_Cost_hp2[0])/1000)
print(f'Total annualized costs ({lifetime} year(s)) /t EUR:', TAC[0]/1000)
print('no. of HE:', sum([z_he[i,j,k][0] for i in range(I) for j in range(J) for k in range(s)]))
print('HP thermal power / kW:', sum([Q_hp[i,j,k][0]+P_hp[i,j,k][0] for i in range(I) for j in range(J) for k in range(s)]))
print('HP elec power / kW:', sum([P_hp[i,j,k][0] for i in range(I) for j in range(J) for k in range(s)]))
print('cold util / kW:', sum(Q_cu[i][0] for i in range(I)),'\n', 
      'hot util / kW:', sum(Q_hu[j][0] for j in range(J)))