Merge branch 'master' into master

.gitignore

@@ -0,0 +1,2 @@
+*.txt~
+*.m~

benchmarks/ACC/reachACC.m

@@ -0,0 +1,59 @@
+%% Reachability analysis of TORA (benchmark 9)
+% Load components and set reachability parameters
+net = Load_nn('controller_5_20.onnx');
+reachStep = 0.01;
+controlPeriod = 0.1;
+output_mat = [0 0 0 0 1 0;1 0 0 -1 0 0; 0 1 0 0 -1 0]; % feedback: relative distance, relative velocity and ego-car velocity
+plant = NonLinearODE(6,1,@dynamicsACC, reachStep, controlPeriod, output_mat);
+% noise = Star(-0.0001, 0.0001);
+% plant.set_taylorTerms(5);
+% plant.set_zonotopeOrder(200);
+% plant.set_polytopeOrder(20);
+% error = 0.1;
+% plant.options.maxError = [error; error; error; error; error; error];
+
+%% Reachability analysis
+tF = 5; % seconds
+time = 0:controlPeriod:5;
+steps = length(time);
+input_ref = [30;1.4];
+% Initial set
+lb = [90; 32; 0; 10; 30; 0];
+ub = [110; 32.2; 0; 0; 30.2; 11];
+init_set = Star(lb,ub);
+% Input set
+lb = 0;
+ub = 0;
+input_set = Star(lb,ub);
+% Store all reachable sets
+% Execute reachabilty analysis
+nncs = NonlinearNNCS(net,plant);
+reachPRM.ref_input = input_ref;
+reachPRM.numSteps = 50;
+reachPRM.init_set = init_set;
+reachPRM.numCores = 1;
+reachPRM.reachMethod = 'approx-star';
+[R,rT] = nncs.reach(reachPRM);
+
+%% Visualize results
+t_gap = 1.4;
+D_default = 10;
+outAll = [];
+safe_dis = [];
+for i=1:length(plant.intermediate_reachSet)
+    outAll = [outAll plant.intermediate_reachSet(i).affineMap(output_mat,[])];
+    safe_dis = [safe_dis plant.intermediate_reachSet(i).affineMap([0 0 0 0 t_gap 0], D_default)];
+end
+times = reachStep:reachStep:tF;
+% save('../../results/ACC_distance','R','rT','outAll','safe_dis','reachStep','error','-v7.3');
+f = figure;
+Star.plotRanges_2D(outAll,2,times,'r');
+hold on;
+Star.plotRanges_2D(safe_dis,1,times,'b');
+title('Relative distance (red) vs. Safe distance (blue)');
+xlabel('Time (s)');
+ylabel('Distance (m)')
+saveas(f,'../../results/ACC_distance_sets.jpg');
+
+
+

benchmarks/Airplane/reach_airplane.m

@@ -0,0 +1,119 @@
+%% Reachability analysis of Airplane Benchmark
+
+%% load the controller
+net = Load_nn('controller_airplane.mat');
+
+% Specify the reach step, has to be smaller than the control period
+reachStep = 0.001;
+%% specify the control period as specified by the benchmark description
+controlPeriod = 0.1;
+
+% define the plant as specified by nnv
+plant = NonLinearODE(12,6,@dynamics, reachStep, controlPeriod, eye(12));
+plant.set_taylorTerms(10)
+plant.set_zonotopeOrder(10);
+% plant.set_polytopeOrder(20);
+%error = 0.01;
+%plant.options.maxError = [error; error;error;error];
+
+%% Reachability analysis
+% Initial set
+lb = [0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0];
+% ub = [0.0;0.0;0.0;1.0;1.0;1.0;1.0;1.0;1.0;0.0;0.0;0.0];
+ub = [0.0;0.0;0.0;0.01;0.01;0.01;0.01;0.01;0.01;0.0;0.0;0.0];
+% ub = [0.0;0.0;0.0;0.001;0.001;0.001;0.001;0.001;0.001;0.0;0.0;0.0];
+init_set = Star(lb,ub);
+% Input set
+lb = [0;0;0;0;0;0];
+ub = [0;0;0;0;0;0];
+input_set = Star(lb,ub);
+% Store all reachable sets
+reachAll = init_set;
+% Execute reachabilty analysis
+num_steps = 13;
+t = tic;
+for i=1:num_steps
+    % Compute plant reachable set
+    init_set = plant.stepReachStar(init_set(1),input_set(1));
+%     init_set = plantReach(plant,init_set,input_set);
+    reachAll = [reachAll init_set];
+    % Compute controller output set
+    input_set = net.reach(init_set,'approx-star');
+end
+timing = toc(t);
+save('../../results/airplaneReach.mat','-v7.3');
+%% Visualize results
+t = tic;
+f1 = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,1,4,'b');
+grid;hold on;
+Star.plotBoxes_2D_noFill(reachAll,1,4,'m');
+title('Airplane x_1 vs. x_4');
+xlabel('x');
+ylabel('u');
+
+f2 = figure;
+Star.plotBoxes_2D_noFill(reachAll,2,5,'b');
+grid;hold on;
+Star.plotBoxes_2D_noFill(reachAll,2,5,'m');
+plot([-0.5 -0.5],[-2 2],'r');
+plot([0.5 0.5],[-2 2],'r');
+title('Airplane x_2 vs x_5');
+xlabel('y');
+ylabel('v');
+saveas(f2,'../../results/reachAirplane_plot2vs5_cP.jpg');
+
+f3 = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,3,6,'b');
+grid;hold on;
+Star.plotBoxes_2D_noFill(reachAll,3,6,'m');
+title('Airplane');
+xlabel('z');
+ylabel('w');
+
+f4 = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,7,10,'b');
+grid;hold on;
+Star.plotBoxes_2D_noFill(reachAll,7,10,'m');
+plot([-1 -1],[-0.2 0.2],'r');
+plot([1 1],[-0.2 0.2],'r');
+title('Airplane x_7 vs. x_{10}');
+xlabel('x_7');
+ylabel('x_{10}');
+saveas(f4,'../../results/reachAirplane_plot7vs10.jpg');
+
+f5 = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,8,11,'b');
+grid;hold on;
+Star.plotBoxes_2D_noFill(reachAll,8,11,'m');
+plot([-1 -1],[-0.2 0.2],'r');
+plot([1 1],[-0.2 0.2],'r');
+title('Airplane x_8 vs. x_{11}');
+xlabel('x_8');
+ylabel('x_{11}');
+saveas(f5,'../../results/reachAirplane_plotvs11.jpg');
+
+f6 = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,9,12,'b');
+grid;hold on;
+Star.plotBoxes_2D_noFill(reachAll,9,12,'m');
+plot([-1 -1],[-0.2 0.2],'r');
+plot([1 1],[-0.2 0.2],'r');
+title('Airplane x_9 vs. x_{12}');
+xlabel('x_9');
+ylabel('x_{12}');
+saveas(f6,'../../results/reachAirplane_plot9vs12.jpg');
+toc(t);
+
+%% Helper function
+function init_set = plantReach(plant,init_set,input_set)
+    nS = length(init_set);
+    nL = length(input_set);
+    ss = [];
+    for k=1:nS
+        for l=1:nL
+            ss =[ss plant.stepReachStar(init_set(k), input_set(l))];
+        end
+    end
+    init_set = ss;
+end

benchmarks/Benchmark10-Unicycle/dynamics10.m

@@ -10,7 +10,7 @@
   %       u(2): steering angle of front wheel
   % disturbance input
   %       w: disturbance with a range (-10^-4,10^4)
-  % Initial state range [9.5, 9.55] × [-4.5, -4.45] × [2.1, 2.11] × [1.5, 1.51]
+  % Initial state range [9.5, 9.55] × [-4.5, -4.45] × [2.1, 2.11] × [1.5, 1.51]
 %
 % The output of the neural network f(x) needs to be normalized 
 % in order to obtain u(1) and u(2).

benchmarks/Benchmark10-Unicycle/reach10.m

@@ -0,0 +1,85 @@
+%% Reachability analysis of the Unicycle (benchmark 10)
+% net = Load_nn('controllerB.onnx');
+net = Load_nn('controllerB_nnv.mat');
+controlPeriod = 0.2;
+% controlPeriod = 0.5;
+reachstep = 0.001;
+plant = NonLinearODE(4,2,@dynamics10, reachstep, controlPeriod, eye(4));
+% noise = Star(-0.0001, 0.0001);
+% plant.set_taylorTerms(2);
+% plant.set_zonotopeOrder(50);
+% plant.set_polytopeOrder(20);
+% error = 0.01;
+% plant.options.maxError = [error; error; error; error];
+tF = 10;
+time = 0:controlPeriod:tF;
+steps = length(time);
+offset = 30;
+offsetM = offset*ones(2,1);
+scale_factor = 1;
+
+%% Reachability analysis
+% Initial set
+lb = [9.5; -4.5; 2.1; 1.5];
+ub = [9.51; -4.49; 2.11; 1.51];
+% ub = [9.55; -4.45; 2.11; 1.51];
+init_set = Star(lb,ub);
+% Input set
+lb = [0;0];
+ub = [0;0];
+input_set = Star(lb,ub);
+% Store all reachable sets
+reachAll = init_set;
+% Execute reachabilty analysis
+steps = 10;
+nI = 1;
+t = tic;
+for i=1:steps
+    % Compute plant reachable set
+%     init_set = plant.stepReachStar(init_set, input_set);
+    init_set = plantReach(plant,init_set(nI),input_set(1));
+    nI = length(init_set);
+    reachAll = [reachAll init_set(nI)];
+    % Compute controller output set
+    inp_set = net.reach(init_set(nI),'approx-star');
+    input_set = [];
+    input_set = inp_set(1).affineMap(eye(2),-offsetM);
+%     for k = 1:length(inp_set)
+%         input_set = [input_set inp_set(k).affineMap(eye(2),-offsetM)];
+%     end
+end
+t = toc(t);
+% save('../../results/reach10','plant','t','reachAll','-v7.3')
+%% Visualize results
+f = figure;
+% Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,1,2,'b');
+hold on;
+Star.plotBoxes_2D_noFill(reachAll,1,2,'b');
+grid;
+title('Benchmark 10 - Unicycle');
+xlabel('x1');
+ylabel('x2');
+% saveas(f,'../../results/reach10_plot1v2.jpg');
+
+f1 = figure;
+% Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,3,4,'b');
+hold on;
+Star.plotBoxes_2D_noFill(reachAll,3,4,'b');
+grid;
+title('Benchmark 10 - Unicycle');
+xlabel('x1');
+ylabel('x2');
+% saveas(f1,'../../results/reach10_plot3v4.jpg');
+
+%% Helper function
+function init_set = plantReach(plant,init_set,input_set)
+    nS = length(init_set);
+    nL = length(input_set);
+    ss = [];
+    for k=1:nS
+        for l=1:nL
+            ss =[ss plant.stepReachStar(init_set(k), input_set(l))];
+        end
+    end
+    init_set = ss;
+end

benchmarks/Benchmark9-Tora/reach9.m

@@ -0,0 +1,86 @@
+%% Reachability analysis of TORA (benchmark 9)
+% Load components and set reachability parameters
+net = Load_nn('controllerTora.onnx');
+reachStep = 0.005;
+% reachStep = 0.01;
+controlPeriod = 1;
+plant = NonLinearODE(4,1,@dynamics9, reachStep, controlPeriod, eye(4));
+% noise = Star(-0.0001, 0.0001);
+plant.set_taylorTerms(10);
+plant.set_zonotopeOrder(100);
+plant.set_polytopeOrder(5);% error = 0.001;
+error = 0.01;
+plant.options.maxError = [error; error; error; error];
+time = 0:controlPeriod:20;
+steps = length(time);
+% Initial set
+lb = [0.69; -0.7; -0.4; 0.5];
+% ub = [0.7; -0.6; -0.3; 0.6];
+% ub = [0.61; -0.69; -0.39; 0.51];
+ub = [0.7; -0.69; -0.39; 0.51];
+% init_set = Box(lb,ub);
+offset = 10;
+scale_factor = 1;
+
+%% Reachability analysis
+init_set = Star(lb,ub);
+% Input set
+lb = 0;
+ub = 0;
+input_set = Star(lb,ub);
+% Store all reachable sets
+reachAll = init_set;
+% Execute reachabilty analysis
+% for i =1:steps
+t = tic;
+for i =1:6
+    % Compute plant reachable set
+    init_set = plantReach(plant,init_set,input_set);
+%     init_set = plant.stepReachStar(init_set, input_set);
+    reachAll = [reachAll init_set];
+    % Compute controller output set
+    inp_set = net.reach(init_set,'approx-star');
+    input_set = [];
+    for k = 1:length(inp_set)
+        input_set = [input_set inp_set(k).affineMap(1,-offset)];
+    end
+end
+disp(' ');
+timing = toc(t);
+save('../../results/reach9','plant','reachAll','timing','-v7.3')
+
+%% Visualize results
+f = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,1,2,'b');
+grid;
+hold on;
+Star.plotBoxes_2D_noFill(reachAll,1,2,'m');
+grid;
+title('Reachable sets Benchmark 9 - TORA')
+xlabel('x1');
+ylabel('x3');
+saveas(f,'../../results/reach9_plot1v3.jpg');
+
+f1 = figure;
+Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,2,4,'b');
+grid;
+hold on;
+Star.plotBoxes_2D_noFill(reachAll,2,4,'m');
+grid;
+title('Reachable sets Benchmark 9 - TORA')
+xlabel('x2');
+ylabel('x4');
+saveas(f1,'../../results/reach9_plot3v4.jpg');
+
+%% Helper function
+function init_set = plantReach(plant,init_set,input_set)
+    nS = length(init_set);
+    nL = length(input_set);
+    ss = [];
+    for k=1:nS
+        for l=1:nL
+            ss =[ss plant.stepReachStar(init_set(k), input_set(l))];
+        end
+    end
+    init_set = ss;
+end

benchmarks/Double_Pendulum/dynamics.m → benchmarks/Double_Pendulum/dynamics_dp.m

@@ -1,4 +1,4 @@
-function [dx] = dynamics(t,x,T)
+function [dx] = dynamics_dp(t,x,T)
 % Ex_double pendulum 
 % constants as per the documentation m=0.5, L=0.5, c= 0, g=1.0
 % description of the pararameters in these equations can be found here