Аппроксимация нелинейных функций

Цель данного руководства — показать, как можно аппроксимировать нелинейные функции с помощью частично целочисленной линейной программы.

В этом руководстве используются следующие пакеты.

using JuMP
import HiGHS
import Plots

Минимизация выпуклой функции (внешняя аппроксимация)

Если аппроксимируемая функция является выпуклой и нужно выполнить минимизацию «вниз» до нее, можно использовать внешнюю аппроксимацию.

Например, — это выпуклая функция:

f(x) = x^2
∇f(x) = 2 * x
plot = Plots.plot(f, -2:0.01:2; ylims = (-0.5, 4), label = false, width = 3)

Так как функция выпуклая, известно, что для любого верно следующее:

for x_k in -2:1:2  ## Совет: попробуйте изменить количество точек 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

Эти касательные плоскости можно использовать в качестве ограничений в модели для построения внешней аппроксимации функции. По мере увеличения количества плоскостей погрешность между истинной функцией и кусочно-линейной внешней аппроксимацией уменьшается.

Вот модель в 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)
    # Совет: попробуйте изменить количество точек 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̄)  # <-- простейшее ограничение только для тестирования.
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return value(y)
end

outer_approximate_x_squared (generic function with 1 method)

Вот несколько значений:

for x̄ in range(; start = -2, stop = 2, length = 15)
    ȳ = outer_approximate_x_squared(x̄)
    Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black)
end
plot

Эта переформулировка не подходит, если нужно максимизировать y.

Максимизация вогнутой функции (внешняя аппроксимация)

Внешняя аппроксимация также подходит, если нужно выполнить максимизацию «вверх» до вогнутой функции.

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  ## Совет: попробуйте изменить количество точек 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

Вот модель в 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)
    # Совет: попробуйте изменить количество точек 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̄)  # <-- простейшее ограничение только для тестирования.
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return value(y)
end

outer_approximate_log (generic function with 1 method)

Вот несколько значений:

for x̄ in range(; start = 0.1, stop = 1.6, length = 15)
    ȳ = outer_approximate_log(x̄)
    Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black)
end
plot

Эта переформулировка не подходит, если нужно минимизировать y.

Минимизация выпуклой функции (внутренняя аппроксимация)

Помимо внешней аппроксимации, можно также построить внутреннюю. Например, истинную функцию можно аппроксимировать с помощью красной кусочно-линейной функции:

f(x) = x^2
x̂ = -2:0.8:2  ## Совет: попробуйте изменить количество точек в 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

Для этого представим переменные решения в виде выпуклой комбинации множества дискретных точек :

Допустимая область выпуклой комбинации фактически допускает любую точку внутри закрашенной области:

I = [1, 2, 3, 4, 5, 6, 1]
Plots.plot!(x̂[I], f.(x̂[I]); fill = (0, 0, "#f004"), width = 0, label = false)
plot

Поэтому эта переформулировка не подходит, если нужно максимизировать .

Вот модель в JuMP:

function inner_approximate_x_squared(x̄)
    f(x) = x^2
    ∇f(x) = 2x
    x̂ = -2:0.8:2  ## Совет: попробуйте изменить количество точек в x̂.
    ŷ = 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] * ŷ[i] for i in 1:n))
    @constraint(model, sum(λ) == 1)
    @objective(model, Min, y)
    @constraint(model, x == x̄)  # <-- простейшее ограничение только для тестирования.
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return value(y)
end

inner_approximate_x_squared (generic function with 1 method)

Вот несколько значений:

for x̄ in range(; start = -2, stop = 2, length = 15)
    ȳ = inner_approximate_x_squared(x̄)
    Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black)
end
plot

Максимизация выпуклой функции (внутренняя аппроксимация)

Внутренняя аппроксимация также подходит, если нужно выполнить максимизацию «вверх» до вогнутой функции.

f(x) = log(x)
x̂ = 0.1:0.5:1.6  ## Совет: попробуйте изменить количество точек в 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

Вот модель в JuMP:

function inner_approximate_log(x̄)
    f(x) = log(x)
    x̂ = 0.1:0.5:1.6  ## Совет: попробуйте изменить количество точек в x̂.
    ŷ = 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] * ŷ[i] for i in 1:n))
    @objective(model, Max, y)
    @constraint(model, x == x̄)  # <-- простейшее ограничение только для тестирования.
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return value(y)
end

inner_approximate_log (generic function with 1 method)

Вот несколько значений:

for x̄ in range(; start = 0.1, stop = 1.6, length = 15)
    ȳ = inner_approximate_log(x̄)
    Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black)
end
plot

Кусочно-линейная аппроксимация

Если модель невыпуклая (или невогнутая), то использовать внешнюю аппроксимацию нельзя, а внутренняя аппроксимация может давать решение, далекое от истинной функции. Например, для мы получаем следующее:

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

Чтобы внутренняя аппроксимация оставалась на красной линии, можно добавить ограничение λ in SOS2(). Множество SOS2 представляет собой особое упорядоченное множество типа 2, которое гарантирует, что не более двух элементов λ могут быть ненулевыми и они должны быть смежными. Это не позволяет модели принимать выпуклую комбинацию точек 1 и 5 и оказаться на нижней границе закрашенной красным области.

Вот модель в JuMP:

function piecewise_linear_sin(x̄)
    f(x) = sin(x)
    # Совет: попробуйте изменить количество точек в x̂.
    x̂ = range(; start = 0, stop = 2π, length = 7)
    ŷ = 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] * ŷ[i] for i in 1:n)
        sum(λ) == 1
        λ in SOS2()  # <-- это новое ограничение
    end)
    @constraint(model, x == x̄)  # <-- простейшее ограничение только для тестирования.
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return value(y)
end

piecewise_linear_sin (generic function with 1 method)

Вот несколько значений:

for x̄ in range(; start = 0, stop = 2π, length = 15)
    ȳ = piecewise_linear_sin(x̄)
    Plots.scatter!(plot, [x̄], [ȳ]; label = false, color = :black)
end
plot