From e19e84589dea1ca4efa3fee6aff1ef9e76f7d767 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 15 Nov 2024 12:25:37 +0100 Subject: [PATCH] Wrote two functional tests (an LP solve and a QP solve) for BQPD --- CMakeLists.txt | 1 + unotest/functional_tests/BQPDSolverTests.cpp | 144 ++++++++++++++++++ unotest/functional_tests/HiGHSSolverTests.cpp | 2 +- 3 files changed, 146 insertions(+), 1 deletion(-) create mode 100644 unotest/functional_tests/BQPDSolverTests.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2d711e8b..34b37a73 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -130,6 +130,7 @@ if(NOT BQPD) message(WARNING "Optional library BQPD was not found.") else() list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp uno/solvers/BQPD/wdotd.f) + list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/BQPDSolverTests.cpp) link_to_uno(bqpd ${BQPD}) endif() diff --git a/unotest/functional_tests/BQPDSolverTests.cpp b/unotest/functional_tests/BQPDSolverTests.cpp new file mode 100644 index 00000000..28c4996e --- /dev/null +++ b/unotest/functional_tests/BQPDSolverTests.cpp @@ -0,0 +1,144 @@ +// Copyright (c) 2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "optimization/Direction.hpp" +#include "linear_algebra/RectangularMatrix.hpp" +#include "linear_algebra/SparseVector.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" +#include "optimization/WarmstartInformation.hpp" +#include "options/Options.hpp" +#include "solvers/BQPD/BQPDSolver.hpp" +#include "tools/Infinity.hpp" + +using namespace uno; + +TEST(BQPDSolver, LP) { + // https://ergo-code.github.io/HiGHS/stable/interfaces/cpp/library/ + // Min f = x_0 + x_1 + 3 + // s.t. x_1 <= 7 + // 5 <= x_0 + 2x_1 <= 15 + // 6 <= 3x_0 + 2x_1 + // 0 <= x_0 <= 4; 1 <= x_1 + + const size_t number_variables = 2; + const size_t number_constraints = 3; + const size_t number_linear_objective_nonzeros = 2; + const size_t number_jacobian_nonzeros = 5; + const size_t number_hessian_nonzeros = 0; + Options options(false); + options["print_subproblem"] = "false"; + BQPDSolver bqpd_solver(number_variables, number_constraints, number_linear_objective_nonzeros, number_jacobian_nonzeros, number_hessian_nonzeros, + BQPDProblemType::LP, options); + + // create the LP + SparseVector linear_objective(number_variables); + linear_objective.insert(0, 1.); + linear_objective.insert(1, 1.); + + const std::vector variables_lower_bounds{0., 1.}; + const std::vector variables_upper_bounds{4., INF}; + const std::vector constraints_lower_bounds{-INF, 5., 6.}; + const std::vector constraints_upper_bounds{7., 15., INF}; + + RectangularMatrix constraint_jacobian(number_constraints, number_variables); + constraint_jacobian[0].insert(1, 1.); + constraint_jacobian[1].insert(0, 1.); + constraint_jacobian[1].insert(1, 2.); + constraint_jacobian[2].insert(0, 3.); + constraint_jacobian[2].insert(1, 2.); + + Direction direction(number_variables, number_constraints); + WarmstartInformation warmstart_information{}; + Vector initial_point{0., 0.}; + + bqpd_solver.solve_LP(number_variables, number_constraints, variables_lower_bounds, variables_upper_bounds, constraints_lower_bounds, + constraints_upper_bounds, linear_objective, constraint_jacobian, initial_point, direction, warmstart_information); + + ASSERT_EQ(direction.status, SubproblemStatus::OPTIMAL); + + const double tolerance = 1e-8; + // check primals + const std::vector primals_reference{0.5, 2.25}; + for (size_t index: Range(number_variables)) { + EXPECT_NEAR(direction.primals[index], primals_reference[index], tolerance); + } + // check duals + const std::vector constraint_duals_reference{0., 0.25, 0.25}; + for (size_t index: Range(number_constraints)) { + EXPECT_NEAR(direction.multipliers.constraints[index], constraint_duals_reference[index], tolerance); + } + const std::vector lower_bound_duals_reference{0., 0.}; + const std::vector upper_bound_duals_reference{0., 0.}; + for (size_t index: Range(number_variables)) { + EXPECT_NEAR(direction.multipliers.lower_bounds[index], lower_bound_duals_reference[index], tolerance); + EXPECT_NEAR(direction.multipliers.upper_bounds[index], upper_bound_duals_reference[index], tolerance); + } +} + +TEST(BQPDSolver, QP) { + // https://doc.cgal.org/latest/QP_solver/index.html#title4 + // Min f = 1/2 * (2 x_0^2 + 8 x_1^2) - 32 x_1 + // s.t. x_0 + x_1 <= 7 + // -x_0 + 2x_1 <= 4 + // 0 <= x_0; 0 <= x_1 <= 4 + + const size_t number_variables = 2; + const size_t number_constraints = 2; + const size_t number_hessian_nonzeros = 2; + const size_t number_linear_objective_nonzeros = 1; + const size_t number_jacobian_nonzeros = 4; + Options options(false); + options["print_subproblem"] = "false"; + options["BQPD_kmax"] = "500"; + BQPDSolver bqpd_solver(number_variables, number_constraints, number_linear_objective_nonzeros, number_jacobian_nonzeros, number_hessian_nonzeros, + BQPDProblemType::QP, options); + + // create the QP + SparseVector linear_objective(number_variables); + linear_objective.insert(1, -32.); + + const std::vector variables_lower_bounds{0., 0.}; + const std::vector variables_upper_bounds{INF, 4.}; + const std::vector constraints_lower_bounds{-INF, -INF}; + const std::vector constraints_upper_bounds{7., 4.}; + + RectangularMatrix constraint_jacobian(number_constraints, number_variables); + constraint_jacobian[0].insert(0, 1.); + constraint_jacobian[0].insert(1, 1.); + constraint_jacobian[1].insert(0, -1.); + constraint_jacobian[1].insert(1, 2.); + + Direction direction(number_variables, number_constraints); + WarmstartInformation warmstart_information{}; + Vector initial_point{0., 0.}; + + SymmetricMatrix hessian(number_variables, number_hessian_nonzeros, false, "CSC"); + hessian.insert(2., 0, 0); + hessian.finalize_column(0); + hessian.insert(8., 1, 1); + hessian.finalize_column(1); + + bqpd_solver.solve_QP(number_variables, number_constraints, variables_lower_bounds, variables_upper_bounds, constraints_lower_bounds, + constraints_upper_bounds, linear_objective, constraint_jacobian, hessian, initial_point, direction, warmstart_information); + + ASSERT_EQ(direction.status, SubproblemStatus::OPTIMAL); + + const double tolerance = 1e-8; + // check primals + const std::vector primals_reference{2., 3.}; + for (size_t index: Range(number_variables)) { + EXPECT_NEAR(direction.primals[index], primals_reference[index], tolerance); + } + // check duals + const std::vector constraint_duals_reference{0., -4.}; + for (size_t index: Range(number_constraints)) { + EXPECT_NEAR(direction.multipliers.constraints[index], constraint_duals_reference[index], tolerance); + } + const std::vector lower_bound_duals_reference{0., 0.}; + const std::vector upper_bound_duals_reference{0., 0.}; + for (size_t index: Range(number_variables)) { + EXPECT_NEAR(direction.multipliers.lower_bounds[index], lower_bound_duals_reference[index], tolerance); + EXPECT_NEAR(direction.multipliers.upper_bounds[index], upper_bound_duals_reference[index], tolerance); + } +} \ No newline at end of file diff --git a/unotest/functional_tests/HiGHSSolverTests.cpp b/unotest/functional_tests/HiGHSSolverTests.cpp index a8c34d3f..de7aea69 100644 --- a/unotest/functional_tests/HiGHSSolverTests.cpp +++ b/unotest/functional_tests/HiGHSSolverTests.cpp @@ -12,7 +12,7 @@ using namespace uno; -TEST(HiGHSSolver, DocumentationLP) { +TEST(HiGHSSolver, LP) { // https://ergo-code.github.io/HiGHS/stable/interfaces/cpp/library/ // Min f = x_0 + x_1 + 3 // s.t. x_1 <= 7