How to write conditional expression using JuMP

Specific Domains Optimization (Mathematical)
1

Hello,

I am using JuMP to solve dynamic programming problem. My constraints are flow and level, which are piecewise and depend on each other. level[r, i] = f(flow[r, i], level[r, i-1]) , flow[r, i] = f(level[r, i], n). The expression should be like as follows. In my case, the only variable is n[r = 1:N, i = 1:numtime]:

for r in 1:N
    for i in 1:numtime
        if i == 1
            # Discharge flow for the 1st time segment from upper reservoir
            flow[r, 1] = (0.685 - 0.19 * n[r, 1] * 0.05 / (upstream[r].initial_level - 
                upstream[r].coef_level)) * 0.5 * n[r, 1] * sqrt(2 * 9.81 * 
                (upstream[r].initial_level - upstream[r].coef_level))
            
            # Water lever of reservoir for the 1st time segment
            level[r, 1] = upstream[r].initial_level + 
                (upstream[r].inflow * raintime[1] * 3600 - flow[r, 1] * 43200) / upstream[r].k
        else
            flow[r, i] == (0.685 - 0.19 * n[r, i] * 0.05 / (level[r, i-1] - 
                upstream[r].coef_level)) * 0.5 * n[r, i] * sqrt(2 * 9.81 * 
                (level[r, i-1] - upstream[r].coef_level))

            level[r, i] == level[r, i-1] + (upstream[r].inflow * raintime[i] * 3600 
                - flow[r, i] * 43200) / upstream[r].k
        end
    end
end

with this written in expression, I want to define the constraints:

# constraint of discharge flow from upstream reservoir
@constraint(
    model, 
    upstream[r].min_flow <= flow[r = 1:N, 1:numtime] <= upstream[r].max_flow
)
# constraint of water flow of upstream reservoir
@constraint(
    model, 
    upstream[r].min_level <= level[r = 1:N, 1:numtime] <= upstream[r].max_level
)

Question: How can I write flow and level using @expression so that the constraints can be defined? Could you please give me any suggestion?

Thanks!

2

Take a read of Please read: make it easier to help you. It’s easier to help if you can provide a reproducible example that someone can copy-paste. If you don’t have that, then simplify the problem as much as possible and provide what you want to do in math. Your question is a little confusing, because there is a lot of other things going on.

What are your decision variables? What is the data? What is the math formulation you’re trying to achieve? I don’t see the piecewise linearity (Piecewise linear function - Wikipedia), only two nonlinear constraints?

You should introduce flow and level as decision variables. Then you’ll need to add the constraints you have using @NLconstraint. You can’t use @constraint or @expression because you have a 1 / level term in there.

1 Like
3

Well, what I am trying to do is as follows:

using JuMP
using AmplNLWriter, Bonmin_jll

struct Reservoir
    initial_level::Float64             
    min_level::Float64              
    max_level::Float64                 
    coef_level::Float64                 
    min_flow::Float64                
    max_flow::Float64              
    d::Float64                           
    k::Float64                           
    inflow::Float64                      
    # states::Dict{Symbol, State}
end

numtime = 6
raintime = [8, 5]
upstream = [
    Reservoir(
        273.00,
        230.00,
        275.17,
        261.00,
        0,
        3140,
        57000,
        2.090104051e7,
        500
    ),
    Reservoir(
        225.00,
        204.00,
        230.80,
        221.00,
        0,
        3140,
        32000,
        2.143410853e8,
        760
    )
)
N = length(upstream)

model = Model(() -> AmplNLWriter.Optimizer(Bonmin_jll.amplexe))

# number of units of spillway hole opening
@variable(
    model, 
    0 <= n[1:N, 1:numtime] <= 60,
    Int
)   
for r in 1:N
    for i in 1:numtime
        if i == 1
            # Discharge flow for the 1st time segment from upper reservoir
            expr_flow[r, 1] = @expression(
                model,
                (0.685 - 0.19 * n[r, 1] * 0.05 / (upstream[r].initial_level - 
                upstream[r].coef_level)) * 0.5 * n[r, 1] * sqrt(2 * 9.81 * 
                (upstream[r].initial_level - upstream[r].coef_level))
            )
            # Water lever of reservoir for the 1st time segment
            expr_level[r, 1] = @expression(
                model,
                upstream[r].initial_level + 
                (upstream[r].inflow * raintime[1] * 3600 - flow[r, 1] * 43200) / upstream[r].k
            )
        else
            expr_flow[r, i] = @expression(
                model,
                (0.685 - 0.19 * n[r, i] * 0.05 / (level[r, i-1] - 
                upstream[r].coef_level)) * 0.5 * n[r, i] * sqrt(2 * 9.81 * 
                (level[r, i-1] - upstream[r].coef_level))
            )
            expr_level[r, i] = @expression(
                model,
                level[r, i-1] + (upstream[r].inflow * raintime[i] * 3600 
                - flow[r, i] * 43200) / upstream[r].k
            )
        end
    end
end

# constraint of discharge flow from upstream reservoir
@constraint(
    model, 
    upstream[r].min_flow <= flow[r = 1:N, 1:numtime] <= upstream[r].max_flow
)
# constraint of water flow of upstream reservoir
@constraint(
    model, 
    upstream[r].min_level <= level[r = 1:N, 1:numtime] <= upstream[r].max_level
)
# objective function 目标函数
@objective(model, Min, sum(flow[r, i] / upstream[r].d for r in 1:N for i in 1:numtime))

# print(model)
optimize!(model)

The problem is that I don’t know how to define level and flow

4

Yes, you’re right. This is not the piecewise expression. The expression is conditional with index of the vector

5

Introduce them as decision variables using @variable. Then re-write your expressions as constraints using @NLconstraint.

6

You mean constraints like this one?

@NLconstraints(
    model,
    begin 
        # Discharge flow for the 1st time segment from upper reservoir
        [r = 1:N, i = 1:1],
        flow[r, 1] == (0.685 - 0.19 * n[r, 1] * 0.05 / (upstream[r].initial_level - 
        upstream[r].coef_level)) * 0.5 * n[r, 1] * sqrt(2 * 9.81 * 
        (upstream[r].initial_level - upstream[r].coef_level))
        
        # Water lever of reservoir for the 1st time segment
        [r = 1:N, i = 1:1],
        level[r, 1] == upstream[r].initial_level + 
        (upstream[r].inflow * raintime[1] * 3600 - flow[r, 1] * 43200) / upstream[r].k

        # Discharge flow 
        [r = 1:N, i = 2:numtime],
        flow[r, i] == (0.685 - 0.19 * n[r, i] * 0.05 / (level[r, i-1] - 
        upstream[r].coef_level)) * 0.5 * n[r, i] * sqrt(2 * 9.81 * 
        (level[r, i-1] - upstream[r].coef_level))
        
        # water lever of reservoir
        [r = 1:N, i = 2:numtime],
        level[r, i] == level[r, i-1] + (upstream[r].inflow * raintime[i] * 3600 
        - flow[r, i] * 43200) / upstream[r].k
    end
)

Infact, I have tried in this way. But I got my problem locally solved. And when I change a little the input(order of 0.1), the zero solution quickly turns to be infeasible. So I am thinking there might be sth wrong with my problem definition

7

You mean constraints like this one?

Yes, exactly like that.

But I got my problem locally solved

Yes. termination_status will be LOCALLY_SOLVED.

when I change a little the input(order of 0.1), the zero solution quickly turns to be infeasible. So I am thinking there might be sth wrong with my problem definition

I can’t offer more suggestions because I can’t run your code, but is there a feasible solution? Did you provide a feasible starting point? You can provide starting points with @variable(model, x, start = 1.2).

8

Did you provide a feasible starting point?

Yes, This is done with set_start_value.() for n, level and flow.

The problem is that my code seems to be only able to get zeros optimal solution. For example, with raintime = [12] , there is a feasible solution, n = [0], which means there is no need to discharge flow for reservoir. But with raintime = [13], the solver can not get feasible solution. I have checked with Excel, when raintime = [13], one possible solution is n = [7].

With this producible code, it’s very nice of you to give me more advice.

using JuMP
using AmplNLWriter, Bonmin_jll

struct Reservoir
    initial_level::Float64              
    min_level::Float64                   
    max_level::Float64                   
    coef_level::Float64                  
    min_flow::Float64                    
    max_flow::Float64                    
    d::Float64                           
    k::Float64                          
    inflow::Float64                     
end

numtime = 1
raintime = [13]

# constant of upstream reservoirs
upstream = [
    Reservoir(
        273.00,
        230.00,
        275.17,
        261.00,
        0,
        3140,
        57000,
        2.090104051e7,
        1000
    ),
]
N = length(upstream)

# optimizer
model = Model(() -> AmplNLWriter.Optimizer(Bonmin_jll.amplexe))

# variable definition
# number of units of spillway hole opening
@variable(
    model, 
    0 <= n[1:N, 1:numtime] <= 60,
    Int
)   
# discharge flow from upstream reservoir
@variable(
    model, 
    upstream[r].min_flow <= flow[r = 1:N, 1:numtime] <= upstream[r].max_flow
)
# water flow of upstream reservoir
@variable(
    model, 
    upstream[r].min_level <= level[r = 1:N, 1:numtime] <= upstream[r].max_level
)

# n0 = 0 * ones(N, numtime)
flow_0 = zeros(N, numtime)
level_0 = upstream[1].initial_level * ones(N, numtime)
n0 = [0]

set_start_value.(n, n0)
set_start_value.(flow, flow_0)
set_start_value.(level, level_0)

@NLconstraints(
    model,
    begin 
        # Discharge flow for the 1st time segment from upper reservoir
        [r = 1:N, i = 1:1],
        flow[r, 1] == (0.685 - 0.19 * n[r, 1] * 0.05 / (upstream[r].initial_level - 
        upstream[r].coef_level)) * 0.5 * n[r, 1] * sqrt(2 * 9.81 * 
        (upstream[r].initial_level - upstream[r].coef_level))
        
        # Water lever of reservoir for the 1st time segment
        [r = 1:N, i = 1:1],
        level[r, 1] == upstream[r].initial_level + 
        (upstream[r].inflow * raintime[1] * 3600 - flow[r, 1] * 43200) / upstream[r].k

        # Discharge flow 
        [r = 1:N, i = 2:numtime],
        flow[r, i] == (0.685 - 0.19 * n[r, i] * 0.05 / (level[r, i-1] - 
        upstream[r].coef_level)) * 0.5 * n[r, i] * sqrt(2 * 9.81 * 
        (level[r, i-1] - upstream[r].coef_level))
        
        # water lever of reservoir
        [r = 1:N, i = 2:numtime],
        level[r, i] == level[r, i-1] + (upstream[r].inflow * raintime[i] * 3600 
        - flow[r, i] * 43200) / upstream[r].k
    end
)


# objective function 目标函数
@objective(model, Min, sum(flow[r, i] / upstream[r].d for r in 1:N for i in 1:numtime))

# print(model)
optimize!(model)

solution_summary(model, verbose=true)
status = termination_status(model)

if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT) && has_values(model)
    if (status == MOI.OPTIMAL)
        println("** Problem solved correctly **")
    else
        println("** Problem returned a (possibly suboptimal) solution **")
    end
    println("- Objective value : ", objective_value(model))
    println("Discharge flow constraint: ---> ", value.(flow))
    println("Water level constraint: ---> ", value.(level))
    println("- Optimal solution is:", value.(n))
    println("- Water level is:", value.(level))
    println("- Discharge flow is:", value.(flow))
else
    println("The model was not solved correctly.")
    println(status)
end
9

Your initial starting point isn’t feasible, so Bonmin encounters an error:

julia> point = Dict(xi => start_value(xi) for xi in all_variables(model))
Dict{VariableRef, Float64} with 3 entries:
  flow[1,1]  => 0.0
  n[1,1]     => 0.0
  level[1,1] => 273.0

julia> primal_feasibility_report(model, point)
Dict{Any, Float64} with 1 entry:
  level[1,1] - (273.0 + (1000.0 * 13.0 * 3600.0 - flow[1,1] * 43200.0) / 2.090104051e7) = 0 => 2.23912

The other problem seems to be the upper bound on the level. You should add a spill variable that allows the reservoir to stay within bounds. Since spilling has a cost, you can penalize it in the objective.

10

Thanks for giving me improving direction for my code!

The upper bound on the level is not satisfied with the initial point, so you recommend to use a spill variable? How can I define the spill variable in my case? Could you please give me an example?

11

Simplifying, you have something like;

level[2] == level[1] + rainfall - flow

But what happens if rainfall is greater than flow? The level increases. What if the level must go above the maximum? The problem is infeasible.

The standard way to fix this is to add a variable which represents spill over the top of the dam:

spill >= 0
level[2] == level[1] + rainfall - flow - spill

Now your problem can never be infeasible, because if it gets too low you will have 0 flow, and if gets too high, you can just spill water.

The other way to model this is to say

level[2] <= level[1] + rainfall - flow

that is, you can’t magically gain water, but you can have less if necessary. In general, more water is good, so this constraint will be tight, except when you need to spill.

12

According to your suggestion, I find a feasible starting point.

julia> # Test the feasibility of starting point
       point = Dict(xi => start_value(xi) for xi in all_variables(model))
Dict{VariableRef, Float64} with 3 entries:
  flow[1,1]  => 36.4898
  n[1,1]     => 7.0
  level[1,1] => 275.164

julia> primal_feasibility_report(model, point)
Dict{Any, Float64}()

With this starting point, the solver still fails to get solution. For now, I have not add a spill variable. But this should not be bother since rainfall is beyond the discharge flow capacity. Start from a feasible starting point, we can always get an optimal solution (local or global) or not?

The model was not solved correctly.
OTHER_ERROR
bonmin: BonHeuristicDiveMIP.cpp:133: virtual int Bonmin::HeuristicDiveMIP::solution(double&, double*): Assertion `isNlpFeasible(minlp, primalTolerance)' failed.
13

I think this is a bug in Bonmin: HeuristicDiveMIP::solution asserts with honor_original_bounds=no from ipopt 3.14.x · Issue #24 · coin-or/Bonmin · GitHub

Adding the spill variable might fix it though.

Otherwise try Juniper.jl as the solver.

14

After adding the spill variable, Bonmin still fails to get solution. So I try Juniper.jl. When I insist on mix-interger problem, Juniper solver gives me an error message.

nl_solver         : MathOptInterface.OptimizerWithAttributes(Ipopt.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("print_level") => 0])
feasibility_pump  : false
log_levels        : [:Options, :Table, :Info]

#Variables: 3
#IntBinVar: 1
Obj Sense: Min

Incumbent using start values: 0.0006401712154631987
Status of relaxation: LOCALLY_SOLVED
Time for relaxation: 0.010133981704711914
Relaxation Obj: 0.0005867002723452641

ONodes   CLevel          Incumbent                   BestBound            Gap    Time   Restarts  GainGap  
============================================================================================================
MethodError: no method matching get_branching_disc_idx!(::Juniper.JuniperProblem, ::Juniper.StepObj, ::Juniper.SolverOptions, ::Vector{Int64}, ::Juniper.GainObj, ::Int32)
Closest candidates are:
 get_branching_disc_idx!(::Any, ::Any, ::Any, ::Any, ::Any) at ~/.julia/packages/Juniper/0Z1vO/src/BnBTree.jl:14
 get_branching_disc_idx!(::Any, ::Any, ::Any, ::Any, ::Any, ::Int64) at ~/.julia/packages/Juniper/0Z1vO/src/BnBTree.jl:14

Stacktrace:
 [1] one_branch_step!(m1::Juniper.JuniperProblem, incumbent::Juniper.Incumbent, opts::Juniper.SolverOptions, step_obj::Juniper.StepObj, disc2var_idx::Vector{Int64}, gains::Juniper.GainObj, counter::Int32)
   @ Juniper ~/.julia/packages/Juniper/0Z1vO/src/BnBTree.jl:367
 [2] solve_sequential(tree::Juniper.BnBTreeObj, last_table_arr::Vector{Any}, time_bnb_solve_start::Float64, fields::Vector{Symbol}, field_chars::Vector{Int32}, time_obj::Juniper.TimeObj)
   @ Juniper ~/.julia/packages/Juniper/0Z1vO/src/BnBTree.jl:468
 [3] solvemip(tree::Juniper.BnBTreeObj)
   @ Juniper ~/.julia/packages/Juniper/0Z1vO/src/BnBTree.jl:743
 [4] optimize!(model::Juniper.Optimizer)
   @ Juniper ~/.julia/packages/Juniper/0Z1vO/src/MOI_wrapper/MOI_wrapper.jl:358
 [5] optimize!
   @ ~/.julia/packages/MathOptInterface/RuRWI/src/Bridges/bridge_optimizer.jl:354 [inlined]
 [6] optimize!
   @ ~/.julia/packages/MathOptInterface/RuRWI/src/MathOptInterface.jl:87 [inlined]
 [7] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Juniper.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
   @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/RuRWI/src/Utilities/cachingoptimizer.jl:316
 [8] optimize!(model::Model; ignore_optimize_hook::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ JuMP ~/.julia/packages/JuMP/Psd1J/src/optimizer_interface.jl:161
 [9] optimize!(model::Model)
   @ JuMP ~/.julia/packages/JuMP/Psd1J/src/optimizer_interface.jl:143
[10] top-level scope
   @ In[186]:189
[11] eval

If I change to a continuous problem, I can get result. The operation above is done by Linux machine.

However, with the Mix-intergered code running on Mac os system, the problem can be solved successfully! Maybe this is a bug in Juniper.jl

15

What is versioninfo() on the linux machine? Is your linux machine 32bit? How did you install Julia?

16

This is the versioninfo():

julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (i686-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-9700 CPU @ 3.00GHz
  WORD_SIZE: 32
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

My linux machine is 64bit.

(base) chenru@chenru-H310M-S2:~/Case_data$ getconf LONG_BIT
64

My Julia is installed by downloading the package julia-1.7.2-linux-i686.tar.gz and decompressing the file

17

WORD_SIZE: 32

You installed the 32-bit version of Julia.

This is a bug in Juniper, but it only happens on the 32-bit version. If you install the 64-bit version things will be fine. (I also wonder if maybe that’s why Bonmin was having trouble.)

Edit: Juniper issue: MethodError on 32-bit linux · Issue #245 · lanl-ansi/Juniper.jl · GitHub

18

Wow, after installing the 64-bit version of Julia, my code runs well. Thanks again for your kind help!