diff --git a/code/nnv/engine/utils/export2vnnlib.m b/code/nnv/engine/utils/export2vnnlib.m index 85d7213627..03e1b29d83 100644 --- a/code/nnv/engine/utils/export2vnnlib.m +++ b/code/nnv/engine/utils/export2vnnlib.m @@ -50,17 +50,30 @@ function export2vnnlib(lb, ub, outsize, property, name) fprintf(fID,"\n"); % 3) Define outputs -nO = length(property.g); % number of output constraints -if nO > 1 +nHS = length(property); % number of halfspaces +if nHS > 1 fprintf(fID,"(assert (or \n"); - for i=1:nO - constraint = halfspaceConstraint2inequality(property.G(i,:), property.g(i)); - fprintf(fID," "+constraint+"\n"); % format used in other vnnlib examples + for j=1:nHS + nO = length(property(j).g); % number of output constraints + if nO > 1 + fprintf(fID," (and "); + for i=1:nO + constraint = halfspaceConstraint2inequality_1(property(j).G(i,:), property(j).g(i)); + fprintf(fID," "+constraint(1:end-1)); % add constrint removing last parenthesis + end + fprintf(fID,")\n"); + else + constraint = halfspaceConstraint2inequality_1(property(j).G, property(j).g); + fprintf(fID," (and "+constraint+"\n"); + end end fprintf(fID,"))"); else - constraint = halfspaceConstraint2inequality_1(property.G, property.g); - fprintf(fID,"(assert "+constraint); + nO = length(property.g); % number of output constraints + for i=1:nO + constraint = halfspaceConstraint2inequality_1(property.G(i,:), property.g(i)); + fprintf(fID,"(assert "+constraint); + end end % close and save file @@ -71,22 +84,22 @@ function export2vnnlib(lb, ub, outsize, property, name) %% Helper functions -function str = halfspaceConstraint2inequality(hRow, hVal) - % Input a halfspace row (G row) and the corresponding g value - % Outputs a string to write in the vnnlib file - - locs = find(hRow ~= 0); % Find indexes that are not zero - if hVal == 0 % Compare two indexes - if hRow(locs(1)) > 0 % - str = "(and (>= Y_"+string(locs(2)-1) + " " + "Y_"+string(locs(1)-1)+ "))"; - else - str = "(and (>= Y_"+string(locs(1)-1) + " " + "Y_"+string(locs(2)-1) + "))"; - end - else % compare index to value - str = "(and (>= Y_"+string(locs(1)-1) + " " + num2str(hVal, '%.16f') + "))"; - end - -end +% function str = halfspaceConstraint2inequality(hRow, hVal) +% % Input a halfspace row (G row) and the corresponding g value +% % Outputs a string to write in the vnnlib file +% +% locs = find(hRow ~= 0); % Find indexes that are not zero +% if hVal == 0 % Compare two indexes +% if hRow(locs(1)) > 0 % +% str = "(and (>= Y_"+string(locs(2)-1) + " " + "Y_"+string(locs(1)-1)+ "))"; +% else +% str = "(and (>= Y_"+string(locs(1)-1) + " " + "Y_"+string(locs(2)-1) + "))"; +% end +% else % compare index to value +% str = "(and (>= Y_"+string(locs(1)-1) + " " + num2str(hVal, '%.16f') + "))"; +% end +% +% end function str = halfspaceConstraint2inequality_1(hRow, hVal) % Input a halfspace row (G row) and the corresponding g value diff --git a/code/nnv/examples/NN/medmnist/generate_vnnlib_files.m b/code/nnv/examples/NN/medmnist/generate_vnnlib_files.m new file mode 100644 index 0000000000..4d9ffc3c4b --- /dev/null +++ b/code/nnv/examples/NN/medmnist/generate_vnnlib_files.m @@ -0,0 +1,74 @@ +%% Let's create some examples for medmnist 2D and 3D + +%% All 2D datasets + +datapath = "data/mat_files/"; +datafiles = ["bloodmnist"; "breastmnist.mat"; "dermamnist.mat"; "octmnist.mat"; "organamnist.mat"; ... + "organcmnist.mat"; "organsmnist.mat"; "pathmnist.mat"; "pneumoniamnist"; "retinamnist.mat"; "tissuemnist.mat"]; + +N = 10; % number of vnnlib files to create +epsilon = [1,2,3]; % {epsilon} pixel color values for every channel + +for k=1:length(datafiles) + % load data + load(datapath + datafiles(k)); + % preprocess dataa + test_images = permute(test_images, [2 3 4 1]); + test_labels = test_labels + 1; + outputSize = length(unique(test_labels)); % number of classes in dataset + % create file name + dataname = split(datafiles(k), '.'); + name = "vnnlib/" + dataname{1} + "_linf_"; + % create vnnlib files + for i=1:N + img = test_images(:,:,:,i); + outputSpec = create_output_spec(outputSize, test_labels(i)); + for j=1:length(epsilon) + [lb,ub] = l_inf_attack(img, epsilon(j), 255, 0); + vnnlibfile = name+string(epsilon(j))+"_"+string(i)+".vnnlib"; + export2vnnlib(lb, ub, outputSize, outputSpec, vnnlibfile); + disp("Created property "+vnnlibfile); + end + end +end + + + + +%% Helper functions + +% Return the bounds an linf attack +function [lb,ub] = l_inf_attack(img, epsilon, max_value, min_value) + imgSize = size(img); + disturbance = epsilon * ones(imgSize, "like", img); % disturbance value + lb = max(img - disturbance, min_value); + ub = min(img + disturbance, max_value); + lb = single(lb); + lb = reshape(lb, [], 1); + ub = single(ub); + ub = reshape(ub, [], 1); +end + +% Define unsafe (not robust) property +function Hs = create_output_spec(outSize, target) + % @Hs: unsafe/not robust region defined as a HalfSpace + % - target: label idx of the given input set + + if target > outSize + error("Target idx must be less than or equal to the output size of the NN."); + end + + % Define HalfSpace Matrix and vector + G = ones(outSize,1); + G = diag(G); + G(target, :) = []; + G = -G; + G(:, target) = 1; + + % Create HalfSapce to define robustness specification + Hs = []; + for i=1:height(G) + Hs = [Hs; HalfSpace(G(i,:), 0)]; + end + +end diff --git a/code/nnv/examples/NN/medmnist/medmnistmodels2onnx.m b/code/nnv/examples/NN/medmnist/medmnistmodels2onnx.m new file mode 100644 index 0000000000..90df7ad2de --- /dev/null +++ b/code/nnv/examples/NN/medmnist/medmnistmodels2onnx.m @@ -0,0 +1,14 @@ +% Transform all models into onnx + +models = dir("models/*.mat"); + +if ~isfolder('onnx') + mkdir('onnx') +end + +for i=1:length(models) + filename = string(models(i).name); + load("models/" + filename); + onnxfile = split(filename, '.'); + exportONNXNetwork(net,['onnx/', onnxfile{1}, '.onnx']); +end diff --git a/code/nnv/examples/NN/xai/single_pixel_range.png b/code/nnv/examples/NN/xai/single_pixel_range.png new file mode 100644 index 0000000000..55f7cdab3e Binary files /dev/null and b/code/nnv/examples/NN/xai/single_pixel_range.png differ diff --git a/code/nnv/examples/NN/xai/verify_baseline.m b/code/nnv/examples/NN/xai/verify_baseline.m new file mode 100644 index 0000000000..65e601988e --- /dev/null +++ b/code/nnv/examples/NN/xai/verify_baseline.m @@ -0,0 +1,83 @@ +%% Verify the importance of pixels using mnist + +%% Part 1. Compute reachability + +% Load the model +net = importNetworkFromONNX("super_resolution.onnx", "InputDataFormats","BCSS", "OutputDataFormats","BC"); +net = matlab2nnv(net); + +% Load data (no download necessary) +digitDatasetPath = fullfile(matlabroot,'toolbox','nnet','nndemos', ... + 'nndatasets','DigitDataset'); +% Images +imds = imageDatastore(digitDatasetPath, ... + 'IncludeSubfolders',true,'LabelSource','foldernames'); + +% Load one image in dataset +[img, fileInfo] = readimage(imds,7010); +target = single(fileInfo.Label); % label = 0 (index 1 for our network) +img = single(img)/255; % change precision + + +%% Verify 1st method +% assume something like brightnening or darkening attack for modifying pixel value, e.g. furthest from current value + +pixel_val_change = 5/255; % 5 pixel color values for range of pixel + +% First, we need to define the reachability options +reachOptions = struct; % initialize +reachOptions.reachMethod = 'approx-star'; % using exact/approx method + +% Reachability analysis +R(28,28) = ImageStar; +for i = 1:size(img,1) + for j = 1:size(img,2) + % Create set with brightening/darkening pixel + lb = img; ub = img; % initialize bounds + if img(i,j) > 122 % darkening + lb(i,j) = 0; + ub(i,j) = pixel_val_change; + else % brightening + lb(i,j) = 1-pixel_val_change; + ub(i,j) = 1; + end + % Create ImageStar + IS = ImageStar(lb,ub); + % Compute reachable set + t = tic; + R(i,j) = net.reach(IS, reachOptions); + toc(t); % track computation time + end +end + +% Visualize results +img_scores = net.evaluate(img); + +ub = zeros(size(img)); +lb = zeros(size(img)); + +for i = 1:size(img,1) + for j = 1:size(img,2) + [l,u] = R(i,j).getRange(1,1,target); + lb(i,j) = l - img_scores(target); + ub(i,j) = u - img_scores(target); + end +end + +diff = ub - lb; + +mapC = hot; % hot map to show the attributions + +% Visualize results +figure; +subplot(2,2,1); +imshow(img); +subplot(2,2,2); +imshow(lb, 'Colormap',winter, 'DisplayRange',[min(lb, [], 'all'), max(lb, [], 'all')]); +colorbar; +subplot(2,2,4); +imshow(ub, 'Colormap',winter, 'DisplayRange',[min(ub, [], 'all'), max(ub, [], 'all')]); +colorbar; +subplot(2,2,3); +imshow(diff, 'Colormap',winter, 'DisplayRange',[min(diff, [], 'all'), max(diff, [], 'all')]); +colorbar; diff --git a/code/nnv/examples/NN/xai/verify_super_resolution.m b/code/nnv/examples/NN/xai/verify_noise_tunnel.m similarity index 57% rename from code/nnv/examples/NN/xai/verify_super_resolution.m rename to code/nnv/examples/NN/xai/verify_noise_tunnel.m index 4716c9c240..cc19809f60 100644 --- a/code/nnv/examples/NN/xai/verify_super_resolution.m +++ b/code/nnv/examples/NN/xai/verify_noise_tunnel.m @@ -18,74 +18,6 @@ target = single(fileInfo.Label); % label = 0 (index 1 for our network) img = single(img)/255; % change precision -% First, we need to define the reachability options -reachOptions = struct; % initialize -reachOptions.reachMethod = 'approx-star'; % using exact/approx method -% does approx method makes sense here? -% use approx for now as exact may compute multiple sets -% at least for now, to show somehing - -max_noise = 0.02; % range -> [-0.02, 0.02] (approx 5 pixel color values) - -% Reachability analysis -R(28,28) = ImageStar; -for i = 1:size(img,1) - for j = 1:size(img,2) - disturbance = zeros(size(img), 'single'); % make a copy - disturbance(i,j) = max_noise; - IS = ImageStar(img, -disturbance, disturbance); - t = tic; - R(i,j) = net.reach(IS, reachOptions); - toc(t); - end -end - -%% Part 2. Verify XAI - -% takes less than a minute to compute this with exact reachability - -img_scores = net.evaluate(img); - -% but what are we verifying exactly? we could compute interval differences -% with original scores (img_scores) to see how things change, but what does -% that mean? - -% Let's do the attribution method (baseline) computation, -% but using the input generation of the noise tunnel - -% We cannot use the exact sets we computed, but we can get the ranges for -% all the outputs and use those for the attribution calculation for each -% pixel - -ub = zeros(size(img)); -lb = zeros(size(img)); - -for i = 1:size(img,1) - for j = 1:size(img,2) - [l,u] = R(i,j).getRange(1,1,target); - lb(i,j) = l - img_scores(target); - ub(i,j) = u - img_scores(target); - end -end - -diff = ub - lb; - -mapC = hot; % hot map to show the attributions - -% Visualize results -figure; -subplot(2,2,1); -imshow(img); -subplot(2,2,2); -imshow(lb, 'Colormap',winter, 'DisplayRange',[min(lb, [], 'all'), max(lb, [], 'all')]); -colorbar; -subplot(2,2,4); -imshow(ub, 'Colormap',winter, 'DisplayRange',[min(ub, [], 'all'), max(ub, [], 'all')]); -colorbar; -subplot(2,2,3); -imshow(diff, 'Colormap',winter, 'DisplayRange',[min(diff, [], 'all'), max(diff, [], 'all')]); -colorbar; - %% Verify 1st method % assume something like brightnening or darkening attack for modifying pixel value, e.g. furthest from current value diff --git a/code/nnv/examples/NN/xai/verify_range.m b/code/nnv/examples/NN/xai/verify_range.m new file mode 100644 index 0000000000..fbd855d49f --- /dev/null +++ b/code/nnv/examples/NN/xai/verify_range.m @@ -0,0 +1,76 @@ +%% Verify the importance of pixels using mnist + +%% Part 1. Compute reachability + +% Load the model +net = importNetworkFromONNX("super_resolution.onnx", "InputDataFormats","BCSS", "OutputDataFormats","BC"); +net = matlab2nnv(net); + +% Load data (no download necessary) +digitDatasetPath = fullfile(matlabroot,'toolbox','nnet','nndemos', ... + 'nndatasets','DigitDataset'); +% Images +imds = imageDatastore(digitDatasetPath, ... + 'IncludeSubfolders',true,'LabelSource','foldernames'); + +% Load one image in dataset +[img, fileInfo] = readimage(imds,7010); +target = single(fileInfo.Label); % label = 0 (index 1 for our network) +img = single(img)/255; % change precision + +% First, we need to define the reachability options +reachOptions = struct; % initialize +% reachOptions.reachMethod = 'exact-star'; % using exact/approx method +reachOptions.reachMethod = 'approx-star'; + +% Reachability analysis +% R(28,28) = ImageStar; +for i = 1:size(img,1) + for j = 1:size(img,2) + lb = img; + ub = img; + lb(i,j) = 0; + ub(i,j) = 1; + IS = ImageStar(lb,ub); + t = tic; + % R{i,j} = net.reach(IS, reachOptions); + R(i,j) = net.reach(IS, reachOptions); + toc(t); + end +end + +%% Part 2. Verify XAI + +% takes less than a minute to compute this with exact reachability + +img_scores = net.evaluate(img); + +ub = zeros(size(img)); +lb = zeros(size(img)); + +for i = 1:size(img,1) + for j = 1:size(img,2) + [l,u] = R(i,j).getRange(1,1,target); + lb(i,j) = l - img_scores(target); + ub(i,j) = u - img_scores(target); + end +end + +diff = ub - lb; + +mapC = hot; % hot map to show the attributions + +% Visualize results +figure; +subplot(2,2,1); +imshow(img); +subplot(2,2,2); +imshow(lb, 'Colormap',winter, 'DisplayRange',[min(lb, [], 'all'), max(lb, [], 'all')]); +colorbar; +subplot(2,2,4); +imshow(ub, 'Colormap',winter, 'DisplayRange',[min(ub, [], 'all'), max(ub, [], 'all')]); +colorbar; +subplot(2,2,3); +imshow(diff, 'Colormap',winter, 'DisplayRange',[min(diff, [], 'all'), max(diff, [], 'all')]); +colorbar; +saveas(gcf, "single_pixel_range.png"); diff --git a/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.h5 b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.h5 new file mode 100644 index 0000000000..6aa242509b Binary files /dev/null and b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.h5 differ diff --git a/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.mat b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.mat new file mode 100644 index 0000000000..1984b024c2 Binary files /dev/null and b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.mat differ diff --git a/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.nnet b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.nnet new file mode 100644 index 0000000000..5f9e900a87 --- /dev/null +++ b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.nnet @@ -0,0 +1,112 @@ +// The contents of this file are licensed under the Creative Commons +// Attribution 4.0 International License: https://creativecommons.org/licenses/by/4.0/ +// Neural Network File Format by Kyle Julian, Stanford 2016 +3, 2, 1, 25 +2, 25, 25,1, 1 +0 +0 +0 +0 +0 +-0.24301, -0.18771 +-0.26350, -0.30369 +0.14823, -0.35473 +-0.32726, 0.23397 +-0.03466, -0.19298 +0.37591, 0.46209 +0.72291, 0.69352 +0.03561, 0.26642 +-0.12026, 0.21470 +-0.16872, -0.30651 +-0.43472, -0.18382 +-0.13991, -0.16939 +-0.40639, 0.06983 +0.27941, -0.20360 +-0.40475, -0.43684 +-0.34817, -0.58015 +0.07824, -0.34545 +-0.36376, -0.39287 +0.65688, 0.18937 +-0.18365, 0.15440 +-0.31044, -0.13218 +-0.07592, -0.43100 +0.42164, -0.02560 +-0.51296, -0.67853 +-0.08458, 0.15043 +-0.00419 +-0.00546 +-0.00949 +-0.00565 +-0.00241 +-0.00529 +0.00291 +-0.00454 +-0.00197 +-0.00673 +-0.00373 +-0.00559 +-0.00445 +-0.00842 +-0.00319 +0.00551 +-0.00701 +-0.00560 +0.00371 +-0.00419 +0.00527 +-0.00635 +-0.00937 +0.00512 +-0.00610 +0.06736, -0.05422, 0.14574, 0.21675, 0.22926, -0.11144, 0.56496, -0.05118, -0.02968, 0.10144, 0.23500, -0.33200, 0.13735, 0.16915, 0.32100, 0.15801, -0.18926, -0.31402, 0.03553, -0.04112, -0.21541, 0.24437, -0.26637, 0.13420, -0.00569 +0.29358, 0.33571, -0.29688, -0.02312, -0.27224, 0.12862, 0.12963, -0.15199, -0.33584, -0.17085, 0.21235, -0.33336, 0.17564, 0.05524, -0.15493, -0.25028, 0.17259, -0.20334, 0.04693, -0.12292, 0.10917, -0.00523, -0.27183, -0.07113, -0.19691 +0.02020, 0.21657, -0.23130, 0.08048, 0.16010, -0.10221, 0.17049, 0.06040, -0.20372, -0.04221, -0.06663, 0.02040, -0.08222, 0.26648, 0.19941, -0.18592, -0.23603, 0.11598, -0.19459, 0.07057, 0.04573, -0.29013, -0.00217, 0.32552, 0.31872 +0.08224, -0.30303, -0.06345, 0.13502, -0.01831, 0.28010, -0.35218, -0.06258, -0.05898, -0.32436, 0.32003, -0.30794, 0.07805, -0.02749, -0.32308, 0.14680, -0.26048, -0.11922, -0.07806, -0.25979, 0.05629, -0.15849, -0.31256, 0.33553, -0.27883 +0.30682, -0.00369, -0.12624, -0.13524, -0.23072, -0.19830, 0.51922, 0.13625, 0.30686, 0.19326, -0.28604, -0.24671, 0.05676, -0.00147, -0.12379, 0.29854, 0.27577, -0.05369, 0.41714, 0.05100, 0.10535, 0.25372, 0.15772, -0.41833, 0.12765 +-0.34396, -0.00131, 0.32711, 0.23731, 0.13850, -0.19803, 0.26637, -0.13583, -0.09831, 0.08030, 0.14047, 0.16551, -0.20448, -0.08249, 0.02232, 0.00606, 0.09334, 0.22577, 0.01190, -0.32693, -0.13444, 0.32171, 0.01188, -0.09272, 0.14908 +-0.28358, -0.23711, -0.00079, -0.20363, 0.26583, -0.25923, -0.01992, 0.30822, 0.19845, -0.02683, 0.02838, 0.07485, 0.06014, 0.13685, -0.18069, 0.31956, -0.29260, 0.01786, -0.24777, -0.13140, -0.13986, -0.16778, 0.29813, -0.27355, 0.13748 +-0.07721, -0.33846, -0.22596, 0.14351, -0.03737, 0.09296, 0.34078, 0.02503, 0.17649, 0.24141, 0.26630, -0.27206, -0.30389, 0.28129, 0.33852, 0.27713, -0.22798, 0.20975, -0.27248, 0.25465, -0.09201, 0.20717, -0.02952, 0.06362, 0.28118 +0.11761, -0.26044, 0.31300, 0.18747, -0.12191, 0.23744, -0.42517, -0.11878, -0.15642, 0.22698, 0.28689, 0.28481, -0.03482, -0.12447, -0.32566, 0.20172, -0.14285, -0.16177, -0.11876, -0.10247, 0.22981, 0.01865, -0.20637, 0.14530, -0.30398 +0.14989, -0.12266, 0.21422, 0.31554, -0.30338, 0.13623, -0.22872, -0.08666, -0.02001, -0.12250, 0.32934, 0.31798, -0.10740, -0.15264, -0.08449, -0.27235, 0.10682, -0.23104, -0.29787, -0.17958, 0.24528, -0.02062, -0.01624, 0.14056, 0.19334 +0.16698, -0.06615, 0.18133, -0.20182, -0.00508, -0.06447, -0.31271, 0.05211, -0.19533, 0.30837, 0.31486, -0.16648, 0.08478, 0.03339, -0.09702, 0.16260, -0.18307, 0.22876, 0.15843, 0.14073, -0.13829, -0.29730, 0.11204, 0.34402, -0.17949 +0.34478, 0.02179, -0.06910, 0.22159, 0.35296, -0.34437, -0.17796, -0.13641, -0.01192, -0.01154, -0.27147, 0.03556, -0.28606, -0.33980, -0.03341, 0.01746, -0.11872, 0.19611, -0.13309, 0.31245, -0.34939, -0.07927, -0.02675, 0.32078, 0.31505 +0.23845, 0.07345, 0.07717, 0.09820, -0.20820, -0.14211, 0.05102, -0.08047, 0.20787, 0.26361, -0.00489, 0.04994, 0.05208, 0.00991, 0.18136, -0.19816, -0.13177, 0.13224, 0.22872, 0.17221, 0.02832, 0.18470, -0.14746, 0.00152, 0.24030 +-0.21038, -0.33311, 0.05829, -0.11248, -0.39031, 0.24018, 0.31533, 0.31591, 0.21842, 0.11106, 0.00419, -0.10125, -0.04118, -0.18811, -0.32428, -0.03926, 0.18926, -0.12291, 0.15375, 0.02808, -0.14635, -0.08132, -0.00497, -0.03953, 0.24568 +0.14025, 0.00668, -0.30074, -0.07440, -0.33066, -0.06600, 0.01166, 0.24954, -0.22035, -0.31518, -0.14028, -0.27341, -0.21291, -0.30773, -0.13274, 0.15092, 0.31991, 0.04984, 0.27492, 0.32448, -0.28637, -0.26806, 0.04600, -0.27472, -0.32435 +-0.34466, -0.28169, 0.16045, -0.29125, -0.30502, 0.21993, 0.00177, 0.00563, -0.15463, 0.06853, -0.13516, 0.02662, 0.18983, -0.22146, -0.14441, 0.06328, -0.05649, 0.30135, -0.23290, -0.21104, 0.01558, -0.07386, 0.31682, -0.33285, -0.10461 +0.26429, 0.11103, -0.28125, 0.20282, 0.09487, 0.00015, 0.17683, -0.19715, 0.17583, 0.14676, -0.26661, -0.21553, -0.23419, 0.12648, 0.02281, -0.34281, -0.14238, -0.11536, -0.29904, 0.24350, -0.06235, 0.04659, -0.08648, 0.21312, 0.23545 +-0.13509, -0.09735, 0.06823, -0.22581, -0.17998, 0.15447, -0.16209, 0.17372, -0.27225, -0.00504, 0.14097, -0.21481, -0.14895, 0.27991, -0.32028, 0.00831, -0.20118, -0.19568, -0.05819, -0.29290, -0.11851, 0.14796, -0.04694, 0.23116, -0.20026 +0.30448, -0.30411, 0.17794, -0.19358, 0.19374, 0.24929, 0.52103, 0.14052, 0.12245, 0.08575, 0.13888, -0.30510, 0.25353, -0.14709, 0.02053, -0.02730, -0.10952, -0.24077, 0.28567, 0.01464, 0.13864, 0.19052, 0.00932, -0.39407, -0.32367 +0.33511, -0.31479, 0.21451, 0.22260, 0.14512, -0.33519, -0.19355, 0.24852, 0.23520, -0.11799, -0.22526, -0.26930, 0.09369, -0.23061, 0.31694, 0.23137, 0.18504, -0.04852, -0.24657, -0.14619, -0.26419, 0.02569, -0.16507, -0.25164, 0.31506 +-0.10569, 0.00200, 0.23895, 0.12868, -0.01618, -0.01120, 0.09951, 0.01289, 0.19420, 0.14726, -0.11491, -0.08230, -0.06727, -0.06584, -0.23567, -0.18129, 0.12943, -0.14550, -0.19815, 0.18831, -0.32903, -0.05185, 0.20706, -0.18634, 0.25237 +0.08708, 0.02897, -0.32331, -0.19755, 0.14713, -0.24031, -0.39770, 0.07247, 0.29314, 0.27825, 0.34283, -0.17500, 0.14072, 0.17569, 0.06979, 0.11654, 0.08414, 0.16877, -0.38219, -0.16694, 0.01645, 0.18350, -0.10728, 0.36223, 0.32663 +0.24512, 0.04184, -0.31656, 0.19614, 0.17455, 0.00337, 0.03893, -0.33625, 0.01757, 0.03453, -0.00278, -0.00553, -0.15259, 0.31461, -0.29216, 0.22329, 0.11153, -0.27728, -0.10089, 0.02872, 0.25985, -0.03283, 0.27549, -0.28317, -0.04284 +0.30643, 0.26365, 0.07507, 0.07706, -0.14626, 0.15628, 0.31865, -0.16474, -0.15439, -0.02804, 0.01072, 0.31345, -0.25433, 0.33175, -0.09756, -0.20805, 0.02763, -0.28118, 0.26005, -0.12078, 0.20447, -0.31157, -0.11742, 0.07141, 0.18978 +-0.34193, -0.20700, 0.07865, -0.01669, 0.04845, -0.06684, -0.49604, -0.02675, 0.07013, -0.15788, 0.13861, -0.15236, -0.06848, -0.11084, 0.33909, 0.27674, -0.29966, 0.22394, 0.20092, -0.04978, 0.17215, 0.07136, -0.24517, 0.19800, -0.30192 +0.00368 +-0.00516 +-0.00375 +0.00510 +0.00306 +-0.00552 +-0.00378 +-0.00580 +0.00648 +-0.00696 +-0.00641 +0.00060 +-0.00563 +-0.00589 +0.00189 +0.00000 +-0.00396 +-0.00536 +0.00366 +-0.00466 +-0.00563 +0.01047 +-0.00733 +0.00197 +0.00208 +-0.39757, -0.47563, -0.13953, 0.32202, -0.25362, 0.32063, 0.01388, -0.12681, 0.43118, -0.33472, 0.09568, 0.00016, 0.45845, -0.01793, -0.00473, -0.30083, -0.16590, 0.32844, -0.36319, -0.33402, 0.45242, 0.34286, 0.05281, -0.12741, 0.40624 +-0.00490 diff --git a/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.onnx b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.onnx new file mode 100644 index 0000000000..25837086f4 Binary files /dev/null and b/code/nnv/examples/NNCS/Single_Pendulum/controller_single_pendulum.onnx differ diff --git a/code/nnv/examples/NNCS/Single_Pendulum/dynamics_sp.m b/code/nnv/examples/NNCS/Single_Pendulum/dynamics_sp.m new file mode 100644 index 0000000000..f427e77d44 --- /dev/null +++ b/code/nnv/examples/NNCS/Single_Pendulum/dynamics_sp.m @@ -0,0 +1,12 @@ +function [dx] = dynamics_sp(x,a) +% Ex_single_pendulum +% constants as per the documentation l=0.5, m=0.5; g= 1; c=0; +% description of the pararameters in these equations can be found here +% https://github.com/amaleki2/benchmark_closedloop_verification/blob/master/AINNC_benchmark.pdf + +l=0.5; m=0.5; g= 1; c=0; + +dx(1,1) = x(2); +dx(2,1) = g/l * sin(x(1)) + (a - c*x(2))/(m*l^2); + +end \ No newline at end of file diff --git a/code/nnv/examples/NNCS/Single_Pendulum/reach.m b/code/nnv/examples/NNCS/Single_Pendulum/reach.m new file mode 100644 index 0000000000..f01405a7b9 --- /dev/null +++ b/code/nnv/examples/NNCS/Single_Pendulum/reach.m @@ -0,0 +1,70 @@ +%% Reachability analysis of Single Pendulum Benchmark + +%% Load Components + +% Load the controller +net = load_NN_from_mat('controller_single_pendulum.mat'); +% Load plant +reachStep = 0.01; +controlPeriod = 0.05; +plant = NonLinearODE(2,1,@dynamics_sp, reachStep, controlPeriod, eye(2)); +% plant.set_tensorOrder(2); + + +%% Reachability analysis + +% Initial set +lb = [0; -0.1]; +ub = [1; 0.1]; +% lb = [1.1999; 0.1999; 0]; +% ub = [1.2; 0.2; 0]; +init_set = Box(lb,ub); +init_partitions = init_set.partition([1 2], [20, 8]); +% Store all reachable sets +outputPartitions = cell(length(init_partitions),1); +num_steps = 20; +reachOptions.reachMethod = 'approx-star'; +for s = 1:length(init_partitions) + % Get initial set for corresponding partition + init_set = init_partitions(s); + init_set = Star(init_set.lb, init_set.ub); + reachAll = init_set; + t = tic; + % Begin NNCS reachability for partition + for i = 1:num_steps + % Compute controller output set + input_set = net.reach(init_set,reachOptions); + input_set = Star.get_hypercube_hull(input_set); + input_set = input_set.toStar; + % Compute plant reachable set + init_set = plant.stepReachStar(init_set, input_set,'lin'); + init_set = Star.get_hypercube_hull(init_set); + init_set = init_set.toStar; + reachAll = [reachAll init_set]; + end + toc(t); + % store results for partition + outputPartitions{s} = reachAll; +end + +%% Visualize results +timeVector = 0:controlPeriod:controlPeriod*num_steps; + +f = figure; +hold on; +rectangle('Position',[0.5,1,1,1],'FaceColor',[1 0 0 0.5],'EdgeColor','r', 'LineWidth',0.1) +for i=1:length(outputPartitions) + Star.plotRanges_2D(outputPartitions{i},1,timeVector,'b'); + hold on; +end +grid; +xlabel('Time (s)'); +ylabel('\theta'); +% xlim([0 0.6]) +% ylim([0.95 1.25]) +% Save figure +if is_codeocean + exportgraphics(f,'/results/logs/singlePendulum.pdf', 'ContentType', 'vector'); +else + exportgraphics(f,'singlePendulum.pdf','ContentType', 'vector'); +end diff --git a/code/nnv/examples/NNCS/Single_Pendulum/reach_parallel.m b/code/nnv/examples/NNCS/Single_Pendulum/reach_parallel.m new file mode 100644 index 0000000000..ce8158db23 --- /dev/null +++ b/code/nnv/examples/NNCS/Single_Pendulum/reach_parallel.m @@ -0,0 +1,70 @@ +%% Reachability analysis of Single Pendulum Benchmark + +%% Load Components + +% Load the controller +net = load_NN_from_mat('controller_single_pendulum.mat'); +% Load plant +reachStep = 0.01; +controlPeriod = 0.05; +plant = NonLinearODE(2,1,@dynamics_sp, reachStep, controlPeriod, eye(2)); +% plant.set_tensorOrder(2); + + +%% Reachability analysis + +% Initial set +lb = [0; -0.1]; +ub = [1; 0.1]; +% lb = [1.1999; 0.1999; 0]; +% ub = [1.2; 0.2; 0]; +init_set = Box(lb,ub); +init_partitions = init_set.partition([1 2], [20, 8]); +% Store all reachable sets +outputPartitions = cell(length(init_partitions),1); +num_steps = 20; +reachOptions.reachMethod = 'approx-star'; +parfor s = 1:length(init_partitions) + % Get initial set for corresponding partition + init_set = init_partitions(s); + init_set = Star(init_set.lb, init_set.ub); + reachAll = init_set; + t = tic; + % Begin NNCS reachability for partition + for i = 1:num_steps + % Compute controller output set + input_set = net.reach(init_set,reachOptions); + input_set = Star.get_hypercube_hull(input_set); + input_set = input_set.toStar; + % Compute plant reachable set + init_set = plant.stepReachStar(init_set, input_set,'lin'); + init_set = Star.get_hypercube_hull(init_set); + init_set = init_set.toStar; + reachAll = [reachAll init_set]; + end + toc(t); + % store results for partition + outputPartitions{s} = reachAll; +end + +%% Visualize results +timeVector = 0:controlPeriod:controlPeriod*num_steps; + +f = figure; +hold on; +rectangle('Position',[0.5,1,1,1],'FaceColor',[1 0 0 0.5],'EdgeColor','r', 'LineWidth',0.1) +for i=1:length(outputPartitions) + Star.plotRanges_2D(outputPartitions{i},1,timeVector,'b'); + hold on; +end +grid; +xlabel('Time (s)'); +ylabel('\theta'); +% xlim([0 0.6]) +% ylim([0.95 1.25]) +% Save figure +if is_codeocean + exportgraphics(f,'/results/logs/singlePendulum.pdf', 'ContentType', 'vector'); +else + exportgraphics(f,'singlePendulum_parallel.pdf','ContentType', 'vector'); +end diff --git a/code/nnv/examples/NNCS/Single_Pendulum/reach_single.m b/code/nnv/examples/NNCS/Single_Pendulum/reach_single.m new file mode 100644 index 0000000000..20318272e6 --- /dev/null +++ b/code/nnv/examples/NNCS/Single_Pendulum/reach_single.m @@ -0,0 +1,57 @@ +%% Reachability analysis of Single Pendulum Benchmark + +%% Load Components + +% Load the controller +net = load_NN_from_mat('controller_single_pendulum.mat'); +% Load plant +reachStep = 0.01; +controlPeriod = 0.05; +plant = NonLinearODE(2,1,@dynamics_sp, reachStep, controlPeriod, eye(2)); +% plant.set_tensorOrder(2); + + +%% Reachability analysis + +% Initial set +lb = [0; -0.1]; +ub = [1; 0.1]; +init_set = Star(lb,ub); +% Store all reachable sets +num_steps = 20; +reachOptions.reachMethod = 'approx-star'; + +reachAll = init_set; +t = tic; +% Begin NNCS reachability for initial set +for i = 1:num_steps + % Compute controller output set + input_set = net.reach(init_set,reachOptions); + input_set = Star.get_hypercube_hull(input_set); + input_set = input_set.toStar; + % Compute plant reachable set + init_set = plant.stepReachStar(init_set, input_set,'lin'); + init_set = Star.get_hypercube_hull(init_set); + init_set = init_set.toStar; + reachAll = [reachAll init_set]; +end +toc(t); + +%% Visualize results +timeVector = 0:controlPeriod:controlPeriod*num_steps; + +f = figure; +hold on; +rectangle('Position',[0.5,1,1,1],'FaceColor',[1 0 0 0.5],'EdgeColor','r', 'LineWidth',0.1) +Star.plotRanges_2D(reachAll,1,timeVector,'b'); +grid; +xlabel('Time (s)'); +ylabel('\theta'); +% xlim([0 0.6]) +% ylim([0.95 1.25]) +% Save figure +if is_codeocean + exportgraphics(f,'/results/logs/singlePendulum.pdf', 'ContentType', 'vector'); +else + exportgraphics(f,'singlePendulum_single.pdf','ContentType', 'vector'); +end diff --git a/code/nnv/examples/NNCS/Single_Pendulum/singlePendulum.pdf b/code/nnv/examples/NNCS/Single_Pendulum/singlePendulum.pdf new file mode 100644 index 0000000000..17bb6014e6 Binary files /dev/null and b/code/nnv/examples/NNCS/Single_Pendulum/singlePendulum.pdf differ diff --git a/code/nnv/examples/NNCS/Single_Pendulum/specifications.txt b/code/nnv/examples/NNCS/Single_Pendulum/specifications.txt new file mode 100644 index 0000000000..4d962f99ca --- /dev/null +++ b/code/nnv/examples/NNCS/Single_Pendulum/specifications.txt @@ -0,0 +1,12 @@ +Initial states: + +x1 = [1.0, 1.2] +x2 = [0.0, 0.2] + +The controller timestep (control period) = 0.05 (seconds). + +Property: + +x1 in [0.0,1.0] for 0.5 <= t <=1 + +Refer to this for more details: https://github.com/amaleki2/benchmark_closedloop_verification/blob/master/AINNC_benchmark.pdf