Solvers integrator step size#788

doc/hybrid_systems.dox

@@ -42,7 +42,7 @@ A <em>hybrid automaton</em> is a tuple \f$H=(Q,E,\gamma,n_q,f_q,I_q,r_{q,e},G_{q
 
  - For each \f$q\in Q\f$,
    - \f$n_q\f$ defines the <em>dimension</em> of the continuous state space \f$\R^{n_q}\f$.
-   - \f$f_q:\R^{n_q}\rightarrow \R^{n_q}\f$ defines the <em>continuous dynamics</em> by \f$\dt{x}=f_q(x)\f$.
+   - \f$f_q:\R^{n_q}\rightarrow \R^{n_q}\f$ defines the <em>continuous dynamics</em> by \f$\dot{x}=f_q(x)\f$.
    - \f$I_q\subset \R^{n_q}\f$ is the <em>invariant</em> which must hold throughout the evolution in mode \f$q\f$.
 
  - For each \f$(q,e)\in\mathrm{dom}(\gamma)\subset Q\times E\f$,
@@ -51,7 +51,7 @@ A <em>hybrid automaton</em> is a tuple \f$H=(Q,E,\gamma,n_q,f_q,I_q,r_{q,e},G_{q
 
 The <em>standard operational semantics</em> of a hybrid automaton \f$H\f$ is as follows.
  - The <em>state</em> of the system is a pair \f$(q,x)\f$ with \f$q\in Q\f$ and \f$x\in\R^{n_q}\f$.
- - The state may evolve according to the differential equation \f$\dt{x}=f_q(x)\f$ as long as \f$x\in I_q\f$.
+ - The state may evolve according to the differential equation \f$\dot{x}=f_q(x)\f$ as long as \f$x\in I_q\f$.
  - The state may evolve according to the discrete transition \f$(q',x')=(\gamma(q,e),r_{q,e}(x))\f$ as long as \f$x\in R_{q,e}\f$.
 
 This definition has some disadvantages:
@@ -79,7 +79,7 @@ More generally, there are many different ways of specifying when continuous evol
 The hybrid automaton interface in %Ariadne is essentially equivalent to the enhanced mathematical definition given above. The dynamic and resets are defined using the VectorFunction class. The guards and invariants are defined using the class ScalarFunction class, with invariants \f$i(x)\leq0\f$ and guards \f$g(x)\geq0\f$.
 A <em>DiscreteMode</em> is a triple \f$(q,f_q,I_q)\f$ where
  - \f$q\f$ is a DiscreteState.
- - \f$f_q\f$ is a VectorFunction giving the dynamic \f$\dt{x}=f_q(x)\f$.
+ - \f$f_q\f$ is a VectorFunction giving the dynamic \f$\dot{x}=f_q(x)\f$.
  - \f$I_q\f$ is a List of ScalarFunction with each \f$i_j:\R^{n_q}\rightarrow \R\f$ giving the invariants \f$i_j(x)\lesssim 0\f$.
 
 A <em>DiscreteTransition</em> is a triple \f$(e,s,t,r,g,u)\f$ where
@@ -100,7 +100,7 @@ A <em>HybridAutomaton</em> is a pair \f$(M,T)\f$ where
 
 A <em>Mode</em> is a tuple \f$(q,f_q)\f$ where
  - \f$q\f$ is a DiscreteLocation.
- - \f$f_q\f$ is a VectorFunction giving the dynamic \f$\dt{x}=f_q(x)\f$.
+ - \f$f_q\f$ is a VectorFunction giving the dynamic \f$\dot{x}=f_q(x)\f$.
 
 A <em>Guard</em> is a tuple \f$(q,e,c,k)\f$ where
  - \f$q\f$ is a DiscreteLocation.

doc/macros.js

@@ -3,6 +3,13 @@ MathJax.Hub.Config({
     Macros: {
 //      jot: "1ex",
 
+      D: "{\\mathrm{D}}",
+      Df: "{\\mathrm{D}f}",
+      Dy: "{\\mathrm{D}y}",
+      d: "{\\mathrm{d}}",
+      dt: "{\\mathrm{d}t}",
+      dx: "{\\mathrm{d}x}",
+      dy: "{\\mathrm{d}y}",
       B: "{\\mathbb{B}}",
       K: "{\\mathbb{K}}",
       S: "{\\mathbb{S}}",
@@ -20,7 +27,6 @@ MathJax.Hub.Config({
       fto: "{\\longrightarrow}",
       interval: ["{[#1]}",1],
       ivl: ["{[#1]}",1],
-      dt: ["{\\dot{#1}}",1],
       unl: ["{\\underline{#1}}",1],
       ovl: ["{\\overline{#1}}",1],
       der: ["{\\dot{#1}}",1],

doc/macros.sty

@@ -2,6 +2,13 @@
 
 \newcommand{\interval}[1]{{[#1]}}
 \newcommand{\ivl}[1]{{[#1]}}
+\newcommand{\D}{\mathrm{D}}
+\newcommand{\Df}{\mathrm{D}f}
+\newcommand{\Dy}{\mathrm{D}y}
+\newcommand{\d}{\mathrm{d}}
+\newcommand{\dt}{\mathrm{d}t}
+\newcommand{\dx}{\mathrm{d}x}
+\newcommand{\dy}{\mathrm{d}y}
 \newcommand{\B}{\mathbb{N}}
 \newcommand{\K}{\mathbb{K}}
 \newcommand{\S}{\mathbb{S}}
@@ -15,7 +22,6 @@
 \newcommand{\Y}{\mathbb{Y}}
 \newcommand{\A}{\mathbb{A}}
 \newcommand{\seq}[1]{\vec{#1}}
-\newcommand{\dt}[1]{\dot{#1}}
 \newcommand{\fto}{\longrightarrow}
 \newcommand{\pfto}{\dashrightarrow}
 

doc/overview.dox

@@ -36,7 +36,7 @@ namespace Ariadne {
 The main functionality of %Ariadne is to perform reachability analysis and verification of hybrid systems.
 A <em>hybrid system</em> is a dynamic system in which continuous evolution is interspersed with discrete jumps.
 The state space of a hybrid system is of the form \f$X=\bigcup_{q\in Q} \{q\}\times \R^{n_q}\f$ with \f$Q\f$ a finite set of <em>locations</em>.
-The evolution of a hybrid system proceeds via a differential equation \f$\dt{x}(t) = f_q(x(t))\f$ in location \f$q\f$
+The evolution of a hybrid system proceeds via a differential equation \f$\dot{x}(t) = f_q(x(t))\f$ in location \f$q\f$
 until a <em>guard condition</em> \f$g_e(x)\geq 0\f$ is satisfied, in which case the <em>event</em> \f$e\f$ occurs,
 the location jumps to a target location \f$q'=\tau(q,e)\f$ and the continuous state jumps to \f$x'=r_{q,e}(x)\f$.
 
@@ -94,8 +94,8 @@ which is a valid interval extension of \f$f\f$ on a box \f$D\subset U\f$.
 
 The main operations which %Ariadne can perform are
  - Solution of differential equations:
-   Given a differential equation of the form \f$\dt{x}=f(x)\f$,
-   computes the flow \f$\phi(x,t)\f$ satisfying \f$\phi(x,0)=x\f$ and \f$\dt{\phi}(x,t) = f(\phi(x,t))\f$.
+   Given a differential equation of the form \f$\dot{x}=f(x)\f$,
+   computes the flow \f$\phi(x,t)\f$ satisfying \f$\phi(x,0)=x\f$ and \f$\dot{\phi}(x,t) = f(\phi(x,t))\f$.
  - Solution of (parameterised) algebraic equations:
    Given an algebraic equation \f$f(a,x)=0\f$, find a
    function \f$h(a)\f$ satisfying \f$f(a,h(a))=0\f$.

doc/python_module_demonstrations.dox

@@ -314,7 +314,7 @@ If \f$f:D_1\times D_2\rightarrow \R^n\f$ with \f$D_2\subset\R^n\f$ then the oper
 \subsection python_differential_equation_solvers_demonstration Solving Differential Equations
 
 We can solve differential equations using the operator \c flow of an \a Integrator class, such as \ref TaylorPicardIntegrator or \ref GradedTaylorSeriesIntegrator.
-If \f$:\R^n\to\R^n\f$, \f$D\subset\R^n\f$ and \f$h>0\f$, then the operator \c flow(f,D,h) computes a function \f$\phi:D\times[0,h]\to\R^n\f$ satisfying \f$\dt{\phi}(x,t)=f(\phi(x,t))\f$ and \f$\phi(x,0)=x\f$.
+If \f$:\R^n\to\R^n\f$, \f$D\subset\R^n\f$ and \f$h>0\f$, then the operator \c flow(f,D,h) computes a function \f$\phi:D\times[0,h]\to\R^n\f$ satisfying \f$\dot{\phi}(x,t)=f(\phi(x,t))\f$ and \f$\phi(x,0)=x\f$.
 
 \skipline integrator=
 \skipline f=

python/source/solver_submodule.cpp

@@ -84,7 +84,9 @@ class IntegratorWrapper
   public:
     IntegratorInterface* clone() const {
         return this->get_override("clone")(); }
-    FlowStepModelType flow_step(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& D, StepSizeType& h) const {
+    FlowStepModelType flow_step(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& D) const {
+        return this->get_override("flow_step")(vf,D); }
+    FlowStepModelType flow_step(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& D, const StepSizeType& h) const {
         return this->get_override("flow_step")(vf,D,h); }
     FlowStepModelType flow_step(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& D, const StepSizeType& h, const UpperBoxType& B) const {
         return this->get_override("flow_step")(vf,D,h,B); }
@@ -135,7 +137,8 @@ Void export_integrators(pybind11::module& module)
     flow_model_class.def("__str__", &__cstr__<FlowModelType>);
 
     pybind11::class_<IntegratorInterface,IntegratorWrapper> integrator_interface_class(module,"IntegratorInterface");
-    integrator_interface_class.def("flow_step",(FlowStepModelType(IntegratorInterface::*)(const ValidatedVectorMultivariateFunction&, const ExactBoxType&, StepSizeType&)const)&IntegratorInterface::flow_step);
+    integrator_interface_class.def("flow_step",(FlowStepModelType(IntegratorInterface::*)(const ValidatedVectorMultivariateFunction&, const ExactBoxType&)const)&IntegratorInterface::flow_step);
+    integrator_interface_class.def("flow_step",(FlowStepModelType(IntegratorInterface::*)(const ValidatedVectorMultivariateFunction&, const ExactBoxType&, const StepSizeType&)const)&IntegratorInterface::flow_step);
     integrator_interface_class.def("flow_step",(FlowStepModelType(IntegratorInterface::*)(const ValidatedVectorMultivariateFunction&,const ExactBoxType&,const StepSizeType&,const UpperBoxType&)const)&IntegratorInterface::flow_step);
     integrator_interface_class.def("__str__", &__cstr__<IntegratorInterface>);
 

source/dynamics/differential_inclusion_evolver.cpp

@@ -212,7 +212,7 @@ auto DifferentialInclusionEvolver::orbit(RealVariablesBox const& initial, Real c
 
     auto initial_box = initial_ranges_to_box(initial);
 
-    StepSizeType hsug(_configuration->maximum_step_size());
+    StepSizeType hsug=_configuration->maximum_step_size();
 
     EulerBounder bounder;
 
@@ -252,7 +252,7 @@ auto DifferentialInclusionEvolver::orbit(RealVariablesBox const& initial, Real c
         UpperBoxType B;
         StepSizeType h;
 
-        CONCLOG_RUN_AT(1,std::tie(h,B)=bounder.compute(_system.function(),domx,_system.inputs(),hsug));
+        CONCLOG_RUN_AT(1,std::tie(h,B)=bounder.compute(_system.function(),domx,_system.inputs(),suggest(hsug)));
         CONCLOG_PRINTLN_AT(2,"flow bounds = "<<B<<" (using h = " << h << ")");
 
         TimeStepType new_t = lower_bound(t+h);

source/dynamics/vector_field_evolver.cpp

@@ -200,8 +200,8 @@ _process_timed_enclosure_step(WorkloadType::Access& workload,
 
     // Compute flow model
     IntegratorInterface const* integrator=this->_integrator.operator->();
-    StepSizeType step_size=maximum_step_size;
-    FlowStepModelType flow_model=integrator->flow_step(dynamic,current_set_bounds,step_size);
+    FlowStepModelType flow_model=integrator->flow_step(dynamic,current_set_bounds,suggest(maximum_step_size));
+    StepSizeType step_size = static_cast<StepSizeType>(flow_model.domain()[flow_model.argument_size()-1u].upper_bound());
     CONCLOG_PRINTLN("step_size = "<<step_size)
     CONCLOG_PRINTLN_AT(1,"flow_model = "<<flow_model)
     FlowStepModelType flow_step_model=partial_evaluate(flow_model,flow_model.domain().size()-1u,step_size);

source/hybrid/hybrid_evolver.cpp

@@ -548,12 +548,11 @@ _compute_flow(EffectiveVectorMultivariateFunction dynamic,
     // We then restrict to the time domain [0,h] since this can make evaluation
     // more accurate, and the time domain might be used explicitly for the domain
     // of the resulting set.
-    StepSizeType step_size=maximum_step_size;
-    ValidatedVectorMultivariateFunctionPatch flow_model=integrator.flow_step(dynamic,initial_box,step_size);
-
+    ValidatedVectorMultivariateFunctionPatch flow_model=integrator.flow_step(dynamic,initial_box,suggest(maximum_step_size));
     CONCLOG_PRINTLN_AT(1,"twosided_flow_model="<<flow_model);
     ExactBoxType flow_domain=flow_model.domain();
-    ARIADNE_ASSERT(step_size==flow_domain[flow_domain.size()-1u].upper_bound());
+    StepSizeType step_size=static_cast<StepSizeType>(flow_domain[flow_domain.size()-1u].upper_bound());
+    ARIADNE_ASSERT(step_size<=maximum_step_size);
     flow_domain[flow_domain.size()-1u]=ExactIntervalType(0,step_size);
     flow_model=restriction(flow_model,flow_domain);
     CONCLOG_PRINTLN_AT(1,"flow_model="<<flow_model);

source/solvers/bounder.cpp

@@ -31,29 +31,29 @@ namespace Ariadne {
 BounderBase::BounderBase(LipschitzTolerance lipschitz, MinimumStepSize minimum_step)
     : _lipschitz_tolerance(lipschitz), _minimum_step_size(minimum_step) { }
 
-Pair<StepSizeType,UpperBoxType> BounderBase::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& hsug) const {
+Pair<StepSizeType,UpperBoxType> BounderBase::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, Suggestion<StepSizeType> const& hsug) const {
     return this->compute(f,D,BoxDomainType(0u),hsug);
 }
 
-Pair<StepSizeType,UpperBoxType> BounderBase::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, StepSizeType const& hsug) const {
+Pair<StepSizeType,UpperBoxType> BounderBase::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, Suggestion<StepSizeType> const& hsug) const {
     return this->compute(f,D,t,BoxDomainType(0u),hsug);
 }
 
 EulerBounder::EulerBounder(LipschitzTolerance lipschitz, MinimumStepSize minimum_step) : BounderBase(lipschitz,minimum_step) { }
 
-Pair<StepSizeType,UpperBoxType> EulerBounder::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, BoxDomainType const& A, StepSizeType const& hsug) const {
+Pair<StepSizeType,UpperBoxType> EulerBounder::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, BoxDomainType const& A, Suggestion<StepSizeType> const& hsug) const {
     ARIADNE_PRECONDITION(f.result_size()==D.dimension());
     ARIADNE_PRECONDITION(f.argument_size()==D.dimension()+A.dimension());
     return this->_compute(f,D,0,A,hsug);
 }
 
-Pair<StepSizeType,UpperBoxType> EulerBounder::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, BoxDomainType const& A, StepSizeType const& hsug) const {
+Pair<StepSizeType,UpperBoxType> EulerBounder::compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, BoxDomainType const& A, Suggestion<StepSizeType> const& hsug) const {
     ARIADNE_PRECONDITION(f.result_size()==D.dimension());
     ARIADNE_PRECONDITION(f.argument_size()==D.dimension()+1u+A.dimension());
     return this->_compute(f,D,t,A,hsug);
 }
 
-Pair<StepSizeType,UpperBoxType> EulerBounder::_compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, BoxDomainType const& A, StepSizeType const& hsug) const {
+Pair<StepSizeType,UpperBoxType> EulerBounder::_compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, BoxDomainType const& A, Suggestion<StepSizeType> const& hsug) const {
     CONCLOG_SCOPE_CREATE;
     const PositiveFloatDP BOX_RADIUS_WIDENING=cast_positive(0.25_exact);
     const PositiveFloatDP NO_WIDENING=cast_positive(1.0_exact);
@@ -63,7 +63,7 @@ Pair<StepSizeType,UpperBoxType> EulerBounder::_compute(ValidatedVectorMultivaria
     const CounterType EXPANSION_STEPS=4;
     const CounterType REFINEMENT_STEPS=4;
 
-    StepSizeType h=hsug;
+    StepSizeType h=hsug.suggestion();
 
     FloatDPUpperBound lipschitz = norm(f.jacobian(Vector<FloatDPBounds>(cast_singleton(product(D,to_time_bounds(t,t+h),A))))).upper();
     StepSizeType hlip = static_cast<StepSizeType>(cast_exact(this->_lipschitz_tolerance.value()/lipschitz));