From 2eac9b9ba6e0694ed6fc54da0120cba9d208382c Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 7 Nov 2023 11:10:47 +1300 Subject: [PATCH 1/5] [docs] add tutorial on piecewise linear --- docs/make.jl | 1 + docs/src/tutorials/linear/piecewise_linear.jl | 183 ++++++++++++++++++ 2 files changed, 184 insertions(+) create mode 100644 docs/src/tutorials/linear/piecewise_linear.jl diff --git a/docs/make.jl b/docs/make.jl index 67f9a0b4136..67fdfd181cb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -291,6 +291,7 @@ const _PAGES = [ "tutorials/linear/multi.md", "tutorials/linear/multi_commodity_network.md", "tutorials/linear/tips_and_tricks.md", + "tutorials/linear/piecewise_linear.md", "tutorials/linear/facility_location.md", "tutorials/linear/finance.md", "tutorials/linear/geographic_clustering.md", diff --git a/docs/src/tutorials/linear/piecewise_linear.jl b/docs/src/tutorials/linear/piecewise_linear.jl new file mode 100644 index 00000000000..e5a1b713e37 --- /dev/null +++ b/docs/src/tutorials/linear/piecewise_linear.jl @@ -0,0 +1,183 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors #src +# This Source Code Form is subject to the terms of the Mozilla Public License #src +# v.2.0. If a copy of the MPL was not distributed with this file, You can #src +# obtain one at https://mozilla.org/MPL/2.0/. #src + +# # Piecewise linear functions + +# The purpose of this tutorial is to explain how to represent piecewise linear +# functions in a JuMP model. + +# This tutorial uses the following packages: + +using JuMP +import Plots + +# ## Minimizing a convex function (outer approximation) + +# If the function you are approximating is convex, and you want to minimize +# "down" onto it, then you can use an outer approximation. + +# For example, $f(x) = x^2$ is a convex function: + +f(x) = x^2 +∇f(x) = 2 * x +plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false) + +# Because $f$ is convex, we know that for any $x_k$, we have: +# $$f(x) \ge f(x_k) + \nabla f(x_k) \cdot (x - x_k)$$ + +for x_k in -2:0.5:2 + g = x -> f(x_k) + ∇f(x_k) * (x - x_k) + Plots.plot!(plot, g, X; color = :gray, label = false) +end +plot + +# We can use these _tangent planes_ as constraints in our model to create an +# outer approximation of the function. As we add more planes, the error between +# the true function and the piecewise linear outer approximation decreases. + +# Here is the model in JuMP: + +model = Model() +@variable(model, -2 <= x <= 2) +@variable(model, y) +@constraint(model, [x_k in -2:0.5:2], y >= f(x_k) + ∇f(x_k) * (x - x_k)) +@objective(model, Min, y) + +# !!! note +# This formulation does not work if we want to maximize `y`. + +# ## Maximizing a concave function (outer approximation) + +# The outer approximation also works if we want to maximize "up" into a concave +# function. + +f(x) = log(x) +∇f(x) = 1 / x +X = 0.1:0.1:2 +plot = Plots.plot(f, X; ylims = (-3, 1), label = false) +for x_k in X + g = x -> f(x_k) + ∇f(x_k) * (x - x_k) + Plots.plot!(plot, g, X; color = :gray, label = false) +end +plot + +# Here is the model in JuMP: + +model = Model() +@variable(model, 0.1 <= x <= 2) +@variable(model, y) +@constraint(model, [x_k in X], y <= f(x_k) + ∇f(x_k) * (x - x_k)) +@objective(model, Max, y) + +# !!! note +# This formulation does not work if we want to minimize `y`. + +# ## Minimizing a convex function (inner approximation) + +# Instead of creating an outer approximation, we can also create an inner +# approximation. For example, given $f(x) = x^2$, we may want to approximate the +# true function by the red piecewise linear function: + +f(x) = x^2 +x̂ = -2:0.8:2 +plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false, linewidth = 3) +Plots.plot!(plot, f, x̂; label = false, color = :red, linewidth = 3) +plot + +# To do so, we represent the decision variables $(x, y)$ by the convex +# combination of a set of discrete points $\{x_k, y_k\}_{k=1}^K$: +# ```math +# \begin{aligned} +# x = \sum\limits_{k=1}^K \lambda_k x_k \\ +# y = \sum\limits_{k=1}^K \lambda_k y_k \\ +# \sum\limits_{k=1}^K \lambda_k = 1 \\ +# \lambda_k \ge 0, k=1,\ldots,k \\ +# \end{aligned} +# ``` + +# The feasible region of the convex combination actually allows any $(x, y)$ +# point inside this shaded region: + +I = [1, 2, 3, 4, 5, 6, 1] +Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false) +plot + +# Thus, this formulation does not work if we want to maximize $y$. + +# Here is the model in JuMP: + +x̂ = -2:0.8:2 +ŷ = f.(x̂) +n = length(x̂) +model = Model() +@variable(model, -2 <= x <= 2) +@variable(model, y) +@variable(model, 0 <= λ[1:n] <= 1) +@constraint(model, x == sum(λ[i] * x̂[i] for i in 1:n)) +@constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) +@constraint(model, sum(λ) == 1) +@objective(model, Min, y) + +# ## Maximizing a convex function (inner approximation) + +# The inner approximation also works if we want to maximize "up" into a concave +# function. + +f(x) = log(x) +x̂ = 0.1:0.5:1.6 +plot = Plots.plot(f, 0.1:0.01:1.6; label = false, linewidth = 3) +Plots.plot!(x̂, f.(x̂), linewidth = 3, color = :red, label = false) +I = [1, 2, 3, 4, 1] +Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false) +plot + +# Here is the model in JuMP: + +ŷ = f.(x̂) +n = length(x̂) +model = Model() +@variable(model, 0.1 <= x <= 1.6) +@variable(model, y) +@variable(model, 0 <= λ[1:n] <= 1) +@constraint(model, sum(λ) == 1) +@constraint(model, x == sum(λ[i] * x̂[i] for i in 1:n)) +@constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) +@objective(model, Max, y) + +# ## Non-convex functions + +# If the model is non-convex (or non-concave), then we cannot use an outer +# approximation, and the inner approximation allows a solution far from the true +# function. For example, for $f(x) = sin(x)$, we have: + +f(x) = sin(x) +plot = Plots.plot(f, 0:0.01:2π; label = false) +x̂ = range(; start = 0, stop = 2π, length = 7) +Plots.plot!(x̂, f.(x̂), linewidth = 3, color = :red, label = false) +I = [1, 5, 6, 7, 3, 2, 1] +Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false) +plot + +# We can force the inner approximation to stay on the red line by adding the +# constraint `λ in SOS2()`. The [`SOS2`](@ref) set is a Special Ordered Set of +# Type 2, and it ensures that at most two elements of `λ` can be non-zero, and +# if they are, that they must be adjacent. This prevents the model from taking +# a convex combination of points 1 and 5 to end up on the lower boundary of the +# shaded red area. + +# Here is the model in JuMP: + +ŷ = f.(x̂) +n = length(x̂) +model = Model(); +@variable(model, 0 <= x <= 2π) +@variable(model, y) +@variable(model, 0 <= λ[1:n] <= 1) +@constraints(model, begin + x == sum(λ[i] * x̂[i] for i in 1:n) + y == sum(λ[i] * ŷ[i] for i in 1:n) + sum(λ) == 1 + λ in SOS2() +end) From a71e62f9fbb8ad90f608ff050ebfee190d29d086 Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 7 Nov 2023 11:54:50 +1300 Subject: [PATCH 2/5] Fix --- docs/src/tutorials/linear/piecewise_linear.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/linear/piecewise_linear.jl b/docs/src/tutorials/linear/piecewise_linear.jl index e5a1b713e37..4d0366032b7 100644 --- a/docs/src/tutorials/linear/piecewise_linear.jl +++ b/docs/src/tutorials/linear/piecewise_linear.jl @@ -29,7 +29,7 @@ plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false) for x_k in -2:0.5:2 g = x -> f(x_k) + ∇f(x_k) * (x - x_k) - Plots.plot!(plot, g, X; color = :gray, label = false) + Plots.plot!(plot, g, -2:0.5:2; color = :gray, label = false) end plot @@ -128,7 +128,7 @@ model = Model() f(x) = log(x) x̂ = 0.1:0.5:1.6 plot = Plots.plot(f, 0.1:0.01:1.6; label = false, linewidth = 3) -Plots.plot!(x̂, f.(x̂), linewidth = 3, color = :red, label = false) +Plots.plot!(x̂, f.(x̂); linewidth = 3, color = :red, label = false) I = [1, 2, 3, 4, 1] Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false) plot @@ -155,7 +155,7 @@ model = Model() f(x) = sin(x) plot = Plots.plot(f, 0:0.01:2π; label = false) x̂ = range(; start = 0, stop = 2π, length = 7) -Plots.plot!(x̂, f.(x̂), linewidth = 3, color = :red, label = false) +Plots.plot!(x̂, f.(x̂); linewidth = 3, color = :red, label = false) I = [1, 5, 6, 7, 3, 2, 1] Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false) plot From e4a6a278c1ba7dd96c57347480df47eedb39f5f8 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 7 Nov 2023 15:08:17 +1300 Subject: [PATCH 3/5] Update piecewise_linear.jl --- docs/src/tutorials/linear/piecewise_linear.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/docs/src/tutorials/linear/piecewise_linear.jl b/docs/src/tutorials/linear/piecewise_linear.jl index 4d0366032b7..758e403b6ef 100644 --- a/docs/src/tutorials/linear/piecewise_linear.jl +++ b/docs/src/tutorials/linear/piecewise_linear.jl @@ -179,5 +179,11 @@ model = Model(); x == sum(λ[i] * x̂[i] for i in 1:n) y == sum(λ[i] * ŷ[i] for i in 1:n) sum(λ) == 1 - λ in SOS2() + λ in SOS2() # <-- this is new end) + +# The feasible region of the linear program is the red line in this plot: + +plot = Plots.plot(f, 0:0.01:2π; label = false) +Plots.plot!(x̂, ŷ; linewidth = 3, color = :red, label = false) +plot From 4c6635fa94325bc6fb08593a94c188f826869f76 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 7 Nov 2023 15:09:56 +1300 Subject: [PATCH 4/5] Update piecewise_linear.jl --- docs/src/tutorials/linear/piecewise_linear.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/tutorials/linear/piecewise_linear.jl b/docs/src/tutorials/linear/piecewise_linear.jl index 758e403b6ef..9acf59fa429 100644 --- a/docs/src/tutorials/linear/piecewise_linear.jl +++ b/docs/src/tutorials/linear/piecewise_linear.jl @@ -3,10 +3,10 @@ # v.2.0. If a copy of the MPL was not distributed with this file, You can #src # obtain one at https://mozilla.org/MPL/2.0/. #src -# # Piecewise linear functions +# # Approximating nonlinear functions -# The purpose of this tutorial is to explain how to represent piecewise linear -# functions in a JuMP model. +# The purpose of this tutorial is to explain how to approximate nonlinear functions +# with a mixed-integer linear program. # This tutorial uses the following packages: @@ -146,7 +146,7 @@ model = Model() @constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) @objective(model, Max, y) -# ## Non-convex functions +# ## Piecewise linear approximation # If the model is non-convex (or non-concave), then we cannot use an outer # approximation, and the inner approximation allows a solution far from the true From c43d0140514391a128166608dbdea1272ff6364c Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 7 Nov 2023 15:45:08 +1300 Subject: [PATCH 5/5] Use HiGHS to solve subproblems --- docs/Project.toml | 2 +- docs/src/tutorials/linear/piecewise_linear.jl | 197 +++++++++++++----- 2 files changed, 142 insertions(+), 57 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index df1640fe8c8..66755e1300b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -50,7 +50,7 @@ Ipopt = "=1.5.0" JSON = "0.21" JSONSchema = "1" Literate = "2.8" -MathOptInterface = "=1.20.1" +MathOptInterface = "=1.22.0" MultiObjectiveAlgorithms = "=1.2.0" PATHSolver = "=1.6.1" Plots = "1" diff --git a/docs/src/tutorials/linear/piecewise_linear.jl b/docs/src/tutorials/linear/piecewise_linear.jl index 9acf59fa429..493604437c0 100644 --- a/docs/src/tutorials/linear/piecewise_linear.jl +++ b/docs/src/tutorials/linear/piecewise_linear.jl @@ -11,6 +11,7 @@ # This tutorial uses the following packages: using JuMP +import HiGHS import Plots # ## Minimizing a convex function (outer approximation) @@ -22,14 +23,14 @@ import Plots f(x) = x^2 ∇f(x) = 2 * x -plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false) +plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false, width = 3) # Because $f$ is convex, we know that for any $x_k$, we have: # $$f(x) \ge f(x_k) + \nabla f(x_k) \cdot (x - x_k)$$ -for x_k in -2:0.5:2 +for x_k in -2:1:2 ## Tip: try changing the number of points x_k g = x -> f(x_k) + ∇f(x_k) * (x - x_k) - Plots.plot!(plot, g, -2:0.5:2; color = :gray, label = false) + Plots.plot!(plot, g, -2:0.01:2; color = :red, label = false, width = 3) end plot @@ -39,11 +40,28 @@ plot # Here is the model in JuMP: -model = Model() -@variable(model, -2 <= x <= 2) -@variable(model, y) -@constraint(model, [x_k in -2:0.5:2], y >= f(x_k) + ∇f(x_k) * (x - x_k)) -@objective(model, Min, y) +function outer_approximate_x_squared(x̄) + f(x) = x^2 + ∇f(x) = 2x + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, -2 <= x <= 2) + @variable(model, y) + ## Tip: try changing the number of points x_k + @constraint(model, [x_k in -2:1:2], y >= f(x_k) + ∇f(x_k) * (x - x_k)) + @objective(model, Min, y) + @constraint(model, x == x̄) # <-- a trivial constraint just for testing. + optimize!(model) + return value(y) +end + +# Here are a few values: + +for x̄ in range(; start = -2, stop = 2, length = 15) + ȳ = outer_approximate_x_squared(x̄) + Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black) +end +plot # !!! note # This formulation does not work if we want to maximize `y`. @@ -55,21 +73,45 @@ model = Model() f(x) = log(x) ∇f(x) = 1 / x -X = 0.1:0.1:2 -plot = Plots.plot(f, X; ylims = (-3, 1), label = false) -for x_k in X +X = 0.1:0.1:1.6 +plot = Plots.plot( + f, + X; + xlims = (0.1, 1.6), + ylims = (-3, log(1.6)), + label = false, + width = 3, +) +for x_k in 0.1:0.5:1.6 ## Tip: try changing the number of points x_k g = x -> f(x_k) + ∇f(x_k) * (x - x_k) - Plots.plot!(plot, g, X; color = :gray, label = false) + Plots.plot!(plot, g, X; color = :red, label = false, width = 3) end plot # Here is the model in JuMP: -model = Model() -@variable(model, 0.1 <= x <= 2) -@variable(model, y) -@constraint(model, [x_k in X], y <= f(x_k) + ∇f(x_k) * (x - x_k)) -@objective(model, Max, y) +function outer_approximate_log(x̄) + f(x) = log(x) + ∇f(x) = 1 / x + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0.1 <= x <= 1.6) + @variable(model, y) + ## Tip: try changing the number of points x_k + @constraint(model, [x_k in 0.1:0.5:2], y <= f(x_k) + ∇f(x_k) * (x - x_k)) + @objective(model, Max, y) + @constraint(model, x == x̄) # <-- a trivial constraint just for testing. + optimize!(model) + return value(y) +end + +# Here are a few values: + +for x̄ in range(; start = 0.1, stop = 1.6, length = 15) + ȳ = outer_approximate_log(x̄) + Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black) +end +plot # !!! note # This formulation does not work if we want to minimize `y`. @@ -81,7 +123,7 @@ model = Model() # true function by the red piecewise linear function: f(x) = x^2 -x̂ = -2:0.8:2 +x̂ = -2:0.8:2 ## Tip: try changing the number of points in x̂ plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false, linewidth = 3) Plots.plot!(plot, f, x̂; label = false, color = :red, linewidth = 3) plot @@ -108,17 +150,33 @@ plot # Here is the model in JuMP: -x̂ = -2:0.8:2 -ŷ = f.(x̂) -n = length(x̂) -model = Model() -@variable(model, -2 <= x <= 2) -@variable(model, y) -@variable(model, 0 <= λ[1:n] <= 1) -@constraint(model, x == sum(λ[i] * x̂[i] for i in 1:n)) -@constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) -@constraint(model, sum(λ) == 1) -@objective(model, Min, y) +function inner_approximate_x_squared(x̄) + f(x) = x^2 + ∇f(x) = 2x + x̂ = -2:0.8:2 ## Tip: try changing the number of points in x̂ + ŷ = f.(x̂) + n = length(x̂) + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, -2 <= x <= 2) + @variable(model, y) + @variable(model, 0 <= λ[1:n] <= 1) + @constraint(model, x == sum(λ[i] * x̂[i] for i in 1:n)) + @constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) + @constraint(model, sum(λ) == 1) + @objective(model, Min, y) + @constraint(model, x == x̄) # <-- a trivial constraint just for testing. + optimize!(model) + return value(y) +end + +# Here are a few values: + +for x̄ in range(; start = -2, stop = 2, length = 15) + ȳ = inner_approximate_x_squared(x̄) + Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black) +end +plot # ## Maximizing a convex function (inner approximation) @@ -126,7 +184,7 @@ model = Model() # function. f(x) = log(x) -x̂ = 0.1:0.5:1.6 +x̂ = 0.1:0.5:1.6 ## Tip: try changing the number of points in x̂ plot = Plots.plot(f, 0.1:0.01:1.6; label = false, linewidth = 3) Plots.plot!(x̂, f.(x̂); linewidth = 3, color = :red, label = false) I = [1, 2, 3, 4, 1] @@ -135,16 +193,32 @@ plot # Here is the model in JuMP: -ŷ = f.(x̂) -n = length(x̂) -model = Model() -@variable(model, 0.1 <= x <= 1.6) -@variable(model, y) -@variable(model, 0 <= λ[1:n] <= 1) -@constraint(model, sum(λ) == 1) -@constraint(model, x == sum(λ[i] * x̂[i] for i in 1:n)) -@constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) -@objective(model, Max, y) +function inner_approximate_log(x̄) + f(x) = log(x) + x̂ = 0.1:0.5:1.6 ## Tip: try changing the number of points in x̂ + ŷ = f.(x̂) + n = length(x̂) + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0.1 <= x <= 1.6) + @variable(model, y) + @variable(model, 0 <= λ[1:n] <= 1) + @constraint(model, sum(λ) == 1) + @constraint(model, x == sum(λ[i] * x̂[i] for i in 1:n)) + @constraint(model, y == sum(λ[i] * ŷ[i] for i in 1:n)) + @objective(model, Max, y) + @constraint(model, x == x̄) # <-- a trivial constraint just for testing. + optimize!(model) + return value(y) +end + +# Here are a few values: + +for x̄ in range(; start = 0.1, stop = 1.6, length = 15) + ȳ = inner_approximate_log(x̄) + Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black) +end +plot # ## Piecewise linear approximation @@ -169,21 +243,32 @@ plot # Here is the model in JuMP: -ŷ = f.(x̂) -n = length(x̂) -model = Model(); -@variable(model, 0 <= x <= 2π) -@variable(model, y) -@variable(model, 0 <= λ[1:n] <= 1) -@constraints(model, begin - x == sum(λ[i] * x̂[i] for i in 1:n) - y == sum(λ[i] * ŷ[i] for i in 1:n) - sum(λ) == 1 - λ in SOS2() # <-- this is new -end) - -# The feasible region of the linear program is the red line in this plot: +function piecewise_linear_sin(x̄) + f(x) = sin(x) + ## Tip: try changing the number of points in x̂ + x̂ = range(; start = 0, stop = 2π, length = 7) + ŷ = f.(x̂) + n = length(x̂) + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0 <= x <= 2π) + @variable(model, y) + @variable(model, 0 <= λ[1:n] <= 1) + @constraints(model, begin + x == sum(λ[i] * x̂[i] for i in 1:n) + y == sum(λ[i] * ŷ[i] for i in 1:n) + sum(λ) == 1 + λ in SOS2() # <-- this is new + end) + @constraint(model, x == x̄) # <-- a trivial constraint just for testing. + optimize!(model) + return value(y) +end + +# Here are a few values: -plot = Plots.plot(f, 0:0.01:2π; label = false) -Plots.plot!(x̂, ŷ; linewidth = 3, color = :red, label = false) +for x̄ in range(; start = 0, stop = 2π, length = 15) + ȳ = piecewise_linear_sin(x̄) + Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black) +end plot