Problem using JuMP with user-defined function

Specific Domains Optimization (Mathematical)
,
1

I am trying to use JuMP to solve a Non-linear optimization problem, but find a problem using the following code:

function Q(H,h_0,C,D,L)
    
    if H <= h_0
        return 0.0
    else
        Δh = H - h_0
        J = Δh/L
        return (J * C^(1.85) * D^(4.87)/10.643)^(1/1.85)*3600
    end
end

using JuMP
using Ipopt

model = Model(Ipopt.Optimizer)

## Parametros
@variables(model,begin
L_1 == 3e3
L_2 == 2e3
D_1 == 250.0e-3
D_2 == 200.0e-3
C_1 == 120.0
C_2 == 120.0
h_1 == 12.0
h_2 == 5.0
end)

## Variáveis, H_t é a independente

@variables(model,begin 
    H_t
    q_1 = Q(H_t,h_1,C_1,D_1,L_1)
    q_2 = Q(H_t,h_2,C_2,D_2,L_2)
    q_t = q_1 + q_2
    H_b = H_m(q_t)
end)


@NLobjective(model, Min, (H_b-H_t)^2)
@constraint(model, c1,
begin
    H_t >=0.0,
    q_1 >=0.0,
    q_2 >=0.0,
    H_b >=0
end)

print(model)
optimize!(model)

@show termination_status(model)
@show primal_status(model)
@show dual_status(model)
@show objective_value(model)

@show value(H_t)
@show value(H_b)
@show value(q_1)
@show value(q_2)
@show value(q_t)

@show shadow_price(c1)

This return the error:

ERROR: LoadError: MethodError: no method matching isless(::VariableRef, ::VariableRef)
Closest candidates are:
  isless(::Any, ::Missing) at missing.jl:88
  isless(::Missing, ::Any) at missing.jl:87
Stacktrace:
 [1] <(x::VariableRef, y::VariableRef)
   @ Base ./operators.jl:279
 [2] <=(x::VariableRef, y::VariableRef)
   @ Base ./operators.jl:328
 [3] Q(H::VariableRef, h_0::VariableRef, C::VariableRef, D::VariableRef, L::VariableRef)
   @ Main ~/Documentos/ITA/6º Semestre/Obrigatórias/HID-32/P1/P1.jl:17
 [4] top-level scope
   @ ~/Documentos/ITA/6º Semestre/Obrigatórias/HID-32/P1/P1.jl:134
in expression starting at /home/vinicius/Documentos/ITA/6º Semestre/Obrigatórias/HID-32/P1/P1.jl:132

There is a way to do this with JuMP?

2

Hi there. The unhelpful error is a shame. I’ll take a look to see if we could improve it.

But it occurs because there are quite a few problems with your syntax (It looks like you’ve used AMPL before; JuMP has quite a few differences).

You should read the JuMP documentation:

Notable differences from AMPL include:

  • Create expressions using @NLexpression and @expression instead of variable = expr
  • Create nonlinear parameters using @NLparameter instead of variable == value
  • You need to register your user-defined functions
  • @constraint should be @constraints, and the c1 name is in the wrong place
  • No commas between lines
  • Add variables bounds in @variables instead of as a @constraint

Here’s how I would write your model. It’s not exactly the same because I don’t know what H_m is.

function Q(H, h_0, C, D, L)
    if H <= h_0
        return 0.0
    end
    Δh = H - h_0
    J = Δh / L
    return (J * C^(1.85) * D^(4.87) / 10.643)^(1 / 1.85) * 3_600
end

using JuMP
using Ipopt

model = Model(Ipopt.Optimizer)
register(model, :Q, 5, Q; autodiff = true)
@variable(model, H_t >= 0)
@NLparameters(model, begin
    L_1 == 3e3
    L_2 == 2e3
    D_1 == 250.0e-3
    D_2 == 200.0e-3
    C_1 == 120.0
    C_2 == 120.0
    h_1 == 12.0
    h_2 == 5.0
end)
@NLexpressions(model, begin
    q_1, Q(H_t, h_1, C_1, D_1, L_1)
    q_2, Q(H_t,h_2,C_2,D_2,L_2)
    q_t, q_1 + q_2
    H_b, q_t  # What is H_m?
end)
@NLconstraints(model, begin
    q_1 >= 0
    q_2 >= 0
    H_b >= 0
end)
@NLobjective(model, Min, (H_b - H_t)^2)
print(model)
optimize!(model)
1 Like