Problem.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                                      *
 *                                                                                               *
 * The Sandals project is distributed under the BSD 2-Clause License.                            *
 *                                                                                               *
 * Davide Stocco                                                               Enrico Bertolazzi *
 * University of Trento                                                     University of Trento *
 * e-mail: davide.stocco@unitn.it                             e-mail: enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SANDALS_PROBLEM_HH
#define SANDALS_PROBLEM_HH
 
#include <Sandals/RungeKutta.hh>
 
namespace Sandals
{
 
  /*\
   |.  ____            _     _
   |. |  _ \ _ __ ___ | |__ | | ___ _ __ ___
   |. | |_) | '__/ _ \| '_ \| |/ _ \ '_ ` _ \
   |. |  __/| | | (_) | |_) | |  __/ | | | | |
   |. |_|   |_|  \___/|_.__/|_|\___|_| |_| |_|
   |.
  \*/

  template <typename Real, Integer N, Integer M, typename Integrator>
  class Problem
  {
  public:
    SANDALS_BASIC_CONSTANTS(Real) 
    const Real SQRT_EPSILON{std::sqrt(EPSILON)};  \
 
    using SystemPtr = typename Implicit<Real, N, M>::Pointer; 
    using IntegratorPtr = std::shared_ptr<Integrator>; 
    using SolutionPtr = std::shared_ptr<Solution<Real, N, M>>; 
 
    using VectorX = typename Integrator::VectorX;
    using MatrixJX = typename Integrator::MatrixJX;
    using VectorF = typename Implicit<Real, N, M>::VectorF;
    using MatrixJF = typename Implicit<Real, N, M>::MatrixJF;
 
  private:
    std::string   m_name{"(undefined)"}; 
    SystemPtr     m_system; 
    IntegratorPtr m_integrator; 
    SolutionPtr   m_solution; 
 
    bool    m_verbose{false};                   
    Real    m_tolerance{EPSILON_HIGH}; 
    Integer m_max_iterations{100};     
 
  public:
    Problem(std::string t_name, SystemPtr t_system, IntegratorPtr t_integrator)
      : m_name(t_name), m_integrator(t_integrator)
    {
      this->m_integrator->system(t_system);
      this->m_solution = std::make_shared<Solution<Real, N, M>>();
    }

    virtual ~Problem() {}

    std::string name() const {return this->m_name;}

    void name(std::string const t_name) {this->m_name = t_name;}

    SystemPtr system() {return this->m_integrator->system();}

    void system(SystemPtr t_system) {this->m_integrator->system(t_system);}

    IntegratorPtr integrator() {return this->m_integrator;}

    void integrator(IntegratorPtr t_integrator) {this->m_integrator = t_integrator;}

    SolutionPtr const solution() const {return this->m_solution;}

    bool verbose_mode() {return this->m_verbose;}

    void verbose_mode(bool t_verbose) {
      this->m_verbose = t_verbose;
      this->m_integrator->verbose_mode(t_verbose);
    }

    void enable_verbose_mode() {this->verbose_mode(true);}

    void disable_verbose_mode() {this->verbose_mode(false);}

    Real tolerance() {return this->m_tolerance;}

    void tolerance(Real const t_tolerance)
    {this->m_tolerance = t_tolerance;}

    Integer & max_iterations() {return this->m_max_iterations;}

    void max_iterations(Integer const t_max_iterations)
    {this->m_max_iterations = t_max_iterations;}

    virtual VectorF b(VectorF const & x_ini, VectorF const & x_end) const = 0;

    virtual MatrixJF Jb_x_ini(VectorF const & x_ini, VectorF const & x_end) const = 0;

    virtual MatrixJF Jb_x_end(VectorF const & x_ini, VectorF const & x_end) const = 0;

    bool single_shooting(VectorX const & t_mesh, VectorF const & ics)
    {
      using ShootingF = Eigen::Matrix<Real, 2*N, 1>;
      using ShootingJF = Eigen::Matrix<Real, 2*N, 2*N>;
 
      #define CMD "Sandals::Problem::single_shooting(...): "
 
      // Temporary variables
      VectorF bcs, x_ini, x_end;
      ShootingF x_sol, x_step, b;
      ShootingJF A(ShootingJF::Zero());
      A.template block<N, N>(0, N) = MatrixJX::Identity();
 
      // Initialize the guess and the solution
      x_sol << ics, ics;
      this->m_solution->clear();
 
      // Solve the boundary value problem using a linearized Newton method
      MatrixJX Jx(MatrixJX::Identity());
      Eigen::FullPivLU<ShootingJF> lu;
      for (Integer i{0}; i < this->m_max_iterations; ++i) {
 
        /* Single shooting method
                  [A]           {x}   =          {f}
         /   -Jx       I     \ /    \   / x_end - x_sol_end \
         |                   | | dx | = |                   |
         \ Jb_x_ini Jb_x_end / \    /   \        -b         /
        */
 
        // Integrate the system and propagate the Jacobian
        if (!this->m_integrator->solve(t_mesh, x_sol.template head<N>(), *this->m_solution, Jx)) {
          SANDALS_ERROR(CMD "failed to integrate the system " << this->m_integrator->name() << ".");
          return false;
        }
 
        // Retrieve the initial and final states
        x_ini = x_sol.template head<N>();
        x_end = this->m_solution->x.col(this->m_solution->t.size() - 1);
 
        // Evaluate the residual of the boundary conditions
        bcs = this->b(x_ini, x_end);
 
        // Print the iteration info
        if (this->m_verbose) {
          std::cout << "Iteration " << i << ": |b| = " << bcs.norm() << std::endl
                    << "  x(" << t_mesh.template head<1>() << ") = " << x_ini.transpose() << std::endl
                    << "  x(" << t_mesh.template tail<1>() << ") = " << x_end.transpose() << std::endl;
        }
 
        // Check if the solution is found (i.e., if the boundary conditions are satisfied)
        if (bcs.norm() < this->m_tolerance) {return true;}
 
        // Build the linear system
        b.template head<N>() = x_end - x_sol.template tail<N>();
        b.template tail<N>() = -bcs;
        A.template block<N, N>(0, 0) = -Jx;
        A.template block<N, N>(N, 0) = this->Jb_x_ini(x_ini, x_end);
        A.template block<N, N>(N, N) = this->Jb_x_end(x_ini, x_end);
 
        // Compute the solution of the linear system
        lu.compute(A);
        SANDALS_ASSERT(lu.rank() == 2*N, CMD "singular Jacobian detected.");
        x_step = lu.solve(b);
 
        // Update the solution
        x_sol += x_step;
      }
 
      // If the loop completes without returning, indicate failure
      if (this->m_verbose) {SANDALS_WARNING(CMD "maximum number of iterations reached.");}
      return false;
 
      #undef CMD
    }
 
  }; // class Problem
 
} // namespace Sandals
 
#endif // SANDALS_PROBLEM_HH