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

**URL:** <https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [April 3, 2022, 5:26pm UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959 "2022-04-03T17:26:29Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [April 3, 2022, 5:26pm UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/1 "2022-04-03T17:26:29Z")

</div>

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

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 3, 2022, 6:56pm UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/2 "2022-04-03T18:56:30Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [April 3, 2022, 7:26pm UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/3 "2022-04-03T19:26:23Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 3, 2022, 7:42pm UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/4 "2022-04-03T19:42:58Z")

</div>

How about something like this:

```julia
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`?

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [April 3, 2022, 8:09pm UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/5 "2022-04-03T20:09:29Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [April 4, 2022, 6:21am UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/6 "2022-04-04T06:21:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 4, 2022, 7:04am UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/7 "2022-04-04T07:04:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [April 4, 2022, 7:44am UTC](https://discourse.julialang.org/t/a-piece-wise-defined-objective-function-can-be-solved-by-gurobi/78959/8 "2022-04-04T07:44:04Z")

</div>

@odow Thanks a lot for the explanation!