Skip to content

Commit

Permalink
Use corresponding QP solver for LCP problem resolution as well
Browse files Browse the repository at this point in the history
  • Loading branch information
olivier-roussel committed Jan 21, 2025
1 parent e9d07de commit 11dfb2c
Show file tree
Hide file tree
Showing 14 changed files with 305 additions and 50 deletions.
47 changes: 29 additions & 18 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,6 @@ set(SOURCE_FILES
${SRC_DIR}/component/solver/QPInverseProblemSolver.cpp
${SRC_DIR}/component/solver/modules/ContactHandler.cpp
${SRC_DIR}/component/solver/modules/ConstraintHandler.cpp
${SRC_DIR}/component/solver/modules/LCPQPSolver.cpp
${SRC_DIR}/component/solver/modules/NLCPSolver.cpp
${SRC_DIR}/component/solver/modules/QPInverseProblem.cpp
${SRC_DIR}/component/solver/modules/QPInverseProblemImpl.cpp
Expand All @@ -139,12 +138,6 @@ endif()
set(DOC_FILES README.md)
file(GLOB_RECURSE EXAMPLE_FILES examples "*.pyscn" "*.py" "*.md" "*.psl" "*.pslx" "*.scn" "*.xml")

option(SOFTROBOTSINVERSE_INSTALL_HEADERS "Install the headers" ON)
if(SOFTROBOTSINVERSE_INSTALL_HEADERS)
add_library(${PROJECT_NAME} SHARED ${SOURCE_FILES} ${DOC_FILES} ${EXAMPLE_FILES} ${HEADER_FILES})
else()
add_library(${PROJECT_NAME} SHARED ${SOURCE_FILES} ${DOC_FILES} ${EXAMPLE_FILES} "${SRC_DIR}/component/config.h.in")
endif()

option(SOFTROBOTSINVERSE_ENABLE_PROXQP "Build proxQP wrapper for inverse problems" OFF)
option(SOFTROBOTSINVERSE_ENABLE_QPOASES "Build qpOASES wrapper for inverse problems" ON)
Expand All @@ -154,15 +147,39 @@ if(NOT SOFTROBOTSINVERSE_ENABLE_PROXQP AND NOT SOFTROBOTSINVERSE_ENABLE_QPOASES)
endif()

if(SOFTROBOTSINVERSE_ENABLE_PROXQP)
find_package(proxsuite REQUIRED)
message("Sofa.SoftRobots.Inverse: proxsuite FOUND. Will build with proxsuite support.")

list(APPEND HEADER_FILES
${SRC_DIR}/component/solver/modules/QPInverseProblemProxQP.h
${SRC_DIR}/component/solver/modules/LCPQPSolverProxQP.h
)
list(APPEND SOURCE_FILES
${SRC_DIR}/component/solver/modules/QPInverseProblemProxQP.cpp
${SRC_DIR}/component/solver/modules/LCPQPSolverProxQP.cpp
)
endif()

if(SOFTROBOTSINVERSE_ENABLE_QPOASES)
list(APPEND HEADER_FILES
${SRC_DIR}/component/solver/modules/QPInverseProblemQPOases.h
${SRC_DIR}/component/solver/modules/LCPQPSolverQPOases.h
)
list(APPEND SOURCE_FILES
${SRC_DIR}/component/solver/modules/QPInverseProblemQPOases.cpp
${SRC_DIR}/component/solver/modules/LCPQPSolverQPOases.cpp
)
endif()

option(SOFTROBOTSINVERSE_INSTALL_HEADERS "Install the headers" ON)
if(SOFTROBOTSINVERSE_INSTALL_HEADERS)
add_library(${PROJECT_NAME} SHARED ${SOURCE_FILES} ${DOC_FILES} ${EXAMPLE_FILES} ${HEADER_FILES})
else()
add_library(${PROJECT_NAME} SHARED ${SOURCE_FILES} ${DOC_FILES} ${EXAMPLE_FILES} "${SRC_DIR}/component/config.h.in")
endif()


if(SOFTROBOTSINVERSE_ENABLE_PROXQP)
find_package(proxsuite REQUIRED)
message("Sofa.SoftRobots.Inverse: proxsuite FOUND. Will build with proxsuite support.")

target_link_libraries(${PROJECT_NAME} PUBLIC proxsuite::proxsuite)
# XXX proxsuite vectorization support via SIMDE and activated by the compilation options '-march=native' or `-mavx2 -mavx512f`
# This will not be available if proxsuite has been fetched
Expand All @@ -174,19 +191,13 @@ endif()
if(SOFTROBOTSINVERSE_ENABLE_QPOASES)
find_package(qpOASES QUIET)

list(APPEND HEADER_FILES
${SRC_DIR}/component/solver/modules/QPInverseProblemQPOases.h
)
list(APPEND SOURCE_FILES
${SRC_DIR}/component/solver/modules/QPInverseProblemQPOases.cpp
)

if(qpOASES_FOUND)
message("Sofa.SoftRobots.Inverse: external qpOASES FOUND. Will build with qpOASES support.")
include_directories(${qpOASES_INCLUDE_DIRS})
target_link_libraries(${PROJECT_NAME} PUBLIC ${qpOASES_LIBRARIES})
set(QPOASES_INCLUDE_DIR ${qpOASES_INCLUDE_DIRS} CACHE INTERNAL "")
else()
message("External package qpOASES not found, using embedded version...")
message("External package qpOASES not found, using embedded version for qpOASES support")
set(OASES_LIBRARY_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/extlibs/qpOASES-3.2.0/")
add_subdirectory(${OASES_LIBRARY_DIRECTORY} extlibs/libqpOASES)
find_package(libqpOASES REQUIRED)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,12 @@ QPInverseProblemSolver::QPInverseProblemSolver()
, m_CP3(nullptr)
{
sofa::helper::OptionsGroup qpSolvers{"qpOASES" , "proxQP"};
#if defined SOFTROBOTSINVERSE_ENABLE_PROXQP && !defined SOFTROBOTSINVERSE_ENABLE_QPOASES
qpSolvers.setSelectedItem(QPSolverImpl::PROXQP);
#else
qpSolvers.setSelectedItem(QPSolverImpl::QPOASES);
#endif

d_qpSolver.setValue(qpSolvers);

d_graph.setWidget("graph");
Expand Down
11 changes: 2 additions & 9 deletions src/SoftRobots.Inverse/component/solver/modules/LCPQPSolver.h
Original file line number Diff line number Diff line change
Expand Up @@ -28,24 +28,17 @@
******************************************************************************/
#pragma once

#include <sofa/type/vector.h>

// LCP solver using qpOASES QP solver
// QP formulation for solving a LCP with a symmetric matrix M.

namespace softrobotsinverse::solver::module {

using sofa::type::vector;

class LCPQPSolver
{

public:
virtual ~LCPQPSolver() = default;

LCPQPSolver(){}
~LCPQPSolver(){}

void solve(int dim, double*q, double**M, double*res);
virtual void solve(int dim, double*q, double**M, double*res) = 0;
};

} // namespace
Expand Down
143 changes: 143 additions & 0 deletions src/SoftRobots.Inverse/component/solver/modules/LCPQPSolverProxQP.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
/******************************************************************************
* SOFA, Simulation Open-Framework Architecture *
* (c) 2006 INRIA, USTL, UJF, CNRS, MGH *
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by *
* the Free Software Foundation; either version 2.1 of the License, or (at *
* your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*******************************************************************************
* Plugin SoftRobots.Inverse *
* *
* This plugin is distributed under the GNU AGPL v3 (Affero General *
* Public License) license. *
* *
* Authors: Christian Duriez, Eulalie Coevoet, Yinoussa Adagolodjo *
* *
* (c) 2023 INRIA *
* *
* Contact information: https://project.inria.fr/softrobot/contact/ *
******************************************************************************/

#include <sofa/helper/logging/Messaging.h>
#include <SoftRobots.Inverse/component/solver/modules/LCPQPSolverProxQP.h>

#include <proxsuite/proxqp/dense/dense.hpp>

#include <optional>

namespace softrobotsinverse::solver::module
{

void LCPQPSolverProxQP::solve(int dim, double*q, double**M, double*res)
{
const int n_in = dim;
const int n_eq = 0;

//////////////////////////////////////////
// Constraints on lambda
//////////////////////////////////////////
Eigen::VectorXd lbox(dim); // l
Eigen::VectorXd ubox(dim); // u
for (int i=0; i<dim; ++i)
{
ubox[i] = 1e99;
lbox[i] = 0.;
}


//////////////////////////////////////////
// Conversion to Eigen type
// TODO: try Eigen Mapping instead
//////////////////////////////////////////
Eigen::MatrixXd H(dim, dim); // Q
Eigen::VectorXd g(dim); // c

for (int i = 0; i < dim; ++i)
{
for (int j = 0; j < dim; ++j)
H(i, j) = M[i][j];
g[i] = q[i];
}


//////////////////////////////////////////
// Inequality constraint matrices
// TODO: try Eigen Mapping instead
//////////////////////////////////////////
Eigen::VectorXd l(n_in); // bl
Eigen::VectorXd u(n_in); // bu
Eigen::MatrixXd C(n_in, n_in); // A

for (int i = 0; i < n_in; ++i)
{
for (int j = 0; j < n_in; ++j)
C(i, j) = M[i][j];

u[i] = 1e99;
l[i] = -q[i];
}


/////////////////////////////////////////
// Solve
/////////////////////////////////////////


proxsuite::proxqp::dense::QP<double> qp(dim, n_eq, n_in, true); // create the QP object with box constraints

const double rho = 1e-14; // default is 1e-6, but has to be set smaller due to low values in Hessian when using data in meters (e.g. from URDF models)

qp.init(H, g, std::nullopt, std::nullopt, C, l, u, lbox, ubox, true, rho); // initialize the model
// qp.init(H, g, A, b, C, l, u, lbox, ubox, true, rho); // initialize the model

qp.settings.eps_rel = 0.;
qp.settings.eps_abs = 1e-12; // default is 1e-5
qp.settings.check_duality_gap = true;
qp.settings.verbose = false;
qp.settings.max_iter = 3000;
// same here, default is 1e-4 for these epsilons, but has to be set smaller due to low values in Hessian when using data in meters (e.g. from URDF models)
qp.settings.eps_primal_inf = 1e-12;
qp.settings.eps_dual_inf = 1e-12;

qp.solve();

switch(qp.results.info.status)
{
case proxsuite::proxqp::QPSolverOutput::PROXQP_SOLVED:
break;
case proxsuite::proxqp::QPSolverOutput::PROXQP_MAX_ITER_REACHED:
msg_warning("QPInverseProblemImpl") << "Solver status: PROXQP_MAX_ITER_REACHED";
break;
case proxsuite::proxqp::QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE:
msg_warning("QPInverseProblemImpl") << "Solver status: PROXQP_PRIMAL_INFEASIBLE";
break;
case proxsuite::proxqp::QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE:
msg_warning("QPInverseProblemImpl") << "Solver status: PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE";
break;
case proxsuite::proxqp::QPSolverOutput::PROXQP_DUAL_INFEASIBLE:
msg_warning("QPInverseProblemImpl") << "Solver status: PROXQP_DUAL_INFEASIBLE";
break;
case proxsuite::proxqp::QPSolverOutput::PROXQP_NOT_RUN:
msg_warning("QPInverseProblemImpl") << "Solver status: PROXQP_NOT_RUN";
break;
default:
msg_error("QPInverseProblemImpl") << "Unknown solver status";
}

// get primal solution
for (int i = 0; i < dim; ++i)
{
res[i] = qp.results.x[i];
}
}

} // namespace
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
/******************************************************************************
* SOFA, Simulation Open-Framework Architecture *
* (c) 2006 INRIA, USTL, UJF, CNRS, MGH *
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by *
* the Free Software Foundation; either version 2.1 of the License, or (at *
* your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*******************************************************************************
* Plugin SoftRobots.Inverse *
* *
* This plugin is distributed under the GNU AGPL v3 (Affero General *
* Public License) license. *
* *
* Authors: Christian Duriez, Eulalie Coevoet, Yinoussa Adagolodjo *
* *
* (c) 2023 INRIA *
* *
* Contact information: https://project.inria.fr/softrobot/contact/ *
******************************************************************************/
#pragma once

#include <SoftRobots.Inverse/component/solver/modules/LCPQPSolver.h>

// LCP solver using proxQP solver
// QP formulation for solving a LCP with a symmetric matrix M.

namespace softrobotsinverse::solver::module {

class LCPQPSolverProxQP: public LCPQPSolver
{

public:

LCPQPSolverProxQP() = default;
virtual ~LCPQPSolverProxQP() = default;

void solve(int dim, double*q, double**M, double*res) override;
};

} // namespace

Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,9 @@
******************************************************************************/

#include <qpOASES.hpp>
#include <sofa/type/vector.h>
#include <sofa/helper/logging/Messaging.h>
#include <SoftRobots.Inverse/component/solver/modules/LCPQPSolver.h>
#include <SoftRobots.Inverse/component/solver/modules/LCPQPSolverQPOases.h>


namespace softrobotsinverse::solver::module
Expand All @@ -45,7 +46,7 @@ using qpOASES::int_t;
using sofa::type::vector;


void LCPQPSolver::solve(int dim, double*q, double**M, double*res)
void LCPQPSolverQPOases::solve(int dim, double*q, double**M, double*res)
{
//////////////////////////////////////////
// Constraints on lambda
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
/******************************************************************************
* SOFA, Simulation Open-Framework Architecture *
* (c) 2006 INRIA, USTL, UJF, CNRS, MGH *
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by *
* the Free Software Foundation; either version 2.1 of the License, or (at *
* your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*******************************************************************************
* Plugin SoftRobots.Inverse *
* *
* This plugin is distributed under the GNU AGPL v3 (Affero General *
* Public License) license. *
* *
* Authors: Christian Duriez, Eulalie Coevoet, Yinoussa Adagolodjo *
* *
* (c) 2023 INRIA *
* *
* Contact information: https://project.inria.fr/softrobot/contact/ *
******************************************************************************/
#pragma once

#include <SoftRobots.Inverse/component/solver/modules/LCPQPSolver.h>

#include <SoftRobots.Inverse/component/solver/modules/LCPQPSolverQPOases.h>

// LCP solver using qpOASES QP solver
// QP formulation for solving a LCP with a symmetric matrix M.

namespace softrobotsinverse::solver::module {

class LCPQPSolverQPOases: public LCPQPSolver
{

public:
LCPQPSolverQPOases() = default;
virtual ~LCPQPSolverQPOases() = default;

void solve(int dim, double*q, double**M, double*res) override;
};

} // namespace

Loading

0 comments on commit 11dfb2c

Please sign in to comment.