From 23a7f002e45f839820869201f014cf946f3b0f15 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 8 Nov 2023 10:06:24 +1300 Subject: [PATCH] [docs] add tutorial on piecewise linear (#3563) --- docs/Project.toml | 2 +- docs/make.jl | 1 + docs/src/tutorials/linear/piecewise_linear.jl | 274 ++++++++++++++++++ 3 files changed, 276 insertions(+), 1 deletion(-) create mode 100644 docs/src/tutorials/linear/piecewise_linear.jl 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/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..493604437c0 --- /dev/null +++ b/docs/src/tutorials/linear/piecewise_linear.jl @@ -0,0 +1,274 @@ +# 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 + +# # Approximating nonlinear functions + +# 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: + +using JuMP +import HiGHS +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, 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: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.01:2; color = :red, label = false, width = 3) +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: + +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`. + +# ## 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: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 = :red, label = false, width = 3) +end +plot + +# Here is the model in JuMP: + +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`. + +# ## 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 ## 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 + +# 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: + +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) + +# 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 ## 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] +Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false) +plot + +# Here is the model in JuMP: + +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 + +# 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: + +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: + +for x̄ in range(; start = 0, stop = 2π, length = 15) + ȳ = piecewise_linear_sin(x̄) + Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black) +end +plot