Abs in JuMP - how (why?) does it work

Specific Domains Optimization (Mathematical)
1

I cannot get my head around what JuMP/Ipopt is doing if I have an absolute value in the constraints. This runs and I cannot find out why:

using JuMP, Ipopt
tf=10

model = Model(Ipopt.Optimizer)

@variable(model, 0<=v) 
@variable(model, x[1:tf]) 
@objective(model, Max, v)
@constraint(model, tryErr[k=1:tf], abs(x[k]).>=v) ###This is equivalent to OR constraints: x[k]>=v OR x[k]<=-v, how does JuMP handle it?
@constraint(model, importantCstr[k=1:tf], abs.(10.0./k.*sin(x[k])).<=0.5)###This one is the same as two constraints: -0.5<=10.0./k.*sin(x[k]) AND 10.0./k.*sin(x[k])<=0.5, OK

optimize!(model)

I found this topic introducing binary variables - but if I have binary variables, then I cannot use Ipopt. So how is JuMP handling the tryErr constraint if the code runs?

The actual problem I am trying to solve can be interpreted as finding a maximum v that pushes all elements of x away from zero while keeping the importantCstr is still satisfied. The problem is nonlinear, so using MILP solvers is out of question.

2

Ipopt is treating abs as a nonlinear function. It computes a derivative by looking at whether the child is positive or negative. The non-smoothness often causes convergence issues for Ipopt.

When you have abs(x) <= 0.5, this constraint is convex, and things should mostly work. Although for some problems Ipopt will not converge because it oscillates around x = 0.

When you have abs(x) >= v, this constraint is non-convex. If Ipopt converges, it may be to a local solution. If I had to guess, it will probably be related to the starting point, and shoot off to whichever “side” of abs is active first.

2 Likes
3

how about using smooth approximations of absolute value function? does it change something?

4

That would fix the derivative/oscillating issue, not the non-convexity of smooth_abs(x) >= v.

5

The “Huber loss” is sometimes used along these lines.

2 Likes
6

I am OK with local solutions (my original problem doesn’t have a unique solution for x anyway, it may have a unique v). I was expecting the code to fail because I couldn’t figure out how the non-convex and non-smooth constraint could be reformulated into a smooth one (the non-smooth and convex is easier). I thought JuMP was doing something smart that I was missing :slight_smile: If ipopt simply treats abs(x) as a nonlinear function, the non-smoothness at zero seems to have little influence as I am pushing away from zero.

I figured out a reformulation like this:

@constraint(model, tryErrMul[k=1:tf], (x[k]-v)*(x[k]+v).>=0.0)

but I will read up on Huber loss, it looks interesting.

2 Likes
7

I would be happy to see what happens when we use

\lambda(x) = \frac{xe^{kx} - xe^{-kx}}{e^{kx} + e^{-kx}}

instead of the abs() function. As @odow states, you will catch the oscillating issue. You also prevent staying near to zero by adding extra contraints also.

Edit: A good discussion on smoothing abs: calculus - Approximate $|x|$ with a smooth function - Mathematics Stack Exchange

Edit 2:

However you lose the convexity.

8

Just playing with the reformulation, here is the code:

using JuMP, Ipopt
tf=10
kVal=10

model = Model(Ipopt.Optimizer)
@variable(model, 0<=v) 
@variable(model, x[1:tf]) 
@objective(model, Max, v)
@expression(model, λ[k=1:tf], (x[k]*exp(x[k]*kVal)-x[k]*exp(-x[k]*kVal))./(exp(x[k]*kVal)+exp(-x[k]*kVal)))
@constraint(model, [k=1:tf], abs.(10.0./k.*sin(x[k])).<=0.5)
@constraint(model, tryErrNew[k=1:tf], λ[k].>=v) 
optimize!(model)

Unfortunately, this reformulation requires me to set the value of kVal that directly affects how well the reformulation approximates the absolute value. For my problem, it is important to make sure that for a given v if x^* solves |x|-v=0, then \lambda(x^*)-v=0 as well. And depending on the chosen kVal, the proposed reformulation \lambda(x) underestimates |x| around zero and the solutions to \lambda(x)-v=0 don’t overlap with the solutions to |x|-v=0. Because I don’t know the actual value of v, I cannot guarantee that the kVal will be large enough.

Besides, the proposed reformulation (or my implementation thereof) seems to mess something up numerically when I rewrite also the convex non-smooth constraint:

using JuMP, Ipopt
tf=10
kVal=10

model = Model(Ipopt.Optimizer)
@variable(model, 0<=v) 
@variable(model, x[1:tf]) 
@objective(model, Max, v)
@expression(model, λ[k=1:tf], (x[k]*exp(x[k]*kVal)-x[k]*exp(-x[k]*kVal))./(exp(x[k]*kVal)+exp(-x[k]*kVal)))
@constraint(model, [k=1:tf], 10.0./k.*sin(x[k]).<=0.5)
@constraint(model, [k=1:tf], 10.0./k.*sin(x[k]).>=-0.5)
@constraint(model, tryErrNew[k=1:tf], λ[k].>=v) 
optimize!(model)

gives
image
whereas

using JuMP, Ipopt
tf=10
kVal=10

model = Model(Ipopt.Optimizer)
@variable(model, 0<=v) 
@variable(model, x[1:tf]) 
@objective(model, Max, v)
@constraint(model, tryErr[k=1:tf], abs(x[k]).>=v) 
@constraint(model, [k=1:tf], 10.0./k.*sin(x[k]).<=0.5)
@constraint(model, [k=1:tf], 10.0./k.*sin(x[k]).>=-0.5)
optimize!(model)

gives:
image

1 Like