diff --git a/.github/workflows/sanitizers.yml b/.github/workflows/sanitizers.yml index fad2fedf..b494c5b0 100644 --- a/.github/workflows/sanitizers.yml +++ b/.github/workflows/sanitizers.yml @@ -8,7 +8,7 @@ concurrency: jobs: build: - timeout-minutes: 10 + timeout-minutes: 15 strategy: fail-fast: false matrix: @@ -16,10 +16,9 @@ jobs: - 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: diff --git a/CMakeLists.txt b/CMakeLists.txt index 39d8867a..43f19ec5 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/sleipnir/optimization/Multistart.hpp b/include/sleipnir/optimization/Multistart.hpp new file mode 100644 index 00000000..e25ad1ea --- /dev/null +++ b/include/sleipnir/optimization/Multistart.hpp @@ -0,0 +1,69 @@ +// Copyright (c) Sleipnir contributors + +#pragma once + +#include +#include +#include +#include + +#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 +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 +MultistartResult Multistart( + function_ref(const DecisionVariables&)> + solve, + std::span initialGuesses) { + std::vector>> futures; + for (const auto& initialGuess : initialGuesses) { + futures.emplace_back(std::async(std::launch::async, solve, initialGuess)); + } + + std::vector> 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 diff --git a/jormungandr/cpp/Docstrings.hpp b/jormungandr/cpp/Docstrings.hpp index e7f5dd29..dedd62f0 100644 --- a/jormungandr/cpp/Docstrings.hpp +++ b/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. diff --git a/jormungandr/optimization/__init__.py b/jormungandr/optimization/__init__.py index b73c0f47..675b17a4 100644 --- a/jormungandr/optimization/__init__.py +++ b/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, + ), + ) diff --git a/jormungandr/test/optimization/nonlinear_problem_test.py b/jormungandr/test/optimization/nonlinear_problem_test.py index 811eca94..cc8242a8 100644 --- a/jormungandr/test/optimization/nonlinear_problem_test.py +++ b/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() diff --git a/test/src/optimization/NonlinearProblemTest.cpp b/test/src/optimization/NonlinearProblemTest.cpp index 2b6c976a..2a1312b6 100644 --- a/test/src/optimization/NonlinearProblemTest.cpp +++ b/test/src/optimization/NonlinearProblemTest.cpp @@ -1,8 +1,11 @@ // Copyright (c) Sleipnir contributors +#include + #include #include #include +#include #include #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 { + 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( + 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;