GCC Code Coverage Report


Directory: ./
File: include/crocoddyl/core/integrator/rk4.hxx
Date: 2025-01-16 08:47:40
Exec Total Coverage
Lines: 0 216 0.0%
Branches: 0 654 0.0%

Line Branch Exec Source
1 ///////////////////////////////////////////////////////////////////////////////
2 // BSD 3-Clause License
3 //
4 // Copyright (C) 2019-2022, LAAS-CNRS, IRI: CSIC-UPC, University of Edinburgh
5 // Heriot-Watt University
6 // Copyright note valid unless otherwise stated in individual files.
7 // All rights reserved.
8 ///////////////////////////////////////////////////////////////////////////////
9
10 #include <iostream>
11
12 #include "crocoddyl/core/integrator/rk4.hpp"
13 #include "crocoddyl/core/utils/exception.hpp"
14
15 namespace crocoddyl {
16
17 template <typename Scalar>
18 IntegratedActionModelRK4Tpl<Scalar>::IntegratedActionModelRK4Tpl(
19 boost::shared_ptr<DifferentialActionModelAbstract> model,
20 boost::shared_ptr<ControlParametrizationModelAbstract> control,
21 const Scalar time_step, const bool with_cost_residual)
22 : Base(model, control, time_step, with_cost_residual) {
23 rk4_c_ = {Scalar(0.), Scalar(0.5), Scalar(0.5), Scalar(1.)};
24 std::cerr
25 << "Deprecated IntegratedActionModelRK4: Use IntegratedActionModelRK"
26 << std::endl;
27 }
28
29 template <typename Scalar>
30 IntegratedActionModelRK4Tpl<Scalar>::IntegratedActionModelRK4Tpl(
31 boost::shared_ptr<DifferentialActionModelAbstract> model,
32 const Scalar time_step, const bool with_cost_residual)
33 : Base(model, time_step, with_cost_residual) {
34 rk4_c_ = {Scalar(0.), Scalar(0.5), Scalar(0.5), Scalar(1.)};
35 std::cerr
36 << "Deprecated IntegratedActionModelRK4: Use IntegratedActionModelRK"
37 << std::endl;
38 }
39
40 template <typename Scalar>
41 IntegratedActionModelRK4Tpl<Scalar>::~IntegratedActionModelRK4Tpl() {}
42
43 template <typename Scalar>
44 void IntegratedActionModelRK4Tpl<Scalar>::calc(
45 const boost::shared_ptr<ActionDataAbstract>& data,
46 const Eigen::Ref<const VectorXs>& x, const Eigen::Ref<const VectorXs>& u) {
47 if (static_cast<std::size_t>(x.size()) != state_->get_nx()) {
48 throw_pretty(
49 "Invalid argument: " << "x has wrong dimension (it should be " +
50 std::to_string(state_->get_nx()) + ")");
51 }
52 if (static_cast<std::size_t>(u.size()) != nu_) {
53 throw_pretty(
54 "Invalid argument: " << "u has wrong dimension (it should be " +
55 std::to_string(nu_) + ")");
56 }
57 const std::size_t nv = state_->get_nv();
58 Data* d = static_cast<Data*>(data.get());
59
60 const boost::shared_ptr<DifferentialActionDataAbstract>& k0_data =
61 d->differential[0];
62 const boost::shared_ptr<ControlParametrizationDataAbstract>& u0_data =
63 d->control[0];
64 control_->calc(u0_data, rk4_c_[0], u);
65 d->ws[0] = u0_data->w;
66 differential_->calc(k0_data, x, d->ws[0]);
67 d->y[0] = x;
68 d->ki[0].head(nv) = d->y[0].tail(nv);
69 d->ki[0].tail(nv) = k0_data->xout;
70 d->integral[0] = k0_data->cost;
71 for (std::size_t i = 1; i < 4; ++i) {
72 const boost::shared_ptr<DifferentialActionDataAbstract>& ki_data =
73 d->differential[i];
74 const boost::shared_ptr<ControlParametrizationDataAbstract>& ui_data =
75 d->control[i];
76 d->dx_rk4[i].noalias() = time_step_ * rk4_c_[i] * d->ki[i - 1];
77 state_->integrate(x, d->dx_rk4[i], d->y[i]);
78 control_->calc(ui_data, rk4_c_[i], u);
79 d->ws[i] = ui_data->w;
80 differential_->calc(ki_data, d->y[i], d->ws[i]);
81 d->ki[i].head(nv) = d->y[i].tail(nv);
82 d->ki[i].tail(nv) = ki_data->xout;
83 d->integral[i] = ki_data->cost;
84 }
85 d->dx =
86 (d->ki[0] + Scalar(2.) * d->ki[1] + Scalar(2.) * d->ki[2] + d->ki[3]) *
87 time_step_ / Scalar(6.);
88 state_->integrate(x, d->dx, d->xnext);
89 d->cost = (d->integral[0] + Scalar(2.) * d->integral[1] +
90 Scalar(2.) * d->integral[2] + d->integral[3]) *
91 time_step_ / Scalar(6.);
92 d->g = k0_data->g;
93 d->h = k0_data->h;
94 if (with_cost_residual_) {
95 d->r = k0_data->r;
96 }
97 }
98
99 template <typename Scalar>
100 void IntegratedActionModelRK4Tpl<Scalar>::calc(
101 const boost::shared_ptr<ActionDataAbstract>& data,
102 const Eigen::Ref<const VectorXs>& x) {
103 if (static_cast<std::size_t>(x.size()) != state_->get_nx()) {
104 throw_pretty(
105 "Invalid argument: " << "x has wrong dimension (it should be " +
106 std::to_string(state_->get_nx()) + ")");
107 }
108 Data* d = static_cast<Data*>(data.get());
109
110 const boost::shared_ptr<DifferentialActionDataAbstract>& k0_data =
111 d->differential[0];
112 differential_->calc(k0_data, x);
113 d->dx.setZero();
114 d->xnext = x;
115 d->cost = k0_data->cost;
116 d->g = k0_data->g;
117 d->h = k0_data->h;
118 if (with_cost_residual_) {
119 d->r = k0_data->r;
120 }
121 }
122
123 template <typename Scalar>
124 void IntegratedActionModelRK4Tpl<Scalar>::calcDiff(
125 const boost::shared_ptr<ActionDataAbstract>& data,
126 const Eigen::Ref<const VectorXs>& x, const Eigen::Ref<const VectorXs>& u) {
127 if (static_cast<std::size_t>(x.size()) != state_->get_nx()) {
128 throw_pretty(
129 "Invalid argument: " << "x has wrong dimension (it should be " +
130 std::to_string(state_->get_nx()) + ")");
131 }
132 if (static_cast<std::size_t>(u.size()) != nu_) {
133 throw_pretty(
134 "Invalid argument: " << "u has wrong dimension (it should be " +
135 std::to_string(nu_) + ")");
136 }
137 const std::size_t nv = state_->get_nv();
138 const std::size_t nu = control_->get_nu();
139 Data* d = static_cast<Data*>(data.get());
140 assert_pretty(
141 MatrixXs(d->dyi_dx[0])
142 .isApprox(MatrixXs::Identity(state_->get_ndx(), state_->get_ndx())),
143 "you have changed dyi_dx[0] values that supposed to be constant.");
144 assert_pretty(
145 MatrixXs(d->dki_dx[0])
146 .topRightCorner(nv, nv)
147 .isApprox(MatrixXs::Identity(nv, nv)),
148 "you have changed dki_dx[0] values that supposed to be constant.");
149
150 for (std::size_t i = 0; i < 4; ++i) {
151 differential_->calcDiff(d->differential[i], d->y[i], d->ws[i]);
152 }
153
154 const boost::shared_ptr<DifferentialActionDataAbstract>& k0_data =
155 d->differential[0];
156 const boost::shared_ptr<ControlParametrizationDataAbstract>& u0_data =
157 d->control[0];
158 d->dki_dx[0].bottomRows(nv) = k0_data->Fx;
159 control_->multiplyByJacobian(
160 u0_data, k0_data->Fu,
161 d->dki_du[0].bottomRows(nv)); // dki_du = dki_dw * dw_du
162
163 d->dli_dx[0] = k0_data->Lx;
164 control_->multiplyJacobianTransposeBy(
165 u0_data, k0_data->Lu,
166 d->dli_du[0]); // dli_du = dli_dw * dw_du
167
168 d->ddli_ddx[0] = k0_data->Lxx;
169 d->ddli_ddw[0] = k0_data->Luu;
170 control_->multiplyByJacobian(
171 u0_data, d->ddli_ddw[0],
172 d->ddli_dwdu[0]); // ddli_dwdu = ddli_ddw * dw_du
173 control_->multiplyJacobianTransposeBy(
174 u0_data, d->ddli_dwdu[0],
175 d->ddli_ddu[0]); // ddli_ddu = dw_du.T * ddli_dwdu
176 d->ddli_dxdw[0] = k0_data->Lxu;
177 control_->multiplyByJacobian(
178 u0_data, d->ddli_dxdw[0],
179 d->ddli_dxdu[0]); // ddli_dxdu = ddli_dxdw * dw_du
180
181 for (std::size_t i = 1; i < 4; ++i) {
182 const boost::shared_ptr<DifferentialActionDataAbstract>& ki_data =
183 d->differential[i];
184 const boost::shared_ptr<ControlParametrizationDataAbstract>& ui_data =
185 d->control[i];
186 d->dyi_dx[i].noalias() = d->dki_dx[i - 1] * rk4_c_[i] * time_step_;
187 d->dyi_du[i].noalias() = d->dki_du[i - 1] * rk4_c_[i] * time_step_;
188 state_->JintegrateTransport(x, d->dx_rk4[i], d->dyi_dx[i], second);
189 state_->Jintegrate(x, d->dx_rk4[i], d->dyi_dx[i], d->dyi_dx[i], first,
190 addto);
191 state_->JintegrateTransport(x, d->dx_rk4[i], d->dyi_du[i],
192 second); // dyi_du = Jintegrate * dyi_du
193
194 // Sparse matrix-matrix multiplication for computing:
195 Eigen::Block<MatrixXs> dkvi_dq = d->dki_dx[i].bottomLeftCorner(nv, nv);
196 Eigen::Block<MatrixXs> dkvi_dv = d->dki_dx[i].bottomRightCorner(nv, nv);
197 Eigen::Block<MatrixXs> dkqi_du = d->dki_du[i].topLeftCorner(nv, nu);
198 Eigen::Block<MatrixXs> dkvi_du = d->dki_du[i].bottomLeftCorner(nv, nu);
199 const Eigen::Block<MatrixXs> dki_dqi = ki_data->Fx.bottomLeftCorner(nv, nv);
200 const Eigen::Block<MatrixXs> dki_dvi =
201 ki_data->Fx.bottomRightCorner(nv, nv);
202 const Eigen::Block<MatrixXs> dqi_dq = d->dyi_dx[i].topLeftCorner(nv, nv);
203 const Eigen::Block<MatrixXs> dqi_dv = d->dyi_dx[i].topRightCorner(nv, nv);
204 const Eigen::Block<MatrixXs> dvi_dq = d->dyi_dx[i].bottomLeftCorner(nv, nv);
205 const Eigen::Block<MatrixXs> dvi_dv =
206 d->dyi_dx[i].bottomRightCorner(nv, nv);
207 const Eigen::Block<MatrixXs> dqi_du = d->dyi_du[i].topLeftCorner(nv, nu);
208 const Eigen::Block<MatrixXs> dvi_du = d->dyi_du[i].bottomLeftCorner(nv, nu);
209 // i. d->dki_dx[i].noalias() = d->dki_dy[i] * d->dyi_dx[i], where dki_dy
210 // is ki_data.Fx
211 d->dki_dx[i].topRows(nv) = d->dyi_dx[i].bottomRows(nv);
212 dkvi_dq.noalias() = dki_dqi * dqi_dq;
213 if (i == 1) {
214 dkvi_dv = time_step_ / Scalar(2.) * dki_dqi;
215 } else {
216 dkvi_dv.noalias() = dki_dqi * dqi_dv;
217 }
218 dkvi_dq.noalias() += dki_dvi * dvi_dq;
219 dkvi_dv.noalias() += dki_dvi * dvi_dv;
220 // ii. d->dki_du[i].noalias() = d->dki_dy[i] * d->dyi_du[i], where dki_dy
221 // is ki_data.Fx
222 dkqi_du = dvi_du;
223 dkvi_du.noalias() = dki_dqi * dqi_du;
224 dkvi_du.noalias() += dki_dvi * dvi_du;
225
226 control_->multiplyByJacobian(ui_data, ki_data->Fu,
227 d->dki_du[i].bottomRows(nv),
228 addto); // dfi_du = dki_dw * dw_du
229
230 d->dli_dx[i].noalias() = ki_data->Lx.transpose() * d->dyi_dx[i];
231 control_->multiplyJacobianTransposeBy(ui_data, ki_data->Lu,
232 d->dli_du[i]); // dli_du = Lu * dw_du
233 d->dli_du[i].noalias() += ki_data->Lx.transpose() * d->dyi_du[i];
234
235 d->Lxx_partialx[i].noalias() = ki_data->Lxx * d->dyi_dx[i];
236 d->ddli_ddx[i].noalias() = d->dyi_dx[i].transpose() * d->Lxx_partialx[i];
237
238 control_->multiplyByJacobian(ui_data, ki_data->Lxu,
239 d->Lxu_i[i]); // Lxu = Lxw * dw_du
240 d->Luu_partialx[i].noalias() = d->Lxu_i[i].transpose() * d->dyi_du[i];
241 d->Lxx_partialu[i].noalias() = ki_data->Lxx * d->dyi_du[i];
242 control_->multiplyByJacobian(
243 ui_data, ki_data->Luu,
244 d->ddli_dwdu[i]); // ddli_dwdu = ddli_ddw * dw_du
245 control_->multiplyJacobianTransposeBy(
246 ui_data, d->ddli_dwdu[i],
247 d->ddli_ddu[i]); // ddli_ddu = dw_du.T * ddli_dwdu
248 d->ddli_ddu[i].noalias() += d->Luu_partialx[i].transpose() +
249 d->Luu_partialx[i] +
250 d->dyi_du[i].transpose() * d->Lxx_partialu[i];
251
252 d->ddli_dxdw[i].noalias() = d->dyi_dx[i].transpose() * ki_data->Lxu;
253 control_->multiplyByJacobian(
254 ui_data, d->ddli_dxdw[i],
255 d->ddli_dxdu[i]); // ddli_dxdu = ddli_dxdw * dw_du
256 d->ddli_dxdu[i].noalias() += d->dyi_dx[i].transpose() * d->Lxx_partialu[i];
257 }
258
259 d->Fx.noalias() = time_step_ / Scalar(6.) *
260 (d->dki_dx[0] + Scalar(2.) * d->dki_dx[1] +
261 Scalar(2.) * d->dki_dx[2] + d->dki_dx[3]);
262 state_->JintegrateTransport(x, d->dx, d->Fx, second);
263 state_->Jintegrate(x, d->dx, d->Fx, d->Fx, first, addto);
264
265 d->Fu.noalias() = time_step_ / Scalar(6.) *
266 (d->dki_du[0] + Scalar(2.) * d->dki_du[1] +
267 Scalar(2.) * d->dki_du[2] + d->dki_du[3]);
268 state_->JintegrateTransport(x, d->dx, d->Fu, second);
269
270 d->Lx.noalias() = time_step_ / Scalar(6.) *
271 (d->dli_dx[0] + Scalar(2.) * d->dli_dx[1] +
272 Scalar(2.) * d->dli_dx[2] + d->dli_dx[3]);
273 d->Lu.noalias() = time_step_ / Scalar(6.) *
274 (d->dli_du[0] + Scalar(2.) * d->dli_du[1] +
275 Scalar(2.) * d->dli_du[2] + d->dli_du[3]);
276
277 d->Lxx.noalias() = time_step_ / Scalar(6.) *
278 (d->ddli_ddx[0] + Scalar(2.) * d->ddli_ddx[1] +
279 Scalar(2.) * d->ddli_ddx[2] + d->ddli_ddx[3]);
280 d->Luu.noalias() = time_step_ / Scalar(6.) *
281 (d->ddli_ddu[0] + Scalar(2.) * d->ddli_ddu[1] +
282 Scalar(2.) * d->ddli_ddu[2] + d->ddli_ddu[3]);
283 d->Lxu.noalias() = time_step_ / Scalar(6.) *
284 (d->ddli_dxdu[0] + Scalar(2.) * d->ddli_dxdu[1] +
285 Scalar(2.) * d->ddli_dxdu[2] + d->ddli_dxdu[3]);
286 d->Gx = k0_data->Gx;
287 d->Hx = k0_data->Hx;
288 d->Gu.resize(differential_->get_ng(), nu_);
289 d->Hu.resize(differential_->get_nh(), nu_);
290 control_->multiplyByJacobian(u0_data, k0_data->Gu, d->Gu);
291 control_->multiplyByJacobian(u0_data, k0_data->Hu, d->Hu);
292 }
293
294 template <typename Scalar>
295 void IntegratedActionModelRK4Tpl<Scalar>::calcDiff(
296 const boost::shared_ptr<ActionDataAbstract>& data,
297 const Eigen::Ref<const VectorXs>& x) {
298 if (static_cast<std::size_t>(x.size()) != state_->get_nx()) {
299 throw_pretty(
300 "Invalid argument: " << "x has wrong dimension (it should be " +
301 std::to_string(state_->get_nx()) + ")");
302 }
303 Data* d = static_cast<Data*>(data.get());
304
305 const boost::shared_ptr<DifferentialActionDataAbstract>& k0_data =
306 d->differential[0];
307 differential_->calcDiff(k0_data, x);
308 d->Lx = k0_data->Lx;
309 d->Lxx = k0_data->Lxx;
310 d->Gx = k0_data->Gx;
311 d->Hx = k0_data->Hx;
312 }
313
314 template <typename Scalar>
315 boost::shared_ptr<ActionDataAbstractTpl<Scalar> >
316 IntegratedActionModelRK4Tpl<Scalar>::createData() {
317 return boost::allocate_shared<Data>(Eigen::aligned_allocator<Data>(), this);
318 }
319
320 template <typename Scalar>
321 bool IntegratedActionModelRK4Tpl<Scalar>::checkData(
322 const boost::shared_ptr<ActionDataAbstract>& data) {
323 boost::shared_ptr<Data> d = boost::dynamic_pointer_cast<Data>(data);
324 if (data != NULL) {
325 return differential_->checkData(d->differential[0]) &&
326 differential_->checkData(d->differential[2]) &&
327 differential_->checkData(d->differential[1]) &&
328 differential_->checkData(d->differential[3]);
329 } else {
330 return false;
331 }
332 }
333
334 template <typename Scalar>
335 void IntegratedActionModelRK4Tpl<Scalar>::quasiStatic(
336 const boost::shared_ptr<ActionDataAbstract>& data, Eigen::Ref<VectorXs> u,
337 const Eigen::Ref<const VectorXs>& x, const std::size_t maxiter,
338 const Scalar tol) {
339 if (static_cast<std::size_t>(u.size()) != nu_) {
340 throw_pretty(
341 "Invalid argument: " << "u has wrong dimension (it should be " +
342 std::to_string(nu_) + ")");
343 }
344 if (static_cast<std::size_t>(x.size()) != state_->get_nx()) {
345 throw_pretty(
346 "Invalid argument: " << "x has wrong dimension (it should be " +
347 std::to_string(state_->get_nx()) + ")");
348 }
349
350 Data* d = static_cast<Data*>(data.get());
351 const boost::shared_ptr<ControlParametrizationDataAbstract>& u0_data =
352 d->control[0];
353 u0_data->w *= 0.;
354 differential_->quasiStatic(d->differential[0], u0_data->w, x, maxiter, tol);
355 control_->params(u0_data, 0., u0_data->w);
356 u = u0_data->u;
357 }
358
359 template <typename Scalar>
360 void IntegratedActionModelRK4Tpl<Scalar>::print(std::ostream& os) const {
361 os << "IntegratedActionModelRK4 {dt=" << time_step_ << ", " << *differential_
362 << "}";
363 }
364
365 } // namespace crocoddyl
366