# Solving a piecewise linear system

**URL:** <https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [December 20, 2018, 4:18pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856 "2018-12-20T16:18:23Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 20, 2018, 4:18pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/1 "2018-12-20T16:18:23Z")

</div>

The help I need here is mostly conceptual, but I suspect that either JuMP.jl or Convex.jl could do this. Any suggestions or pointers to texts/tutorials would be appreciated, I just don’t know how to reformulate this.

I need to solve a system in n+1 unknowns s\_0, s\_1, \dots, s\_n (the first one is handled specially), described as follows.

Let x^+ = \max(0, x) as usual, and apply elementwise for vectors. Let \mathbf{s} denote the vector [s\_1, \dots, s\_n] (note: s\_0 is not in there).

Let a \> 0, b \> 0, c\_0 \> 0, \kappa \> 0 be constants, and \mathbf{c} and \mathbf{f} n-element vectors. Then the system is

a s\_0 = c\_0 + \kappa \langle \mathbf{s}^+ , f\rangle

and for all i = 1, \dots, n,

a s\_0 + b s\_i = c\_i + \langle (\mathbf{s} - s\_i)^+, f \rangle

where \langle \cdot, \cdot \rangle is the dot product.

I read a suggestion that such problems can be cast as LP problems, but I have never encountered this before so I don’t know how to start.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [December 20, 2018, 5:28pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/2 "2018-12-20T17:28:43Z")

</div>

How big is n?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 20, 2018, 6:09pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/3 "2018-12-20T18:09:40Z")

</div>

1 write as min 0 under some constraints

2 rewrite the nonlinear operations as linear constraints, eg max(0,x) becomes y, with constraint y \>= x, y \>= 0

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 20, 2018, 6:20pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/4 "2018-12-20T18:20:57Z")

</div>

At first glance this looks like an MILP. You could use [GitHub - rdeits/ConditionalJuMP.jl: Automatic transformation of implications and complementarity into mixed-integer models in Julia](https://github.com/rdeits/ConditionalJuMP.jl) to set it up. There might be a way to formulate it as an LP, but I’m not seeing it straight away. What @antoine-levitt proposed constitutes a relaxation of the problem, I’d think?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 20, 2018, 6:51pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/5 "2018-12-20T18:51:49Z")

</div>

If the I stands for Integer in MILP, I don’t think it is one, all variables are real. What @antoine-levitt suggested may be the best solution; I was thinking of a Gauss-Seidel scheme but this should be better. n is around 10^3 to 10^5, more is always better but I can live with less if I need to.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 20, 2018, 6:56pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/6 "2018-12-20T18:56:33Z")

</div>

Ah sorry of course what I said doesn’t work because y isn’t constrained to be as small as possible (I was thinking of a case where you minimize a max).

Edit : maybe minimize y then? But there must be a cleaner way

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [December 20, 2018, 7:14pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/7 "2018-12-20T19:14:34Z")

</div>

IIUC, @antoine-levitt 's solution would work if you use the Simplex method. Proof: We have n+1 original variables, and the same number of equality constraints, and we will introduce n^2 new variables and 2n^2 new inequality constraints. Since a basic feasible solution will have at least as many constraints active as the number of variables, all equality constraints will be active naturally and at least n^2 of the inequality constraints will be active. The latter means that the max relationship will hold for each variable.

The main problem here is that if n is large, you will have an LP with many constraints and variables, so it may be pretty hard to solve. 10^10 variables and constraints is a pretty large LP but _maybe_ CPLEX can do it.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 20, 2018, 9:16pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/8 "2018-12-20T21:16:04Z")

</div>

I tried setting up the relaxation proposed by @antoine-levitt and @mohamed82008 in JuMP and solving it using a simplex method, see below. But the relaxation really doesn’t appear to be tight.

```julia
using JuMP
using Random
using LinearAlgebra

# problem data
Random.seed!(1)
n = 10
a = rand()
b = rand()
κ = rand()
c₀ = rand()
c = randn(n)
f = randn(n)

using Gurobi
model = Model(solver=GurobiSolver(Method=0 #=Primal simplex=#))
# using Clp
# model = Model(solver=ClpSolver(SolveType=1 #=Primal simplex=#))
@variable model s₀
@variable model s[1 : n]
@variable model s⁺[1 : n] >= 0
@constraint model s⁺ .>= s
@variable model δ[1 : n, 1 : n] >= 0 # (s - s[i])⁺
for i = 1 : n
    @constraint model δ[:, i] .>= s - s[i]
end
@constraint model a * s₀ == c₀ + κ * s⁺ ⋅ f
for i = 1 : n
    @constraint model a * s₀ + b * s[i] == c[i] + δ[:, i] ⋅ f
end

solve(model)

using Test
@test getvalue.(s⁺) ≈ max.(getvalue.(s), 0) atol=1e-6 # passes
for i = 1 : n
    @test getvalue.(δ[:, i]) ≈ max.(getvalue.(s) .- getvalue(s[i]), 0) atol=1e-6 # fails
end

```

Am I doing something wrong?

> [@Tamas\_Papp](#):
>
> If the I stands for Integer in MILP, I don’t think it is one, all variables are real.

Yes, integer. All variables in the original nonlinear program are real, but I was proposing to reformulate it as an MILP, with integer variables used to indicate which piece of each piecewise function you’re on.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 20, 2018, 9:21pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/9 "2018-12-20T21:21:20Z")

</div>

Yes, my “solution” was very wrong. I think if you optimize over the sum of the extra variables it will work, but that looks like a hack, there should be a better way.

From what I understand (but I could be wrong) LP solvers are usually very well optimized, and it takes some effort to beat them, even if the LP reformulation is slightly suboptimal.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 20, 2018, 9:25pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/10 "2018-12-20T21:25:30Z")

</div>

Adding

```julia
@objective model Min sum(s⁺) + sum(δ)

```

does not make the tests pass though.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 20, 2018, 9:29pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/11 "2018-12-20T21:29:45Z")

</div>

> [@tkoolen](#):
>
> At first glance this looks like an MILP. You could use [https://github.com/rdeits/ConditionalJuMP.jl](https://github.com/rdeits/ConditionalJuMP.jl) to set it up.

Actually, [GitHub - joehuchette/PiecewiseLinearOpt.jl: Solve optimization problems containing piecewise linear functions](https://github.com/joehuchette/PiecewiseLinearOpt.jl) is probably a better fit.

---

<div class="post-metadata">

**Author:** ![Azamat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/azamat/32/6892_2.png) [@Azamat](https://discourse.julialang.org/u/Azamat)\
**Post date:** [December 20, 2018, 9:39pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/12 "2018-12-20T21:39:19Z")

</div>

> [@Tamas\_Papp](#):
>
> as0+bsi=ci+⟨(s−si)+,f⟩

Just to clarify, you are subtracting element from the vector in (\mathbf{s} − s\_i). Did you mean elementwise difference? I.e. `s .- s_i`?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 20, 2018, 9:45pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/13 "2018-12-20T21:45:03Z")

</div>

Ah yes, there’s no particular reason it should. That’ll teach me to actually think before typing. I’ll shut up now 🙂

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 20, 2018, 9:54pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/14 "2018-12-20T21:54:52Z")

</div>

The following works for small `n` at least, if you’re OK with adding artificial bounds on `s` (needed for these reformulations to work):

```julia
using JuMP
using Random
using LinearAlgebra
using PiecewiseLinearOpt

# problem data
Random.seed!(1)
n = 5
a = rand()
b = rand()
κ = rand()
c₀ = rand()
c = randn(n)
f = randn(n)

# using Gurobi
# model = Model(solver=GurobiSolver())
using Cbc
model = Model(solver=CbcSolver())
@variable model s₀
@variable model s[1 : n]
s⁺ = piecewiselinear.(Ref(model), s, Ref([-10, 0, 10]), x -> max(0, x)) # 10 is arbitrary
δ = Matrix{Variable}(undef, n, n)
for i = 1 : n
    δ[:, i] = piecewiselinear.(Ref(model), s - s[i], Ref([-10, 0, 10]), x -> max(0, x)) # 10 is arbitrary
end
@constraint model a * s₀ == c₀ + κ * s⁺ ⋅ f
for i = 1 : n
    @constraint model a * s₀ + b * s[i] == c[i] + δ[:, i] ⋅ f
end

solve(model)

```

Edit: this requires the master branch of PiecewiseLinearOpt on Julia 1.0.

---

<div class="post-metadata">

**Author:** ![Azamat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/azamat/32/6892_2.png) [@Azamat](https://discourse.julialang.org/u/Azamat)\
**Post date:** [December 20, 2018, 10:44pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/15 "2018-12-20T22:44:18Z")

</div>

Here is my reformulation of your problem as a linear program:

 ![58%20PM](https://global.discourse-cdn.com/julialang/original/3X/f/a/fa4d1b851036f77c9c8edf11780807745a5a4993.png)

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 20, 2018, 11:24pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/16 "2018-12-20T23:24:28Z")

</div>

What values should be chosen for \alpha and \beta?

---

<div class="post-metadata">

**Author:** ![Azamat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/azamat/32/6892_2.png) [@Azamat](https://discourse.julialang.org/u/Azamat)\
**Post date:** [December 21, 2018, 12:40am UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/17 "2018-12-21T00:40:18Z")

</div>

Actually, after thinking a bit more about this, I don’t think my recasting is completely equivalent.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 21, 2018, 7:26am UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/18 "2018-12-21T07:26:23Z")

</div>

> [@Azamat](#):
>
> Did you mean elementwise difference? I.e. `s .- s_i` ?

Yes. Thanks for the clarification.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [December 21, 2018, 11:57am UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/20 "2018-12-21T11:57:29Z")

</div>

This is weird. Either my understanding is off, not sure how, or the solver is doing something funny behind the scene.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 21, 2018, 1:32pm UTC](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856/21 "2018-12-21T13:32:43Z")

</div>

I have figured it out: it is trivial to show (eg proof by contradiction) that c\_i \< c\_j \Rightarrow s\_i \< s\_j, so we have the ordering and can just write the linear system because we know which s\_k - s\_j is positive, and solve the whole system given s\_0, so it becomes a rootfinding problem in a single variable.

Which is great, because I need to solve it millions of times 😉

Thanks for all the suggestions, I learned a lot from them.

[Next page](https://discourse.julialang.org/t/solving-a-piecewise-linear-system/18856.md?page=2)