Numerical error in conversion Star.getZono() causes 'lin-adaptive' algorithm for nonlinearODE reachability to fail

[mldiego]

Functions to reproduce error:

function vanderpol_error_nnv()
%% Reachability analysis of the nonlinear system Vanderpol  
reachstep = 0.05; % step size to compute reach sets
final_time = 30.0; % Time horizon
Initial_radius = 0.01; % Uncertainty in dynamics.
model = NonLinearODE(2,1,@vanderpol, reachstep, final_time,eye(2));
model.options.alg = 'lin-adaptive';

% Initial states
x0 = double([-1;-1]);
lb = x0 - Initial_radius;
ub = x0 + Initial_radius;
init_set = Star(lb,ub);
input_set = Star(0,0); % No inputs, but need to define it
I = init_set.getZono; % This conversion carries some numerical errors od 1e-18 which leads to errors in adaptive algorithm (why? not sure)
U = input_set.getZono;
R0 = zonotope(I.c, I.V);
U = zonotope(U.c, U.V);

model.set_R0(R0);
model.set_U(U);
model.set_timeStep(reachstep);
model.set_tFinal(final_time);

% Compute reach set using CORA
options = model.options;
params = model.params;
Z0 = zonotope([-1;-1],[0.01 0;0 0.01]);
params.R0 = Z0;
sys = nonlinearSys(model.dynamics_func, model.dim, model.nI); % CORA nonlinearSys class
Rc = reach(sys, params, options); % This works

mattt = double(single(init_set.V)); % This works (converting so single first eliminates the small numerical errors from onversion)
R0 = zonotope(mattt);
params.R0 = R0;
sys = nonlinearSys(model.dynamics_func, model.dim, model.nI); % CORA nonlinearSys class
Rc2 = reach(sys, params, options);

mattt = double(init_set.V);
R0 = zonotope(mattt);
params.R0 = R0;
sys = nonlinearSys(model.dynamics_func, model.dim, model.nI); % CORA nonlinearSys class
Rc3 = reach(sys, params, options); % This returns an error

% Compute reachability analysis
t = tic;
R_nnv = model.stepReachStar(init_set,input_set); % Same error as above
toc(t);
end

Vanderpol dynamics function to use:

function dx = vanderpol(x,u)
  dx(1,1) = x(2);
  dx(2,1) = (x(1)*x(1)-1)*x(2)-x(1);
end

Error log

Error using interval
Wrong input arguments for constructor of class interval!
   Information: Lower limit larger than upper limit.

Error in expmtie_adaptive (line 91)
        Asum = interval(Asum_neg,Asum_pos);

Error in linearSys/initReach_adaptive (line 37)
[sys,linOptions] = expmtie_adaptive(sys,linOptions);

Error in nonlinearSys/linReach_adaptive (line 94)
        [Rlin,options,linOptions] = initReach_adaptive(linSys,Rdelta,options,linOptions);

Error in nonlinearSys/reach_adaptive (line 69)
    [Rnext.ti,Rnext.tp,options] = linReach_adaptive(obj,options,options.R);

Error in contDynamics/reach (line 54)
            [timeInt,timePoint,res,tVec,options] = reach_adaptive(obj,params,options);

Error in vanderpol_error_nnv (line 43)
Rc3 = reach(sys, params, options); % This returns an error

Numerical error:

>> R0.Z - Z0.Z

ans =

   1.0e-17 *

                   0   0.867361737988404                   0
                   0                   0   0.867361737988404

We need to discuss how we want to deal with these numerical errors when converting from Star to Box to Zono to use these methods. For now, we could add some warnings if the adaptive algorithms are used.

This is not a one-of-a-kind error; I have been able to reproduce this error with other examples. A quick fix could be to overapproximate the values of ub and lb of the Box conversion with single precision and once the zonotope is generated convert the center and generators of the zonotope to double precision. This seems to fix this problem, at least in the examples where I have encountered this error.