Implement multistart (#500)

Fixes #359.

.github/workflows/sanitizers.yml

@@ -8,18 +8,17 @@ concurrency:
 
 jobs:
   build:
-    timeout-minutes: 10
+    timeout-minutes: 15
     strategy:
       fail-fast: false
       matrix:
         include:
           - name: asan
             cmake-flags: -DCMAKE_BUILD_TYPE=Asan -DBUILD_PYTHON=OFF
             ctest-env:
-          # tsan takes a long time, and we don't currently use threading
-          #- name: tsan
-          #  cmake-flags: -DCMAKE_BUILD_TYPE=Tsan -DBUILD_PYTHON=OFF
-          #  ctest-env: TSAN_OPTIONS=second_deadlock_stack=1
+          - name: tsan
+            cmake-flags: -DCMAKE_BUILD_TYPE=Tsan -DBUILD_PYTHON=OFF
+            ctest-env: TSAN_OPTIONS=second_deadlock_stack=1
           - name: ubsan
             cmake-flags: -DCMAKE_BUILD_TYPE=Ubsan -DBUILD_PYTHON=OFF
             ctest-env:

CMakeLists.txt

@@ -94,6 +94,11 @@ option(USE_SYSTEM_EIGEN "Use system eigen" OFF)
 option(USE_SYSTEM_FMT "Use system fmt" OFF)
 option(USE_SYSTEM_PYBIND "Use system pybind" OFF)
 
+# Required for std::async()
+find_package(Threads REQUIRED)
+
+target_link_libraries(Sleipnir PUBLIC Threads::Threads)
+
 # Eigen dependency
 if(NOT USE_SYSTEM_EIGEN)
     fetchcontent_declare(
@@ -325,7 +330,10 @@ if(BUILD_PYTHON)
             ${fmt_SOURCE_DIR}/include
             ${CMAKE_CURRENT_SOURCE_DIR}/jormungandr/cpp
     )
-    target_link_libraries(_jormungandr PUBLIC pybind11::module Eigen3::Eigen)
+    target_link_libraries(
+        _jormungandr
+        PUBLIC pybind11::module Threads::Threads Eigen3::Eigen
+    )
 
     install(
         TARGETS _jormungandr

include/sleipnir/optimization/Multistart.hpp

@@ -0,0 +1,69 @@
+// Copyright (c) Sleipnir contributors
+
+#pragma once
+
+#include <algorithm>
+#include <future>
+#include <span>
+#include <vector>
+
+#include "sleipnir/optimization/SolverStatus.hpp"
+#include "sleipnir/util/FunctionRef.hpp"
+
+namespace sleipnir {
+
+/**
+ * The result of a multistart solve.
+ *
+ * @tparam DecisionVariables The type containing the decision variable initial
+ *   guess.
+ */
+template <typename DecisionVariables>
+struct MultistartResult {
+  SolverStatus status;
+  DecisionVariables variables;
+};
+
+/**
+ * Solves an optimization problem from different starting points in parallel,
+ * then returns the solution with the lowest cost.
+ *
+ * Each solve is performed on a separate thread. Solutions from successful
+ * solves are always preferred over solutions from unsuccessful solves, and cost
+ * (lower is better) is the tiebreaker between successful solves.
+ *
+ * @tparam DecisionVariables The type containing the decision variable initial
+ *   guess.
+ * @param solve A user-provided function that takes a decision variable initial
+ *   guess and returns a MultistartResult.
+ * @param initialGuesses A list of decision variable initial guesses to try.
+ */
+template <typename DecisionVariables>
+MultistartResult<DecisionVariables> Multistart(
+    function_ref<MultistartResult<DecisionVariables>(const DecisionVariables&)>
+        solve,
+    std::span<const DecisionVariables> initialGuesses) {
+  std::vector<std::future<MultistartResult<DecisionVariables>>> futures;
+  for (const auto& initialGuess : initialGuesses) {
+    futures.emplace_back(std::async(std::launch::async, solve, initialGuess));
+  }
+
+  std::vector<MultistartResult<DecisionVariables>> results;
+  for (auto& future : futures) {
+    results.emplace_back(future.get());
+  }
+
+  return *std::min_element(
+      results.cbegin(), results.cend(), [](const auto& a, const auto& b) {
+        // Prioritize successful solve
+        if (a.status.exitCondition == SolverExitCondition::kSuccess &&
+            b.status.exitCondition != SolverExitCondition::kSuccess) {
+          return true;
+        }
+
+        // Otherwise prioritize solution with lower cost
+        return a.status.cost < b.status.cost;
+      });
+}
+
+}  // namespace sleipnir

jormungandr/cpp/Docstrings.hpp

@@ -293,6 +293,35 @@ Parameter ``lhs``:
 Parameter ``rhs``:
     Right-hand side.)doc";
 
+static const char *__doc_sleipnir_Multistart =
+R"doc(Solves an optimization problem from different starting points in
+parallel, then returns the solution with the lowest cost.
+
+Each solve is performed on a separate thread. Solutions from
+successful solves are always preferred over solutions from
+unsuccessful solves, and cost (lower is better) is the tiebreaker
+between successful solves.
+
+Template parameter ``DecisionVariables``:
+    The type containing the decision variable initial guess.
+
+Parameter ``solve``:
+    A user-provided function that takes a decision variable initial
+    guess and returns a MultistartResult.
+
+Parameter ``initialGuesses``:
+    A list of decision variable initial guesses to try.)doc";
+
+static const char *__doc_sleipnir_MultistartResult =
+R"doc(The result of a multistart solve.
+
+Template parameter ``DecisionVariables``:
+    The type containing the decision variable initial guess.)doc";
+
+static const char *__doc_sleipnir_MultistartResult_status = R"doc()doc";
+
+static const char *__doc_sleipnir_MultistartResult_variables = R"doc()doc";
+
 static const char *__doc_sleipnir_OCPSolver =
 R"doc(This class allows the user to pose and solve a constrained optimal
 control problem (OCP) in a variety of ways.

jormungandr/optimization/__init__.py

@@ -1 +1,39 @@
+import concurrent.futures
+
 from .._jormungandr.optimization import *
+
+
+def multistart(solve, initial_guesses):
+    """
+    Solves an optimization problem from different starting points in parallel,
+    then returns the solution with the lowest cost.
+
+    Each solve is performed on a separate thread. Solutions from successful
+    solves are always preferred over solutions from unsuccessful solves, and
+    cost (lower is better) is the tiebreaker between successful solves.
+
+    Parameter ``solve``:
+        A user-provided function that takes a decision variable initial guess
+        and returns a MultistartResult.
+
+    Parameter ``initial_guesses``:
+        A list of decision variable initial guesses to try.
+    """
+    with concurrent.futures.ThreadPoolExecutor(
+        max_workers=len(initial_guesses)
+    ) as executor:
+        futures = [
+            executor.submit(solve, initial_guess) for initial_guess in initial_guesses
+        ]
+        results = [
+            future.result() for future in concurrent.futures.as_completed(futures)
+        ]
+
+    # Prioritize successful solve, otherwise prioritize solution with lower cost
+    return min(
+        results,
+        key=lambda x: (
+            int(x[0].exit_condition != SolverExitCondition.SUCCESS),
+            x[0].cost,
+        ),
+    )

jormungandr/test/optimization/nonlinear_problem_test.py

@@ -2,7 +2,11 @@
 
 import jormungandr.autodiff as autodiff
 from jormungandr.autodiff import ExpressionType
-from jormungandr.optimization import OptimizationProblem, SolverExitCondition
+from jormungandr.optimization import (
+    OptimizationProblem,
+    SolverExitCondition,
+    multistart,
+)
 
 import numpy as np
 import pytest
@@ -91,6 +95,47 @@ def test_rosenbrock_with_disk_constraint():
             assert y.value() == pytest.approx(1.0, abs=1e-1)
 
 
+def mishras_bird_function_solve(input):
+    problem = OptimizationProblem()
+
+    x = problem.decision_variable()
+    x.set_value(input.x)
+    y = problem.decision_variable()
+    y.set_value(input.y)
+
+    # https://en.wikipedia.org/wiki/Test_functions_for_optimization#Test_functions_for_constrained_optimization
+    problem.minimize(
+        autodiff.sin(y) * autodiff.exp((1 - autodiff.cos(x)) ** 2)
+        + autodiff.cos(x) * autodiff.exp((1 - autodiff.sin(y)) ** 2)
+        + (x - y) ** 2
+    )
+
+    problem.subject_to((x + 5) ** 2 + (y + 5) ** 2 < 25)
+
+    return (problem.solve(), DecisionVariables(x.value(), y.value()))
+
+
+class DecisionVariables:
+    def __init__(self, x, y):
+        self.x = x
+        self.y = y
+
+
+def test_mishras_bird_function():
+    status, variables = multistart(
+        mishras_bird_function_solve,
+        [DecisionVariables(-3, -8), DecisionVariables(-3, -1.5)],
+    )
+
+    assert status.cost_function_type == ExpressionType.NONLINEAR
+    assert status.equality_constraint_type == ExpressionType.NONE
+    assert status.inequality_constraint_type == ExpressionType.QUADRATIC
+    assert status.exit_condition == SolverExitCondition.SUCCESS
+
+    assert variables.x == pytest.approx(-3.130246803458174, abs=1e-15)
+    assert variables.y == pytest.approx(-1.5821421769364057, abs=1e-15)
+
+
 def test_narrow_feasible_region():
     problem = OptimizationProblem()
 

test/src/optimization/NonlinearProblemTest.cpp

@@ -1,8 +1,11 @@
 // Copyright (c) Sleipnir contributors
 
+#include <vector>
+
 #include <catch2/catch_approx.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <fmt/core.h>
+#include <sleipnir/optimization/Multistart.hpp>
 #include <sleipnir/optimization/OptimizationProblem.hpp>
 
 #include "CatchStringConverters.hpp"
@@ -105,6 +108,47 @@ TEST_CASE("NonlinearProblem - Rosenbrock with disk constraint",
   }
 }
 
+TEST_CASE("NonlinearProblem - Mishra's Bird function", "[NonlinearProblem]") {
+  struct DecisionVariables {
+    double x;
+    double y;
+  };
+
+  auto Solve = [](const DecisionVariables& input)
+      -> sleipnir::MultistartResult<DecisionVariables> {
+    sleipnir::OptimizationProblem problem;
+
+    auto x = problem.DecisionVariable();
+    x.SetValue(input.x);
+    auto y = problem.DecisionVariable();
+    y.SetValue(input.y);
+
+    // https://en.wikipedia.org/wiki/Test_functions_for_optimization#Test_functions_for_constrained_optimization
+    problem.Minimize(sleipnir::sin(y) *
+                         sleipnir::exp(sleipnir::pow(1 - sleipnir::cos(x), 2)) +
+                     sleipnir::cos(x) *
+                         sleipnir::exp(sleipnir::pow(1 - sleipnir::sin(y), 2)) +
+                     sleipnir::pow(x - y, 2));
+
+    problem.SubjectTo(sleipnir::pow(x + 5, 2) + sleipnir::pow(y + 5, 2) < 25);
+
+    return {problem.Solve(), DecisionVariables{x.Value(), y.Value()}};
+  };
+
+  auto [status, variables] = sleipnir::Multistart<DecisionVariables>(
+      Solve,
+      std::vector{DecisionVariables{-3, -8}, DecisionVariables{-3, -1.5}});
+
+  CHECK(status.costFunctionType == sleipnir::ExpressionType::kNonlinear);
+  CHECK(status.equalityConstraintType == sleipnir::ExpressionType::kNone);
+  CHECK(status.inequalityConstraintType ==
+        sleipnir::ExpressionType::kQuadratic);
+  CHECK(status.exitCondition == sleipnir::SolverExitCondition::kSuccess);
+
+  CHECK(variables.x == Catch::Approx(-3.130246803458174).margin(1e-15));
+  CHECK(variables.y == Catch::Approx(-1.5821421769364057).margin(1e-15));
+}
+
 TEST_CASE("NonlinearProblem - Narrow feasible region", "[NonlinearProblem]") {
   sleipnir::OptimizationProblem problem;