Program Listing for File ODESolver.h
↰ Return to documentation for file (src/ompl/control/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 ¤t, 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