This microsite is associated to the paper
Online interval depth localization of an underwater robot with ballast

This paper work presents an efficient online method to simulate a dynamical system with interval uncertainties. These uncertainties, can be either on the initial state vector, on the time-dependent inputs, or on the evolution function. Compared to other techniques used for the guaranteed integration of differential inclusion, the presented approach is online and requires a fixed number of operations at each sampling time. An illustration related to the depth reachability of an underwater robot with ballast is considered to illustrate the efficiency of the approach.

Interval Runge-Kutta method

from roblib import *
from codac import *

def uniplot(T,X,col): plot(T,X,col)

def simu_RK_interval(V0,t1,A,B):
    T,Vmin,Vmax=[],[],[]
    vmin=V0.lb()
    vmax=V0.ub()
    for t in arange(0,t1,dt):
        Vmin.append(vmin)
        Vmax.append(vmax)
        T.append(t)
        def fmin(v):  return (A-B*v*abs(v)).lb()
        def fmax(v):  return (A-B*v*abs(v)).ub()
        vmin=vmin+dt*fmin(vmin+(dt/2)*fmin(vmin))
        vmax=vmax+dt*fmax(vmax+(dt/2)*fmax(vmax))
    uniplot(T,Vmin,"red")
    uniplot(T,Vmax,"red")
    return


figure()
simu_RK_interval(Interval(v0-dv,v0+dv),t1,Interval(a-da,a+da),Interval(b-db,b+db))

Sinking body problem

def psiplus(v0,t,a,b):
    vbar=sqrt(a/b)
    def expc(x):
        if abs(x)<1E-9: return 1.0
        return ((exp(x)-1)/x)
    ee=exp(2*b*vbar*t)
    return (ee*(vbar+v0)+v0-vbar)/(1+ee+2*v0*b*t*expc(2*b*vbar*t))


def psipsiplus(v0,T,a,B):
    if type(B)==float : B=Interval(B)
    def psipsiplus_b(v0,T,a,B):
        if type(T)==float : T=Interval(T)
        Vbar=sqrt(a/B)
        def expc(x):
            X=Interval(x,x)
            if abs(x)<0.01:
                return 1+Interval(0.5)*X+0.00002*Interval(-1,1) # 0.00002 : upper bound for the Taylor remainder
            return (exp(X)-1)/X
        def Expc(X): return expc(X.lb())|expc(X.ub())
        EE=exp(2*sqrt(a*B)*T)
        return (EE*(Vbar+v0)+v0-Vbar)/(1+EE+2*v0*B*T*Expc(2*sqrt(a*B)*T))
    return psipsiplus_b(v0,T,a,Interval(B.ub()))|psipsiplus_b(v0,T,a,Interval(B.lb()))


def psiminus(v0,t,a,b):  return (-sqrt(-a/b))*tan(np.arctan(v0/(-sqrt(-a/b)))-b*(-sqrt(-a/b))*t)



def t2(v0,a,b): return -np.sign(a)*(np.arctan(v0*(np.sqrt(b/np.abs(a)))))/(np.sqrt(np.abs(a)*b))

def phi(v0,t,a,b):
    s=np.sign(a)
    if (a*v0>=0): return s*psiplus(s*v0,t,s*a,b)
    else:
        if t=0): return s*psipsiplus(s*v0,t,s*a,B)
    else:
        T2=Interval(t2(v0,a,B.lb()))|Interval(t2(v0,a,B.ub()))
        if not(t in T2):
            return Interval(phi(v0,t,a,B.lb())) | Interval(phi(v0,t,a,B.ub()))
        else: return s*psipsiplus(0.0,t-T2,s*a,B)


def phiphi_v0TaB(v0,T,a,B):
    return phiphi_v0taB(v0,T.lb(),a,B)|phiphi_v0taB(v0,T.ub(),a,B)


def phiphi(V0,T,A,B):
    if B.lb()<=0: print("Error in Psi : B should be positive")
    return phiphi_v0TaB(V0.lb(),T,A.lb(),B)|phiphi_v0TaB(V0.ub(),T,A.ub(),B)




def simu_phi(v0,t1,a,b):
    dv,da,db=0.1,0.1,0.1
    dt1=0.1
    V0=Interval(v0-dv,v0+dv)
    A=Interval(a-da,a+da)
    B=Interval(b-db,b+db)
    for t in arange(0,t1,dt1):
        T1=Interval(t,t+dt1)
        V=phiphi(V0,T1,A,B)
        draw_box_border(T1.lb(),T1.ub(),V.lb(),V.ub(),'black',1)
        T1=Interval(t+dt1)
        V=phiphi(V0,T1,A,B)
        draw_box_border(T1.lb(),T1.ub(),V.lb(),V.ub(),'blue',3)
    draw_box_border(0,0,V0.lb(),V0.ub(),'magenta',3)
    return


def drawpsiminus_b(v0,t,a,b1):
    P,B=[],[]
    figure()
    for b in arange(0.001,b1,0.01):
        P.append(psiminus(v0,t,a,b))
        B.append(b)
    uniplot(B,P,"red")


dt=0.001
simu_phi(1,5,1,2)
simu_phi(-1,5,-1,2)
simu_phi(1,5,-1,1)
simu_phi(-1,5,1,1)

Reachability of the float

from roblib import *
from codac import *
import numpy as np
from reachball import phiphi

dt = 0.1
tmax=20 #2
kmax=int(tmax/dt)
g0,l0,beta0,cx0=9.81,1,0.1,0.9
S0=Interval([0,0.1])
V0=Interval([0,0.1])
D0=Interval([0,0.1])
S,V,D=S0,V0,D0
figure()
for k in range(0,kmax):
  T=Interval([k*dt,(k+1)*dt])
  _U=exp(-T)
  draw_box_border(T.lb(),T.ub(),_U.lb(),_U.ub(),'black',1)
  draw_box_border(T.ub(),T.ub(),_U.lb(),_U.ub(),'blue',2)
  _S=S+_U*Interval(0,dt)
  S=S+_U*dt
  T1=Interval((k+1)*dt)
  draw_box_border(T.lb(),T.ub(),_S.lb(),_S.ub(),'black',1)
  draw_box_border(T1.lb(),T1.ub(),S.lb(),S.ub(),'blue',2)  
  _A=g0*(1-1/(1+beta0*_S))
  _B=((0.5*cx0)/((1+beta0*_S)*l0))
  _V=phiphi(V,Interval(0,dt),_A,_B)
  V=phiphi(V,Interval(dt),_A,_B)
  draw_box_border(T.lb(),T.ub(),_V.lb(),_V.ub(),'black',1)
  draw_box_border(T1.lb(),T1.ub(),V.lb(),V.ub(),'blue',2)

  _D=D+_V*Interval(0,dt)
  D=D+_V*Interval(dt)
  draw_box_border(T.lb(),T.ub(),_D.lb(),_D.ub(),'black',1)
  draw_box_border(T1.lb(),T1.ub(),D.lb(),D.ub(),'blue',2)