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

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
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) / 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-L11)^2 + w2 * (y*Ldf-L21)^2) * Dv_dL * 3600
return (w1 * α(x*Ldf)  + w2*α(y*Ldf))*Dindf + 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), value(y*Ldf))
``````

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.

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.

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

``````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) / 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
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 - L11)^2 + w2 * (y * Ldf - L21)^2) * Dv_dL * 3600,
)
@NLexpression(
model,
α_x,
ifelse(
x * Ldf < Lnom,
a * (x * Ldf - Lnom)^2 + α_0,
b * (x * Ldf - Lnom)^2 + α_0,
),
)
@NLexpression(
model,
α_y,
ifelse(
y * Ldf < Lnom,
a * (y * Ldf - Lnom)^2 + α_0,
b * (y * Ldf - Lnom)^2 + α_0,
),
)
@NLobjective(model, Min, Dindf * (w1 * α_x  + w2 * α_y) + f0dL)
optimize!(model)
return value(x*Ldf), value(y*Ldf)
end
main()
``````

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

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

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

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.

@odow Thanks a lot for the explanation!

1 Like