A piece-wise defined objective function can be solved by Gurobi?

Specific Domains Optimization (Mathematical)
1

Hi, I have an optimization with a piece-wise defined objective as follows:

using JuMP, Ipopt
using Gurobi

T_fc = [0, 185, 348, 515, 658, 830]
R_fc = [0.1803, 0.1978, 0.2033, 0.2092, 0.2157, 0.2254]

const R0 = R_fc[1]
const Imin, Inom, Imax = 0.2, 0.7, 2
Lmin, Lnom, Lmax = 0.8305, 2.38, 3.1

const FT = 0.1803 * (1+5.39*0.1)
Rfc1 = 0.1803 + 0.01
Rfc2 = 0.1803
L11, L21 = 2.1, 2.1
Dindf = [250, 250]
Ldf = [5, 4]

function MoM(T, R)
    det_R = [r-R[i] for (i, r) in enumerate(R[2:end])]
    det_w = [t-T[j] for (j, t) in enumerate(T[2:end])]
    me_R = (R[end] - R[1]) / T[end]

    A = R[end]*(1 - sum(det_w.^2) / (sum(det_w)^2))
    B = sum((det_R - me_R * det_w).^2)

    # Based on the equations from reference, calculate the estimates for alpha and beta
    beta_est = B / A
    alpha_est = me_R / beta_est
    return alpha_est, beta_est
end
C = 20
α_est, β_est = MoM(T_fc, R_fc)
const α_0, β_0 = α_est/C, β_est*C

K_dL = 30 * 5.39 # constant to tune load varying casued R aging rate -> KdL
Dv_dL = K_dL * 5.93e-7 * R0 / (Lmax - Lmin)

# weight factors
w1 = (FT - Rfc2) / (2FT - Rfc1 - Rfc2)
w2 = 1 - w1

function α(x)
    a = 0.0019727939
    b = 0.0078887
    # α_0 = 0.006226
    if x< Lnom
        return a * (x-Lnom)^2 + α_0
    else
        return b * (x-Lnom)^2 + α_0
    end
end
# objective function
function obj(x, y)
    f0dL = (w1 * (x*Ldf[1]-L11)^2 + w2 * (y*Ldf[1]-L21)^2) * Dv_dL * 3600
    return (w1 * α(x*Ldf[1])  + w2*α(y*Ldf[1]))*Dindf[1] + f0dL
end

# model = Model(Ipopt.Optimizer)
model = Model(Gurobi.Optimizer)

lb = maximum(Lmin ./ Ldf)
ub = minimum(Lmax ./ Ldf)
@variable(model, lb<=x<=ub, start=0.5)
@variable(model, lb<=y<=ub, start=0.5)

# @NLconstraint(model, x+y==1)
@constraint(model, x+y==1)

register(model, :obj, 2, obj; autodiff=true)
# @NLobjective(model, Min, obj(x,y))
@objective(model, Min, obj(x,y))

optimize!(model)
println(value(x*Ldf[1]), value(y*Ldf[1]))

The Ipopt solver works fine, but when I tried with Gurobi, I received an error as:

MethodError: no method matching isless(::AffExpr, ::Float64)
Closest candidates are:
  isless(::T, ::T) where T<:Union{Float16, Float32, Float64} at ~/.julia/juliaup/julia-1.7.2+0~x64/share/julia/base/float.jl:460
  isless(::Union{StatsBase.PValue, StatsBase.TestStat}, ::AbstractFloat) at ~/.julia/packages/StatsBase/pJqvO/src/statmodels.jl:99
  isless(::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at ~/.julia/packages/StatsBase/pJqvO/src/statmodels.jl:90
  ...

I have checked the Gurobi website, it seems there is no example for defining such an objective function (piece-wise and nonlinear). Can this problem be solved by Gurobi?
Thank you very much for checking my problem.

2

You cannot use Gurobi for this directly because you cannot use @objective with a user defined function, and Gurobi does not support the @NL macros.

3

Ok, I see. Thank you.
I need to solve such optimization problems repeatedly and the Ipopt
the solver takes some time (maybe I didn’t write it efficiently).

4

How about something like this:

using JuMP
import Ipopt

function MoM(T, R)
    det_R = [r - R[i] for (i, r) in enumerate(R[2:end])]
    det_w = [t - T[j] for (j, t) in enumerate(T[2:end])]
    me_R = (R[end] - R[1]) / T[end]
    A = R[end] * (1 - sum(det_w.^2) / sum(det_w)^2)
    B = sum((det_R - me_R * det_w).^2)
    # Based on the equations from reference, calculate the estimates for alpha 
    # and beta
    beta_est = B / A
    alpha_est = me_R / beta_est
    return alpha_est, beta_est
end

function main()
    T_fc = [0, 185, 348, 515, 658, 830]
    R_fc = [0.1803, 0.1978, 0.2033, 0.2092, 0.2157, 0.2254]
    R0 = R_fc[1]
    Imin, Inom, Imax = 0.2, 0.7, 2
    Lmin, Lnom, Lmax = 0.8305, 2.38, 3.1
    FT = 0.1803 * (1+5.39*0.1)
    Rfc1 = 0.1803 + 0.01
    Rfc2 = 0.1803
    L11, L21 = 2.1, 2.1
    Dindf = [250, 250]
    Ldf = [5, 4]
    C = 20
    α_est, β_est = MoM(T_fc, R_fc)
    α_0, β_0 = α_est / C, β_est * C
    K_dL = 30 * 5.39 # constant to tune load varying casued R aging rate -> KdL
    Dv_dL = K_dL * 5.93e-7 * R0 / (Lmax - Lmin)
    w1 = (FT - Rfc2) / (2FT - Rfc1 - Rfc2)
    w2 = 1 - w1
    a, b = 0.0019727939, 0.0078887
    lb = maximum(Lmin ./ Ldf)
    ub = minimum(Lmax ./ Ldf)
    # Model
    model = Model(Ipopt.Optimizer)
    @variable(model, lb <= x <= ub, start = 0.5)
    @variable(model, lb <= y <= ub, start = 0.5)
    @constraint(model, x + y == 1)
    @NLexpression(
        model,
        f0dL,
        (w1 * (x * Ldf[1] - L11)^2 + w2 * (y * Ldf[1] - L21)^2) * Dv_dL * 3600,
    )
    @NLexpression(
        model,
        α_x,
        ifelse(
            x * Ldf[1] < Lnom,
            a * (x * Ldf[1] - Lnom)^2 + α_0,
            b * (x * Ldf[1] - Lnom)^2 + α_0,
        ),
    )
    @NLexpression(
        model,
        α_y,
        ifelse(
            y * Ldf[1] < Lnom,
            a * (y * Ldf[1] - Lnom)^2 + α_0,
            b * (y * Ldf[1] - Lnom)^2 + α_0,
        ),
    )
    @NLobjective(model, Min, Dindf[1] * (w1 * α_x  + w2 * α_y) + f0dL)
    optimize!(model)
    return value(x*Ldf[1]), value(y*Ldf[1])
end
main()

I think you have a few typos, for example it should probably be β_0 in the ifelse?

5

Thank you very much for providing this demo.
It helps me a lot. I am going to reorganize my program.
I will recheck the possible typos as you mentioned, thanks!

1 Like
6

I have tried your solution, it was about 10 s fast than my previous Ipopt program.
So it is better to put variables as a local within a function and use @NLexpression
to configure the objective function?
Thanks again @odow

7

If you use a user-defined function we disable second derivative information. So using the macros is preferred.

Avoiding global variables is good practice for Julia.

8

@odow Thanks a lot for the explanation!

1 Like