# How to write conditional expression using JuMP

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 * 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!

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

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 * 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`

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

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

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 * 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

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

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 = ` , there is a feasible solution, n = , which means there is no need to discharge flow for reservoir. But with `raintime = `, the solver can not get feasible solution. I have checked with Excel, when `raintime = `, one possible solution is n = .

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 = 

# 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.initial_level * ones(N, numtime)
n0 = 

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 * 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
``````

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.

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?

Simplifying, you have something like;

``````level == level + 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 == level + 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 <= level + 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.

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

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.

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:
 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
 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
 solvemip(tree::Juniper.BnBTreeObj)
@ Juniper ~/.julia/packages/Juniper/0Z1vO/src/BnBTree.jl:743
 optimize!(model::Juniper.Optimizer)
@ Juniper ~/.julia/packages/Juniper/0Z1vO/src/MOI_wrapper/MOI_wrapper.jl:358
 optimize!
@ ~/.julia/packages/MathOptInterface/RuRWI/src/Bridges/bridge_optimizer.jl:354 [inlined]
 optimize!
@ ~/.julia/packages/MathOptInterface/RuRWI/src/MathOptInterface.jl:87 [inlined]
 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
 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
 optimize!(model::Model)
@ JuMP ~/.julia/packages/JuMP/Psd1J/src/optimizer_interface.jl:143
 top-level scope
@ In:189
 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`

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

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

`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

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