diff --git a/.github/workflows/ros-build-test.yml b/.github/workflows/ros-build-test.yml
index 2580bdc93..0f3ac9433 100644
--- a/.github/workflows/ros-build-test.yml
+++ b/.github/workflows/ros-build-test.yml
@@ -30,7 +30,7 @@ jobs:
- name: System deps
run: |
apt-get update
- apt-get install -y git ninja-build liburdfdom-dev liboctomap-dev libassimp-dev
+ apt-get install -y git ninja-build liburdfdom-dev liboctomap-dev libassimp-dev checkinstall wget rsync
- uses: actions/checkout@v2
with:
@@ -46,14 +46,33 @@ jobs:
git clone --recurse-submodules https://github.com/leggedrobotics/pinocchio.git src/pinocchio
git clone --recurse-submodules https://github.com/leggedrobotics/hpp-fcl.git src/hpp-fcl
git clone https://github.com/leggedrobotics/ocs2_robotic_assets.git src/ocs2_robotic_assets
+
+ - name: Install RaiSim
+ run: |
+ git clone --depth 1 https://github.com/raisimTech/raisimLib.git -b v1.1.01 src/raisim_tech
+ cd src/raisim_tech
+ mkdir build
+ cd build
+ cmake .. -DPYTHON_EXECUTABLE=$(python3 -c "import sys; print(sys.executable)")
+ make -j4 && checkinstall
+
+ - name: Install ONNX Runtime
+ run: |
+ wget https://github.com/microsoft/onnxruntime/releases/download/v1.7.0/onnxruntime-linux-x64-1.7.0.tgz -P tmp/microsoft
+ tar xf tmp/microsoft/onnxruntime-linux-x64-1.7.0.tgz -C tmp/microsoft
+ rsync -a tmp/microsoft/onnxruntime-linux-x64-1.7.0/include/ /usr/local/include/onnxruntime
+ rsync -a tmp/microsoft/onnxruntime-linux-x64-1.7.0/lib/ /usr/local/lib
+ rsync -a src/ocs2/ocs2_mpcnet/ocs2_mpcnet_core/misc/onnxruntime/cmake/ /usr/local/share/onnxruntime
+
- name: Build (${{ matrix.build_type }})
shell: bash
run: |
source /opt/ros/noetic/setup.bash
- catkin_make_isolated --use-ninja --merge --only-pkg-with-deps ocs2 --ignore-pkg ocs2_raisim --cmake-args -DCMAKE_BUILD_TYPE=${{ matrix.build_type }}
+ catkin_make_isolated --use-ninja --merge --only-pkg-with-deps ocs2 --cmake-args -DCMAKE_BUILD_TYPE=${{ matrix.build_type }}
- name: Test
shell: bash
run: |
source devel_isolated/setup.bash
catkin_make_isolated --use-ninja --merge --only-pkg-with-deps ocs2 --catkin-make-args run_tests
+ catkin_test_results
diff --git a/.gitignore b/.gitignore
index 15f4690f1..46f83bb33 100644
--- a/.gitignore
+++ b/.gitignore
@@ -25,3 +25,4 @@ ocs2_ddp/test/ddp_test_generated/
*.out
*.synctex.gz
.vscode/
+runs/
diff --git a/jenkins-pipeline b/jenkins-pipeline
index 9c91ed1ea..dbe1a5c14 100644
--- a/jenkins-pipeline
+++ b/jenkins-pipeline
@@ -1,5 +1,5 @@
library 'continuous_integration_pipeline'
-ciPipeline("--ros-distro noetic --publish-doxygen --recipes raisimlib\
+ciPipeline("--ros-distro noetic --publish-doxygen --recipes onnxruntime raisimlib\
--dependencies 'git@github.com:leggedrobotics/hpp-fcl.git;master;git'\
'git@github.com:leggedrobotics/pinocchio.git;master;git'\
'git@github.com:leggedrobotics/ocs2_robotic_assets.git;main;git'\
diff --git a/ocs2/package.xml b/ocs2/package.xml
index 834b4090e..fcd5dd285 100644
--- a/ocs2/package.xml
+++ b/ocs2/package.xml
@@ -16,14 +16,17 @@
ocs2_oc
ocs2_qp_solver
ocs2_ddp
+ ocs2_slp
ocs2_sqp
ocs2_ros_interfaces
ocs2_python_interface
ocs2_pinocchio
ocs2_robotic_tools
+ ocs2_perceptive
ocs2_robotic_examples
ocs2_thirdparty
ocs2_raisim
+ ocs2_mpcnet
diff --git a/ocs2_core/cmake/ocs2_cxx_flags.cmake b/ocs2_core/cmake/ocs2_cxx_flags.cmake
index 1735d016f..7f38eb333 100644
--- a/ocs2_core/cmake/ocs2_cxx_flags.cmake
+++ b/ocs2_core/cmake/ocs2_cxx_flags.cmake
@@ -22,5 +22,5 @@ list(APPEND OCS2_CXX_FLAGS
)
# Cpp standard version
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
diff --git a/ocs2_core/include/ocs2_core/Types.h b/ocs2_core/include/ocs2_core/Types.h
index d09cc26c0..3c161783f 100644
--- a/ocs2_core/include/ocs2_core/Types.h
+++ b/ocs2_core/include/ocs2_core/Types.h
@@ -71,17 +71,8 @@ using matrix_array2_t = std::vector;
/** Array of arrays of dynamic matrix trajectory type. */
using matrix_array3_t = std::vector;
-/** Eigen scalar type. */
-using eigen_scalar_t = Eigen::Matrix;
-/** Eigen scalar trajectory type. */
-using eigen_scalar_array_t = std::vector;
-/** Array of eigen scalar trajectory type. */
-using eigen_scalar_array2_t = std::vector;
-/** Array of arrays of eigen scalar trajectory type. */
-using eigen_scalar_array3_t = std::vector;
-
/**
- * Defines the linear approximation
+ * Defines the linear approximation of a scalar function
* f(x,u) = dfdx' dx + dfdu' du + f
*/
struct ScalarFunctionLinearApproximation {
@@ -134,7 +125,7 @@ std::ostream& operator<<(std::ostream& out, const ScalarFunctionLinearApproximat
*
* @param[in] stateDim: Number of states.
* @param[in] inputDim: Number of inputs.
- * @param[in] data: Given quadratic approximation.
+ * @param[in] data: Given linear approximation.
* @param[in] dataName: The name of the data which appears in the output error message.
* @return The description of the error. If there was no error it would be empty;
*/
@@ -148,7 +139,7 @@ inline ScalarFunctionLinearApproximation operator*(scalar_t scalar, ScalarFuncti
}
/**
- * Defines the quadratic approximation
+ * Defines the quadratic approximation of a scalar function
* f(x,u) = 1/2 dx' dfdxx dx + du' dfdux dx + 1/2 du' dfduu du + dfdx' dx + dfdu' du + f
*/
struct ScalarFunctionQuadraticApproximation {
diff --git a/ocs2_core/include/ocs2_core/constraint/StateConstraintCollection.h b/ocs2_core/include/ocs2_core/constraint/StateConstraintCollection.h
index 3bb8a0d53..79b112fce 100644
--- a/ocs2_core/include/ocs2_core/constraint/StateConstraintCollection.h
+++ b/ocs2_core/include/ocs2_core/constraint/StateConstraintCollection.h
@@ -49,11 +49,14 @@ class StateConstraintCollection : public Collection {
~StateConstraintCollection() override = default;
StateConstraintCollection* clone() const override;
- /** Returns the number of active constraints at given time. */
- virtual size_t getNumConstraints(scalar_t time) const;
+ /** Returns the number of active constraints at a given time. */
+ size_t getNumConstraints(scalar_t time) const;
- /** Get the constraint vector value */
- virtual vector_t getValue(scalar_t time, const vector_t& state, const PreComputation& preComp) const;
+ /** Returns the number of active constraints at a given time for each term. If a term is inactive, its size is zero. */
+ size_array_t getTermsSize(scalar_t time) const;
+
+ /** Get an array of all constraints. If a term is inactive, the corresponding element is a vector of size zero. */
+ virtual vector_array_t getValue(scalar_t time, const vector_t& state, const PreComputation& preComp) const;
/** Get the constraint linear approximation */
virtual VectorFunctionLinearApproximation getLinearApproximation(scalar_t time, const vector_t& state,
@@ -68,4 +71,4 @@ class StateConstraintCollection : public Collection {
StateConstraintCollection(const StateConstraintCollection& other);
};
-} // namespace ocs2
\ No newline at end of file
+} // namespace ocs2
diff --git a/ocs2_core/include/ocs2_core/constraint/StateInputConstraintCollection.h b/ocs2_core/include/ocs2_core/constraint/StateInputConstraintCollection.h
index f102c593b..aecee543f 100644
--- a/ocs2_core/include/ocs2_core/constraint/StateInputConstraintCollection.h
+++ b/ocs2_core/include/ocs2_core/constraint/StateInputConstraintCollection.h
@@ -49,11 +49,14 @@ class StateInputConstraintCollection : public Collection {
~StateInputConstraintCollection() override = default;
StateInputConstraintCollection* clone() const override;
- /** Returns the number of active constraints at given time. */
- virtual size_t getNumConstraints(scalar_t time) const;
+ /** Returns the number of active constraints at a given time. */
+ size_t getNumConstraints(scalar_t time) const;
- /** Get the constraint vector value */
- virtual vector_t getValue(scalar_t time, const vector_t& state, const vector_t& input, const PreComputation& preComp) const;
+ /** Returns the number of active constraints at a given time for each term. If a term is inactive, its size is zero. */
+ size_array_t getTermsSize(scalar_t time) const;
+
+ /** Get an array of all constraints. If a term is inactive, the corresponding element is a vector of size zero. */
+ virtual vector_array_t getValue(scalar_t time, const vector_t& state, const vector_t& input, const PreComputation& preComp) const;
/** Get the constraint linear approximation */
virtual VectorFunctionLinearApproximation getLinearApproximation(scalar_t time, const vector_t& state, const vector_t& input,
@@ -68,4 +71,4 @@ class StateInputConstraintCollection : public Collection {
StateInputConstraintCollection(const StateInputConstraintCollection& other);
};
-} // namespace ocs2
\ No newline at end of file
+} // namespace ocs2
diff --git a/ocs2_core/include/ocs2_core/control/ControllerType.h b/ocs2_core/include/ocs2_core/control/ControllerType.h
index 65aa1962f..bbae196ee 100644
--- a/ocs2_core/include/ocs2_core/control/ControllerType.h
+++ b/ocs2_core/include/ocs2_core/control/ControllerType.h
@@ -34,6 +34,6 @@ namespace ocs2 {
/**
* Enum class for specifying controller type
*/
-enum class ControllerType { UNKNOWN, FEEDFORWARD, LINEAR, NEURAL_NETWORK };
+enum class ControllerType { UNKNOWN, FEEDFORWARD, LINEAR, ONNX, BEHAVIORAL };
} // namespace ocs2
diff --git a/ocs2_core/include/ocs2_core/integration/TrapezoidalIntegration.h b/ocs2_core/include/ocs2_core/integration/TrapezoidalIntegration.h
index 6adb50176..37088b127 100644
--- a/ocs2_core/include/ocs2_core/integration/TrapezoidalIntegration.h
+++ b/ocs2_core/include/ocs2_core/integration/TrapezoidalIntegration.h
@@ -29,22 +29,32 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#pragma once
+#include
#include
namespace ocs2 {
-template
-SCALAR_T trapezoidalIntegration(const std::vector& timeTrajectory, const std::vector& valueTrajectory) {
+/**
+ * Compute the trapezoidal integration of a trajectory of VALUE_T given the time stamp timeTrajectory and initial value initialValue.
+ *
+ * @note It requires that the VALUE_T has overwrite operator+(VALUE_T, VALUE_T) and define VALUE_T::operator*(SCALAR_T)
+ */
+template
+VALUE_T trapezoidalIntegration(const std::vector& timeTrajectory, const std::vector& valueTrajectory,
+ VALUE_T initialValue) {
+ assert(timeTrajectory.size() == valueTrajectory.size());
+
if (timeTrajectory.size() < 2) {
- return 0.0;
+ return initialValue;
}
- SCALAR_T areaUnderCurve = 0.0;
for (size_t k = 1; k < timeTrajectory.size(); k++) {
- areaUnderCurve += 0.5 * (valueTrajectory[k] + valueTrajectory[k - 1]) * (timeTrajectory[k] - timeTrajectory[k - 1]);
+ auto temp = valueTrajectory[k - 1] + valueTrajectory[k];
+ temp *= (0.5 * (timeTrajectory[k] - timeTrajectory[k - 1]));
+ initialValue += temp;
} // end of k loop
- return areaUnderCurve;
+ return initialValue;
}
} // namespace ocs2
diff --git a/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateConstraint.h b/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateConstraint.h
index 0c94db368..7e96eb387 100644
--- a/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateConstraint.h
+++ b/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateConstraint.h
@@ -47,9 +47,10 @@ class LoopshapingStateConstraint final : public StateConstraintCollection {
LoopshapingStateConstraint* clone() const override { return new LoopshapingStateConstraint(*this); }
- vector_t getValue(scalar_t time, const vector_t& state, const PreComputation& preComp) const override;
+ vector_array_t getValue(scalar_t time, const vector_t& state, const PreComputation& preComp) const override;
VectorFunctionLinearApproximation getLinearApproximation(scalar_t time, const vector_t& state,
const PreComputation& preComp) const override;
+
VectorFunctionQuadraticApproximation getQuadraticApproximation(scalar_t time, const vector_t& state,
const PreComputation& preComp) const override;
diff --git a/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateInputConstraint.h b/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateInputConstraint.h
index 7314c1182..14830ed3c 100644
--- a/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateInputConstraint.h
+++ b/ocs2_core/include/ocs2_core/loopshaping/constraint/LoopshapingStateInputConstraint.h
@@ -41,7 +41,7 @@ class LoopshapingStateInputConstraint : public StateInputConstraintCollection {
public:
~LoopshapingStateInputConstraint() override = default;
- vector_t getValue(scalar_t time, const vector_t& state, const vector_t& input, const PreComputation& preComp) const override;
+ vector_array_t getValue(scalar_t time, const vector_t& state, const vector_t& input, const PreComputation& preComp) const override;
protected:
LoopshapingStateInputConstraint(const StateInputConstraintCollection& systemConstraint,
diff --git a/ocs2_core/include/ocs2_core/misc/LinearAlgebra.h b/ocs2_core/include/ocs2_core/misc/LinearAlgebra.h
index b16175052..29c23a8d4 100644
--- a/ocs2_core/include/ocs2_core/misc/LinearAlgebra.h
+++ b/ocs2_core/include/ocs2_core/misc/LinearAlgebra.h
@@ -38,41 +38,18 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
namespace ocs2 {
namespace LinearAlgebra {
-// forward declarations
-void makePsdEigenvalue(matrix_t& squareMatrix, scalar_t minEigenvalue);
-
-void makePsdCholesky(matrix_t& A, scalar_t minEigenvalue);
-
-void computeConstraintProjection(const matrix_t& Dm, const matrix_t& RmInvUmUmT, matrix_t& DmDagger, matrix_t& DmDaggerTRmDmDaggerUUT,
- matrix_t& RmInvConstrainedUUT);
-
/**
* Set the eigenvalues of a triangular matrix to a minimum magnitude (maintaining the sign).
*/
-inline void setTriangularMinimumEigenvalues(matrix_t& Lr, scalar_t minEigenValue = numeric_traits::weakEpsilon()) {
- for (Eigen::Index i = 0; i < Lr.rows(); ++i) {
- scalar_t& eigenValue = Lr(i, i); // diagonal element is the eigenvalue
- if (eigenValue < 0.0) {
- eigenValue = std::min(-minEigenValue, eigenValue);
- } else {
- eigenValue = std::max(minEigenValue, eigenValue);
- }
- }
-}
+void setTriangularMinimumEigenvalues(matrix_t& Lr, scalar_t minEigenValue = numeric_traits::weakEpsilon());
/**
* Makes the input matrix PSD using a eigenvalue decomposition.
*
- * @tparam Derived type.
- * @param squareMatrix: The matrix to become PSD.
+ * @param [in, out] squareMatrix: The matrix to become PSD.
* @param [in] minEigenvalue: minimum eigenvalue.
*/
-template
-void makePsdEigenvalue(Eigen::MatrixBase& squareMatrix, scalar_t minEigenvalue = numeric_traits::limitEpsilon()) {
- matrix_t mat = squareMatrix;
- makePsdEigenvalue(mat, minEigenvalue);
- squareMatrix = mat;
-}
+void makePsdEigenvalue(matrix_t& squareMatrix, scalar_t minEigenvalue = numeric_traits::limitEpsilon());
/**
* Makes the input matrix PSD based on Gershgorin circle theorem. If the input matrix is positive definite and diagonally dominant,
@@ -89,20 +66,10 @@ void makePsdEigenvalue(Eigen::MatrixBase& squareMatrix, scalar_t minEig
* (1) Aii < minEigenvalue + Ri ==> minEigenvalue < lambda < minEigenvalue + 2 Ri
* (2) Aii > minEigenvalue + Ri ==> minEigenvalue < Aii - Ri < lambda < Aii + Ri
*
- * @tparam Derived type.
- * @param squareMatrix: The matrix to become PSD.
+ * @param [in, out] squareMatrix: The matrix to become PSD.
* @param [in] minEigenvalue: minimum eigenvalue.
*/
-template
-void makePsdGershgorin(Eigen::MatrixBase& squareMatrix, scalar_t minEigenvalue = numeric_traits::limitEpsilon()) {
- assert(squareMatrix.rows() == squareMatrix.cols());
- squareMatrix = 0.5 * (squareMatrix + squareMatrix.transpose()).eval();
- for (size_t i = 0; i < squareMatrix.rows(); i++) {
- // Gershgorin radius: since the matrix is symmetric we use column sum instead of row sum
- auto Ri = squareMatrix.col(i).cwiseAbs().sum() - std::abs(squareMatrix(i, i));
- squareMatrix(i, i) = std::max(squareMatrix(i, i), Ri + minEigenvalue);
- }
-}
+void makePsdGershgorin(matrix_t& squareMatrix, scalar_t minEigenvalue = numeric_traits::limitEpsilon());
/**
* Makes the input matrix PSD based on modified Cholesky decomposition.
@@ -116,32 +83,19 @@ void makePsdGershgorin(Eigen::MatrixBase& squareMatrix, scalar_t minEig
* References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with Limited memory, SIAM J. Sci. Comput.
* 21(1), pp. 24-45, 1999
*
- * @tparam Derived type.
- * @param A: The matrix to become PSD.
+ * @param [in, out] A: The matrix to become PSD.
* @param [in] minEigenvalue: minimum eigenvalue.
*/
-template
-void makePsdCholesky(Eigen::MatrixBase& A, scalar_t minEigenvalue = numeric_traits::limitEpsilon()) {
- matrix_t mat = A;
- makePsdCholesky(mat, minEigenvalue);
- A = mat;
-}
+void makePsdCholesky(matrix_t& A, scalar_t minEigenvalue = numeric_traits::limitEpsilon());
/**
* Computes the U*U^T decomposition associated to the inverse of the input matrix, where U is an upper triangular
* matrix. Note that the U*U^T decomposition is different from the Cholesky decomposition (U^T*U).
*
- * @tparam Derived type.
* @param [in] Am: A symmetric square positive definite matrix
* @param [out] AmInvUmUmT: The upper-triangular matrix associated to the UUT decomposition of inv(Am) matrix.
*/
-template
-void computeInverseMatrixUUT(const Derived& Am, Derived& AmInvUmUmT) {
- // Am = Lm Lm^T --> inv(Am) = inv(Lm^T) inv(Lm) where Lm^T is upper triangular
- Eigen::LLT lltOfA(Am);
- AmInvUmUmT.setIdentity(Am.rows(), Am.cols()); // for dynamic size matrices
- lltOfA.matrixU().solveInPlace(AmInvUmUmT);
-}
+void computeInverseMatrixUUT(const matrix_t& Am, matrix_t& AmInvUmUmT);
/**
* Computes constraint projection for linear constraints C*x + D*u - e = 0, with the weighting inv(Rm)
@@ -154,11 +108,36 @@ void computeInverseMatrixUUT(const Derived& Am, Derived& AmInvUmUmT) {
* @param [out] RmInvConstrainedUUT: The VVT decomposition of (I-DmDagger*Dm) * inv(Rm) * (I-DmDagger*Dm)^T where V is of
* the dimension n_u*(n_u-n_c) with n_u = Rm.rows() and n_c = Dm.rows().
*/
-template
-void computeConstraintProjection(const matrix_t& Dm, const DerivedInputMatrix& RmInvUmUmT, matrix_t& DmDagger,
- matrix_t& DmDaggerTRmDmDaggerUUT, matrix_t& RmInvConstrainedUUT) {
- computeConstraintProjection(Dm, matrix_t(RmInvUmUmT), DmDagger, DmDaggerTRmDmDaggerUUT, RmInvConstrainedUUT);
-}
+void computeConstraintProjection(const matrix_t& Dm, const matrix_t& RmInvUmUmT, matrix_t& DmDagger, matrix_t& DmDaggerTRmDmDaggerUUT,
+ matrix_t& RmInvConstrainedUUT);
+
+/**
+ * Returns the linear projection
+ * u = Pu * \tilde{u} + Px * x + Pe
+ *
+ * s.t. C*x + D*u + e = 0 is satisfied for any \tilde{u}
+ *
+ * Implementation based on the QR decomposition
+ *
+ * @param [in] constraint : C = dfdx, D = dfdu, e = f;
+ * @return Projection terms Px = dfdx, Pu = dfdu, Pe = f (first) and left pseudo-inverse of D^T (second);
+ */
+std::pair qrConstraintProjection(const VectorFunctionLinearApproximation& constraint);
+
+/**
+ * Returns the linear projection
+ * u = Pu * \tilde{u} + Px * x + Pe
+ *
+ * s.t. C*x + D*u + e = 0 is satisfied for any \tilde{u}
+ *
+ * Implementation based on the LU decomposition
+ *
+ * @param [in] constraint : C = dfdx, D = dfdu, e = f;
+ * @param [in] extractPseudoInverse : If true, left pseudo-inverse of D^T is returned. If false, an empty matrix is returned;
+ * @return Projection terms Px = dfdx, Pu = dfdu, Pe = f (first) and left pseudo-inverse of D^T (second);
+ */
+std::pair luConstraintProjection(const VectorFunctionLinearApproximation& constraint,
+ bool extractPseudoInverse = false);
/** Computes the rank of a matrix */
template
diff --git a/ocs2_core/include/ocs2_core/misc/Numerics.h b/ocs2_core/include/ocs2_core/misc/Numerics.h
index 5003791ce..d2e45d2f4 100644
--- a/ocs2_core/include/ocs2_core/misc/Numerics.h
+++ b/ocs2_core/include/ocs2_core/misc/Numerics.h
@@ -42,22 +42,39 @@ namespace numerics {
*
* @tparam T1: data type of x.
* @tparam T2: data type of y.
- * @param [in] x: a floating-point number.
- * @param [in] y: a floating-point number.
+ * @tparam T3: data type of prec.
+ * @param [in] x: First floating-point number.
+ * @param [in] y: Second floating-point number.
+ * @param [in] prec: The comparison precision.
* @return bool: true if x=y.
*/
-template
-typename std::enable_if::type>::value ||
- std::is_floating_point::type>::value,
- bool>::type
-almost_eq(T1&& x, T2&& y) {
- using TypeResult = typename std::conditional::type>::value,
- typename std::remove_reference::type, typename std::remove_reference::type>::type;
+template
+bool almost_eq(T1&& x, T2&& y, T3&& prec) {
+ static_assert(std::is_floating_point::type>::value, "First argument is not floating point!");
+ static_assert(std::is_floating_point::type>::value, "Second argument is not floating point!");
+ static_assert(std::is_floating_point::type>::value, "prec is not floating point!");
// the machine epsilon has to be scaled to the magnitude of the values used
- // and multiplied by the desired precision in ULPs (units in the last place)
- return std::abs(x - y) <= std::numeric_limits::epsilon() * std::abs(x + y)
- // unless the result is subnormal
- || std::abs(x - y) < std::numeric_limits::min();
+ // and multiplied by the desired precision unless the result is subnormal
+ using Type = const std::remove_reference_t;
+ const auto absDiff = std::abs(x - static_cast(y));
+ const auto magnitude = std::min(std::abs(x), std::abs(static_cast(y)));
+ return absDiff <= static_cast(prec) * magnitude || absDiff < std::numeric_limits::min();
+}
+
+/**
+ * Almost equal which uses machine epsilon to compare floating-point values for equality.
+ * refer to: https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon
+ *
+ * @tparam T1: data type of x.
+ * @tparam T2: data type of y.
+ * @param [in] x: First floating-point number.
+ * @param [in] y: Second floating-point number.
+ * @return bool: true if x=y.
+ */
+template
+bool almost_eq(T1&& x, T2&& y) {
+ const auto prec = std::numeric_limits>::epsilon();
+ return almost_eq(x, y, prec);
}
/**
@@ -65,15 +82,28 @@ almost_eq(T1&& x, T2&& y) {
*
* @tparam T1: data type of x.
* @tparam T2: data type of y.
- * @param [in] x: a floating-point number.
- * @param [in] y: a floating-point number.
+ * @tparam T3: data type of prec.
+ * @param [in] x: First floating-point number.
+ * @param [in] y: Second floating-point number.
+ * @param [in] prec: The comparison precision.
* @return bool: true if x<=y.
*/
-template
-typename std::enable_if::type>::value ||
- std::is_floating_point::type>::value,
- bool>::type
-almost_le(T1&& x, T2&& y) {
+template
+bool almost_le(T1&& x, T2&& y, T3&& prec) {
+ return x < y || almost_eq(x, y, prec);
+}
+
+/**
+ * Almost less-equal which uses machine epsilon to compare floating-point values for equality.
+ *
+ * @tparam T1: data type of x.
+ * @tparam T2: data type of y.
+ * @param [in] x: First floating-point number.
+ * @param [in] y: Second floating-point number.
+ * @return bool: true if x<=y.
+ */
+template
+bool almost_le(T1&& x, T2&& y) {
return x < y || almost_eq(x, y);
}
@@ -82,15 +112,28 @@ almost_le(T1&& x, T2&& y) {
*
* @tparam T1: data type of x.
* @tparam T2: data type of y.
- * @param [in] x: a floating-point number.
- * @param [in] y: a floating-point number.
+ * @tparam T3: data type of prec.
+ * @param [in] x: First floating-point number.
+ * @param [in] y: Second floating-point number.
+ * @param [in] prec: The comparison precision.
+ * @return bool: true if x>=y.
+ */
+template
+bool almost_ge(T1&& x, T2&& y, T3&& prec) {
+ return x > y || almost_eq(x, y, prec);
+}
+
+/**
+ * Almost greater-equal which uses machine epsilon to compare floating-point values for equality.
+ *
+ * @tparam T1: data type of x.
+ * @tparam T2: data type of y.
+ * @param [in] x: First floating-point number.
+ * @param [in] y: Second floating-point number.
* @return bool: true if x>=y.
*/
template
-typename std::enable_if::type>::value ||
- std::is_floating_point::type>::value,
- bool>::type
-almost_ge(T1&& x, T2&& y) {
+bool almost_ge(T1&& x, T2&& y) {
return x > y || almost_eq(x, y);
}
diff --git a/ocs2_core/include/ocs2_core/model_data/Metrics.h b/ocs2_core/include/ocs2_core/model_data/Metrics.h
index 2d205d377..dc01f54d5 100644
--- a/ocs2_core/include/ocs2_core/model_data/Metrics.h
+++ b/ocs2_core/include/ocs2_core/model_data/Metrics.h
@@ -53,23 +53,33 @@ struct LagrangianMetricsConstRef {
};
/**
- * The collection of cost, equality constraints, and LagrangianMetrics structure for all possible constraint terms (handled by
- * Lagrangian method) in a particular time point.
- * cost : The total cost in a particular time point.
- * stateEqConstraint : A vector of all active state equality constraints.
- * stateInputEqConstraint : A vector of all active state-input equality constraints.
- * stateEqLagrangian : An array of state equality constraint terms handled by Lagrangian method.
- * stateIneqLagrangian : An array of state inequality constraint terms handled by Lagrangian method.
- * stateInputEqLagrangian : An array of state-input equality constraint terms handled by Lagrangian method.
- * stateInputIneqLagrangian : An array of state-input inequality constraint terms handled by Lagrangian method.
+ * The collection of cost, dynamics violation, constraints, and LagrangianMetrics structure for all possible constraint
+ * terms (handled by Lagrangian method) in a point of time.
+ * cost : The total cost in a particular time point.
+ * dynamicsViolation : The vector of dynamics violation.
+ * stateEqConstraint : An array of all state equality constraints.
+ * stateInputEqConstraint : An array of all state-input equality constraints.
+ * stateIneqConstraint : An array of all state inequality constraints.
+ * stateInputIneqConstraint : An array of all state-input inequality constraints.
+ * stateEqLagrangian : An array of state equality constraint terms handled by Lagrangian method.
+ * stateIneqLagrangian : An array of state inequality constraint terms handled by Lagrangian method.
+ * stateInputEqLagrangian : An array of state-input equality constraint terms handled by Lagrangian method.
+ * stateInputIneqLagrangian : An array of state-input inequality constraint terms handled by Lagrangian method.
*/
-struct MetricsCollection {
+struct Metrics {
// Cost
scalar_t cost;
+ // Dynamics violation
+ vector_t dynamicsViolation;
+
// Equality constraints
- vector_t stateEqConstraint;
- vector_t stateInputEqConstraint;
+ vector_array_t stateEqConstraint;
+ vector_array_t stateInputEqConstraint;
+
+ // Inequality constraints
+ vector_array_t stateIneqConstraint;
+ vector_array_t stateInputIneqConstraint;
// Lagrangians
std::vector stateEqLagrangian;
@@ -77,49 +87,62 @@ struct MetricsCollection {
std::vector stateInputEqLagrangian;
std::vector stateInputIneqLagrangian;
- /** Exchanges the values of MetricsCollection */
- void swap(MetricsCollection& other) {
- // Cost
- std::swap(cost, other.cost);
- // Equality constraints
- stateEqConstraint.swap(other.stateEqConstraint);
- stateInputEqConstraint.swap(other.stateInputEqConstraint);
- // Lagrangians
- stateEqLagrangian.swap(other.stateEqLagrangian);
- stateIneqLagrangian.swap(other.stateIneqLagrangian);
- stateInputEqLagrangian.swap(other.stateInputEqLagrangian);
- stateInputIneqLagrangian.swap(other.stateInputIneqLagrangian);
- }
+ /** Exchanges the values of Metrics */
+ void swap(Metrics& other);
- /** Clears the value of the MetricsCollection */
- void clear() {
- // Cost
- cost = 0.0;
- // Equality constraints
- stateEqConstraint = vector_t();
- stateInputEqConstraint = vector_t();
- // Lagrangians
- stateEqLagrangian.clear();
- stateIneqLagrangian.clear();
- stateInputEqLagrangian.clear();
- stateInputIneqLagrangian.clear();
- }
+ /** Clears the value of the Metrics */
+ void clear();
+
+ /** Returns true if *this is approximately equal to other, within the precision determined by prec. */
+ bool isApprox(const Metrics& other, scalar_t prec = 1e-8) const;
};
-/** Sums penalties of an array of LagrangianMetrics */
+/** Sums penalties of an array of LagrangianMetrics. */
inline scalar_t sumPenalties(const std::vector& metricsArray) {
scalar_t s = 0.0;
std::for_each(metricsArray.begin(), metricsArray.end(), [&s](const LagrangianMetrics& m) { s += m.penalty; });
return s;
}
-/** Computes the sum of squared norm of constraints of an array of LagrangianMetrics */
+/** Computes the sum of squared norm of constraints of an array of LagrangianMetrics. */
inline scalar_t constraintsSquaredNorm(const std::vector& metricsArray) {
scalar_t s = 0.0;
std::for_each(metricsArray.begin(), metricsArray.end(), [&s](const LagrangianMetrics& m) { s += m.constraint.squaredNorm(); });
return s;
}
+/** Computes the sum of squared norm of a vector of equality constraints violation. */
+inline scalar_t getEqConstraintsSSE(const vector_t& eqConstraint) {
+ if (eqConstraint.size() == 0) {
+ return 0.0;
+ } else {
+ return eqConstraint.squaredNorm();
+ }
+}
+
+/** Computes the sum of squared norm of a vector of inequality constraints violation. */
+inline scalar_t getIneqConstraintsSSE(const vector_t& ineqConstraint) {
+ if (ineqConstraint.size() == 0) {
+ return 0.0;
+ } else {
+ return ineqConstraint.cwiseMin(0.0).squaredNorm();
+ }
+}
+
+/** Computes the sum of squared norm of an array of equality constraints violation. */
+inline scalar_t getEqConstraintsSSE(const vector_array_t& eqConstraint) {
+ scalar_t s = 0.0;
+ std::for_each(eqConstraint.begin(), eqConstraint.end(), [&s](const vector_t& v) { s += getEqConstraintsSSE(v); });
+ return s;
+}
+
+/** Computes the sum of squared norm of an array of inequality constraints violation. */
+inline scalar_t getIneqConstraintsSSE(const vector_array_t& ineqConstraint) {
+ scalar_t s = 0.0;
+ std::for_each(ineqConstraint.begin(), ineqConstraint.end(), [&s](const vector_t& v) { s += getIneqConstraintsSSE(v); });
+ return s;
+}
+
/**
* Serializes an array of LagrangianMetrics structures associated to an array of constraint terms.
*
@@ -128,14 +151,40 @@ inline scalar_t constraintsSquaredNorm(const std::vector& met
*/
vector_t toVector(const std::vector& termsLagrangianMetrics);
+/**
+ * Serializes an array of constraint terms.
+ *
+ * @ param [in] constraintArray : An array of constraint terms.
+ * @return Serialized vector of the format : (..., constraintArray[i], ...).
+ */
+vector_t toVector(const vector_array_t& constraintArray);
+
/**
* Gets the size of constraint terms.
*
+ * @ param [in] constraintArray : An array of constraint terms.
+ * @return An array of constraint terms size. It has the same size as the input array.
+ */
+size_array_t getSizes(const vector_array_t& constraintArray);
+
+/**
+ * Gets the size of constraint Lagrangian terms.
+ *
* @ param [in] termsLagrangianMetrics : LagrangianMetrics associated to an array of constraint terms.
* @return An array of constraint terms size. It has the same size as the input array.
*/
size_array_t getSizes(const std::vector& termsLagrangianMetrics);
+/**
+ * Deserializes the vector to an array of constraint terms.
+ *
+ * @param [in] termsSize : An array of constraint terms size. It as the same size as the output array.
+ * @param [in] vec : Serialized array of constraint terms of the format :
+ * (..., constraintArray[i], ...)
+ * @return An array of constraint terms.
+ */
+vector_array_t toConstraintArray(const size_array_t& termsSize, const vector_t& vec);
+
/**
* Deserializes the vector to an array of LagrangianMetrics structures based on size of constraint terms.
*
@@ -161,13 +210,13 @@ namespace LinearInterpolation {
LagrangianMetrics interpolate(const index_alpha_t& indexAlpha, const std::vector& dataArray);
/**
- * Linearly interpolates a trajectory of MetricsCollection.
+ * Linearly interpolates a trajectory of Metrics.
*
* @param [in] indexAlpha : index and interpolation coefficient (alpha) pair.
- * @param [in] dataArray : A trajectory of MetricsCollection.
- * @return The interpolated MetricsCollection at indexAlpha.
+ * @param [in] dataArray : A trajectory of Metrics.
+ * @return The interpolated Metrics at indexAlpha.
*/
-MetricsCollection interpolate(const index_alpha_t& indexAlpha, const std::vector& dataArray);
+Metrics interpolate(const index_alpha_t& indexAlpha, const std::vector& dataArray);
} // namespace LinearInterpolation
} // namespace ocs2
diff --git a/ocs2_core/include/ocs2_core/penalties/augmented/ModifiedRelaxedBarrierPenalty.h b/ocs2_core/include/ocs2_core/penalties/augmented/ModifiedRelaxedBarrierPenalty.h
index 239e20b53..a5bb74645 100644
--- a/ocs2_core/include/ocs2_core/penalties/augmented/ModifiedRelaxedBarrierPenalty.h
+++ b/ocs2_core/include/ocs2_core/penalties/augmented/ModifiedRelaxedBarrierPenalty.h
@@ -76,7 +76,7 @@ class ModifiedRelaxedBarrierPenalty final : public AugmentedPenaltyBase {
/** Factory function */
static std::unique_ptr create(Config config) {
- return std::unique_ptr(new ModifiedRelaxedBarrierPenalty(std::move(config)));
+ return std::make_unique(std::move(config));
}
~ModifiedRelaxedBarrierPenalty() override = default;
diff --git a/ocs2_core/include/ocs2_core/penalties/augmented/QuadraticPenalty.h b/ocs2_core/include/ocs2_core/penalties/augmented/QuadraticPenalty.h
index 317516ab5..7a8b0578c 100644
--- a/ocs2_core/include/ocs2_core/penalties/augmented/QuadraticPenalty.h
+++ b/ocs2_core/include/ocs2_core/penalties/augmented/QuadraticPenalty.h
@@ -72,9 +72,7 @@ class QuadraticPenalty final : public AugmentedPenaltyBase {
explicit QuadraticPenalty(Config config) : config_(std::move(config)) {}
/** Factory function */
- static std::unique_ptr create(Config config) {
- return std::unique_ptr(new QuadraticPenalty(std::move(config)));
- }
+ static std::unique_ptr create(Config config) { return std::make_unique(std::move(config)); }
~QuadraticPenalty() override = default;
QuadraticPenalty* clone() const override { return new QuadraticPenalty(*this); }
diff --git a/ocs2_core/include/ocs2_core/penalties/augmented/SlacknessSquaredHingePenalty.h b/ocs2_core/include/ocs2_core/penalties/augmented/SlacknessSquaredHingePenalty.h
index f4ca9da26..09cadba0f 100644
--- a/ocs2_core/include/ocs2_core/penalties/augmented/SlacknessSquaredHingePenalty.h
+++ b/ocs2_core/include/ocs2_core/penalties/augmented/SlacknessSquaredHingePenalty.h
@@ -75,7 +75,7 @@ class SlacknessSquaredHingePenalty final : public AugmentedPenaltyBase {
/** Factory function */
static std::unique_ptr create(Config config) {
- return std::unique_ptr(new SlacknessSquaredHingePenalty(std::move(config)));
+ return std::make_unique(std::move(config));
}
~SlacknessSquaredHingePenalty() override = default;
diff --git a/ocs2_core/include/ocs2_core/penalties/augmented/SmoothAbsolutePenalty.h b/ocs2_core/include/ocs2_core/penalties/augmented/SmoothAbsolutePenalty.h
index 836ba28b5..d6f1e6b0f 100644
--- a/ocs2_core/include/ocs2_core/penalties/augmented/SmoothAbsolutePenalty.h
+++ b/ocs2_core/include/ocs2_core/penalties/augmented/SmoothAbsolutePenalty.h
@@ -76,9 +76,7 @@ class SmoothAbsolutePenalty final : public AugmentedPenaltyBase {
explicit SmoothAbsolutePenalty(Config config) : config_(std::move(config)) {}
/** Factory function */
- static std::unique_ptr create(Config config) {
- return std::unique_ptr(new SmoothAbsolutePenalty(std::move(config)));
- }
+ static std::unique_ptr create(Config config) { return std::make_unique(std::move(config)); }
~SmoothAbsolutePenalty() override = default;
SmoothAbsolutePenalty* clone() const override { return new SmoothAbsolutePenalty(*this); }
diff --git a/ocs2_core/include/ocs2_core/thread_support/BufferedValue.h b/ocs2_core/include/ocs2_core/thread_support/BufferedValue.h
index 9fc28dca7..8c1464b78 100644
--- a/ocs2_core/include/ocs2_core/thread_support/BufferedValue.h
+++ b/ocs2_core/include/ocs2_core/thread_support/BufferedValue.h
@@ -58,10 +58,10 @@ class BufferedValue {
T& get() { return activeValue_; }
/** Copy a new value into the buffer. */
- void setBuffer(const T& value) { buffer_.reset(std::unique_ptr(new T(value))); }
+ void setBuffer(const T& value) { buffer_.reset(std::make_unique(value)); }
/** Move a new value into the buffer. */
- void setBuffer(T&& value) { buffer_.reset(std::unique_ptr(new T(std::move(value)))); }
+ void setBuffer(T&& value) { buffer_.reset(std::make_unique(std::move(value))); }
/**
* Replaces the active value with the value in the buffer.
@@ -87,4 +87,4 @@ class BufferedValue {
Synchronized buffer_;
};
-} // namespace ocs2
\ No newline at end of file
+} // namespace ocs2
diff --git a/ocs2_core/include/ocs2_core/thread_support/ThreadPool.h b/ocs2_core/include/ocs2_core/thread_support/ThreadPool.h
index 1ae4324dd..106e9fc21 100644
--- a/ocs2_core/include/ocs2_core/thread_support/ThreadPool.h
+++ b/ocs2_core/include/ocs2_core/thread_support/ThreadPool.h
@@ -141,7 +141,7 @@ struct ThreadPool::Task final : public ThreadPool::TaskBase {
/**************************************************************************************************/
template
std::future::type> ThreadPool::run(Functor taskFunction) {
- std::unique_ptr> taskPtr(new Task(std::move(taskFunction)));
+ auto taskPtr = std::make_unique>(std::move(taskFunction));
auto future = taskPtr->packagedTask.get_future();
if (workerThreads_.empty()) {
diff --git a/ocs2_core/src/Types.cpp b/ocs2_core/src/Types.cpp
index c6d4e4583..3645fa031 100644
--- a/ocs2_core/src/Types.cpp
+++ b/ocs2_core/src/Types.cpp
@@ -446,6 +446,8 @@ VectorFunctionQuadraticApproximation& VectorFunctionQuadraticApproximation::setZ
}
} else {
dfdu = matrix_t();
+ dfdux.resize(nv);
+ dfduu.resize(nv);
for (size_t i = 0; i < nv; i++) {
dfdux[i] = matrix_t();
dfduu[i] = matrix_t();
diff --git a/ocs2_core/src/augmented_lagrangian/AugmentedLagrangian.cpp b/ocs2_core/src/augmented_lagrangian/AugmentedLagrangian.cpp
index 01db3586e..d24cb8746 100644
--- a/ocs2_core/src/augmented_lagrangian/AugmentedLagrangian.cpp
+++ b/ocs2_core/src/augmented_lagrangian/AugmentedLagrangian.cpp
@@ -36,7 +36,7 @@ namespace ocs2 {
/******************************************************************************************************/
std::unique_ptr create(std::unique_ptr constraintPtr,
std::vector> penaltyPtrArray) {
- return std::unique_ptr(new StateAugmentedLagrangian(std::move(constraintPtr), std::move(penaltyPtrArray)));
+ return std::make_unique(std::move(constraintPtr), std::move(penaltyPtrArray));
}
/******************************************************************************************************/
@@ -44,7 +44,7 @@ std::unique_ptr create(std::unique_ptr create(std::unique_ptr constraintPtr,
std::unique_ptr penaltyPtr) {
- return std::unique_ptr(new StateAugmentedLagrangian(std::move(constraintPtr), std::move(penaltyPtr)));
+ return std::make_unique(std::move(constraintPtr), std::move(penaltyPtr));
}
/******************************************************************************************************/
@@ -52,8 +52,7 @@ std::unique_ptr create(std::unique_ptr create(std::unique_ptr constraintPtr,
std::vector> penaltyPtrArray) {
- return std::unique_ptr(
- new StateInputAugmentedLagrangian(std::move(constraintPtr), std::move(penaltyPtrArray)));
+ return std::make_unique(std::move(constraintPtr), std::move(penaltyPtrArray));
}
/******************************************************************************************************/
@@ -61,7 +60,7 @@ std::unique_ptr create(std::unique_ptr create(std::unique_ptr constraintPtr,
std::unique_ptr penaltyPtr) {
- return std::unique_ptr(new StateInputAugmentedLagrangian(std::move(constraintPtr), std::move(penaltyPtr)));
+ return std::make_unique(std::move(constraintPtr), std::move(penaltyPtr));
}
} // namespace ocs2
diff --git a/ocs2_core/src/augmented_lagrangian/StateAugmentedLagrangianCollection.cpp b/ocs2_core/src/augmented_lagrangian/StateAugmentedLagrangianCollection.cpp
index 4189cffcf..b002109f9 100644
--- a/ocs2_core/src/augmented_lagrangian/StateAugmentedLagrangianCollection.cpp
+++ b/ocs2_core/src/augmented_lagrangian/StateAugmentedLagrangianCollection.cpp
@@ -90,10 +90,18 @@ ScalarFunctionQuadraticApproximation StateAugmentedLagrangianCollection::getQuad
// accumulate terms
for (size_t i = firstActiveInd + 1; i < terms_.size(); i++) {
if (terms_[i]->isActive(time)) {
- penalty += terms_[i]->getQuadraticApproximation(time, state, termsMultiplier[i], preComp);
+ const auto termPenalty = terms_[i]->getQuadraticApproximation(time, state, termsMultiplier[i], preComp);
+ penalty.f += termPenalty.f;
+ penalty.dfdx += termPenalty.dfdx;
+ penalty.dfdxx += termPenalty.dfdxx;
}
}
+ // make sure that input derivatives are empty
+ penalty.dfdu = vector_t();
+ penalty.dfduu = matrix_t();
+ penalty.dfdux = matrix_t();
+
return penalty;
}
diff --git a/ocs2_core/src/constraint/StateConstraintCollection.cpp b/ocs2_core/src/constraint/StateConstraintCollection.cpp
index 77086ea07..94dc0634a 100644
--- a/ocs2_core/src/constraint/StateConstraintCollection.cpp
+++ b/ocs2_core/src/constraint/StateConstraintCollection.cpp
@@ -48,34 +48,37 @@ StateConstraintCollection* StateConstraintCollection::clone() const {
/******************************************************************************************************/
size_t StateConstraintCollection::getNumConstraints(scalar_t time) const {
size_t numConstraints = 0;
-
- // accumulate number of constraints for each constraintTerm
for (const auto& constraintTerm : this->terms_) {
if (constraintTerm->isActive(time)) {
numConstraints += constraintTerm->getNumConstraints(time);
}
}
-
return numConstraints;
}
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
-vector_t StateConstraintCollection::getValue(scalar_t time, const vector_t& state, const PreComputation& preComp) const {
- vector_t constraintValues;
- constraintValues.resize(getNumConstraints(time));
-
- // append vectors of constraint values from each constraintTerm
- size_t i = 0;
- for (const auto& constraintTerm : this->terms_) {
- if (constraintTerm->isActive(time)) {
- const auto constraintTermValues = constraintTerm->getValue(time, state, preComp);
- constraintValues.segment(i, constraintTermValues.rows()) = constraintTermValues;
- i += constraintTermValues.rows();
+size_array_t StateConstraintCollection::getTermsSize(scalar_t time) const {
+ size_array_t termsSize(this->terms_.size(), 0);
+ for (size_t i = 0; i < this->terms_.size(); ++i) {
+ if (this->terms_[i]->isActive(time)) {
+ termsSize[i] = this->terms_[i]->getNumConstraints(time);
}
}
+ return termsSize;
+}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+vector_array_t StateConstraintCollection::getValue(scalar_t time, const vector_t& state, const PreComputation& preComp) const {
+ vector_array_t constraintValues(this->terms_.size());
+ for (size_t i = 0; i < this->terms_.size(); ++i) {
+ if (this->terms_[i]->isActive(time)) {
+ constraintValues[i] = this->terms_[i]->getValue(time, state, preComp);
+ }
+ } // end of i loop
return constraintValues;
}
diff --git a/ocs2_core/src/constraint/StateInputConstraintCollection.cpp b/ocs2_core/src/constraint/StateInputConstraintCollection.cpp
index 44264a06b..ae3daa870 100644
--- a/ocs2_core/src/constraint/StateInputConstraintCollection.cpp
+++ b/ocs2_core/src/constraint/StateInputConstraintCollection.cpp
@@ -49,35 +49,38 @@ StateInputConstraintCollection* StateInputConstraintCollection::clone() const {
/******************************************************************************************************/
size_t StateInputConstraintCollection::getNumConstraints(scalar_t time) const {
size_t numConstraints = 0;
-
- // accumulate number of constraints for each constraintTerm
for (const auto& constraintTerm : this->terms_) {
if (constraintTerm->isActive(time)) {
numConstraints += constraintTerm->getNumConstraints(time);
}
}
-
return numConstraints;
}
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
-vector_t StateInputConstraintCollection::getValue(scalar_t time, const vector_t& state, const vector_t& input,
- const PreComputation& preComp) const {
- vector_t constraintValues;
- constraintValues.resize(getNumConstraints(time));
-
- // append vectors of constraint values from each constraintTerm
- size_t i = 0;
- for (const auto& constraintTerm : this->terms_) {
- if (constraintTerm->isActive(time)) {
- const auto constraintTermValues = constraintTerm->getValue(time, state, input, preComp);
- constraintValues.segment(i, constraintTermValues.rows()) = constraintTermValues;
- i += constraintTermValues.rows();
+size_array_t StateInputConstraintCollection::getTermsSize(scalar_t time) const {
+ size_array_t termsSize(this->terms_.size(), 0);
+ for (size_t i = 0; i < this->terms_.size(); ++i) {
+ if (this->terms_[i]->isActive(time)) {
+ termsSize[i] = this->terms_[i]->getNumConstraints(time);
}
}
+ return termsSize;
+}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+vector_array_t StateInputConstraintCollection::getValue(scalar_t time, const vector_t& state, const vector_t& input,
+ const PreComputation& preComp) const {
+ vector_array_t constraintValues(this->terms_.size());
+ for (size_t i = 0; i < this->terms_.size(); ++i) {
+ if (this->terms_[i]->isActive(time)) {
+ constraintValues[i] = this->terms_[i]->getValue(time, state, input, preComp);
+ }
+ } // end of i loop
return constraintValues;
}
diff --git a/ocs2_core/src/integration/Integrator.cpp b/ocs2_core/src/integration/Integrator.cpp
index 01efb48d3..e7f1ab696 100644
--- a/ocs2_core/src/integration/Integrator.cpp
+++ b/ocs2_core/src/integration/Integrator.cpp
@@ -80,24 +80,24 @@ IntegratorType fromString(const std::string& name) {
std::unique_ptr newIntegrator(IntegratorType integratorType, const std::shared_ptr& eventHandlerPtr) {
switch (integratorType) {
case (IntegratorType::EULER):
- return std::unique_ptr(new IntegratorEuler(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
case (IntegratorType::ODE45):
- return std::unique_ptr(new ODE45(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
case (IntegratorType::ODE45_OCS2):
- return std::unique_ptr(new RungeKuttaDormandPrince5(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
case (IntegratorType::ADAMS_BASHFORTH):
- return std::unique_ptr(new IntegratorAdamsBashforth<1>(eventHandlerPtr));
+ return std::make_unique>(eventHandlerPtr);
case (IntegratorType::BULIRSCH_STOER):
- return std::unique_ptr(new IntegratorBulirschStoer(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
case (IntegratorType::MODIFIED_MIDPOINT):
- return std::unique_ptr(new IntegratorModifiedMidpoint(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
case (IntegratorType::RK4):
- return std::unique_ptr(new IntegratorRK4(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
case (IntegratorType::RK5_VARIABLE):
- return std::unique_ptr(new IntegratorRK5Variable(eventHandlerPtr));
+ return std::make_unique(eventHandlerPtr);
#if (BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 > 55)
case (IntegratorType::ADAMS_BASHFORTH_MOULTON):
- return std::unique_ptr(new IntegratorAdamsBashforthMoulton<1>(eventHandlerPtr));
+ return std::make_unique>(eventHandlerPtr);
#endif
default:
throw std::runtime_error("Integrator of type " + integrator_type::toString(integratorType) + " not supported.");
diff --git a/ocs2_core/src/loopshaping/augmented_lagrangian/LoopshapingAugmentedLagrangian.cpp b/ocs2_core/src/loopshaping/augmented_lagrangian/LoopshapingAugmentedLagrangian.cpp
index 223f8d7c3..9ec9286f5 100644
--- a/ocs2_core/src/loopshaping/augmented_lagrangian/LoopshapingAugmentedLagrangian.cpp
+++ b/ocs2_core/src/loopshaping/augmented_lagrangian/LoopshapingAugmentedLagrangian.cpp
@@ -41,8 +41,7 @@ namespace LoopshapingAugmentedLagrangian {
/******************************************************************************************************/
std::unique_ptr create(const StateAugmentedLagrangianCollection& systemAugmentedLagrangian,
std::shared_ptr loopshapingDefinition) {
- return std::unique_ptr(
- new LoopshapingStateAugmentedLagrangian(systemAugmentedLagrangian, std::move(loopshapingDefinition)));
+ return std::make_unique(systemAugmentedLagrangian, std::move(loopshapingDefinition));
}
/******************************************************************************************************/
@@ -52,11 +51,9 @@ std::unique_ptr create(const StateInput
std::shared_ptr loopshapingDefinition) {
switch (loopshapingDefinition->getType()) {
case LoopshapingType::outputpattern:
- return std::unique_ptr(
- new LoopshapingAugmentedLagrangianOutputPattern(systemAugmentedLagrangian, std::move(loopshapingDefinition)));
+ return std::make_unique(systemAugmentedLagrangian, std::move(loopshapingDefinition));
case LoopshapingType::eliminatepattern:
- return std::unique_ptr(
- new LoopshapingAugmentedLagrangianEliminatePattern(systemAugmentedLagrangian, std::move(loopshapingDefinition)));
+ return std::make_unique(systemAugmentedLagrangian, std::move(loopshapingDefinition));
default:
throw std::runtime_error("[LoopshapingAugmentedLagrangian::create] invalid loopshaping type");
}
diff --git a/ocs2_core/src/loopshaping/constraint/LoopshapingConstraint.cpp b/ocs2_core/src/loopshaping/constraint/LoopshapingConstraint.cpp
index 09dc48011..e99e8dc48 100644
--- a/ocs2_core/src/loopshaping/constraint/LoopshapingConstraint.cpp
+++ b/ocs2_core/src/loopshaping/constraint/LoopshapingConstraint.cpp
@@ -41,7 +41,7 @@ namespace LoopshapingConstraint {
/******************************************************************************************************/
std::unique_ptr create(const StateConstraintCollection& systemConstraint,
std::shared_ptr loopshapingDefinition) {
- return std::unique_ptr(new LoopshapingStateConstraint(systemConstraint, loopshapingDefinition));
+ return std::make_unique(systemConstraint, std::move(loopshapingDefinition));
}
/******************************************************************************************************/
@@ -51,11 +51,9 @@ std::unique_ptr create(const StateInputConstrain
std::shared_ptr loopshapingDefinition) {
switch (loopshapingDefinition->getType()) {
case LoopshapingType::outputpattern:
- return std::unique_ptr(
- new LoopshapingConstraintOutputPattern(systemConstraint, std::move(loopshapingDefinition)));
+ return std::make_unique(systemConstraint, std::move(loopshapingDefinition));
case LoopshapingType::eliminatepattern:
- return std::unique_ptr(
- new LoopshapingConstraintEliminatePattern(systemConstraint, std::move(loopshapingDefinition)));
+ return std::make_unique(systemConstraint, std::move(loopshapingDefinition));
default:
throw std::runtime_error("[LoopshapingConstraint::create] invalid loopshaping type");
}
diff --git a/ocs2_core/src/loopshaping/constraint/LoopshapingStateConstraint.cpp b/ocs2_core/src/loopshaping/constraint/LoopshapingStateConstraint.cpp
index 1df9ded24..af67c3761 100644
--- a/ocs2_core/src/loopshaping/constraint/LoopshapingStateConstraint.cpp
+++ b/ocs2_core/src/loopshaping/constraint/LoopshapingStateConstraint.cpp
@@ -38,9 +38,9 @@ namespace ocs2 {
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
-vector_t LoopshapingStateConstraint::getValue(scalar_t t, const vector_t& x, const PreComputation& preComp) const {
+vector_array_t LoopshapingStateConstraint::getValue(scalar_t t, const vector_t& x, const PreComputation& preComp) const {
if (this->empty()) {
- return vector_t::Zero(0);
+ return vector_array_t();
}
const LoopshapingPreComputation& preCompLS = cast(preComp);
diff --git a/ocs2_core/src/loopshaping/constraint/LoopshapingStateInputConstraint.cpp b/ocs2_core/src/loopshaping/constraint/LoopshapingStateInputConstraint.cpp
index 14125e7f2..9ef11c9b8 100644
--- a/ocs2_core/src/loopshaping/constraint/LoopshapingStateInputConstraint.cpp
+++ b/ocs2_core/src/loopshaping/constraint/LoopshapingStateInputConstraint.cpp
@@ -38,9 +38,10 @@ namespace ocs2 {
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
-vector_t LoopshapingStateInputConstraint::getValue(scalar_t t, const vector_t& x, const vector_t& u, const PreComputation& preComp) const {
+vector_array_t LoopshapingStateInputConstraint::getValue(scalar_t t, const vector_t& x, const vector_t& u,
+ const PreComputation& preComp) const {
if (this->empty()) {
- return vector_t::Zero(0);
+ return vector_array_t();
}
const LoopshapingPreComputation& preCompLS = cast(preComp);
diff --git a/ocs2_core/src/loopshaping/cost/LoopshapingCost.cpp b/ocs2_core/src/loopshaping/cost/LoopshapingCost.cpp
index 44a10dc89..df5e301b7 100644
--- a/ocs2_core/src/loopshaping/cost/LoopshapingCost.cpp
+++ b/ocs2_core/src/loopshaping/cost/LoopshapingCost.cpp
@@ -41,7 +41,7 @@ namespace LoopshapingCost {
/******************************************************************************************************/
std::unique_ptr create(const StateCostCollection& systemCost,
std::shared_ptr loopshapingDefinition) {
- return std::unique_ptr(new LoopshapingStateCost(systemCost, loopshapingDefinition));
+ return std::make_unique(systemCost, std::move(loopshapingDefinition));
}
/******************************************************************************************************/
@@ -51,9 +51,9 @@ std::unique_ptr create(const StateInputCostCollection&
std::shared_ptr loopshapingDefinition) {
switch (loopshapingDefinition->getType()) {
case LoopshapingType::outputpattern:
- return std::unique_ptr(new LoopshapingCostOutputPattern(systemCost, std::move(loopshapingDefinition)));
+ return std::make_unique(systemCost, std::move(loopshapingDefinition));
case LoopshapingType::eliminatepattern:
- return std::unique_ptr(new LoopshapingCostEliminatePattern(systemCost, std::move(loopshapingDefinition)));
+ return std::make_unique(systemCost, std::move(loopshapingDefinition));
default:
throw std::runtime_error("[LoopshapingStateInputCost::create] invalid loopshaping type");
}
diff --git a/ocs2_core/src/loopshaping/dynamics/LoopshapingDynamics.cpp b/ocs2_core/src/loopshaping/dynamics/LoopshapingDynamics.cpp
index d57cc1fe4..3ebb53890 100644
--- a/ocs2_core/src/loopshaping/dynamics/LoopshapingDynamics.cpp
+++ b/ocs2_core/src/loopshaping/dynamics/LoopshapingDynamics.cpp
@@ -125,11 +125,9 @@ std::unique_ptr LoopshapingDynamics::create(const SystemDyn
switch (loopshapingDefinition->getType()) {
case LoopshapingType::outputpattern:
- return std::unique_ptr(
- new LoopshapingDynamicsOutputPattern(systemDynamics, std::move(loopshapingDefinition), preComputation));
+ return std::make_unique(systemDynamics, std::move(loopshapingDefinition), preComputation);
case LoopshapingType::eliminatepattern:
- return std::unique_ptr(
- new LoopshapingDynamicsEliminatePattern(systemDynamics, std::move(loopshapingDefinition), preComputation));
+ return std::make_unique(systemDynamics, std::move(loopshapingDefinition), preComputation);
default:
throw std::runtime_error("[LoopshapingDynamics::create] invalid loopshaping type");
}
diff --git a/ocs2_core/src/loopshaping/soft_constraint/LoopshapingSoftConstraint.cpp b/ocs2_core/src/loopshaping/soft_constraint/LoopshapingSoftConstraint.cpp
index 31d97b7ca..6755ad0f0 100644
--- a/ocs2_core/src/loopshaping/soft_constraint/LoopshapingSoftConstraint.cpp
+++ b/ocs2_core/src/loopshaping/soft_constraint/LoopshapingSoftConstraint.cpp
@@ -42,7 +42,7 @@ namespace LoopshapingSoftConstraint {
std::unique_ptr create(const StateCostCollection& systemSoftConstraint,
std::shared_ptr loopshapingDefinition) {
// State-only soft constraint wrapper is identical to LoopshapingStateCost wrapper.
- return std::unique_ptr(new LoopshapingStateCost(systemSoftConstraint, loopshapingDefinition));
+ return std::make_unique(systemSoftConstraint, std::move(loopshapingDefinition));
}
/******************************************************************************************************/
@@ -52,11 +52,9 @@ std::unique_ptr create(const StateInputCostCollection&
std::shared_ptr loopshapingDefinition) {
switch (loopshapingDefinition->getType()) {
case LoopshapingType::outputpattern:
- return std::unique_ptr(
- new LoopshapingSoftConstraintOutputPattern(systemSoftConstraint, std::move(loopshapingDefinition)));
+ return std::make_unique(systemSoftConstraint, std::move(loopshapingDefinition));
case LoopshapingType::eliminatepattern:
- return std::unique_ptr(
- new LoopshapingSoftConstraintEliminatePattern(systemSoftConstraint, std::move(loopshapingDefinition)));
+ return std::make_unique(systemSoftConstraint, std::move(loopshapingDefinition));
default:
throw std::runtime_error("[LoopshapingStateInputSoftConstraint::create] invalid loopshaping type");
}
diff --git a/ocs2_core/src/misc/LinearAlgebra.cpp b/ocs2_core/src/misc/LinearAlgebra.cpp
index 40a2fb6b0..a51a9031d 100644
--- a/ocs2_core/src/misc/LinearAlgebra.cpp
+++ b/ocs2_core/src/misc/LinearAlgebra.cpp
@@ -32,6 +32,20 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
namespace ocs2 {
namespace LinearAlgebra {
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+void setTriangularMinimumEigenvalues(matrix_t& Lr, scalar_t minEigenValue) {
+ for (Eigen::Index i = 0; i < Lr.rows(); ++i) {
+ scalar_t& eigenValue = Lr(i, i); // diagonal element is the eigenvalue
+ if (eigenValue < 0.0) {
+ eigenValue = std::min(-minEigenValue, eigenValue);
+ } else {
+ eigenValue = std::max(minEigenValue, eigenValue);
+ }
+ }
+}
+
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
@@ -57,6 +71,19 @@ void makePsdEigenvalue(matrix_t& squareMatrix, scalar_t minEigenvalue) {
}
}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+void makePsdGershgorin(matrix_t& squareMatrix, scalar_t minEigenvalue) {
+ assert(squareMatrix.rows() == squareMatrix.cols());
+ squareMatrix = 0.5 * (squareMatrix + squareMatrix.transpose()).eval();
+ for (size_t i = 0; i < squareMatrix.rows(); i++) {
+ // Gershgorin radius: since the matrix is symmetric we use column sum instead of row sum
+ auto Ri = squareMatrix.col(i).cwiseAbs().sum() - std::abs(squareMatrix(i, i));
+ squareMatrix(i, i) = std::max(squareMatrix(i, i), Ri + minEigenvalue);
+ }
+}
+
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
@@ -86,6 +113,16 @@ void makePsdCholesky(matrix_t& A, scalar_t minEigenvalue) {
A.diagonal().array() += minEigenvalue;
}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+void computeInverseMatrixUUT(const matrix_t& Am, matrix_t& AmInvUmUmT) {
+ // Am = Lm Lm^T --> inv(Am) = inv(Lm^T) inv(Lm) where Lm^T is upper triangular
+ Eigen::LLT lltOfA(Am);
+ AmInvUmUmT.setIdentity(Am.rows(), Am.cols()); // for dynamic size matrices
+ lltOfA.matrixU().solveInPlace(AmInvUmUmT);
+}
+
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
@@ -117,6 +154,50 @@ void computeConstraintProjection(const matrix_t& Dm, const matrix_t& RmInvUmUmT,
RmInvConstrainedUUT.noalias() = RmInvUmUmT * QRof_RmInvUmUmTT_DmT_Qu;
}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+std::pair qrConstraintProjection(const VectorFunctionLinearApproximation& constraint) {
+ // Constraint Projectors are based on the QR decomposition
+ const auto numConstraints = constraint.dfdu.rows();
+ const auto numInputs = constraint.dfdu.cols();
+ const Eigen::HouseholderQR QRof_DT(constraint.dfdu.transpose());
+
+ const matrix_t Q = QRof_DT.householderQ();
+ const auto Q1 = Q.leftCols(numConstraints);
+
+ const auto R = QRof_DT.matrixQR().topRows(numConstraints).triangularView();
+ const matrix_t pseudoInverse = R.solve(Q1.transpose()); // left pseudo-inverse of D^T
+
+ VectorFunctionLinearApproximation projectionTerms;
+ projectionTerms.dfdu = Q.rightCols(numInputs - numConstraints);
+ projectionTerms.dfdx.noalias() = -pseudoInverse.transpose() * constraint.dfdx;
+ projectionTerms.f.noalias() = -pseudoInverse.transpose() * constraint.f;
+
+ return std::make_pair(std::move(projectionTerms), std::move(pseudoInverse));
+}
+
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+std::pair luConstraintProjection(const VectorFunctionLinearApproximation& constraint,
+ bool extractPseudoInverse) {
+ // Constraint Projectors are based on the LU decomposition
+ const Eigen::FullPivLU lu(constraint.dfdu);
+
+ VectorFunctionLinearApproximation projectionTerms;
+ projectionTerms.dfdu = lu.kernel();
+ projectionTerms.dfdx.noalias() = -lu.solve(constraint.dfdx);
+ projectionTerms.f.noalias() = -lu.solve(constraint.f);
+
+ matrix_t pseudoInverse;
+ if (extractPseudoInverse) {
+ pseudoInverse = lu.solve(matrix_t::Identity(constraint.f.size(), constraint.f.size())).transpose(); // left pseudo-inverse of D^T
+ }
+
+ return std::make_pair(std::move(projectionTerms), std::move(pseudoInverse));
+}
+
// Explicit instantiations for dynamic sized matrices
template int rank(const matrix_t& A);
template Eigen::VectorXcd eigenvalues(const matrix_t& A);
diff --git a/ocs2_core/src/model_data/Metrics.cpp b/ocs2_core/src/model_data/Metrics.cpp
index 8d1066934..39e440aa7 100644
--- a/ocs2_core/src/model_data/Metrics.cpp
+++ b/ocs2_core/src/model_data/Metrics.cpp
@@ -29,8 +29,67 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include "ocs2_core/model_data/Metrics.h"
+#include
+
namespace ocs2 {
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+void Metrics::swap(Metrics& other) {
+ // Cost
+ std::swap(cost, other.cost);
+ // Dynamics violation
+ dynamicsViolation.swap(other.dynamicsViolation);
+ // Equality constraints
+ stateEqConstraint.swap(other.stateEqConstraint);
+ stateInputEqConstraint.swap(other.stateInputEqConstraint);
+ // Inequality constraints
+ stateIneqConstraint.swap(other.stateIneqConstraint);
+ stateInputIneqConstraint.swap(other.stateInputIneqConstraint);
+ // Lagrangians
+ stateEqLagrangian.swap(other.stateEqLagrangian);
+ stateIneqLagrangian.swap(other.stateIneqLagrangian);
+ stateInputEqLagrangian.swap(other.stateInputEqLagrangian);
+ stateInputIneqLagrangian.swap(other.stateInputIneqLagrangian);
+}
+
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+void Metrics::clear() {
+ // Cost
+ cost = 0.0;
+ // Dynamics violation
+ dynamicsViolation = vector_t();
+ // Equality constraints
+ stateEqConstraint.clear();
+ stateInputEqConstraint.clear();
+ // Inequality constraints
+ stateIneqConstraint.clear();
+ stateInputIneqConstraint.clear();
+ // Lagrangians
+ stateEqLagrangian.clear();
+ stateIneqLagrangian.clear();
+ stateInputEqLagrangian.clear();
+ stateInputIneqLagrangian.clear();
+}
+
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+bool Metrics::isApprox(const Metrics& other, scalar_t prec) const {
+ return numerics::almost_eq(this->cost, other.cost, prec) && this->dynamicsViolation.isApprox(other.dynamicsViolation, prec) &&
+ toVector(this->stateEqConstraint).isApprox(toVector(other.stateEqConstraint), prec) &&
+ toVector(this->stateInputEqConstraint).isApprox(toVector(other.stateInputEqConstraint), prec) &&
+ toVector(this->stateIneqConstraint).isApprox(toVector(other.stateIneqConstraint), prec) &&
+ toVector(this->stateInputIneqConstraint).isApprox(toVector(other.stateInputIneqConstraint), prec) &&
+ toVector(this->stateEqLagrangian).isApprox(toVector(other.stateEqLagrangian), prec) &&
+ toVector(this->stateIneqLagrangian).isApprox(toVector(other.stateIneqLagrangian), prec) &&
+ toVector(this->stateInputEqLagrangian).isApprox(toVector(other.stateInputEqLagrangian), prec) &&
+ toVector(this->stateInputIneqLagrangian).isApprox(toVector(other.stateInputIneqLagrangian), prec);
+}
+
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
@@ -50,6 +109,39 @@ vector_t toVector(const std::vector& termsLagrangianMetrics)
return vec;
}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+vector_t toVector(const vector_array_t& constraintArray) {
+ size_t n = 0;
+ std::for_each(constraintArray.begin(), constraintArray.end(), [&n](const vector_t& v) { n += v.size(); });
+
+ vector_t vec(n);
+ size_t head = 0;
+ for (const auto& v : constraintArray) {
+ vec.segment(head, v.size()) = v;
+ head += v.size();
+ } // end of i loop
+
+ return vec;
+}
+
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+vector_array_t toConstraintArray(const size_array_t& termsSize, const vector_t& vec) {
+ vector_array_t constraintArray;
+ constraintArray.reserve(termsSize.size());
+
+ size_t head = 0;
+ for (const auto& l : termsSize) {
+ constraintArray.emplace_back(vec.segment(head, l));
+ head += l;
+ } // end of i loop
+
+ return constraintArray;
+}
+
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
@@ -66,6 +158,16 @@ std::vector toLagrangianMetrics(const size_array_t& termsSize
return lagrangianMetrics;
}
+/******************************************************************************************************/
+/******************************************************************************************************/
+/******************************************************************************************************/
+size_array_t getSizes(const vector_array_t& constraintArray) {
+ size_array_t s(constraintArray.size());
+ std::transform(constraintArray.begin(), constraintArray.end(), s.begin(),
+ [](const vector_t& v) { return static_cast(v.size()); });
+ return s;
+}
+
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
@@ -99,71 +201,103 @@ LagrangianMetrics interpolate(const index_alpha_t& indexAlpha, const std::vector
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
-MetricsCollection interpolate(const index_alpha_t& indexAlpha, const std::vector& dataArray) {
+Metrics interpolate(const index_alpha_t& indexAlpha, const std::vector& dataArray) {
// number of terms
const auto ind = indexAlpha.second > 0.5 ? indexAlpha.first : indexAlpha.first + 1;
- const size_t mumStateEq = dataArray[ind].stateEqLagrangian.size();
- const size_t mumStateIneq = dataArray[ind].stateIneqLagrangian.size();
- const size_t mumStateInputEq = dataArray[ind].stateInputEqLagrangian.size();
- const size_t mumStateInputIneq = dataArray[ind].stateInputIneqLagrangian.size();
+ const size_t numStateEqConst = dataArray[ind].stateEqConstraint.size();
+ const size_t numStateInputEqCost = dataArray[ind].stateInputEqConstraint.size();
+ const size_t numStateIneqConst = dataArray[ind].stateIneqConstraint.size();
+ const size_t numStateInputIneqCost = dataArray[ind].stateInputIneqConstraint.size();
+ const size_t numStateEqLag = dataArray[ind].stateEqLagrangian.size();
+ const size_t numStateIneqLag = dataArray[ind].stateIneqLagrangian.size();
+ const size_t numStateInputEqLag = dataArray[ind].stateInputEqLagrangian.size();
+ const size_t numStateInputIneqLag = dataArray[ind].stateInputIneqLagrangian.size();
- MetricsCollection out;
+ Metrics out;
// cost
- out.cost = interpolate(indexAlpha, dataArray,
- [](const std::vector& array, size_t t) -> const scalar_t& { return array[t].cost; });
-
- // constraints
- out.stateEqConstraint = interpolate(indexAlpha, dataArray, [](const std::vector& array, size_t t) -> const vector_t& {
- return array[t].stateEqConstraint;
- });
- out.stateInputEqConstraint =
- interpolate(indexAlpha, dataArray,
- [](const std::vector& array, size_t t) -> const vector_t& { return array[t].stateInputEqConstraint; });
+ out.cost =
+ interpolate(indexAlpha, dataArray, [](const std::vector& array, size_t t) -> const scalar_t& { return array[t].cost; });
+
+ // dynamics violation
+ out.dynamicsViolation = interpolate(
+ indexAlpha, dataArray, [](const std::vector& array, size_t t) -> const vector_t& { return array[t].dynamicsViolation; });
+
+ // equality constraints
+ out.stateEqConstraint.reserve(numStateEqConst);
+ for (size_t i = 0; i < numStateEqConst; ++i) {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ return array[t].stateEqConstraint[i];
+ });
+ out.stateEqConstraint.emplace_back(std::move(constraint));
+ }
+ out.stateInputEqConstraint.reserve(numStateInputEqCost);
+ for (size_t i = 0; i < numStateInputEqCost; ++i) {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ return array[t].stateInputEqConstraint[i];
+ });
+ out.stateInputEqConstraint.emplace_back(std::move(constraint));
+ }
+
+ // inequality constraints
+ out.stateIneqConstraint.reserve(numStateIneqConst);
+ for (size_t i = 0; i < numStateIneqConst; ++i) {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ return array[t].stateIneqConstraint[i];
+ });
+ out.stateIneqConstraint.emplace_back(std::move(constraint));
+ }
+ out.stateInputIneqConstraint.reserve(numStateInputIneqCost);
+ for (size_t i = 0; i < numStateInputIneqCost; ++i) {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ return array[t].stateInputIneqConstraint[i];
+ });
+ out.stateInputIneqConstraint.emplace_back(std::move(constraint));
+ }
// state equality Lagrangian
- out.stateEqLagrangian.reserve(mumStateEq);
- for (size_t i = 0; i < mumStateEq; i++) {
- auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
+ out.stateEqLagrangian.reserve(numStateEqLag);
+ for (size_t i = 0; i < numStateEqLag; ++i) {
+ auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
return array[t].stateEqLagrangian[i].penalty;
});
- auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
return array[t].stateEqLagrangian[i].constraint;
});
out.stateEqLagrangian.emplace_back(penalty, std::move(constraint));
} // end of i loop
// state inequality Lagrangian
- out.stateIneqLagrangian.reserve(mumStateIneq);
- for (size_t i = 0; i < mumStateIneq; i++) {
- auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
+ out.stateIneqLagrangian.reserve(numStateIneqLag);
+ for (size_t i = 0; i < numStateIneqLag; ++i) {
+ auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
return array[t].stateIneqLagrangian[i].penalty;
});
- auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
return array[t].stateIneqLagrangian[i].constraint;
});
out.stateIneqLagrangian.emplace_back(penalty, std::move(constraint));
} // end of i loop
// state-input equality Lagrangian
- out.stateInputEqLagrangian.reserve(mumStateInputEq);
- for (size_t i = 0; i < mumStateInputEq; i++) {
- auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
+ out.stateInputEqLagrangian.reserve(numStateInputEqLag);
+ for (size_t i = 0; i < numStateInputEqLag; ++i) {
+ auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
return array[t].stateInputEqLagrangian[i].penalty;
});
- auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
return array[t].stateInputEqLagrangian[i].constraint;
});
out.stateInputEqLagrangian.emplace_back(penalty, std::move(constraint));
} // end of i loop
// state-input inequality Lagrangian
- out.stateInputIneqLagrangian.reserve(mumStateInputIneq);
- for (size_t i = 0; i < mumStateInputIneq; i++) {
- auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
+ out.stateInputIneqLagrangian.reserve(numStateInputIneqLag);
+ for (size_t i = 0; i < numStateInputIneqLag; ++i) {
+ auto penalty = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const scalar_t& {
return array[t].stateInputIneqLagrangian[i].penalty;
});
- auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
+ auto constraint = interpolate(indexAlpha, dataArray, [i](const std::vector& array, size_t t) -> const vector_t& {
return array[t].stateInputIneqLagrangian[i].constraint;
});
out.stateInputIneqLagrangian.emplace_back(penalty, std::move(constraint));
diff --git a/ocs2_core/src/penalties/MultidimensionalPenalty.cpp b/ocs2_core/src/penalties/MultidimensionalPenalty.cpp
index b15175abc..1eaf51b3c 100644
--- a/ocs2_core/src/penalties/MultidimensionalPenalty.cpp
+++ b/ocs2_core/src/penalties/MultidimensionalPenalty.cpp
@@ -63,7 +63,7 @@ class PenaltyBaseWrapper final : public augmented::AugmentedPenaltyBase {
};
std::unique_ptr createWrapper(std::unique_ptr penaltyPtr) {
- return std::unique_ptr(new PenaltyBaseWrapper(std::move(penaltyPtr)));
+ return std::make_unique(std::move(penaltyPtr));
}
scalar_t getMultiplier(const vector_t* l, size_t ind) {
diff --git a/ocs2_core/test/constraint/testConstraintCollection.cpp b/ocs2_core/test/constraint/testConstraintCollection.cpp
index 418275b98..aef7c3fe5 100644
--- a/ocs2_core/test/constraint/testConstraintCollection.cpp
+++ b/ocs2_core/test/constraint/testConstraintCollection.cpp
@@ -27,6 +27,8 @@ OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
******************************************************************************/
+#include
+
#include
#include
@@ -37,7 +39,7 @@ TEST(TestConstraintCollection, add) {
ocs2::StateInputConstraintCollection constraintCollection;
// Add after construction
- std::unique_ptr constraintTerm(new TestEmptyConstraint());
+ auto constraintTerm = std::make_unique();
ASSERT_NO_THROW({ constraintCollection.add("Constraint1", std::move(constraintTerm)); });
}
@@ -48,7 +50,7 @@ TEST(TestConstraintCollection, numberOfConstraints) {
EXPECT_EQ(constraintCollection.getNumConstraints(0.0), 0);
// Add Linear inequality constraint term, which has 2 constraints
- std::unique_ptr constraintTerm(new TestDummyConstraint());
+ auto constraintTerm = std::make_unique();
const size_t addedConstraints = constraintTerm->getNumConstraints(0.0);
constraintCollection.add("Constraint1", std::move(constraintTerm));
@@ -56,6 +58,37 @@ TEST(TestConstraintCollection, numberOfConstraints) {
EXPECT_EQ(constraintCollection.getNumConstraints(0.0), addedConstraints);
}
+TEST(TestConstraintCollection, termsSize) {
+ ocs2::StateInputConstraintCollection constraintCollection;
+
+ // Initially we have zero constraints for all types
+ auto termsSize = constraintCollection.getTermsSize(0.0);
+ EXPECT_EQ(termsSize.size(), 0);
+ EXPECT_EQ(std::accumulate(termsSize.begin(), termsSize.end(), 0), 0);
+
+ // Add 2 Linear inequality constraint term, which has 2 constraints
+ constraintCollection.add("Constraint1", std::make_unique());
+ constraintCollection.add("Constraint2", std::make_unique());
+ constraintCollection.add("Constraint3", std::make_unique());
+ auto& constraint1 = constraintCollection.get("Constraint1");
+ const size_t constraint1Size = constraint1.getNumConstraints(0.0);
+
+ // Check the right constraint size
+ termsSize = constraintCollection.getTermsSize(0.0);
+ EXPECT_EQ(termsSize.size(), 3);
+ if (termsSize.size() == 3) {
+ EXPECT_EQ(termsSize[0], constraint1Size);
+ EXPECT_EQ(termsSize[1], constraint1Size);
+ EXPECT_EQ(termsSize[2], constraint1Size);
+ }
+
+ // Deactivate constraint1
+ constraint1.setActivity(false);
+ termsSize = constraintCollection.getTermsSize(0.0);
+ EXPECT_EQ(termsSize.size(), 3);
+ EXPECT_EQ(std::accumulate(termsSize.begin(), termsSize.end(), 0), 2 * constraint1Size);
+}
+
TEST(TestConstraintCollection, activatingConstraints) {
ocs2::StateInputConstraintCollection constraintCollection;
@@ -63,7 +96,7 @@ TEST(TestConstraintCollection, activatingConstraints) {
EXPECT_EQ(constraintCollection.getNumConstraints(0.0), 0);
// Add Linear inequality constraint term, which has 2 constraints
- std::unique_ptr constraintTerm(new TestDummyConstraint());
+ auto constraintTerm = std::make_unique();
const size_t addedConstraints = constraintTerm->getNumConstraints(0.0);
constraintCollection.add("Constraint1", std::move(constraintTerm));
@@ -78,7 +111,7 @@ TEST(TestConstraintCollection, activatingConstraints) {
TEST(TestConstraintCollection, clone) {
ocs2::StateInputConstraintCollection constraintCollection;
- std::unique_ptr constraintTerm(new TestDummyConstraint());
+ auto constraintTerm = std::make_unique();
const size_t addedConstraints = constraintTerm->getNumConstraints(0.0);
constraintCollection.add("Constraint1", std::move(constraintTerm));
@@ -103,17 +136,16 @@ TEST(TestConstraintCollection, getValue) {
EXPECT_EQ(constraintValues.size(), 0);
// Add Linear inequality constraint term, which has 2 constraints, twice
- std::unique_ptr constraintTerm1(new TestDummyConstraint());
- std::unique_ptr constraintTerm2(new TestDummyConstraint());
+ auto constraintTerm1 = std::make_unique();
+ auto constraintTerm2 = std::make_unique();
constraintCollection.add("Constraint1", std::move(constraintTerm1));
constraintCollection.add("Constraint2", std::move(constraintTerm2));
-
+ const ocs2::vector_t expectedValue = (ocs2::vector_t(2) << 1.0, 2.0).finished();
constraintValues = constraintCollection.getValue(t, x, u, ocs2::PreComputation());
- ASSERT_EQ(constraintValues.rows(), 4);
- EXPECT_EQ(constraintValues[0], 1.0);
- EXPECT_EQ(constraintValues[1], 2.0);
- EXPECT_EQ(constraintValues[2], 1.0);
- EXPECT_EQ(constraintValues[3], 2.0);
+
+ ASSERT_EQ(constraintValues.size(), 2);
+ EXPECT_TRUE(constraintValues[0].isApprox(expectedValue));
+ EXPECT_TRUE(constraintValues[1].isApprox(expectedValue));
}
TEST(TestConstraintCollection, getLinearApproximation) {
@@ -128,8 +160,8 @@ TEST(TestConstraintCollection, getLinearApproximation) {
x.setZero();
// Add Linear inequality constraint term, which has 2 constraints, twice
- std::unique_ptr