Hi all
I am trying to construct conditioinal statements using the @NLexpression macro but it seems expressions and even variables can’t be used in this context. Is there a way to do this? My aim is to do the following:
if NLexpression satisfies condition
NLconstraint1
else
NLconstraint2
end
Attaching a MWP to show the issue.
model = Model()
@variable(model,x)
@expression(model,exp,2*x)
if x > 1
println("yes")
end
if exp > 1
println("yes")
end
Not a JuMP expert but x is a variable reference to a variable that doesn’t have a value because the model has not been solved yet. After you optimize it the decision variables will have values you can query with value(x).
If you need to use if-else conditions within your mathematical model, you can use a big-M type of solution (see e.g. optimization - if-else condition for the objective variable using big M notation - Operations Research Stack Exchange) or maybe look at something like hdavid16/DisjunctiveProgramming.jl: A JuMP extension for Generalized Disjunctive Programming (github.com)
Hi @AshwinSandeep1, welcome to the forum.
As @slwu89 says, you can’t write exactly what you’ve tried to do, because when you are building the model we don’t know what the value of x is.
One option is to use ifelse.
using JuMP, Ipopt
model = Model(Ipopt.Optimizer)
@variable(model, x)
@expression(model, expr, 2*x)
@NLconstraint(model, ifelse(expr < 0, x, x^2) <= 0)
Do you have a larger reproducible example with the formulation you’re trying to model? People may have suggestions for how to write it.
Thank you @slwu89 and @odow for your quick responses. What @odow suggested is exactly what I wanted! I think using the DisjunctiveProgramming package may have been possible but would have required more work on my side. Just for your info, this is the code, I tried to implement originally:
## Define SOC
@NLexpression(EP, SOC_interior[t in INTERIOR_SUBPERIODS], EP[:vS][end,t-1]/EP[:eTotalCapEnergy][end])
@NLexpression(EP, SOC_start[t in START_SUBPERIODS], EP[:vS][end,t+hours_per_subperiod-1]/EP[:eTotalCapEnergy][end])
## Define the constraint on discharging power
@NLconstraint(EP, TEGSeffpower_interior[t in INTERIOR_SUBPERIODS], ifelse(EP[:SOC_interior][t] >= 0.6717, EP[:vP][end,t] <= EP[:eTotalCap][end],
ifelse(EP[:SOC_interior][t] < 0.6717 && EP[:SOC_interior][t] >= 0.3671,EP[:vP][end,t] <= EP[:eTotalCap][end]*(0.6564*(EP[:SOC_interior][t])+0.5582),
ifelse(EP[:SOC_interior][t] < 0.3671 && EP[:SOC_interior][t] >= 0.1372,EP[:vP][end,t] <= EP[:eTotalCap][end]*(1.7581*(EP[:SOC_interior][t])+0.1546),
ifelse(EP[:SOC_interior][t] < 0.1372 && EP[:SOC_interior][t] >= 0.0195,EP[:vP][end,t] <= EP[:eTotalCap][end]*(3.1561*(EP[:SOC_interior][t])-0.0379),
EP[:vP][end,t] == 0)))))
This gave an error which I didn’t understand. So, I replaced @NLconstraint with @NLexpression and it worked. But when I analysed the results, I realised this was not really constraining the problem. Now, I realise that I made a tiny but significant syntax error where I didn’t put the “<=0” at the end of the constraint. Now, it works as expected!
Thank you again.
This isn’t doing what you think it’s doing.
I think you need instead:
@NLexpression(
EP,
SOC_interior[t in INTERIOR_SUBPERIODS],
vS[end, t-1] / eTotalCapEnergy[end],
)
@NLconstraint(
EP,
TEGSeffpower_interior[t in INTERIOR_SUBPERIODS],
vP[end, t] <= eTotalCap[end] * ifelse(
SOC_interior[t] >= 0.6717,
1,
ifelse(
0.3671 <= SOC_interior[t] < 0.6717,
0.6564 * SOC_interior[t] + 0.5582,
ifelse(
0.1372 <= SOC_interior[t] < 0.3671,
1.7581 * SOC_interior[t] + 0.1546,
ifelse(
0.0195 <= SOC_interior[t] < 0.1372,
3.1561 * SOC_interior[t] - 0.0379,
0.0,
),
),
),
),
)
Note that ifelse can only be used to compute coefficients. It can’t turn on and off constraints.
Oh wow! I would have never thought of multiplying a variable by an ifelse statement. Thank you so much for this. Would you be able to tell me what my code is doing (I tested your version and mine and the results are very different) or direct me to some documentation about using ifelse. I couldn’t find much online.
ifelse(a, x, y) is the expression “if a == true then x else y”.
From the REPL
help?> ifelse
search: ifelse
ifelse(condition::Bool, x, y)
Return x if condition is true, otherwise return y. This differs from ? or if in that it is an ordinary function,
so all the arguments are evaluated first. In some cases, using ifelse instead of an if statement can eliminate the
branch in generated code and provide higher performance in tight loops.
Examples
≡≡≡≡≡≡≡≡≡≡
julia> ifelse(1 > 2, 1, 2)
2
Your code is:
@NLconstraint(
EP,
TEGSeffpower_interior[t in INTERIOR_SUBPERIODS],
ifelse(
EP[:SOC_interior][t] >= 0.6717,
EP[:vP][end,t] <= EP[:eTotalCap][end],
... other stuff ...
) <= 0
)
which says:
(if EP[:SOC_interior][t] >= 0.6717
EP[:vP][end,t] <= EP[:eTotalCap][end]
else
... other stuff
end) <= 0
or in other words, if EP[:SOC_interior][t] >= 0.6717, then
EP[:vP][end,t] <= EP[:eTotalCap][end] must be false. This is the opposite of
what you want.
Hi @odow. Thank you again. I didn’t quite understand why EP[:vP][end,t] <= EP[:eTotalCap][end] must be false in my code. But, I will take your word for it since when I simulated my version, the optimization stops after some iterations saying the problem is infeasible. But, with your version, it keeps running. This actually leads me to 2 questions. I hope you can help.
Cannot call restoration phase at a point that is almost feasible for the restoration NLP. Abort in line search due to no other fall back. Step computation in the restoration phase failed.. And there is another version of the problem which simply keeps running! I am using the Ipopt solver (with ma57 linear solver) with other settings left at default. The primal and dual infeasibility errors drop below their tolerances in some iterations individually but not together (which i believe is the requirement for convergence). Do you have any suggestions on what I could do to help the problem converge? I even tried settings all tolerance values very high (e.g. from the default 1e-08 to 1) and it still did not converge! @expression(EP, SOC_interior[t in INTERIOR_SUBPERIODS], EP[:vS][end,t-1]/1000)
@constraint(
EP,
TEGSeffpower_interior[t in INTERIOR_SUBPERIODS],
EP[:vP][end, t] <= 1000 * ifelse(
EP[:SOC_interior][t] >= 0.6717,
1,
ifelse(
EP[:SOC_interior][t] < 0.6717 && EP[:SOC_interior][t] >= 0.3671,
0.6564 * EP[:SOC_interior][t] + 0.5582,
ifelse(
EP[:SOC_interior][t] < 0.3671 && EP[:SOC_interior][t] >= 0.1372,
1.7581 * EP[:SOC_interior][t] + 0.1546,
ifelse(
EP[:SOC_interior][t] < 0.1372 && EP[:SOC_interior][t] >= 0.0195,
3.1561 * EP[:SOC_interior][t] - 0.0379,
0.0,
),
),
),
),
)
As you can see, SOC is now linear and so I would expect the constraint to be linear as well in principle. However, running this gives me an error saying this syntax is not valid. So, I am guessing the use of ifelse statements in this way is only valid for NL constraints. Do you agree? Does this mean I have no option of reformulating the NLP as an LP?
Thank you
Ashwin
Your constraint is @NLconstraint(EP, ifelse(condition, a <= b, ...other stuff...) <= 0)
If condition is true, then the constraint is equivalent to (a <= b) <= 0, and since (a <= b) evaluates to either 0 (if false) or 1 (if true), then the constraint says that a <= b must be false.
The bigger problem is that Ipopt assumes that the functions are smooth and twice differentiable. It will accept non-differentiable functions like ifelse, but they can cause problems.
Do you always want to maximize vP? If so, instead of constructing the piecewise linear function directly, consider using a formulation like this:
for t in INTERIOR_SUBPERIODS
@NLconstraints(model, begin
vP[end, t] <= eTotalCap[end]
vP[end, t] <= eTotalCap[end] * (0.6564 * SOC_interior[t] + 0.5582)
vP[end, t] <= eTotalCap[end] * (1.7581 * SOC_interior[t] + 0.1546)
vP[end, t] <= eTotalCap[end] * (3.1561 * SOC_interior[t] - 0.0379)
end)
end
You might need to tweak it at the bottom end for the < 0.0195 case, but hopefully you get the idea.
I am guessing the use of
ifelsestatements in this way is only valid for NL constraints
Correct. The epigraph formulation (above) could be written as an LP with @constraint though.
I think I understand what you mean. Attaching a figure just to make sure we are on the same page. But, I have also shown a scenario where this approach would not work (unfortunately, this is a scenario of interest):
There is another scenario wherein I want to replace the PL fits with NL fits where this approach wouldn’t be applicable either due to the curvature of the NL functions. I do have another strategy for modeling which I will describe in another response for clarity.
Just want to thank you again for taking the time to offer great suggestions! This discussion is proving very useful.
Yes, that’s what I mean, and yes, this only works if the curvature is the “right” way. If it’s the “wrong” way, then that means your problem is non-convex, and so all bets are off.
I do have another strategy for modeling which I will describe in another response for clarity.
Sure. It’d help to give a bit more detail on the model. There are quite a few “it depends” answers for what you should do.
One option would be to try fitting some smooth curve instead of piecewise linear.
To give more context about the model: I am trying to model the discharge and charge behaviour of an energy storage system like a battery. The variable vP which you have seen in the code is the discharge power of the system (there is an equivalent charging power vCHARGE as well). In earlier work, these variables were modeled as being either constant with SOC (state of charge) or linearly decreasing/increasing. These were straightforward linear models which solved instantly with Gurobi. However, since this was a poor representation of the atual behaviour, the next step was to use piecewise-linear and piecewise-nonlinear fits to model the behaviour. This is the background of the discussion.
As you said, there is the option of using a single average linear fit (y=mx+c) and a non-linear fit instead of going the piecewise way. One question is whether the non-convex scenario would cause numerical issues even if I use this single fit without ifelse statements?
The other strategy I have is to use linear interpolation of the raw data to compute the power at an SOC using the code below:
## Import and condition discharge power data
discharge_df = DataFrame(CSV.File(joinpath(dirname(dirname(dirname(@__FILE__))),"data","fitting_data","discharge_tau10_ff1.csv")))
push!(discharge_df,(0.0,0.0))
discharge_df = reverse(discharge_df)
push!(discharge_df,(1.0,1.0))
## Define SOC
@NLexpression(EP, SOC_interior[t in INTERIOR_SUBPERIODS], EP[:vS][end,t-1]/(EP[:eTotalCapEnergy][end]))
## Define the discharge interpolation function
function interp_func_dis(SOC)
interp_linear_dis = linear_interpolation(discharge_df[:,1], discharge_df[:,2])
return interp_linear_dis(SOC)
end
register(EP, :interp_func_dis, 1, interp_func_dis; autodiff = true)
# Define the constraint on discharging power
@NLconstraint(EP, TEGSeffpower_interior[t in INTERIOR_SUBPERIODS], EP[:vP][end,t] <= EP[:eTotalCap][end] * interp_func_dis(EP[:SOC_interior][t]))
This again keeps running without convergence which I now understand based on your arguments is because the interp function is not twice differentiable. Also, when I try linearising SOC and repalce the NLconstraint with a linear constraint, it gives a method error which I guess is because this syntax is again probably only valid for NLconstraints. So here, my question would be whether you see another way of implementing the linear interpolation method such that it converges?
If the rest of your model is linear (or quadratic) and you can use Gurobi, then use a special ordered set of type 2: Tips and tricks · JuMP
Something like this:
using JuMP Gurobi
x = [0.0000, 0.0195, 0.1372, 0.3671, 0.6717, 1.0000]
y = [0.0000, 0.0000, 0.3951, 0.7970, 1.0000, 1.0000]
model = Model(Gurobi.Optimizer)
for t in INTERIOR_SUBPERIODS
λ = @variable(model, [1:6], lower_bound = 0, upper_bound = 1)
@constraint(model, sum(λ) == 1)
@constraint(model, λ in SOS2())
@constraint(model, SOC_interior[t] == λ' * x)
@constraint(model, vP[end, t] <= 1000 * λ' * y)
end
I think I understand how this works. So the exact code would be:
## Define SOC
@expression(EP, SOC_interior[t in INTERIOR_SUBPERIODS], EP[:vS][end,t-1])
x = [0.0000, 0.0195, 0.1372, 0.3671, 0.6717, 1.0000]
y = [0.0000, 0.0238, 0.3958, 0.7991, 1.0000, 1.0000]
for t in INTERIOR_SUBPERIODS
λ = @variable(EP, [1:6], lower_bound = 0, upper_bound = 1)
@constraint(EP, sum(λ) == 1)
@constraint(EP, λ in SOS2())
@constraint(EP, EP[:SOC_interior][t] == λ' * x)
@constraint(EP, EP[:vP][end, t] <= 1000 * λ' * y)
end
where the x and y vectors indicate the transition points between each line of the PL model. Based on the value of SOC at a particular time step, the algorithm decides which 2 consective lambda values to be non-zero and thus which line to use. However, the caveats of this method are:
Do you agree?
Yes that looks about right. And all your understanding is correct.
This is very interesting. It is a shame that it only works for linear expressions and constraints with Gurobi. Do you know from experience whether it would be possible to manually implement the algorithm behind SOS2 such that I can use NLexpression and NLconstraint with Ipopt? My first thought would be no since some ifelse statements would be required in the background to decide which lambdas should be non-zero.
Do you know from experience whether it would be possible to manually implement the algorithm behind SOS2 such that I can use NLexpression and NLconstraint with Ipopt?
In general, no. Implementing SOS constraints requires integer variables. (Gurobi implements SOS constraints by branching during the solution process.)
The code with the ifelse above is exactly what you’d write, but Ipopt struggles with this because of the non-differentiability.
Excluding this function, what is the final formulation for the model? Is it nonlinear, quadratic, continuous/discrete, non-convex?
Modeling usually requires some approximations. You’ll probably find that if you add enough points to the piecewise linear function you’ll get a reasonable solution.
Excluding this work ,the model is linear.
Using SOS2 seems a reasonable way to proceed. It’s just that ideally I would have liked to keep SOC an @NLexpression. Simplifying it to an @expression is a compromise but if there is no other option, it is good enough.
Why? What is the true function you are trying to model?
You might be interested to read this paper
Non-linear charge-based battery storage optimization model with bi-variate cubic spline constraints
It uses a cubic splines implementation where “the resulting function is differentiable in Julia/JuMP.”