Couenne no results - seems stuck

Specific Domains Optimization (Mathematical)
1

Hi all,

I am using Couenne for my MINLP problem. I have a problem using Couenne when I add a nonlinear term to the objective. I have tried different options in the couenne.opt file based on different posts here, still no results (Clp0000I Optimal - objective value 0). Here is a simple version of my code with only the problematic objective function:

using JuMP
using AmplNLWriter

m = Model(with_optimizer(AmplNLWriter.Optimizer, "path to couenne.exe"))
#m = Model(with_optimizer(Ipopt.Optimizer))


    
    min_velocity = 0.6
    max_velocity = 1.5
    nu = 1e-6
    g = 9.81
    rho = 1000
    pump_efficiency = 0.6
    Cost_per_kWh = 0.24 #Dollars/kWh
    pipe_roughness = 0.00015
    
    Di_standard = 1e-3*[0; 32; ;40 ;  50;  63;  75; 90; 110; 125; 140; 160; 180; 200; 225; 250]
    Price_per_Unit_pipe = [0; 1.08; 1.72; 2.65; 4.22; 5.95; 8.56; 12.74; 16.25; 20.40; 26.51; 33.68; 41.56; 52.71; 64.60]
    QP = [0.06, 0.04]
    Dis = [20 100; 100 20]

    NP = size(QP,1) # number of units
    NDi = size(Di_standard,1) # number of available standard diameters

    @variable(m, x[1:NP, 1:NP] >= 0)
    @variable(m, Di[1:NP, 1:NP] >= 0)
    @variable(m, QFW[1:NP] >= 0)

    @variable(m, v[1:NP, 1:NP] >= 0)
    @variable(m, 0 <= Re[1:NP, 1:NP])
    @variable(m, 0 <= f[1:NP, 1:NP])
    @variable(m, 0 <= hf[1:NP, 1:NP])

    @variable(m, Di_idx[1:NP, 1:NP, 1:NDi], Bin)
    @variable(m, x_idx[1:NP, 1:NP], Bin)

    @constraint(m, 0 .<= QP - x'*QP)

    @constraint(m, 1 .>= sum(x, dims=2))
    @constraint(m, 1 .== sum(Di_idx, dims=3))

    @constraint(m, [i = 1:NP, j = 1:NP], Di[i,j] == sum(Di_standard[n]*Di_idx[i,j,n] for n = 1:NDi))

    @constraint(m, [i = 1:NP, j = 1:NP], QFW[i] == QP[i] - sum(x[j,i]*QP[j] for j = 1:NP))
    @constraint(m, sum(QFW[i] for i =1:NP)/sum(QP[i] for  i = 1:NP) == 0.2)

    @NLconstraint(m, [i = 1:NP, j = 1:NP], v[i,j] == 4*x[i,j]*QP[i]/(pi*Di[i,j]*Di[i,j]))

    #Big-M formulation

    @constraint(m, x .<= x_idx)
    @constraint(m, x .>= 0)

    @variable(m, Di_aux[1:NP, 1:NP])
    @constraint(m, Di - Di_aux .<= (1-x_idx))
    @constraint(m, 0 .<= Di - Di_aux)
    @constraint(m, 0 .<= Di_aux )
    @constraint(m, Di_aux .<= x_idx)


    M1 = 10
    @variable(m, v_aux[1:NP, 1:NP])
    @constraint(m, v - v_aux .<=  M1*(1-x_idx))
    @constraint(m, 0 .<= v - v_aux)
    @constraint(m, min_velocity*x_idx .<= v_aux )
    @constraint(m, v_aux.<= max_velocity*x_idx)


    @NLconstraint(m, [i = 1:NP, j = 1:NP], Re[i,j] == Di[i,j]*v[i,j]/nu)
    @NLconstraint(m, [i = 1:NP, j = 1:NP], f[i,j] == 0.316/(Re[i,j]^0.25)*x_idx[i,j])
    @variable(m, f_aux[1:NP, 1:NP])
    @constraint(m, f- f_aux .<= (1-x_idx))
    @constraint(m, 0 .<= f - f_aux)
    @constraint(m, 0 .<= f_aux)
    @constraint(m, f_aux .<= x_idx)

    M2 = 10
    @NLconstraint(m, [i = 1:NP, j = 1:NP], hf[i,j] == (f[i,j]/Di[i,j])*Dis[i,j]*(v[i,j]^2/(2*g))) # Darcy equation
    @variable(m, hf_aux[1:NP, 1:NP])
    @constraint(m, hf- hf_aux .<= M2*(1-x_idx))
    @constraint(m, 0 .<= hf - hf_aux)
    @constraint(m, 0 .<= hf_aux)
    @constraint(m, hf_aux .<= M2*x_idx)

    @NLexpression(m, energy_cost[i = 1:NP, j = 1:NP], rho*x[i,j]*QP[i]*hf[i,j]*g/pump_efficiency*Cost_per_kWh/1000)


    @NLobjective(m, Min, sum(energy_cost[i,j] for i = 1:NP, j = 1:NP)/sum(QP[i] for i = 1:NP, j = 1:NP))

    JuMP.optimize!(m)
2

Please post the output, including what options you have tried, how long you ran it for, etc.

In addition, try to simplify the model as much as possible while retaining the problematic behavior (e.g., by removing constants and variables and constraints if possible).

You may want to take a read of Please read: make it easier to help you.

3

Thanks for the reply. I have tried different gaps “allowable_fraction_gap 1.0 or 5.0” for now. The code runs for some time for such a simple case and when I use a time_limit of 1000 I get NaN results.

Cbc0010I After 52100 nodes, 25548 on tree, 1e+050 best solution, best possible 0

even when I remove a constraint that should give me the solution of all variables equal to zero, COuenne cannot do it:
#@constraint(m, sum(QFW[i] for i =1:NP)/sum(QP[i] for i = 1:NP) == 0.2)

Here is a more simplified version of the code:

min_velocity = 0.6
max_velocity = 1.5

Di_standard = 1e-3*[0; 32; ;40 ;  50;  63;  75; 90; 110; 125; 140; 160; 180; 200; 225; 250]
Price_per_Unit_pipe = [0; 1.08; 1.72; 2.65; 4.22; 5.95; 8.56; 12.74; 16.25; 20.40; 26.51; 33.68; 41.56; 52.71; 64.60]
QP = [0.06, 0.04]
Dis = [20 100; 100 20]

NP = size(QP,1) # number of units
NDi = size(Di_standard,1) # number of available standard diameters

@variable(m, x[1:NP, 1:NP] >= 0)
@variable(m, Di[1:NP, 1:NP] >= 0)
@variable(m, QFW[1:NP] >= 0)

@variable(m, v[1:NP, 1:NP] >= 0)
@variable(m, 0 <= Re[1:NP, 1:NP])
@variable(m, 0 <= f[1:NP, 1:NP])
@variable(m, 0 <= hf[1:NP, 1:NP])

@variable(m, Di_idx[1:NP, 1:NP, 1:NDi], Bin)
@variable(m, x_idx[1:NP, 1:NP], Bin)

@constraint(m, 0 .<= QP - x'*QP)

@constraint(m, 1 .>= sum(x, dims=2))
@constraint(m, 1 .== sum(Di_idx, dims=3))

@constraint(m, [i = 1:NP, j = 1:NP], Di[i,j] == sum(Di_standard[n]*Di_idx[i,j,n] for n = 1:NDi))

@constraint(m, [i = 1:NP, j = 1:NP], QFW[i] == QP[i] - sum(x[j,i]*QP[j] for j = 1:NP))
@constraint(m, sum(QFW[i] for i =1:NP)/sum(QP[i] for  i = 1:NP) == 0.2)

@NLconstraint(m, [i = 1:NP, j = 1:NP], v[i,j] == 4*x[i,j]*QP[i]/(pi*Di[i,j]*Di[i,j]))

#Big-M formulation

@constraint(m, x .<= x_idx)
@constraint(m, x .>= 0)

@variable(m, Di_aux[1:NP, 1:NP])
@constraint(m, Di - Di_aux .<= (1-x_idx))
@constraint(m, 0 .<= Di - Di_aux)
@constraint(m, 0 .<= Di_aux )
@constraint(m, Di_aux .<= x_idx)


M1 = 10
@variable(m, v_aux[1:NP, 1:NP])
@constraint(m, v - v_aux .<=  M1*(1-x_idx))
@constraint(m, 0 .<= v - v_aux)
@constraint(m, min_velocity*x_idx .<= v_aux )
@constraint(m, v_aux.<= max_velocity*x_idx)


@NLconstraint(m, [i = 1:NP, j = 1:NP], Re[i,j] == Di[i,j]*v[i,j]/10e-6)
@NLconstraint(m, [i = 1:NP, j = 1:NP], f[i,j] == 0.316/(Re[i,j]^0.25)*x_idx[i,j])
@variable(m, f_aux[1:NP, 1:NP])
@constraint(m, f- f_aux .<= (1-x_idx))
@constraint(m, 0 .<= f - f_aux)
@constraint(m, 0 .<= f_aux)
@constraint(m, f_aux .<= x_idx)

M2 = 10
@NLconstraint(m, [i = 1:NP, j = 1:NP], hf[i,j] == (f[i,j]/Di[i,j])*Dis[i,j]*(v[i,j]^2/(2*9.81))) 
@variable(m, hf_aux[1:NP, 1:NP])
@constraint(m, hf- hf_aux .<= M2*(1-x_idx))
@constraint(m, 0 .<= hf - hf_aux)
@constraint(m, 0 .<= hf_aux)
@constraint(m, hf_aux .<= M2*x_idx)

@NLexpression(m, energy_cost[i = 1:NP, j = 1:NP], x[i,j]*QP[i]*hf[i,j])


@NLobjective(m, Min, sum(energy_cost[i,j] for i = 1:NP, j = 1:NP)/sum(QP[i] for i = 1:NP, j = 1:NP))

JuMP.optimize!(m)