Add OuterApproximation algorithm to LocalImprovementSearch

src/plugins/duality_handlers.jl

@@ -217,7 +217,7 @@ function get_dual_solution(node::Node, lagrange::LagrangianDuality)
         return conic_obj, conic_dual
     end
     L_star, λ_star =
-        LocalImprovementSearch.minimize(lagrange.method, λ_star) do x
+        LocalImprovementSearch.minimize(lagrange.method, λ_star, conic_obj) do x
             L_k = _solve_primal_problem(node.subproblem, x, h_expr, h_k)
             return L_k === nothing ? nothing : (s * L_k, s * h_k)
         end

src/plugins/local_improvement_search.jl

@@ -5,29 +5,22 @@
 
 module LocalImprovementSearch
 
+import JuMP
+
 _norm(x) = sqrt(sum(xi^2 for xi in x))
 
 abstract type AbstractSearchMethod end
 
-function minimize(f::Function, x₀::Vector{Float64})
-    return minimize(f, BFGS(100), x₀)
-end
-
-###
-### BFGS
-###
-struct BFGS <: AbstractSearchMethod
-    evaluation_limit::Int
-end
-
 """
-    minimize(f::Function, x₀::Vector{Float64})
+    minimize(
+        f::Function,
+        [method::AbstractSearchMethod = BFGS(100)],
+        x₀::Vector{Float64},
+        lower_bound::Float64 = -Inf,
+    )
 
 Minimizes a convex function `f` using first-order information.
 
-The algorithm is a modified version of BFGS, with a specialized back-tracking
-inexact line-search.
-
 Compared to off-the-shelf implementations, it has a number of features tailored
 to this purpose:
 
@@ -44,8 +37,32 @@ to this purpose:
    * `(f, Δf)::Tuple{Float64,Vector{Float64}`:  a tuple of the function
      evaluation and first-order gradient information.
  * `x₀::Vector{Float64}`: a feasible starting point.
+
+## Default method
+
+The default algorithm is a modified version of BFGS, with a specialized
+back-tracking inexact line-search.
 """
-function minimize(f::F, bfgs::BFGS, x₀::Vector{Float64}) where {F<:Function}
+function minnimize end
+
+function minimize(f::Function, x₀::Vector{Float64}, lower_bound::Float64 = -Inf)
+    return minimize(f, BFGS(100), x₀, lower_bound)
+end
+
+###
+### BFGS
+###
+
+struct BFGS <: AbstractSearchMethod
+    evaluation_limit::Int
+end
+
+function minimize(
+    f::F,
+    bfgs::BFGS,
+    x₀::Vector{Float64},
+    lower_bound::Float64 = -Inf,
+) where {F<:Function}
     # Initial estimte for the Hessian matrix in BFGS
     B = zeros(length(x₀), length(x₀))
     for i in 1:size(B, 1)
@@ -129,4 +146,59 @@ function _line_search(
     return 0.0, fₖ, ∇fₖ
 end
 
+###
+### Cutting plane
+###
+
+struct OuterApproximation{O} <: AbstractSearchMethod
+    optimizer::O
+end
+
+function minimize(
+    f::F,
+    method::OuterApproximation,
+    x₀::Vector{Float64},
+    lower_bound::Float64,
+) where {F<:Function}
+    model = JuMP.Model(method.optimizer)
+    JuMP.set_silent(model)
+    n = length(x₀)
+    JuMP.@variable(model, x[i in 1:n], start = x₀[i])
+    JuMP.@variable(model, θ >= lower_bound)
+    JuMP.@objective(model, Min, θ)
+    xₖ = x₀
+    fₖ, ∇fₖ = f(xₖ)::Tuple{Float64,Vector{Float64}}
+    upper_bound = fₖ
+    JuMP.@constraint(model, θ >= fₖ + ∇fₖ' * (x - xₖ))
+    evals = Ref(0)
+    d_step = Inf
+    while d_step > 1e-8 && evals[] < 20
+        JuMP.optimize!(model)
+        lower_bound, xₖ₊₁ = JuMP.value(θ), JuMP.value.(x)
+        ret = f(xₖ₊₁)
+        while ret === nothing
+            # point is infeasible
+            xₖ₊₁ = 0.5 * (xₖ + xₖ₊₁)
+            ret = f(xₖ₊₁)
+        end
+        fₖ₊₁, ∇fₖ₊₁ = ret::Tuple{Float64,Vector{Float64}}
+        evals[] += 1
+        upper_bound = fₖ₊₁
+        JuMP.@constraint(model, θ >= fₖ₊₁ + ∇fₖ₊₁' * (x - xₖ₊₁))
+        d = xₖ₊₁ - xₖ
+        d_step = _norm(d)
+        if sign(∇fₖ' * d) != sign(∇fₖ₊₁' * d)
+            # There is a kink between the x
+            xₖ₊₂ = 0.5 * (xₖ + xₖ₊₁)
+            fₖ₊₂, ∇fₖ₊₂ = f(xₖ₊₂)::Tuple{Float64,Vector{Float64}}
+            evals[] += 1
+            JuMP.@constraint(model, θ >= fₖ₊₂ + ∇fₖ₊₂' * (x - xₖ₊₂))
+            fₖ, xₖ = fₖ₊₂, xₖ₊₂
+        else
+            fₖ, xₖ = fₖ₊₁, xₖ₊₁
+        end
+    end
+    return fₖ, xₖ
+end
+
 end

test/plugins/local_improvement_search.jl

@@ -1,7 +1,13 @@
+#  Copyright (c) 2017-23, Oscar Dowson and SDDP.jl contributors.
+#  This Source Code Form is subject to the terms of the Mozilla Public
+#  License, v. 2.0. If a copy of the MPL was not distributed with this
+#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
 module TestLocalImprovementSearch
 
 using Test
-import SDDP.LocalImprovementSearch
+import SDDP: LocalImprovementSearch
+import HiGHS
 
 function runtests()
     for name in names(@__MODULE__, all = true)
@@ -66,6 +72,61 @@ function test_piecewise()
     return
 end
 
+function test_x_squared_outer_approximation()
+    calls = 0
+    solver = LocalImprovementSearch.OuterApproximation(HiGHS.Optimizer)
+    f, x = LocalImprovementSearch.minimize(solver, [0.0], 0.0) do x
+        f = 2.1 + (x[1] - 1.1)^2
+        f′ = [2 * (x[1] - 1.1)]
+        calls += 1
+        return f, f′
+    end
+    @info "OA(squared) = $(calls)"
+    @test isapprox(f, 2.1, atol = 1e-6)
+    @test isapprox(x, [1.1], atol = 1e-4)
+    return
+end
+
+function test_exp_outer_approximation()
+    calls = 0
+    solver = LocalImprovementSearch.OuterApproximation(HiGHS.Optimizer)
+    f, x = LocalImprovementSearch.minimize(solver, [1.0], 0.0) do x
+        calls += 1
+        if x[1] < 0.1 || x[1] > 20
+            return nothing
+        end
+        return exp(x[1]), [exp(x[1])]
+    end
+    @info "OA(exp) = $(calls)"
+    @test isapprox(f, exp(0.1), atol = 1e-2)
+    @test isapprox(x, [0.1], atol = 1e-2)
+    return
+end
+
+function test_piecewise_outer_approximation()
+    calls = 0
+    solver = LocalImprovementSearch.OuterApproximation(HiGHS.Optimizer)
+    f, x = LocalImprovementSearch.minimize(solver, [0.05], -1.0) do x
+        calls += 1
+        if x[1] < 0.0
+            return nothing
+        elseif 0.0 <= x[1] < 0.1
+            return -0.1 - 1 * (x[1] - 0.0), [-1.0]
+        elseif 0.1 <= x[1] < 0.4
+            return -0.2 - 0.8 * (x[1] - 0.1), [-0.8]
+        elseif 0.4 <= x[1] <= 1.0
+            return -0.44 + 0.1 * (x[1] - 0.4), [0.1]
+        else
+            @assert 1.0 <= x[1]
+            return nothing
+        end
+    end
+    @info "OA(piecewise) = $(calls)"
+    @test isapprox(f, -0.44, atol = 1e-3)
+    @test isapprox(x, [0.4], atol = 1e-3)
+    return
+end
+
 end
 
 TestLocalImprovementSearch.runtests()