Gekko large MINLP semi-continuous variable issue

Try reducing the problem by simplifying constraints such as

soe_low = 0.1
soe_high = 0.9
new_eq(soe[t] >= soe_low)
new_eq(soe[t] <= soe_high)

to be replaced by:

soe = new_array(lb=soe_low,ub=soe_high)

This reduces the number of degrees of freedom. Equations can also be rearranged to avoid divide-by-zero if there are variables in the denominator.

#new_eq(soe[t] == e_stock[t] / Ebat_rated / soh)
new_eq(Ebat_rated * soh * soe[t] == e_stock[t])

Another thing to try is an if2() statement that uses an MPCC instead of a binary variable. Mixed integer optimization problems can be challenging to solve. It appears that there are additional constraints after IPOPT solves the NLP (as you mentioned):

for t in range(n_step):
    if t <= 20 or t >= 76: # Reduction of the binary problem
        new_eq(y_gen[t] == 1)

Your approach to change the options for the MINLP solver is good. It can take a long time to solve, especially with many binary variables. You may also want to increase the max iteration limit. APOPT (NLP) required 753 iterations (barely under the limit of 1000). Most optimization problems are typically solved in 5-50 iterations so this problem is particularly challenging. Here is something that may be closer to solving successfully.

import functools

from gekko import GEKKO
import pandas as pd
import matplotlib.pyplot as plt

QUARTER_HOUR = pd.Timedelta(hours=0.25)
MinutesPerHour = 60


##########################################################
#####################  MODEL DEFINITION  #################
##########################################################

gekko_model = GEKKO(remote=True)

##########################################################
########################  PARAMETER ######################
##########################################################

# Design PV + inverter + transfo
Pmax_inverter_pv = 600
eff_ondul_pv = 0.98
eff_transfo_pv = 0.99

# Design battery + inverter + transfo
n_bat = 3
p_bat_min = -750
p_bat_max = 750
Smax_inverter_bat = p_bat_max
Ebat_rated = 980
eff_ondul_bat = 0.98
eff_transfo_bat = 0.99
soe_low = 0.1
soe_high = 0.9
eff_bat = 0.91 ** 0.5
tanphi_min = 5


# Design Diesel Generator
n_gen = 2
Sgen_max = 1000
Pgen_min = 150
Pgen_max = 800
eff_transfo_gen = 0.99
conso_idle = 1
conso_rated = 10

# Auxiliaries
eff_transfo_aux = 0.99
hours_aux = 10
p_aux = 10
p_aux_hv = p_aux / eff_transfo_aux

# Max production power
p_mpp_dc = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.40, 6.67, 18.66, 34.52, 53.56, 73.02, 91.67, 96.52, 178,
            203.54, 201.52, 294.74, 310.06, 340.17, 370.87, 423.51,
            447.55, 410.88, 437.4, 463.51, 493.55, 499.92, 522.68, 509.18,
            502.81, 522.82, 541.00, 515.43, 514.88, 465.46, 370.32, 523.98,
            324.65, 376.04, 401.06, 276.15, 192.13, 166.11, 152.92, 211.46,
            99.90, 166.58, 139.01, 71.05, 80.44, 46.71, 27.31, 38.20, 32.59,
            16.52, 5.11, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

p_mpp_ac_hv = [p_dc * eff_ondul_pv * eff_transfo_pv for p_dc in p_mpp_dc]

# Timestep
n_step = 96
dt = 0.25

# Load
p_load = [300] * n_step
q_load = [150] * n_step

# initial conditions
soh = 1
soe_init = 0.2
soh_gen = 1
e_stock_init = soe_init * soh * Ebat_rated
p_gen_lv_init = 300

# Energy reserve for auxiliaries
soe_reserve_aux = soe_low + hours_aux * p_aux / (eff_bat * eff_ondul_bat
                                                 * eff_transfo_bat
                                                 * eff_transfo_aux
                                                 * Ebat_rated * soh)

# Optimization constants
ZERO = 0
UB = 10_000

##########################################################
#####  DECLARATION & INITIALIZATION of VARIABLES  ########
##########################################################

# Shortcut for creation of variables and arrays
new_var = functools.partial(gekko_model.Var,
                            value=ZERO, lb=ZERO, ub=UB, integer=False)
new_array = functools.partial(gekko_model.Array,
                              new_var, n_step)

# Definition of variables
smoothing_penalty = new_array()
p_pv_dc = new_array()
p_pv_lv = new_array()
p_pv_curtail = new_array()
p_pv_hv = new_array()
e_stock = new_array()
p_bat_dc_charge = new_array()
p_bat_dc_discharge = new_array()
p_bat_dc = new_array(lb=-UB)
soe = new_array(lb=soe_low,ub=soe_high)
s_bat_lv = new_array()
p_bat_lv = new_array(lb=-UB)
p_bat_lv_charge = new_array()
p_bat_lv_discharge = new_array()
q_bat_lv = new_array(lb=-UB)
q_bat_lv_charge = new_array()
q_bat_lv_discharge = new_array()
p_bat_hv = new_array(lb=-UB)
p_bat_hv_charge = new_array()
p_bat_hv_discharge = new_array()
q_bat_hv = new_array(lb=-UB)
q_bat_hv_charge = new_array()
q_bat_hv_discharge = new_array()
y_gen = new_array(ub=1, integer=True)
s_gen_lv = new_array()
p_gen_lv = new_array()
q_gen_lv = new_array()
p_gen_hv = new_array()
q_gen_hv = new_array()
p_enr = new_array()
conso_gen = new_array()
conso_gen_tot = new_var()
penalty_tot = new_var(lb=-UB)
curtail_penalty = new_var()
obj = new_var(lb=-UB)

##########################################################
###############  OPTIMIZATION CONSTRAINTS  ###############
##########################################################

# Creation of new equation
new_eq = gekko_model.Equation

# ---------------------------
# PV field constraints ------
# ---------------------------

for t in range(n_step):
    new_eq(p_mpp_dc[t] == p_pv_dc[t] + p_pv_curtail[t])
    new_eq(p_pv_dc[t] * eff_ondul_pv == p_pv_lv[t])
    new_eq(p_pv_lv[t] <= Pmax_inverter_pv)
    new_eq(p_pv_lv[t] * eff_transfo_pv == p_pv_hv[t])

# -----------------------------
# Battery constraints ---------
# -----------------------------

for t in range(n_step):

    if t == 0:
        new_eq(e_stock[t]
               == e_stock_init +
               dt * (p_bat_dc_charge[t] * eff_bat
                     - p_bat_dc_discharge[t] / eff_bat))
    else:
        new_eq(e_stock[t]
               == e_stock[t - 1] +
               dt * (p_bat_dc_charge[t] * eff_bat
                     - p_bat_dc_discharge[t] / eff_bat))

    new_eq(p_bat_dc[t] == p_bat_dc_discharge[t] - p_bat_dc_charge[t])
    # new_eq(p_bat_dc_discharge[t] * p_bat_dc_charge[t] <= 0)

    new_eq(Ebat_rated * soh * soe[t] == e_stock[t])
    if t == n_step - 1:
        new_eq(soe[t] >= soe_reserve_aux)

    new_eq(e_stock[t] <= soh * Ebat_rated)
    new_eq(p_bat_dc[t] <= p_bat_max)
    new_eq(p_bat_dc[t] >= p_bat_min)

    new_eq(p_bat_lv_discharge[t] == p_bat_dc_discharge[t] * eff_ondul_bat)
    new_eq(p_bat_lv_charge[t] * eff_ondul_bat == p_bat_dc_charge[t])
    new_eq(p_bat_lv[t] == p_bat_lv_discharge[t] - p_bat_lv_charge[t])
    # new_eq(p_bat_lv_discharge[t] * p_bat_lv_charge[t] <= 0)

    new_eq(s_bat_lv[t] <= Smax_inverter_bat)
    new_eq(s_bat_lv[t] ** 2 == p_bat_lv[t] ** 2 + q_bat_lv[t] ** 2)

    new_eq(q_bat_lv[t] == q_bat_lv_discharge[t] - q_bat_lv_charge[t])
    # new_eq(q_bat_lv_discharge[t] * q_bat_lv_charge[t] <= 0)

    new_eq(q_bat_lv_discharge[t] <= tanphi_min * p_bat_lv_discharge[t])
    new_eq(q_bat_lv_discharge[t] >= - tanphi_min * p_bat_lv_discharge[t])
    new_eq(q_bat_lv_charge[t] <= tanphi_min * p_bat_lv_charge[t])
    new_eq(q_bat_lv_charge[t] >= - tanphi_min * p_bat_lv_charge[t])

    new_eq(p_bat_hv_discharge[t] == p_bat_lv_discharge[t] * eff_transfo_bat)
    new_eq(p_bat_hv_charge[t] * eff_transfo_bat == p_bat_lv_charge[t])
    new_eq(p_bat_hv[t] == p_bat_hv_discharge[t] - p_bat_hv_charge[t])
    # new_eq(p_bat_hv_discharge[t] * p_bat_hv_charge[t] <= 0)

    new_eq(q_bat_hv_discharge[t] == q_bat_lv_discharge[t] * eff_transfo_bat)
    new_eq(q_bat_hv_charge[t] * eff_transfo_bat == q_bat_lv_charge[t])
    new_eq(q_bat_hv[t] == q_bat_hv_discharge[t] - q_bat_hv_charge[t])
    # new_eq(q_bat_hv_discharge[t] * q_bat_hv_charge[t] <= 0)


# -------------------------------
# Diesel Generator constraints --
# -------------------------------

for t in range(n_step):

    new_eq(p_gen_lv[t] <= y_gen[t] * Pgen_max * soh_gen)
    new_eq(p_gen_lv[t] >= y_gen[t] * Pgen_min)
    new_eq(s_gen_lv[t] <= Sgen_max * y_gen[t])

    new_eq(s_gen_lv[t] ** 2 == p_gen_lv[t] ** 2 + q_gen_lv[t] ** 2)

    if t == 0:
        new_eq(conso_gen[t]
               == conso_idle
               + ((conso_rated - conso_idle)
                  * (p_gen_lv[t] + p_gen_lv_init)
                  * dt
                  / (2 * Sgen_max)))
    else:
        new_eq(conso_gen[t]
               == conso_idle
               + ((conso_rated - conso_idle)
                  * ((p_gen_lv[t] + p_gen_lv[t - 1]) +
                     (q_gen_lv[t] + q_gen_lv[t - 1]))
                  * dt
                  / (2 * Sgen_max)))

    new_eq(p_gen_lv[t] * eff_transfo_gen == p_gen_hv[t])
    new_eq(q_gen_lv[t] * eff_transfo_gen == q_gen_hv[t])

# --------------------
# Other constraints --
# --------------------

for t in range(n_step):
    new_eq(p_load[t] + p_aux_hv == p_pv_hv[t] + p_gen_hv[t] + p_bat_hv[t])
    new_eq(q_load[t] == q_bat_hv[t] + q_gen_hv[t])

    # Penalty
    if t == 0:
        new_eq(smoothing_penalty[t] == 0)
    else:
        new_eq(smoothing_penalty[t]
               == ((p_gen_lv[t] - p_gen_lv[t - 1]) ** 2
                   + (q_gen_lv[t] - q_gen_lv[t - 1]) ** 2))

##########################################################
#################  OBJECTIVE FUNCTION ####################
##########################################################

array_sum = gekko_model.sum

new_eq(conso_gen_tot == array_sum(conso_gen))
new_eq(curtail_penalty == 5e-4 * array_sum(p_pv_curtail))
new_eq(penalty_tot == 1e-4 * array_sum(smoothing_penalty))
new_eq(obj == - conso_gen_tot - curtail_penalty - penalty_tot)

gekko_model.Maximize(obj)

##########################################################
##################### SOLVE SEQUENCE #####################
##########################################################

gekko_model.options.IMODE = 3
gekko_model.options.MAX_TIME = 900

#for t in range(n_step):
#    if t <= 20 or t >= 76: # Reduction of the binary problem
#        new_eq(y_gen[t] == 1)

# Initialize, check solvability
print("Problem setup OK, Initializing ...")
gekko_model.options.COLDSTART = 2
gekko_model.options.SOLVER = 3
gekko_model.solve(disp=False)

# Solve with IPOPT
gekko_model.options.SOLVER = 3
gekko_model.options.COLDSTART = 0
gekko_model.options.MAX_ITER = 10000
print('Initialization OK --> Solving NLP ...')
gekko_model.solve(disp=True)

# Solve with APOPT and continuous variables
gekko_model.options.SOLVER = 1
gekko_model.solver_options = [
    f"minlp_maximum_iterations 1000",
    f"minlp_as_nlp 1"
]
print('Solving NLP OK --> Solving NLP with APOPT ...')

#for t in range(n_step):
#    if t <= 20 or t >= 76: # Reduction of the binary problem
#        new_eq(y_gen[t] == 1)

gekko_model.solve(disp=True)


# Solve with APOPT and binary variables
gekko_model.options.SOLVER = 1
gekko_model.solver_options = [
    f"minlp_maximum_iterations 1000",
    f"minlp_gap_tol 1e-2",
    f"minlp_integer_tol 1e-2",
    f"minlp_max_iter_with_int_sol 500",
    f"objective_convergence_tolerance 1e-4",
    f"constraint_convergence_tolerance 1e-4",
    f"nlp_maximum_iterations 500",
    f"minlp_branch_method 1",
    f"minlp_as_nlp 0"
]
print('Solving NLP OK --> Solving MINLP with APOPT ...')

#for t in range(n_step):
#    if t <= 20 or t >= 76: # Reduction of the binary problem
#        new_eq(y_gen[t] == 1)

gekko_model.solve(disp=True)


##########################################################
########### SAVING VARIABLES IN DATAFRAME ###############
##########################################################

var_dict = {
    "__time__": {
        "time": n_step,
    },
    "__list__": {
        "p_mpp_dc": p_mpp_dc,
        "p_mpp_ac_hv": p_mpp_ac_hv,
        "Pload": p_load,
        "Qload": q_load,
    },
    "__array_1D__": {
        "p_pv_dc": p_pv_dc,
        "p_pv_curtail": p_pv_curtail,
        "p_pv_ac_lv": p_pv_lv,
        "p_pv_ac_hv": p_pv_hv,
        "p_bat_dc": p_bat_dc,
        "p_bat_ac_lv": p_bat_lv,
        "q_bat_ac_lv": q_bat_lv,
        "s_bat_ac_lv": s_bat_lv,
        "p_bat_ac_hv": p_bat_hv,
        "q_bat_ac_hv": q_bat_hv,
        "soe": soe,
        "e_stock": e_stock,
        "y_gen": y_gen,
        "p_gen_ac_lv": p_gen_lv,
        "q_gen_ac_lv": q_gen_lv,
        "s_gen_ac_lv": s_gen_lv,
        "p_gen_ac_hv": p_gen_hv,
        "q_gen_ac_hv": q_gen_hv,
        "conso_gen": conso_gen,
        "p_bat_dc_charge": p_bat_dc_charge,
        "p_bat_dc_discharge": p_bat_dc_discharge,
        "p_bat_lv_charge": p_bat_lv_charge,
        "p_bat_lv_discharge": p_bat_lv_discharge,
        "q_bat_lv_charge": q_bat_lv_charge,
        "q_bat_lv_discharge": q_bat_lv_discharge,
        "p_bat_hv_charge": p_bat_hv_charge,
        "p_bat_hv_discharge": p_bat_hv_discharge,
        "q_bat_hv_charge": q_bat_hv_charge,
        "q_bat_hv_discharge": q_bat_hv_discharge,
    },
    "__array_nD__": {
    },
    "__scalar__": {
        "conso_gen_tot": conso_gen_tot
    },
}

t0 = index_start = pd.to_datetime('2019-08-15 00:00:00')
n_step = var_dict["__time__"]["time"]
time_index = pd.DatetimeIndex(
    [t0 + n * QUARTER_HOUR for n in range(n_step)],
    freq=QUARTER_HOUR,
)
result = pd.DataFrame(index=time_index, columns=[])

for name, value in var_dict["__list__"].items():
    result[name] = value

for name, value in var_dict["__array_1D__"].items():
    result[name] = [inner_list[0] for inner_list in value]

for name, value in var_dict["__array_nD__"].items():
    for dim in range(1, value.shape[1] + 1):
        inner_list = list(value[:, dim - 1])
        for p in range(n_step):
            inner_list[p] = inner_list[p].VALUE[0]

        result[f"{name}_{dim}"] = inner_list

for name, gekko_value in var_dict["__scalar__"].items():
    result[name] = gekko_value.value[0]

##########################################################
#################### PLOT RESULTS ########################
##########################################################

# Tracé des courbes principales
fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True)
ax1.plot(result[result.columns[2]], label="Pload")
ax1.plot(result[result.columns[7]], label="p_pv_ac_hv")
ax1.plot(result[result.columns[12]], label="p_bat_ac_hv")
ax1.plot(result[result.columns[20]], label="p_gen_ac_hv")
ax2.plot(result[result.columns[14]], label="soe")
ax2.plot(result[result.columns[16]], label="y_gen")
ax3.plot(result[result.columns[3]], label="Qload")
ax3.plot(result[result.columns[13]], label="q_bat_ac_hv")
ax3.plot(result[result.columns[21]], label="q_gen_ac_hv")

ax1.get_shared_x_axes().join(ax1, ax2, ax3)
ax1.set_xticklabels([])
box = ax1.get_position()
ax1.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax1.legend(loc='center left', bbox_to_anchor=(1, 0.5))
box = ax2.get_position()
ax2.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax2.legend(loc='center left', bbox_to_anchor=(1, 0.5))
box = ax3.get_position()
ax3.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax3.legend(loc='center left', bbox_to_anchor=(1, 0.5))

plt.show()

Thank you very much to @John Hedengren for their advice. The only thing is, this is being developed for the purpose of an Energy Management System providing setpoints for a model predictive controller every fifteen minutes. So, it has to be solved relatively quickly.

The solution I found is the following. I added a quadratic penalization to the objective function in the following form:

y_gen[t] * (1 - y_gen[t]))

This penalizes non-binary values. Then, I solved my problem a first time with IPOPT. I recovered the result, post-treated it, and solved it again with IPOPT but with a "binary" variable which is parameterized (array of 0 or 1 coming from the post-treated previous result).

I think this is nice way to solve my problem quickly even though the penalization is neither linear nor convex. Each solve iteration takes less than a minute.