Update reach files

benchmarks/Airplane/reach_airplane.m

@@ -1,24 +1,27 @@
-%% Reachability analysis of Double Pendulum Benchmark
+%% 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.1;
+reachStep = 0.001;
 %% specify the control period as specified by the benchmark description
-controlPeriod = 0.05;
+controlPeriod = 0.1;
 
 % define the plant as specified by nnv
 plant = NonLinearODE(12,6,@dynamics, reachStep, controlPeriod, eye(12));
-%plant.set_zonotopeOrder(50);
+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;0.0;0.0; 0.0;1.0;1.0;1.0;1.0;1.0;1.0];
+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];
@@ -27,34 +30,90 @@
 % Store all reachable sets
 reachAll = init_set;
 % Execute reachabilty analysis
-% for i =1:steps
-num_steps =10;
+num_steps = 13;
+t = tic;
 for i=1:num_steps
     % Compute plant reachable set
-    init_set = plant.stepReachStar(init_set, input_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');
 
-%times = reachStep:reachStep:(num_steps*controlPeriod);
-%Star.plotRanges_2D(plant.intermediate_reachSet,1,times,'r')
-
-% times = 0:controlPeriod:(num_steps*controlPeriod);
-% Star.plotRanges_2D(reachAll,1,times,'r')
-% 
-% 
-% f = figure;
-% Star.plotBoxes_2D_noFill(reachAll,1,2,'b');
-% grid;
-% title('Single Pendulum reachable sets');
-% xlabel('x1');
-% ylabel('x2');
-% 
-% f1 = figure;
-% Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,1,2,'b');
-% grid;
-% title('Single Pendulum reachable sets');
-% xlabel('x1');
-% ylabel('x2');
+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/Airplane/specifications.txt

@@ -3,20 +3,21 @@ Initial states:
 x1 = 0
 x2 = 0
 x3 = 0
-x4 = 0
-x5 = 0
-x6 = 0
+x4 = [0.0, 1.0]
+x5 = [0.0, 1.0]
+x6 = [0.0, 1.0]
 x7 = [0.0, 1.0]
 x8 = [0.0, 1.0]
 x9 = [0.0, 1.0]
-x10 = [0.0, 1.0]
-x11 = [0.0, 1.0]
-x12 = [0.0, 1.0]
+x10 = 0
+x11 = 0
+x12 = 0
 
-t = 20 seconds
+t = 20 control steps = 2 seconds
+control step = 0.1
 
 Property:
 
-x2 should be in [−0.5, 0.5] and x10,x11,x12 should be in [-1.0,1.0]
+x2 should be in [−0.5, 0.5] and x7,x8,x9 should be in [-1.0,1.0]
 
 Refer to this for more details: https://github.com/amaleki2/benchmark_closedloop_verification/blob/master/AINNC_benchmark.pdf

benchmarks/Benchmark10-Unicycle/dynamics10.m

@@ -1,4 +1,4 @@
-function [dx] = dynamics10(t,x,u,w)
+function [dx] = dynamics10(t,x,u)
 % Ex_car_model
   %vehicleODE Bicycle model of a vehicle with
   % states
@@ -15,6 +15,6 @@
 dx(1,1) = x(4) * cos(x(3));
 dx(2,1) = x(4) * sin(x(3));
 dx(3,1) = u(2);
-dx(4,1) = u(1) + w;
-
+dx(4,1) = u(1); % No noise
+% dx(4,1) = u(1) + w; % noise
 end

benchmarks/Benchmark10-Unicycle/reach10.m

@@ -2,32 +2,27 @@
 % net = Load_nn('controllerB.onnx');
 net = Load_nn('controllerB_nnv.mat');
 controlPeriod = 0.2;
-reachstep = 0.01;
+% 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_taylorTerms(2);
+% plant.set_zonotopeOrder(50);
 % plant.set_polytopeOrder(20);
-error = 0.01;
-plant.options.maxError = [error; error; error; error];
+% error = 0.01;
+% plant.options.maxError = [error; error; error; error];
 tF = 10;
 time = 0:controlPeriod:tF;
 steps = length(time);
-offset = 20;
+offset = 30;
 offsetM = offset*ones(2,1);
 scale_factor = 1;
 
-%% Simulate system
-% There in an offset in the outputs, cannot simulate it like this
-% nncs = NonlinearNNCS(net,plant);
-% [simTrace, controlTrace] = nncs.evaluate(controlPeriod,40,[9.5;-4.5;2.2;1.5],[]);
-% figure;
-% plot(simTrace(1,:),simTrace(2,:));
-
 %% Reachability analysis
 % Initial set
 lb = [9.5; -4.5; 2.1; 1.5];
-ub = [9.55; -4.45; 2.11; 1.51];
+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];
@@ -36,25 +31,55 @@
 % Store all reachable sets
 reachAll = init_set;
 % Execute reachabilty analysis
-% for i =1:steps
+steps = 10;
+nI = 1;
 t = tic;
-for i=1:30
+for i=1:steps
     % Compute plant reachable set
-    init_set = plant.stepReachStar(init_set, input_set);
-    reachAll = [reachAll init_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
-    input_set = net.reach(init_set,'approx-star');
-    input_set = input_set.affineMap(eye(2),-offsetM);
+    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','-v7.3')
-
+% save('../../results/reach10','plant','t','reachAll','-v7.3')
+%% Visualize results
 f = figure;
-Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,1,2,'b');
+% Star.plotBoxes_2D_noFill(plant.intermediate_reachSet,1,2,'b');
 hold on;
 Star.plotBoxes_2D_noFill(reachAll,1,2,'b');
 grid;
-title('Benchmark 10 reachable sets');
+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(f,'../../results/reach10_plot.jpg');
+% 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