5 #include "pinocchio/algorithm/jacobian.hpp" 6 #include "pinocchio/algorithm/rnea.hpp" 7 #include "pinocchio/algorithm/rnea-derivatives.hpp" 8 #include "pinocchio/algorithm/kinematics-derivatives.hpp" 9 #include "pinocchio/algorithm/contact-dynamics.hpp" 10 #include "pinocchio/algorithm/joint-configuration.hpp" 11 #include "pinocchio/parsers/sample-models.hpp" 15 #include <boost/test/unit_test.hpp> 16 #include <boost/utility/binary.hpp> 18 BOOST_AUTO_TEST_SUITE ( BOOST_TEST_MODULE )
20 using namespace Eigen;
37 const std::string RF =
"rleg6_joint";
44 Motion::Vector6 gamma_RF; gamma_RF.setZero();
46 gamma_RF +=
data.a[RF_id].toVector();
50 Force::Vector6 contact_force_ref =
data.lambda_c;
56 rnea(model,data_check,q,v,ddq_ref,fext);
58 BOOST_CHECK(data_check.
tau.isApprox(tau));
60 BOOST_CHECK(data_check.a[RF_id].toVector().isApprox(-gamma_RF));
67 const double eps = 1e-8;
70 MatrixXd dddq_dtau(model.
nv,model.
nv);
73 for(
int k = 0; k < model.
nv; ++k)
79 Force::Vector6 contact_force_plus = data_fd.lambda_c;
81 dddq_dtau.col(k) = (ddq_plus - ddq_ref)/eps;
82 dlambda_dtau.col(k) = (contact_force_plus - contact_force_ref)/eps;
87 MatrixXd
A(model.
nv+6,model.
nv+6);
88 data.M.transpose().triangularView<Eigen::Upper>() =
data.M.triangularView<Eigen::Upper>();
89 A.topLeftCorner(model.
nv,model.
nv) =
data.M;
90 A.bottomLeftCorner(6, model.
nv) = J_RF;
91 A.topRightCorner(model.
nv, 6) = J_RF.transpose();
92 A.bottomRightCorner(6,6).setZero();
94 MatrixXd Ainv = A.inverse();
95 BOOST_CHECK(Ainv.topRows(model.
nv).leftCols(model.
nv).isApprox(dddq_dtau,std::sqrt(eps)));
96 BOOST_CHECK(Ainv.bottomRows(6).leftCols(model.
nv).isApprox(-dlambda_dtau,std::sqrt(eps)));
99 MatrixXd dddq_dv(model.
nv,model.
nv);
102 for(
int k = 0; k < model.
nv; ++k)
108 Force::Vector6 contact_force_plus = data_fd.lambda_c;
110 dddq_dv.col(k) = (ddq_plus - ddq_ref)/eps;
111 dlambda_dv.col(k) = (contact_force_plus - contact_force_ref)/eps;
117 MatrixXd dddq_dv_anal = -Ainv.topRows(model.
nv).leftCols(model.
nv) * data_check.
dtau_dv;
118 MatrixXd dlambda_dv_anal = -Ainv.bottomRows(6).leftCols(model.
nv) * data_check.
dtau_dv;
120 BOOST_CHECK(dddq_dv_anal.isApprox(dddq_dv,std::sqrt(eps)));
121 BOOST_CHECK(dlambda_dv_anal.isApprox(-dlambda_dv,std::sqrt(eps)));
123 MatrixXd dddq_dq(model.
nv,model.
nv);
126 for(
int k = 0; k < model.
nv; ++k)
135 Force::Vector6 contact_force_plus = data_fd.lambda_c;
137 dddq_dq.col(k) = (ddq_plus - ddq_ref)/eps;
138 dlambda_dq.col(k) = (contact_force_plus - contact_force_ref)/eps;
144 Data::Matrix6x v_partial_dq(6,model.
nv), a_partial_dq(6,model.
nv), a_partial_dv(6,model.
nv), a_partial_da(6,model.
nv);
145 v_partial_dq.setZero(); a_partial_dq.setZero(); a_partial_dv.setZero(); a_partial_da.setZero();
146 Data data_kin(model);
150 MatrixXd dddq_dq_anal = -Ainv.topRows(model.
nv).leftCols(model.
nv) * data_check.
dtau_dq;
151 dddq_dq_anal -= Ainv.topRows(model.
nv).rightCols(6) * a_partial_dq;
153 MatrixXd dlambda_dq_anal = Ainv.bottomRows(6).leftCols(model.
nv) * data_check.
dtau_dq;
154 dlambda_dq_anal += Ainv.bottomRows(6).rightCols(6) * a_partial_dq;
156 BOOST_CHECK(dddq_dq_anal.isApprox(dddq_dq,std::sqrt(eps)));
157 BOOST_CHECK(dlambda_dq_anal.isApprox(dlambda_dq,std::sqrt(eps)));
161 template<
typename ConfigVectorType,
typename TangentVectorType1,
typename TangentVectorType2>
163 const Eigen::MatrixBase<ConfigVectorType> & q,
164 const Eigen::MatrixBase<TangentVectorType1> &
v,
165 const Eigen::MatrixBase<TangentVectorType2> &
tau,
172 Motion::Vector6 gamma;
174 gamma = data.a[
id].toVector();
179 res.head(model.
nv) = data.
ddq;
197 const std::string RF =
"rleg6_joint";
203 Motion::Vector6 gamma_RF; gamma_RF.setZero();
207 Force::Vector6 contact_force_ref = x_ref.tail(6);
213 rnea(model,data_check,q,v,ddq_ref,fext);
215 BOOST_CHECK(data_check.
tau.isApprox(tau));
217 BOOST_CHECK(data_check.a[RF_id].toVector().isZero());
225 const double eps = 1e-8;
228 MatrixXd dddq_dtau(model.
nv,model.
nv);
231 for(
int k = 0; k < model.
nv; ++k)
237 Force::Vector6 contact_force_plus = x_plus.tail(6);
239 dddq_dtau.col(k) = (ddq_plus - ddq_ref)/eps;
240 dlambda_dtau.col(k) = (contact_force_plus - contact_force_ref)/eps;
245 MatrixXd
A(model.
nv+6,model.
nv+6);
246 data.M.transpose().triangularView<Eigen::Upper>() =
data.M.triangularView<Eigen::Upper>();
247 A.topLeftCorner(model.
nv,model.
nv) =
data.M;
248 A.bottomLeftCorner(6, model.
nv) = J_RF;
249 A.topRightCorner(model.
nv, 6) = J_RF.transpose();
250 A.bottomRightCorner(6,6).setZero();
252 MatrixXd Ainv = A.inverse();
253 BOOST_CHECK(Ainv.topRows(model.
nv).leftCols(model.
nv).isApprox(dddq_dtau,std::sqrt(eps)));
254 BOOST_CHECK(Ainv.bottomRows(6).leftCols(model.
nv).isApprox(-dlambda_dtau,std::sqrt(eps)));
257 MatrixXd dddq_dv(model.
nv,model.
nv);
260 for(
int k = 0; k < model.
nv; ++k)
266 Force::Vector6 contact_force_plus = x_plus.tail(6);
268 dddq_dv.col(k) = (ddq_plus - ddq_ref)/eps;
269 dlambda_dv.col(k) = (contact_force_plus - contact_force_ref)/eps;
276 Data::Matrix6x v_partial_dq(6,model.
nv), a_partial_dq(6,model.
nv), a_partial_dv(6,model.
nv), a_partial_da(6,model.
nv);
277 v_partial_dq.setZero(); a_partial_dq.setZero(); a_partial_dv.setZero(); a_partial_da.setZero();
278 Data data_kin(model);
282 MatrixXd dddq_dv_anal = -Ainv.topRows(model.
nv).leftCols(model.
nv) * data_check.
dtau_dv;
283 dddq_dv_anal -= Ainv.topRows(model.
nv).rightCols(6) * a_partial_dv;
284 MatrixXd dlambda_dv_anal = -Ainv.bottomRows(6).leftCols(model.
nv) * data_check.
dtau_dv;
285 dlambda_dv_anal -= Ainv.bottomRows(6).rightCols(6) * a_partial_dv;
287 BOOST_CHECK(dddq_dv_anal.isApprox(dddq_dv,std::sqrt(eps)));
288 BOOST_CHECK(dlambda_dv_anal.isApprox(-dlambda_dv,std::sqrt(eps)));
291 MatrixXd dddq_dq(model.
nv,model.
nv);
294 for(
int k = 0; k < model.
nv; ++k)
302 Force::Vector6 contact_force_plus = x_plus.tail(6);
304 dddq_dq.col(k) = (ddq_plus - ddq_ref)/eps;
305 dlambda_dq.col(k) = (contact_force_plus - contact_force_ref)/eps;
311 v_partial_dq.setZero(); a_partial_dq.setZero(); a_partial_dv.setZero(); a_partial_da.setZero();
315 MatrixXd dddq_dq_anal = -Ainv.topRows(model.
nv).leftCols(model.
nv) * data_check.
dtau_dq;
316 dddq_dq_anal -= Ainv.topRows(model.
nv).rightCols(6) * a_partial_dq;
318 BOOST_CHECK(dddq_dq_anal.isApprox(dddq_dq,std::sqrt(eps)));
320 MatrixXd dlambda_dq_anal = Ainv.bottomRows(6).leftCols(model.
nv) * data_check.
dtau_dq;
321 dlambda_dq_anal += Ainv.bottomRows(6).rightCols(6) * a_partial_dq;
323 BOOST_CHECK(dlambda_dq_anal.isApprox(dlambda_dq,std::sqrt(eps)));
327 BOOST_AUTO_TEST_SUITE_END ()
TangentVectorType tau
Vector of joint torques (dim model.nv).
MatrixXs dtau_dv
Partial derivative of the joint torque vector with respect to the joint velocity. ...
def forwardDynamics(model, data, args)
void rnea(const int num_threads, ModelPoolTpl< Scalar, Options, JointCollectionTpl > &pool, const Eigen::MatrixBase< ConfigVectorPool > &q, const Eigen::MatrixBase< TangentVectorPool1 > &v, const Eigen::MatrixBase< TangentVectorPool2 > &a, const Eigen::MatrixBase< TangentVectorPool3 > &tau)
The Recursive Newton-Euler algorithm. It computes the inverse dynamics, aka the joint torques accordi...
void integrate(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, const Eigen::MatrixBase< ConfigVectorType > &q, const Eigen::MatrixBase< TangentVectorType > &v, const Eigen::MatrixBase< ReturnType > &qout)
Integrate a configuration vector for the specified model for a tangent vector during one unit time...
int njoints
Number of joints.
JointIndex id(const JointModelTpl< Scalar, Options, JointCollectionTpl > &jmodel)
Visit a JointModelTpl through JointIdVisitor to get the index of the joint in the kinematic chain...
MatrixXs dtau_dq
Partial derivative of the joint torque vector with respect to the joint configuration.
void getJointJacobian(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, const DataTpl< Scalar, Options, JointCollectionTpl > &data, const typename ModelTpl< Scalar, Options, JointCollectionTpl >::JointIndex jointId, const ReferenceFrame rf, const Eigen::MatrixBase< Matrix6Like > &J)
Computes the Jacobian of a specific joint frame expressed either in the world (rf = WORLD) frame or i...
void forwardKinematics(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, DataTpl< Scalar, Options, JointCollectionTpl > &data, const Eigen::MatrixBase< ConfigVectorType > &q)
Update the joint placements according to the current joint configuration.
The LOCAL frame convention corresponds to the frame directly attached to the moving part (Joint...
bp::tuple computeRNEADerivatives(const Model &model, Data &data, const Eigen::VectorXd &q, const Eigen::VectorXd &v, const Eigen::VectorXd &a)
#define PINOCCHIO_ALIGNED_STD_VECTOR(Type)
pinocchio::JointIndex JointIndex
void computeForwardKinematicsDerivatives(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, DataTpl< Scalar, Options, JointCollectionTpl > &data, const Eigen::MatrixBase< ConfigVectorType > &q, const Eigen::MatrixBase< TangentVectorType1 > &v, const Eigen::MatrixBase< TangentVectorType2 > &a)
Computes all the terms required to compute the derivatives of the placement, spatial velocity and acc...
void getJointAccelerationDerivatives(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, const DataTpl< Scalar, Options, JointCollectionTpl > &data, const Model::JointIndex jointId, const ReferenceFrame rf, const Eigen::MatrixBase< Matrix6xOut1 > &v_partial_dq, const Eigen::MatrixBase< Matrix6xOut2 > &a_partial_dq, const Eigen::MatrixBase< Matrix6xOut3 > &a_partial_dv, const Eigen::MatrixBase< Matrix6xOut4 > &a_partial_da)
Computes the partial derivaties of the spatial acceleration of a given with respect to the joint conf...
TangentVectorType ddq
The joint accelerations computed from ABA.
VectorXs lambda_c
Lagrange Multipliers corresponding to the contact forces in pinocchio::forwardDynamics.
Main pinocchio namespace.
JointIndex getJointId(const std::string &name) const
Return the index of a joint given by its name.
const DataTpl< Scalar, Options, JointCollectionTpl >::Matrix6x & computeJointJacobians(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, DataTpl< Scalar, Options, JointCollectionTpl > &data, const Eigen::MatrixBase< ConfigVectorType > &q)
Computes the full model Jacobian, i.e. the stack of all motion subspace expressed in the world frame...
int nv
Dimension of the velocity vector space.
void normalize(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, const Eigen::MatrixBase< ConfigVectorType > &qout)
Normalize a configuration vector.
void humanoidRandom(ModelTpl< Scalar, Options, JointCollectionTpl > &model, bool usingFF=true)
Create a humanoid kinematic tree with 6-DOF limbs and random joint placements.
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > VectorXd
Eigen::Matrix< Scalar, 6, Eigen::Dynamic, Options > Matrix6x
The 6d jacobian type (temporary)
int nq
Dimension of the configuration vector representation.
VectorXs TangentVectorType
Dense vectorized version of a joint tangent vector (e.g. velocity, acceleration, etc). It also handles the notion of co-tangent vector (e.g. torque, etc).