ODESolver.h

/*********************************************************************
* Software License Agreement (BSD License)
*
*  Copyright (c) 2011, Rice University
*  All rights reserved.
*
*  Redistribution and use in source and binary forms, with or without
*  modification, are permitted provided that the following conditions
*  are met:
*
*   * Redistributions of source code must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*   * Redistributions in binary form must reproduce the above
*     copyright notice, this list of conditions and the following
*     disclaimer in the documentation and/or other materials provided
*     with the distribution.
*   * Neither the name of the Rice University nor the names of its
*     contributors may be used to endorse or promote products derived
*     from this software without specific prior written permission.
*
*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
*  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
*  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
*  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
*  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
*  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
*  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
*  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
*  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
*  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
*  POSSIBILITY OF SUCH DAMAGE.
*********************************************************************/
 
/* Author: Ryan Luna */
 
#ifndef OMPL_CONTROL_ODESOLVER_
#define OMPL_CONTROL_ODESOLVER_
 
#include "ompl/control/Control.h"
#include "ompl/control/SpaceInformation.h"
#include "ompl/control/StatePropagator.h"
#include "ompl/util/Console.h"
#include "ompl/util/ClassForward.h"
 
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
#include <boost/numeric/odeint/stepper/generation/make_controlled.hpp>
#include <boost/numeric/odeint/integrate/integrate_const.hpp>
#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
namespace odeint = boost::numeric::odeint;
#include <functional>
#include <cassert>
#include <utility>
#include <vector>
 
namespace ompl
{
    namespace control
    {
        OMPL_CLASS_FORWARD(ODESolver);
 
 
        class ODESolver
        {
        public:
            using StateType = std::vector<double>;
 
            using ODE = std::function<void(const StateType &, const Control *, StateType &)>;
 
            using PostPropagationEvent =
                std::function<void(const base::State *, const Control *, double, base::State *)>;
 
            ODESolver(SpaceInformationPtr si, ODE ode, double intStep)
              : si_(std::move(si)), ode_(std::move(ode)), intStep_(intStep)
            {
            }
 
            virtual ~ODESolver() = default;
 
            void setODE(const ODE &ode)
            {
                ode_ = ode;
            }
 
            double getIntegrationStepSize() const
            {
                return intStep_;
            }
 
            void setIntegrationStepSize(double intStep)
            {
                intStep_ = intStep;
            }
 
            const SpaceInformationPtr &getSpaceInformation() const
            {
                return si_;
            }
 
            static StatePropagatorPtr getStatePropagator(ODESolverPtr solver,
                                                         const PostPropagationEvent &postEvent = nullptr)
            {
                class ODESolverStatePropagator : public StatePropagator
                {
                public:
                    ODESolverStatePropagator(const ODESolverPtr& solver, PostPropagationEvent pe)
                      : StatePropagator(solver->si_), solver_(solver), postEvent_(std::move(pe))
                    {
                        if (solver == nullptr)
                            OMPL_ERROR("ODESolverPtr does not reference a valid ODESolver object");
                    }
 
                    void propagate(const base::State *state, const Control *control, double duration,
                                   base::State *result) const override
                    {
                        ODESolver::StateType reals;
                        si_->getStateSpace()->copyToReals(reals, state);
                        solver_->solve(reals, control, duration);
                        si_->getStateSpace()->copyFromReals(result, reals);
 
                        if (postEvent_)
                            postEvent_(state, control, duration, result);
                    }
 
                protected:
                    ODESolverPtr solver_;
                    ODESolver::PostPropagationEvent postEvent_;
                };
                return std::make_shared<ODESolverStatePropagator>(std::move(solver), postEvent);
            }
 
        protected:
            virtual void solve(StateType &state, const Control *control, double duration) const = 0;
 
            const SpaceInformationPtr si_;
 
            ODE ode_;
 
            double intStep_;
 
            // Functor used by the boost::numeric::odeint stepper object
            struct ODEFunctor
            {
                ODEFunctor(ODE o, const Control *ctrl) : ode(std::move(o)), control(ctrl)
                {
                }
 
                // boost::numeric::odeint will callback to this method during integration to evaluate the system
                void operator()(const StateType &current, StateType &output, double /*time*/)
                {
                    ode(current, control, output);
                }
 
                ODE ode;
                const Control *control;
            };
        };
 
        template <class Solver = odeint::runge_kutta4<ODESolver::StateType>>
        class ODEBasicSolver : public ODESolver
        {
        public:
            ODEBasicSolver(const SpaceInformationPtr &si, const ODESolver::ODE &ode, double intStep = 1e-2)
              : ODESolver(si, ode, intStep)
            {
            }
 
        protected:
            void solve(StateType &state, const Control *control, double duration) const override
            {
                Solver solver;
                ODESolver::ODEFunctor odefunc(ode_, control);
                odeint::integrate_const(solver, odefunc, state, 0.0, duration, intStep_);
            }
        };
 
        template <class Solver = odeint::runge_kutta_cash_karp54<ODESolver::StateType>>
        class ODEErrorSolver : public ODESolver
        {
        public:
            ODEErrorSolver(const SpaceInformationPtr &si, const ODESolver::ODE &ode, double intStep = 1e-2)
              : ODESolver(si, ode, intStep)
            {
            }
 
            ODESolver::StateType getError()
            {
                return error_;
            }
 
        protected:
            void solve(StateType &state, const Control *control, double duration) const override
            {
                ODESolver::ODEFunctor odefunc(ode_, control);
 
                if (error_.size() != state.size())
                    error_.assign(state.size(), 0.0);
 
                Solver solver;
                solver.adjust_size(state);
 
                double time = 0.0;
                while (time < duration + std::numeric_limits<float>::epsilon())
                {
                    solver.do_step(odefunc, state, time, intStep_, error_);
                    time += intStep_;
                }
            }
 
            mutable ODESolver::StateType error_;
        };
 
        template <class Solver = odeint::runge_kutta_cash_karp54<ODESolver::StateType>>
        class ODEAdaptiveSolver : public ODESolver
        {
        public:
            ODEAdaptiveSolver(const SpaceInformationPtr &si, const ODESolver::ODE &ode, double intStep = 1e-2)
              : ODESolver(si, ode, intStep), maxError_(1e-6), maxEpsilonError_(1e-7)
            {
            }
 
            double getMaximumError() const
            {
                return maxError_;
            }
 
            void setMaximumError(double error)
            {
                maxError_ = error;
            }
 
            double getMaximumEpsilonError() const
            {
                return maxEpsilonError_;
            }
 
            void setMaximumEpsilonError(double error)
            {
                maxEpsilonError_ = error;
            }
 
        protected:
            void solve(StateType &state, const Control *control, double duration) const override
            {
                ODESolver::ODEFunctor odefunc(ode_, control);
                auto solver = make_controlled(1.0e-6, 1.0e-6, Solver());
                odeint::integrate_adaptive(solver, odefunc, state, 0.0, duration, intStep_);
            }
 
            double maxError_;
 
            double maxEpsilonError_;
        };
    }
}
 
#endif