GCC Code Coverage Report


Directory: ./
File: include/crocoddyl/multibody/actuations/floating-base-thrusters.hpp
Date: 2025-03-26 19:23:43
Exec Total Coverage
Lines: 93 132 70.5%
Branches: 77 226 34.1%

Line Branch Exec Source
1 ///////////////////////////////////////////////////////////////////////////////
2 // BSD 3-Clause License
3 //
4 // Copyright (C) 2014-2025, Heriot-Watt University
5 // Copyright note valid unless otherwise stated in individual files.
6 // All rights reserved.
7 ///////////////////////////////////////////////////////////////////////////////
8
9 #ifndef CROCODDYL_MULTIBODY_ACTUATIONS_PROPELLERS_HPP_
10 #define CROCODDYL_MULTIBODY_ACTUATIONS_PROPELLERS_HPP_
11
12 #include <pinocchio/multibody/fwd.hpp>
13 #include <pinocchio/spatial/se3.hpp>
14 #include <vector>
15
16 #include "crocoddyl/core/actuation-base.hpp"
17 #include "crocoddyl/core/utils/conversions.hpp"
18 #include "crocoddyl/multibody/states/multibody.hpp"
19
20 namespace crocoddyl {
21
22 enum ThrusterType { CW = 0, CCW };
23
24 template <typename _Scalar>
25 struct ThrusterTpl {
26 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
27
28 typedef _Scalar Scalar;
29 typedef pinocchio::SE3Tpl<Scalar> SE3;
30 typedef ThrusterTpl<Scalar> Thruster;
31
32 /**
33 * @brief Initialize the thruster in a give pose from the root joint.
34 *
35 * @param pose[in] Pose from root joint
36 * @param ctorque[in] Coefficient of generated torque per thrust
37 * @param type[in] Type of thruster (clockwise or counterclockwise,
38 * default clockwise)
39 * @param[in] min_thrust[in] Minimum thrust (default 0.)
40 * @param[in] max_thrust[in] Maximum thrust (default inf number))
41 */
42 6648 ThrusterTpl(const SE3& pose, const Scalar ctorque,
43 const ThrusterType type = CW,
44 const Scalar min_thrust = Scalar(0.),
45 const Scalar max_thrust = std::numeric_limits<Scalar>::infinity())
46 6648 : pose(pose),
47 6648 ctorque(ctorque),
48 6648 type(type),
49 6648 min_thrust(min_thrust),
50 6648 max_thrust(max_thrust) {}
51
52 /**
53 * @brief Initialize the thruster in a pose in the origin of the root joint.
54 *
55 * @param pose[in] Pose from root joint
56 * @param ctorque[in] Coefficient of generated torque per thrust
57 * @param type[in] Type of thruster (clockwise or counterclockwise,
58 * default clockwise)
59 * @param[in] min_thrust[in] Minimum thrust (default 0.)
60 * @param[in] max_thrust[in] Maximum thrust (default inf number))
61 */
62 ThrusterTpl(const Scalar ctorque, const ThrusterType type = CW,
63 const Scalar min_thrust = Scalar(0.),
64 const Scalar max_thrust = std::numeric_limits<Scalar>::infinity())
65 : pose(SE3::Identity()),
66 ctorque(ctorque),
67 type(type),
68 min_thrust(min_thrust),
69 max_thrust(max_thrust) {}
70 12499 ThrusterTpl(const ThrusterTpl<Scalar>& clone)
71 12545 : pose(clone.pose),
72 12545 ctorque(clone.ctorque),
73 12545 type(clone.type),
74 12545 min_thrust(clone.min_thrust),
75 12545 max_thrust(clone.max_thrust) {}
76
77 template <typename NewScalar>
78 20 ThrusterTpl<NewScalar> cast() const {
79 typedef ThrusterTpl<NewScalar> ReturnType;
80
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
20 ReturnType ret(
81 20 pose.template cast<NewScalar>(), scalar_cast<NewScalar>(ctorque), type,
82 20 scalar_cast<NewScalar>(min_thrust), scalar_cast<NewScalar>(max_thrust));
83 20 return ret;
84 }
85
86 ThrusterTpl& operator=(const ThrusterTpl<Scalar>& other) {
87 if (this != &other) {
88 pose = other.pose;
89 ctorque = other.ctorque;
90 type = other.type;
91 min_thrust = other.min_thrust;
92 max_thrust = other.max_thrust;
93 }
94 return *this;
95 }
96
97 template <typename OtherScalar>
98 bool operator==(const ThrusterTpl<OtherScalar>& other) const {
99 return (pose == other.pose && ctorque == other.ctorque &&
100 type == other.type && min_thrust == other.min_thrust &&
101 max_thrust == other.max_thrust);
102 }
103
104 friend std::ostream& operator<<(std::ostream& os,
105 const ThrusterTpl<Scalar>& X) {
106 os << " pose:" << std::endl
107 << X.pose << " ctorque: " << X.ctorque << std::endl
108 << " type: " << X.type << std::endl
109 << "min_thrust: " << X.min_thrust << std::endl
110 << "max_thrust: " << X.max_thrust << std::endl;
111 return os;
112 }
113
114 SE3 pose; //!< Thruster pose
115 Scalar ctorque; //!< Coefficient of generated torque per thrust
116 ThrusterType type; //!< Type of thruster (CW and CCW for clockwise and
117 //!< counterclockwise, respectively)
118 Scalar min_thrust; //!< Minimum thrust
119 Scalar max_thrust; //!< Minimum thrust
120 };
121
122 /**
123 * @brief Actuation models for floating base systems actuated with thrusters
124 *
125 * This actuation model models floating base robots equipped with thrusters,
126 * e.g., multicopters or marine robots equipped with manipulators. It control
127 * inputs are the thrusters' thrust (i.e., forces) and joint torques.
128 *
129 * Both actuation and Jacobians are computed analytically by `calc` and
130 * `calcDiff`, respectively.
131 *
132 * We assume the robot velocity to zero for easily related square thruster
133 * velocities with thrust and torque generated. This approach is similarly
134 * implemented in M. Geisert and N. Mansard, "Trajectory generation for
135 * quadrotor based systems using numerical optimal control", (ICRA). See Section
136 * III.C.
137 *
138 * \sa `ActuationModelAbstractTpl`, `calc()`, `calcDiff()`, `createData()`
139 */
140 template <typename _Scalar>
141 class ActuationModelFloatingBaseThrustersTpl
142 : public ActuationModelAbstractTpl<_Scalar> {
143 public:
144 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
145
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
20 CROCODDYL_DERIVED_CAST(ActuationModelBase,
146 ActuationModelFloatingBaseThrustersTpl)
147
148 typedef _Scalar Scalar;
149 typedef MathBaseTpl<Scalar> MathBase;
150 typedef ActuationModelAbstractTpl<Scalar> Base;
151 typedef ActuationDataAbstractTpl<Scalar> Data;
152 typedef StateMultibodyTpl<Scalar> StateMultibody;
153 typedef ThrusterTpl<Scalar> Thruster;
154 typedef typename MathBase::Vector3s Vector3s;
155 typedef typename MathBase::VectorXs VectorXs;
156 typedef typename MathBase::MatrixXs MatrixXs;
157
158 /**
159 * @brief Initialize the floating base actuation model equipped with
160 * thrusters
161 *
162 * @param[in] state State of the dynamical system
163 * @param[in] thrusters Vector of thrusters
164 */
165 194 ActuationModelFloatingBaseThrustersTpl(std::shared_ptr<StateMultibody> state,
166 const std::vector<Thruster>& thrusters)
167 : Base(state,
168 194 state->get_nv() -
169 194 state->get_pinocchio()
170 194 ->joints[(
171
3/6
✓ Branch 5 taken 194 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 194 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 194 times.
✗ Branch 11 not taken.
388 state->get_pinocchio()->existJointName("root_joint")
172
4/12
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 194 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 194 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 194 times.
✗ Branch 14 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
388 ? state->get_pinocchio()->getJointId("root_joint")
173 : 0)]
174
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
194 .nv() +
175 194 thrusters.size()),
176 194 thrusters_(thrusters),
177 194 n_thrusters_(thrusters.size()),
178
2/4
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 194 times.
✗ Branch 6 not taken.
194 W_thrust_(state_->get_nv(), nu_),
179
5/10
✓ Branch 3 taken 194 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 194 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 194 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 194 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 194 times.
✗ Branch 16 not taken.
776 update_data_(true) {
180
4/8
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 194 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 194 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 194 times.
194 if (!state->get_pinocchio()->existJointName("root_joint")) {
181 throw_pretty(
182 "Invalid argument: "
183 << "the first joint has to be a root one (e.g., free-flyer joint)");
184 }
185 // Update the joint actuation part
186
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
194 W_thrust_.setZero();
187
2/2
✓ Branch 0 taken 83 times.
✓ Branch 1 taken 111 times.
194 if (nu_ > n_thrusters_) {
188
1/2
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
83 W_thrust_.bottomRightCorner(nu_ - n_thrusters_, nu_ - n_thrusters_)
189 .diagonal()
190
2/4
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 83 times.
✗ Branch 5 not taken.
83 .setOnes();
191 }
192 // Update the floating base actuation part
193
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
194 set_thrusters(thrusters_);
194 194 }
195 402 virtual ~ActuationModelFloatingBaseThrustersTpl() = default;
196
197 /**
198 * @brief Compute the actuation signal and actuation set from its thrust
199 * and joint torque inputs \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
200 *
201 * @param[in] data Floating base thrusters actuation data
202 * @param[in] x State point \f$\mathbf{x}\in\mathbb{R}^{ndx}\f$
203 * @param[in] u Joint-torque input \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
204 */
205 8870 virtual void calc(const std::shared_ptr<Data>& data,
206 const Eigen::Ref<const VectorXs>&,
207 const Eigen::Ref<const VectorXs>& u) override {
208
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8870 times.
8870 if (static_cast<std::size_t>(u.size()) != nu_) {
209 throw_pretty(
210 "Invalid argument: " << "u has wrong dimension (it should be " +
211 std::to_string(nu_) + ")");
212 }
213
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8870 times.
8870 if (update_data_) {
214 updateData(data);
215 }
216
2/4
✓ Branch 4 taken 8870 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8870 times.
✗ Branch 8 not taken.
8870 data->tau.noalias() = data->dtau_du * u;
217 8870 }
218
219 /**
220 * @brief Compute the Jacobians of the floating base thruster actuation
221 * function
222 *
223 * @param[in] data Floating base thrusters actuation data
224 * @param[in] x State point \f$\mathbf{x}\in\mathbb{R}^{ndx}\f$
225 * @param[in] u Joint-torque input \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
226 */
227 #ifndef NDEBUG
228 1604 virtual void calcDiff(const std::shared_ptr<Data>& data,
229 const Eigen::Ref<const VectorXs>&,
230 const Eigen::Ref<const VectorXs>&) override {
231 #else
232 virtual void calcDiff(const std::shared_ptr<Data>&,
233 const Eigen::Ref<const VectorXs>&,
234 const Eigen::Ref<const VectorXs>&) override {
235 #endif
236 // The derivatives has constant values which were set in createData.
237
2/12
✓ Branch 4 taken 1604 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1604 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
1604 assert_pretty(MatrixXs(data->dtau_du).isApprox(W_thrust_),
238 "dtau_du has wrong value");
239 1604 }
240
241 2314 virtual void commands(const std::shared_ptr<Data>& data,
242 const Eigen::Ref<const VectorXs>&,
243 const Eigen::Ref<const VectorXs>& tau) override {
244
2/4
✓ Branch 4 taken 2314 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2314 times.
✗ Branch 8 not taken.
2314 data->u.noalias() = data->Mtau * tau;
245 2314 }
246
247 #ifndef NDEBUG
248 630 virtual void torqueTransform(const std::shared_ptr<Data>& data,
249 const Eigen::Ref<const VectorXs>&,
250 const Eigen::Ref<const VectorXs>&) override {
251 #else
252 virtual void torqueTransform(const std::shared_ptr<Data>&,
253 const Eigen::Ref<const VectorXs>&,
254 const Eigen::Ref<const VectorXs>&) override {
255 #endif
256 // The torque transform has constant values which were set in createData.
257
2/12
✓ Branch 4 taken 630 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 630 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
630 assert_pretty(MatrixXs(data->Mtau).isApprox(Mtau_), "Mtau has wrong value");
258 630 }
259
260 /**
261 * @brief Create the floating base thruster actuation data
262 *
263 * @return the actuation data
264 */
265 6518 virtual std::shared_ptr<Data> createData() override {
266
1/2
✓ Branch 1 taken 6518 times.
✗ Branch 2 not taken.
6518 std::shared_ptr<Data> data =
267 6518 std::allocate_shared<Data>(Eigen::aligned_allocator<Data>(), this);
268
1/2
✓ Branch 1 taken 6518 times.
✗ Branch 2 not taken.
6518 updateData(data);
269 6518 return data;
270 }
271
272 template <typename NewScalar>
273 5 ActuationModelFloatingBaseThrustersTpl<NewScalar> cast() const {
274 typedef ActuationModelFloatingBaseThrustersTpl<NewScalar> ReturnType;
275 typedef StateMultibodyTpl<NewScalar> StateType;
276 typedef ThrusterTpl<NewScalar> ThrusterType;
277
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
5 std::vector<ThrusterType> thrusters = vector_cast<NewScalar>(thrusters_);
278
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
10 ReturnType ret(
279 5 std::static_pointer_cast<StateType>(state_->template cast<NewScalar>()),
280 thrusters);
281 10 return ret;
282 5 }
283
284 /**
285 * @brief Return the vector of thrusters
286 */
287 const std::vector<Thruster>& get_thrusters() const { return thrusters_; }
288
289 /**
290 * @brief Return the number of thrusters
291 */
292 std::size_t get_nthrusters() const { return n_thrusters_; }
293
294 /**
295 * @brief Modify the vector of thrusters
296 *
297 * Since we don't want to allocate data, we request to pass the same
298 * number of thrusters.
299 *
300 * @param[in] thrusters Vector of thrusters
301 */
302 194 void set_thrusters(const std::vector<Thruster>& thrusters) {
303
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 194 times.
194 if (static_cast<std::size_t>(thrusters.size()) != n_thrusters_) {
304 throw_pretty("Invalid argument: "
305 << "the number of thrusters is wrong (it should be " +
306 std::to_string(n_thrusters_) + ")");
307 }
308 194 thrusters_ = thrusters;
309 // Update the mapping matrix from thrusters thrust to body force/moments
310
2/2
✓ Branch 0 taken 776 times.
✓ Branch 1 taken 194 times.
970 for (std::size_t i = 0; i < n_thrusters_; ++i) {
311 776 const Thruster& p = thrusters_[i];
312
4/8
✓ Branch 1 taken 776 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 776 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 776 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 776 times.
✗ Branch 11 not taken.
776 const Vector3s& f_z = p.pose.rotation() * Vector3s::UnitZ();
313
3/6
✓ Branch 1 taken 776 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 776 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 776 times.
✗ Branch 8 not taken.
776 W_thrust_.template topRows<3>().col(i) += f_z;
314
4/8
✓ Branch 1 taken 776 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 776 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 776 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 776 times.
✗ Branch 11 not taken.
776 W_thrust_.template middleRows<3>(3).col(i).noalias() +=
315
2/4
✓ Branch 1 taken 776 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 776 times.
✗ Branch 5 not taken.
776 p.pose.translation().cross(f_z);
316
2/3
✓ Branch 0 taken 388 times.
✓ Branch 1 taken 388 times.
✗ Branch 2 not taken.
776 switch (p.type) {
317 388 case CW:
318
4/8
✓ Branch 1 taken 388 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 388 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 388 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 388 times.
✗ Branch 11 not taken.
388 W_thrust_.template middleRows<3>(3).col(i) += p.ctorque * f_z;
319 388 break;
320 388 case CCW:
321
4/8
✓ Branch 1 taken 388 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 388 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 388 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 388 times.
✗ Branch 11 not taken.
388 W_thrust_.template middleRows<3>(3).col(i) -= p.ctorque * f_z;
322 388 break;
323 }
324 }
325 // Compute the torque transform matrix from generalized torques to joint
326 // torque inputs
327
1/2
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
194 Mtau_ = pseudoInverse(W_thrust_);
328
2/4
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 194 times.
✗ Branch 6 not taken.
194 S_.noalias() = W_thrust_ * Mtau_;
329 194 update_data_ = true;
330 194 }
331
332 const MatrixXs& get_Wthrust() const { return W_thrust_; }
333
334 const MatrixXs& get_S() const { return S_; }
335
336 void print(std::ostream& os) const override {
337 os << "ActuationModelFloatingBaseThrusters {nu=" << nu_
338 << ", nthrusters=" << n_thrusters_ << ", thrusters=" << std::endl;
339 for (std::size_t i = 0; i < n_thrusters_; ++i) {
340 os << std::to_string(i) << ": " << thrusters_[i];
341 }
342 os << "}";
343 }
344
345 protected:
346 std::vector<Thruster> thrusters_; //!< Vector of thrusters
347 std::size_t n_thrusters_; //!< Number of thrusters
348 MatrixXs W_thrust_; //!< Matrix from thrusters thrusts to body wrench
349 MatrixXs Mtau_; //!< Constaint torque transform from generalized torques to
350 //!< joint torque inputs
351 MatrixXs S_; //!< Selection matrix for under-actuation part
352
353 bool update_data_;
354 using Base::nu_;
355 using Base::state_;
356
357 private:
358 6518 void updateData(const std::shared_ptr<Data>& data) {
359 6518 data->dtau_du = W_thrust_;
360 6518 data->Mtau = Mtau_;
361 6518 const std::size_t nv = state_->get_nv();
362
2/2
✓ Branch 0 taken 41050 times.
✓ Branch 1 taken 6518 times.
47568 for (std::size_t k = 0; k < nv; ++k) {
363
3/3
✓ Branch 2 taken 13026 times.
✓ Branch 3 taken 28004 times.
✓ Branch 4 taken 20 times.
41050 if (fabs(S_(k, k)) < std::numeric_limits<Scalar>::epsilon()) {
364 13036 data->tau_set[k] = false;
365 } else {
366 28014 data->tau_set[k] = true;
367 }
368 }
369 6518 update_data_ = false;
370 6518 }
371 };
372
373 } // namespace crocoddyl
374
375 extern template class CROCODDYL_EXPLICIT_INSTANTIATION_DECLARATION_DLLAPI
376 crocoddyl::ActuationModelFloatingBaseThrustersTpl<double>;
377 extern template class CROCODDYL_EXPLICIT_INSTANTIATION_DECLARATION_DLLAPI
378 crocoddyl::ActuationModelFloatingBaseThrustersTpl<float>;
379
380 #endif // CROCODDYL_MULTIBODY_ACTUATIONS_PROPELLERS_HPP_
381