integrator_impl.hpp
1 /*
2  * This file is part of CasADi.
3  *
4  * CasADi -- A symbolic framework for dynamic optimization.
5  * Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
6  * KU Leuven. All rights reserved.
7  * Copyright (C) 2011-2014 Greg Horn
8  *
9  * CasADi is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 3 of the License, or (at your option) any later version.
13  *
14  * CasADi is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with CasADi; if not, write to the Free Software
21  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22  *
23  */
24 
25 
26 #ifndef CASADI_INTEGRATOR_IMPL_HPP
27 #define CASADI_INTEGRATOR_IMPL_HPP
28 
29 #include "integrator.hpp"
30 #include "oracle_function.hpp"
31 #include "plugin_interface.hpp"
32 #include "casadi_enum.hpp"
33 
35 
36 namespace casadi {
37 
41 struct CASADI_EXPORT IntegratorMemory : public OracleMemory {
42  // Work vectors, forward problem
43  double *q, *x, *z, *p, *u, *e, *edot, *old_e, *xdot, *zdot;
44  // Work vectors, backward problem
45  double *adj_x, *adj_z, *adj_p, *adj_q;
46  // Temporary work vectors of length max(nx + nz, nrx, nrz)
47  double *tmp1, *tmp2;
48  // Current control interval
49  casadi_int k;
50  // Current time
51  double t;
52  // Next time to be visited by the integrator
53  double t_next;
54  // Time not to be exceeded by during integrator integration
55  double t_stop;
56  // Time at the beginning of the current control interval
57  double t_start;
58  // Next output time
59  double t_next_out;
60  // Next stop time due to step change in input
61  double t_step;
62  // Which events have been triggered
63  casadi_int *event_triggered;
64  // Do we need to reset the solver?
65  bool reset_solver;
66  // Number of root-finding iterations for a single event
67  casadi_int event_iter;
68  // Number of events encountered thus far
69  casadi_int num_events;
70  // Index of event last triggered
71  casadi_int event_index;
72 };
73 
75 struct CASADI_EXPORT SpForwardMem {
76  const bvec_t** arg;
77  bvec_t** res;
78  casadi_int* iw;
79  bvec_t* w;
80 };
81 
83 struct CASADI_EXPORT SpReverseMem {
84  bvec_t** arg;
85  bvec_t** res;
86  casadi_int* iw;
87  bvec_t* w;
88 };
89 
97 class CASADI_EXPORT
98 Integrator : public OracleFunction, public PluginInterface<Integrator> {
99  public:
103  Integrator(const std::string& name, const Function& oracle,
104  double t0, const std::vector<double>& tout);
105 
109  ~Integrator() override=0;
110 
112 
115  size_t get_n_in() override { return INTEGRATOR_NUM_IN;}
116  size_t get_n_out() override { return INTEGRATOR_NUM_OUT;}
118 
120 
123  Sparsity get_sparsity_in(casadi_int i) override;
124  Sparsity get_sparsity_out(casadi_int i) override;
126 
128 
131  std::string get_name_in(casadi_int i) override { return integrator_in(i);}
132  std::string get_name_out(casadi_int i) override { return integrator_out(i);}
134 
138  int init_mem(void* mem) const override;
139 
141 
144  static const Options options_;
145  const Options& get_options() const override { return options_;}
147 
151  void init(const Dict& opts) override;
152 
156  void set_work(void* mem, const double**& arg, double**& res,
157  casadi_int*& iw, double*& w) const override;
158 
160  virtual Function create_advanced(const Dict& opts);
161 
162  virtual MX algebraic_state_init(const MX& x0, const MX& z0) const { return z0; }
163  virtual MX algebraic_state_output(const MX& Z) const { return Z; }
164 
165  // Set the quadrature states
166  void set_q(IntegratorMemory* m, const double* q) const;
167 
168  // Set the differential states
169  void set_x(IntegratorMemory* m, const double* x) const;
170 
171  // Set the algebraic variables
172  void set_z(IntegratorMemory* m, const double* z) const;
173 
174  // Set the parameters
175  void set_p(IntegratorMemory* m, const double* p) const;
176 
177  // Set the controls
178  void set_u(IntegratorMemory* m, const double* u) const;
179 
180  // Get the quadrature states
181  void get_q(IntegratorMemory* m, double* q) const;
182 
183  // Get the differential states
184  void get_x(IntegratorMemory* m, double* x) const;
185 
186  // Get the algebraic variables
187  void get_z(IntegratorMemory* m, double* z) const;
188 
192  virtual void reset(IntegratorMemory* mem, bool first_call) const {}
193 
197  casadi_int next_stop(casadi_int k, const double* u) const;
198 
202  int calc_edot(IntegratorMemory* m) const;
203 
207  int predict_events(IntegratorMemory* m) const;
208 
212  int trigger_event(IntegratorMemory* m, casadi_int* ind) const;
213 
217  int advance(IntegratorMemory* m) const;
218 
222  virtual int advance_noevent(IntegratorMemory* mem) const = 0;
223 
227  virtual void resetB(IntegratorMemory* mem) const = 0;
228 
232  casadi_int next_stopB(casadi_int k, const double* u) const;
233 
237  virtual void impulseB(IntegratorMemory* mem,
238  const double* adj_x, const double* adj_z, const double* adj_q) const = 0;
239 
243  virtual void retreat(IntegratorMemory* mem, const double* u,
244  double* adj_x, double* adj_p, double* adj_u) const = 0;
245 
249  int eval(const double** arg, double** res, casadi_int* iw, double* w, void* mem) const override;
250 
254  virtual void print_stats(IntegratorMemory* mem) const {}
255 
257  int fdae_sp_forward(SpForwardMem* m, const bvec_t* x,
258  const bvec_t* p, const bvec_t* u, bvec_t* ode, bvec_t* alg) const;
259 
261  int fquad_sp_forward(SpForwardMem* m, const bvec_t* x, const bvec_t* z,
262  const bvec_t* p, const bvec_t* u, bvec_t* quad) const;
263 
265  int bdae_sp_forward(SpForwardMem* m, const bvec_t* x, const bvec_t* z,
266  const bvec_t* p, const bvec_t* u, const bvec_t* adj_ode, const bvec_t* adj_quad,
267  bvec_t* adj_x, bvec_t* adj_z) const;
268 
270  int bquad_sp_forward(SpForwardMem* m, const bvec_t* x, const bvec_t* z,
271  const bvec_t* p, const bvec_t* u, const bvec_t* adj_ode, const bvec_t* adj_alg,
272  const bvec_t* adj_quad, bvec_t* adj_p, bvec_t* adj_u) const;
273 
277  int sp_forward(const bvec_t** arg, bvec_t** res,
278  casadi_int* iw, bvec_t* w, void* mem) const override;
279 
281  int fdae_sp_reverse(SpReverseMem* m, bvec_t* x,
282  bvec_t* p, bvec_t* u, bvec_t* ode, bvec_t* alg) const;
283 
285  int fquad_sp_reverse(SpReverseMem* m, bvec_t* x, bvec_t* z,
286  bvec_t* p, bvec_t* u, bvec_t* quad) const;
287 
289  int bdae_sp_reverse(SpReverseMem* m, bvec_t* x, bvec_t* z,
290  bvec_t* p, bvec_t* u, bvec_t* adj_ode, bvec_t* adj_quad,
291  bvec_t* adj_x, bvec_t* adj_z) const;
292 
294  int bquad_sp_reverse(SpReverseMem* m, bvec_t* x, bvec_t* z,
295  bvec_t* p, bvec_t* u, bvec_t* adj_ode, bvec_t* adj_alg, bvec_t* adj_quad,
296  bvec_t* adj_p, bvec_t* adj_u) const;
297 
301  int sp_reverse(bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w, void* mem) const override;
302 
305  bool has_spfwd() const override { return true;}
306  bool has_sprev() const override { return true;}
308 
310 
313  Function get_forward(casadi_int nfwd, const std::string& name,
314  const std::vector<std::string>& inames,
315  const std::vector<std::string>& onames,
316  const Dict& opts) const override;
317  bool has_forward(casadi_int nfwd) const override { return true;}
319 
321 
324  Function get_reverse(casadi_int nadj, const std::string& name,
325  const std::vector<std::string>& inames,
326  const std::vector<std::string>& onames,
327  const Dict& opts) const override;
328  bool has_reverse(casadi_int nadj) const override { return ne_ == 0;}
330 
334  virtual Dict getDerivativeOptions(bool fwd) const;
335 
337 
340  template<typename MatType> Function get_forward_dae(const std::string& name) const;
341  Function augmented_dae() const;
343 
345  static bool all_zero(const double* v, casadi_int n);
346 
348  Sparsity sp_jac_aug(const Sparsity& J, const Sparsity& J1) const;
349 
351  Sparsity sp_jac_dae();
352 
354  Sparsity sp_jac_rdae();
355 
357  Sparsity sp_jac_dae_, sp_jac_rdae_;
358 
360  inline casadi_int nt() const { return tout_.size();}
361 
363 
366  enum DaeOut { DAE_ODE, DAE_ALG, DAE_NUM_OUT};
367  static std::vector<std::string> dae_out() { return {"ode", "alg"}; }
368  enum QuadOut { QUAD_QUAD, QUAD_NUM_OUT};
369  static std::vector<std::string> quad_out() { return {"quad"}; }
370  enum BDynIn { BDYN_T, BDYN_X, BDYN_Z, BDYN_P, BDYN_U,
371  BDYN_OUT_ODE, BDYN_OUT_ALG, BDYN_OUT_QUAD, BDYN_OUT_ZERO,
372  BDYN_ADJ_ODE, BDYN_ADJ_ALG, BDYN_ADJ_QUAD, BDYN_ADJ_ZERO, BDYN_NUM_IN};
373  static std::string bdyn_in(casadi_int i);
374  static std::vector<std::string> bdyn_in();
375  enum BDynOut { BDYN_ADJ_T, BDYN_ADJ_X, BDYN_ADJ_Z, BDYN_ADJ_P, BDYN_ADJ_U, BDYN_NUM_OUT};
376  static std::string bdyn_out(casadi_int i);
377  static std::vector<std::string> bdyn_out();
378  enum DAEBOut { BDAE_ADJ_X, BDAE_ADJ_Z, BDAE_NUM_OUT};
379  static std::vector<std::string> bdae_out() { return {"adj_x", "adj_z"}; }
380  enum QuadBOut { BQUAD_ADJ_P, BQUAD_ADJ_U, BQUAD_NUM_OUT};
381  static std::vector<std::string> bquad_out() { return {"adj_p", "adj_u"}; }
383 
385  double t0_;
386 
388  std::vector<double> tout_;
389 
391  casadi_int nfwd_, nadj_;
392 
394  Function rdae_;
395 
397  casadi_int nx_, nz_, nq_, nx1_, nz1_, nq1_;
398 
400  casadi_int nrx_, nrz_, nrq_, nuq_, nrx1_, nrz1_, nrq1_, nuq1_;
401 
403  casadi_int np_, nrp_, np1_, nrp1_;
404 
406  casadi_int nu_, nu1_;
407 
409  casadi_int ne_;
410 
412  casadi_int ntmp_;
413 
414  // Nominal values for states
415  std::vector<double> nom_x_, nom_z_;
416 
418  Dict augmented_options_;
419 
421  Dict opts_;
422 
424  bool print_stats_;
425 
427  Function transition_;
428 
430  casadi_int max_event_iter_;
431 
433  casadi_int max_events_;
434 
436  double event_tol_;
437 
439  double event_acceptable_tol_;
440 
441  // Creator function for internal class
442  typedef Integrator* (*Creator)(const std::string& name, const Function& oracle,
443  double t0, const std::vector<double>& tout);
444 
445  // No static functions exposed
446  struct Exposed{ };
447 
449  static std::map<std::string, Plugin> solvers_;
450 
451 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
452  static std::mutex mutex_solvers_;
453 #endif // CASADI_WITH_THREADSAFE_SYMBOLICS
454 
456  static const std::string infix_;
457 
459  template<typename XType>
460  static Function map2oracle(const std::string& name, const std::map<std::string, XType>& d);
461 
465  void serialize_body(SerializingStream &s) const override;
469  void serialize_type(SerializingStream &s) const override;
470 
474  static ProtoFunction* deserialize(DeserializingStream& s);
475 
479  std::string serialize_base_function() const override { return "Integrator"; }
480 
482  static bool grid_in(casadi_int i);
483 
485  static bool grid_out(casadi_int i);
486 
488  static casadi_int adjmap_out(casadi_int i);
489 
490  protected:
494  explicit Integrator(DeserializingStream& s);
495 };
496 
498 enum StepIn {
500  STEP_T,
502  STEP_H,
504  STEP_X0,
506  STEP_V0,
508  STEP_P,
510  STEP_U,
512  STEP_NUM_IN
513 };
514 
516 enum StepOut {
518  STEP_XF,
520  STEP_VF,
522  STEP_QF,
524  STEP_NUM_OUT
525 };
526 
528 enum BStepIn {
529  BSTEP_T,
530  BSTEP_H,
531  BSTEP_X0,
532  BSTEP_V0,
533  BSTEP_P,
534  BSTEP_U,
535  BSTEP_OUT_XF,
536  BSTEP_OUT_VF,
537  BSTEP_OUT_QF,
538  BSTEP_ADJ_XF,
539  BSTEP_ADJ_VF,
540  BSTEP_ADJ_QF,
541  BSTEP_NUM_IN
542 };
543 
545 enum BStepOut {
546  BSTEP_ADJ_T,
547  BSTEP_ADJ_H,
548  BSTEP_ADJ_X0,
549  BSTEP_ADJ_V0,
550  BSTEP_ADJ_P,
551  BSTEP_ADJ_U,
552  BSTEP_NUM_OUT
553 };
554 
555 struct CASADI_EXPORT FixedStepMemory : public IntegratorMemory {
557  double *v, *v_prev, *q_prev;
558 
560  double *rv, *adj_u, *adj_p_prev, *adj_u_prev;
561 
563  double *x_tape, *v_tape;
564 };
565 
566 class CASADI_EXPORT FixedStepIntegrator : public Integrator {
567  public:
568 
570  explicit FixedStepIntegrator(const std::string& name, const Function& dae,
571  double t0, const std::vector<double>& tout);
572 
574  ~FixedStepIntegrator() override;
575 
577 
580  static const Options options_;
581  const Options& get_options() const override { return options_;}
583 
585  void init(const Dict& opts) override;
586 
590  void set_work(void* mem, const double**& arg, double**& res,
591  casadi_int*& iw, double*& w) const override;
592 
594  Function create_advanced(const Dict& opts) override;
595 
599  void* alloc_mem() const override { return new FixedStepMemory();}
600 
604  int init_mem(void* mem) const override;
605 
609  void free_mem(void *mem) const override { delete static_cast<FixedStepMemory*>(mem);}
610 
612  virtual void setup_step() = 0;
613 
617  void reset(IntegratorMemory* mem, bool first_call) const override;
618 
622  int advance_noevent(IntegratorMemory* mem) const override;
623 
625  void resetB(IntegratorMemory* mem) const override;
626 
628  void impulseB(IntegratorMemory* mem,
629  const double* adj_x, const double* adj_z, const double* adj_q) const override;
630 
634  void retreat(IntegratorMemory* mem, const double* u,
635  double* adj_x, double* adj_p, double* adj_u) const override;
636 
638  void stepF(FixedStepMemory* m, double t, double h,
639  const double* x0, const double* v0, double* xf, double* vf, double* qf) const;
640 
642  void stepB(FixedStepMemory* m, double t, double h,
643  const double* x0, const double* xf, const double* vf,
644  const double* adj_xf, const double* rv0,
645  double* adj_x0, double* adj_p, double* adj_u) const;
646 
647  // Target number of finite elements
648  casadi_int nk_target_;
649 
650  // Number of steps per control interval
651  std::vector<casadi_int> disc_;
652 
654  casadi_int nv_, nv1_, nrv_, nrv1_;
655 
659  void serialize_body(SerializingStream &s) const override;
660 
661 protected:
665  explicit FixedStepIntegrator(DeserializingStream& s);
666 };
667 
668 class CASADI_EXPORT ImplicitFixedStepIntegrator : public FixedStepIntegrator {
669  public:
670 
672  explicit ImplicitFixedStepIntegrator(const std::string& name, const Function& dae,
673  double t0, const std::vector<double>& tout);
674 
676  ~ImplicitFixedStepIntegrator() override;
677 
679 
682  static const Options options_;
683  const Options& get_options() const override { return options_;}
685 
687  void init(const Dict& opts) override;
688 
692  void serialize_body(SerializingStream &s) const override;
693 
694 protected:
698  explicit ImplicitFixedStepIntegrator(DeserializingStream& s);
699 };
700 
701 } // namespace casadi
703 
704 #endif // CASADI_INTEGRATOR_IMPL_HPP
casadi::integrator_in
CASADI_EXPORT std::vector< std::string > integrator_in()
Get input scheme of integrators.
casadi::integrator_out
CASADI_EXPORT std::vector< std::string > integrator_out()
Get integrator output scheme of integrators.
casadi
The casadi namespace.
Definition: archiver.hpp:32
casadi::Dict
GenericType::Dict Dict
C++ equivalent of Python's dict or MATLAB's struct.
Definition: generic_type.hpp:258