5 #include "pinocchio/autodiff/casadi.hpp" 7 #include "pinocchio/algorithm/kinematics.hpp" 8 #include "pinocchio/algorithm/frames.hpp" 9 #include "pinocchio/algorithm/jacobian.hpp" 10 #include "pinocchio/algorithm/crba.hpp" 11 #include "pinocchio/algorithm/rnea.hpp" 12 #include "pinocchio/algorithm/aba.hpp" 13 #include "pinocchio/algorithm/joint-configuration.hpp" 15 #include "pinocchio/parsers/sample-models.hpp" 17 #include <boost/test/unit_test.hpp> 18 #include <boost/utility/binary.hpp> 20 BOOST_AUTO_TEST_SUITE(BOOST_TEST_MODULE)
25 typedef casadi::SX ADScalar;
35 model.lowerPositionLimit.head<3>().
fill(-1.);
36 model.upperPositionLimit.head<3>().
fill(1.);
39 typedef Model::ConfigVectorType ConfigVector;
40 typedef Model::TangentVectorType TangentVector;
41 ConfigVector
q(model.nq);
43 TangentVector
v(TangentVector::Random(model.nv));
45 typedef ADModel::ConfigVectorType ConfigVectorAD;
46 typedef ADModel::TangentVectorType TangentVectorAD;
47 ADModel ad_model = model.cast<ADScalar>();
48 ADData ad_data(ad_model);
51 Data::Matrix6x jacobian_local(6,model.nv), jacobian_world(6,model.nv);
52 jacobian_local.setZero(); jacobian_world.setZero();
54 BOOST_CHECK(jacobian_local.isZero() && jacobian_world.isZero());
61 ConfigVectorAD q_ad(model.nq);
62 for(Eigen::DenseIndex k = 0; k < model.nq; ++k)
66 std::cout <<
"q =\n " << q_ad << std::endl;
69 TangentVectorAD v_ad(model.nv);
70 for(Eigen::DenseIndex k = 0; k < model.nv; ++k)
74 std::cout <<
"v =\n " << v_ad << std::endl;
78 MotionAD & v_local = ad_data.v[(size_t)joint_id];
79 MotionAD v_world = ad_data.oMi[(size_t)joint_id].
act(v_local);
81 casadi::SX cs_v_local(6,1), cs_v_world(6,1);
82 for(Eigen::DenseIndex k = 0; k < 6; ++k)
84 cs_v_local(k) = v_local.toVector()[k];
85 cs_v_world(k) = v_world.toVector()[k];
87 std::cout <<
"v_local = " << cs_v_local << std::endl;
88 std::cout <<
"v_world = " << cs_v_world << std::endl;
90 casadi::Function eval_velocity_local(
"eval_velocity_local",
91 casadi::SXVector {cs_q, cs_v},
92 casadi::SXVector {cs_v_local});
93 std::cout <<
"eval_velocity_local = " << eval_velocity_local << std::endl;
95 casadi::Function eval_velocity_world(
"eval_velocity_world",
96 casadi::SXVector {cs_q, cs_v},
97 casadi::SXVector {cs_v_world});
98 std::cout <<
"eval_velocity_world = " << eval_velocity_world << std::endl;
100 casadi::SX
dv_dv_local = jacobian(cs_v_local, cs_v);
101 casadi::Function eval_jacobian_local(
"eval_jacobian_local",
102 casadi::SXVector {cs_q,cs_v},
103 casadi::SXVector {dv_dv_local});
104 std::cout <<
"eval_jacobian_local = " << eval_jacobian_local << std::endl;
106 casadi::SX dv_dv_world = jacobian(cs_v_world, cs_v);
107 casadi::Function eval_jacobian_world(
"eval_jacobian_world",
108 casadi::SXVector {cs_q,cs_v},
109 casadi::SXVector {dv_dv_world});
111 std::vector<double> q_vec((
size_t)model.nq);
112 Eigen::Map<ConfigVector>(q_vec.data(),model.nq,1) = q;
114 std::vector<double> v_vec((
size_t)model.nv);
115 Eigen::Map<TangentVector>(v_vec.data(),model.nv,1) = v;
117 casadi::DMVector v_local_res = eval_velocity_local(casadi::DMVector {q_vec,v_vec});
118 casadi::DMVector J_local_res = eval_jacobian_local(casadi::DMVector {q_vec,v_vec});
119 std::cout <<
"J_local_res:" << J_local_res << std::endl;
121 std::vector<double> v_local_vec(
static_cast< std::vector<double>
>(v_local_res[0]));
122 BOOST_CHECK((jacobian_local*v).isApprox(Eigen::Map<pinocchio::Motion::Vector6>(v_local_vec.data())));
124 casadi::DMVector v_world_res = eval_velocity_world(casadi::DMVector {q_vec,v_vec});
125 casadi::DMVector J_world_res = eval_jacobian_world(casadi::DMVector {q_vec,v_vec});
127 std::vector<double> v_world_vec(
static_cast< std::vector<double>
>(v_world_res[0]));
128 BOOST_CHECK((jacobian_world*v).isApprox(Eigen::Map<pinocchio::Motion::Vector6>(v_world_vec.data())));
132 std::vector<double> J_local_vec(
static_cast< std::vector<double>
>(J_local_res[0]));
133 J_local_mat = Eigen::Map<Data::Matrix6x>(J_local_vec.data(),6,model.nv);
134 BOOST_CHECK(jacobian_local.isApprox(J_local_mat));
136 std::vector<double> J_world_vec(
static_cast< std::vector<double>
>(J_world_res[0]));
137 J_world_mat = Eigen::Map<Data::Matrix6x>(J_world_vec.data(),6,model.nv);
138 BOOST_CHECK(jacobian_world.isApprox(J_world_mat));
144 typedef casadi::SX ADScalar;
154 model.lowerPositionLimit.head<3>().
fill(-1.);
155 model.upperPositionLimit.head<3>().
fill(1.);
158 typedef Model::ConfigVectorType ConfigVector;
159 typedef Model::TangentVectorType TangentVector;
160 ConfigVector
q(model.nq);
162 TangentVector
v(TangentVector::Random(model.nv));
163 TangentVector
a(TangentVector::Random(model.nv));
167 typedef ADModel::ConfigVectorType ConfigVectorAD;
168 typedef ADModel::TangentVectorType TangentVectorAD;
169 ADModel ad_model = model.cast<ADScalar>();
170 ADData ad_data(ad_model);
173 ConfigVectorAD q_ad(model.nq);
177 TangentVectorAD v_ad(model.nv);
181 TangentVectorAD a_ad(model.nv);
193 typedef casadi::SX ADScalar;
203 model.lowerPositionLimit.head<3>().
fill(-1.);
204 model.upperPositionLimit.head<3>().
fill(1.);
207 typedef Model::ConfigVectorType ConfigVector;
208 typedef Model::TangentVectorType TangentVector;
209 ConfigVector
q(model.nq);
211 TangentVector
v(TangentVector::Random(model.nv));
212 TangentVector
a(TangentVector::Random(model.nv));
214 typedef ADModel::ConfigVectorType ConfigVectorAD;
215 typedef ADModel::TangentVectorType TangentVectorAD;
216 ADModel ad_model = model.cast<ADScalar>();
217 ADData ad_data(ad_model);
222 ConfigVectorAD q_ad(model.nq);
223 q_ad = Eigen::Map<ConfigVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_q).data(),model.nq,1);
226 TangentVectorAD v_ad(model.nv);
227 v_ad = Eigen::Map<TangentVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_v).data(),model.nv,1);
230 TangentVectorAD a_ad(model.nv);
231 a_ad = Eigen::Map<TangentVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_a).data(),model.nv,1);
233 rnea(ad_model,ad_data,q_ad,v_ad,a_ad);
234 casadi::SX tau_ad(model.nv,1);
237 for(Eigen::DenseIndex k = 0; k < model.nv; ++k)
238 tau_ad(k) = ad_data.tau[k];
239 casadi::Function eval_rnea(
"eval_rnea",
240 casadi::SXVector {cs_q, cs_v, cs_a},
241 casadi::SXVector {tau_ad});
243 std::vector<double> q_vec((
size_t)model.nq);
244 Eigen::Map<ConfigVector>(q_vec.data(),model.nq,1) = q;
246 std::vector<double> v_vec((
size_t)model.nv);
247 Eigen::Map<TangentVector>(v_vec.data(),model.nv,1) = v;
249 std::vector<double> a_vec((
size_t)model.nv);
250 Eigen::Map<TangentVector>(a_vec.data(),model.nv,1) = a;
251 casadi::DM tau_res = eval_rnea(casadi::DMVector {q_vec,v_vec,a_vec})[0];
252 std::cout <<
"tau_res = " << tau_res << std::endl;
253 Data::TangentVectorType tau_vec = Eigen::Map<Data::TangentVectorType>(
static_cast< std::vector<double>
>(tau_res).data(),model.nv,1);
255 BOOST_CHECK(data.tau.isApprox(tau_vec));
261 typedef casadi::SX ADScalar;
271 model.lowerPositionLimit.head<3>().
fill(-1.);
272 model.upperPositionLimit.head<3>().
fill(1.);
275 typedef Model::ConfigVectorType ConfigVector;
276 typedef Model::TangentVectorType TangentVector;
277 ConfigVector
q(model.nq);
279 TangentVector
v(TangentVector::Random(model.nv));
280 TangentVector
a(TangentVector::Random(model.nv));
282 typedef ADModel::ConfigVectorType ConfigVectorAD;
283 typedef ADModel::TangentVectorType TangentVectorAD;
284 ADModel ad_model = model.cast<ADScalar>();
285 ADData ad_data(ad_model);
288 data.M.triangularView<Eigen::StrictlyLower>()
289 = data.M.transpose().triangularView<Eigen::StrictlyLower>();
293 ConfigVectorAD q_ad(model.nq);
294 q_ad = Eigen::Map<ConfigVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_q).data(),model.nq,1);
297 TangentVectorAD v_ad(model.nv);
298 v_ad = Eigen::Map<TangentVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_v).data(),model.nv,1);
301 TangentVectorAD a_ad(model.nv);
302 a_ad = Eigen::Map<TangentVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_a).data(),model.nv,1);
305 rnea(ad_model,ad_data,q_ad,v_ad,a_ad);
306 casadi::SX cs_tau(model.nv,1);
309 for(Eigen::DenseIndex k = 0; k < model.nv; ++k)
310 cs_tau(k) = ad_data.tau[k];
311 casadi::Function eval_rnea(
"eval_rnea",
312 casadi::SXVector {cs_q, cs_v, cs_a},
313 casadi::SXVector {cs_tau});
315 crba(ad_model,ad_data,q_ad);
316 ad_data.M.triangularView<Eigen::StrictlyLower>()
317 = ad_data.M.transpose().triangularView<Eigen::StrictlyLower>();
318 casadi::SX M_ad(model.nv,model.nv);
319 for(Eigen::DenseIndex j = 0; j < model.nv; ++j)
321 for(Eigen::DenseIndex
i = 0;
i < model.nv; ++
i)
323 M_ad(
i,j) = ad_data.M(
i,j);
327 std::vector<double> q_vec((
size_t)model.nq);
328 Eigen::Map<ConfigVector>(q_vec.data(),model.nq,1) = q;
330 std::vector<double> v_vec((
size_t)model.nv);
331 Eigen::Map<TangentVector>(v_vec.data(),model.nv,1) = v;
333 std::vector<double> a_vec((
size_t)model.nv);
334 Eigen::Map<TangentVector>(a_vec.data(),model.nv,1) = a;
336 casadi::Function eval_crba(
"eval_crba",
337 casadi::SXVector {cs_q},
338 casadi::SXVector {M_ad});
339 casadi::DM M_res = eval_crba(casadi::DMVector {q_vec})[0];
340 Data::MatrixXs M_mat = Eigen::Map<Data::MatrixXs>(
static_cast< std::vector<double>
>(M_res).data(),
343 BOOST_CHECK(data.M.isApprox(M_mat));
345 casadi::SX
dtau_da = jacobian(cs_tau,cs_a);
346 casadi::Function eval_dtau_da(
"eval_dtau_da",
347 casadi::SXVector {cs_q,cs_v,cs_a},
348 casadi::SXVector {dtau_da});
349 casadi::DM dtau_da_res = eval_dtau_da(casadi::DMVector {q_vec, v_vec, a_vec})[0];
350 Data::MatrixXs dtau_da_mat = Eigen::Map<Data::MatrixXs>(
static_cast< std::vector<double>
>(dtau_da_res).data(),
352 BOOST_CHECK(data.M.isApprox(dtau_da_mat));
358 typedef casadi::SX ADScalar;
368 model.lowerPositionLimit.head<3>().
fill(-1.);
369 model.upperPositionLimit.head<3>().
fill(1.);
372 typedef Model::ConfigVectorType ConfigVector;
373 typedef Model::TangentVectorType TangentVector;
374 ConfigVector
q(model.nq);
376 TangentVector
v(TangentVector::Random(model.nv));
377 TangentVector
tau(TangentVector::Random(model.nv));
379 typedef ADModel::ConfigVectorType ConfigVectorAD;
380 typedef ADModel::TangentVectorType TangentVectorAD;
381 ADModel ad_model = model.cast<ADScalar>();
382 ADData ad_data(ad_model);
387 ConfigVectorAD q_ad(model.nq);
388 q_ad = Eigen::Map<ConfigVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_q).data(),model.nq,1);
391 TangentVectorAD v_ad(model.nv);
392 v_ad = Eigen::Map<TangentVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_v).data(),model.nv,1);
395 TangentVectorAD tau_ad(model.nv);
396 tau_ad = Eigen::Map<TangentVectorAD>(
static_cast< std::vector<ADScalar>
>(cs_tau).data(),model.nv,1);
399 aba(ad_model,ad_data,q_ad,v_ad,tau_ad);
400 casadi::SX cs_ddq(model.nv,1);
401 for(Eigen::DenseIndex k = 0; k < model.nv; ++k)
402 cs_ddq(k) = ad_data.ddq[k];
403 casadi::Function eval_aba(
"eval_aba",
404 casadi::SXVector {cs_q, cs_v, cs_tau},
405 casadi::SXVector {cs_ddq});
407 std::vector<double> q_vec((
size_t)model.nq);
408 Eigen::Map<ConfigVector>(q_vec.data(),model.nq,1) = q;
410 std::vector<double> v_vec((
size_t)model.nv);
411 Eigen::Map<TangentVector>(v_vec.data(),model.nv,1) = v;
413 std::vector<double> tau_vec((
size_t)model.nv);
414 Eigen::Map<TangentVector>(tau_vec.data(),model.nv,1) = tau;
416 casadi::DM ddq_res = eval_aba(casadi::DMVector {q_vec, v_vec, tau_vec})[0];
417 Data::TangentVectorType ddq_mat = Eigen::Map<Data::TangentVectorType>(
static_cast< std::vector<double>
>(ddq_res).data(),
420 BOOST_CHECK(ddq_mat.isApprox(data.ddq));
423 BOOST_AUTO_TEST_SUITE_END()
const DataTpl< Scalar, Options, JointCollectionTpl >::TangentVectorType & aba(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 > &tau)
The Articulated-Body algorithm. It computes the forward dynamics, aka the joint accelerations given t...
The WORLD frame convention corresponds to the frame concident with the Universe/Inertial frame but mo...
BOOST_AUTO_TEST_CASE(test_jacobian)
const DataTpl< Scalar, Options, JointCollectionTpl >::MatrixXs & crba(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, DataTpl< Scalar, Options, JointCollectionTpl > &data, const Eigen::MatrixBase< ConfigVectorType > &q)
Computes the upper triangular part of the joint space inertia matrix M by using the Composite Rigid B...
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 randomConfiguration(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, const Eigen::MatrixBase< ConfigVectorIn1 > &lowerLimits, const Eigen::MatrixBase< ConfigVectorIn2 > &upperLimits, const Eigen::MatrixBase< ReturnType > &qout)
Generate a configuration vector uniformly sampled among provided limits.
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...
static void act(const Eigen::MatrixBase< Mat > &iV, const ForceDense< ForceDerived > &f, Eigen::MatrixBase< MatRet > const &jF)
Action of a motion set on a force object. The input motion set is represented by a 6xN matrix whose e...
void sym(const Eigen::MatrixBase< MatrixDerived > &eig_mat, std::string const &name)
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.
void updateGlobalPlacements(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, DataTpl< Scalar, Options, JointCollectionTpl > &data)
Update the global placement of the joints oMi according to the relative placements of the joints...
The LOCAL frame convention corresponds to the frame directly attached to the moving part (Joint...
Eigen::Matrix< double, 6, Eigen::Dynamic > Matrix6x
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...
void updateFramePlacements(const ModelTpl< Scalar, Options, JointCollectionTpl > &model, DataTpl< Scalar, Options, JointCollectionTpl > &data)
Updates the position of each frame contained in the model.
void copy(::casadi::Matrix< Scalar > const &src, Eigen::MatrixBase< MT > &dst)
void humanoidRandom(ModelTpl< Scalar, Options, JointCollectionTpl > &model, bool usingFF=true)
Create a humanoid kinematic tree with 6-DOF limbs and random joint placements.
void fill(Eigen::Ref< MatType > mat, const typename MatType::Scalar &value)