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...
bool isApprox(const Shape &s1, const Shape &s2, const FCL_REAL tol)
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.