I want to do reachability analysis of a non linear system using CORA. Can anyone guide me about invariant and reset conditions?

%% i have tried to put the lower limit of 'b' that it should not go below zero in my invariant. But it seems it is not working

%WPP Hybrid Automata

% System Dynamics ---------------------------------------------------------

% parameter

j=3410432;

R=35;

wr=2.58;

% b w bd d V t

% differential equations

f1 = @(x,u) [...

-5;

0;

0;

0;

0;

1;

];

m1 = nonlinearSys('m1',f1);

syms b w bd d V t;

inv = levelSet([0-b; 0-w;w-0.2314],[b;w;bd;d;V;t],'<=');

resetA = [1 0 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1];

resetc = [0;0.2314;0;0;0;0];

reset = struct('A', resetA, 'c', resetc);

guard = levelSet(b-0,[b;w;bd;d;V;t],'==');

trans = transition(guard, reset, 2);

loc(1) = location('loc1',inv,trans,m1);

%----------------------------------------------------------------------

f2 = @(x,u) [...

0;

(0.5*3.14*R^3*x(5)^2*1.225*(0.00006384*(x(2)*R/x(5))^4- 0.00176*(x(2)*R/x(5))^3+0.01452*(x(2)*R/x(5))^2- 0.03162*(x(2)*R/x(5))+0.02551))/j;

0;

0;

0;

1;

];

m2 = nonlinearSys('m2',f2);

inv = levelSet([b-0;-w+0;-w+0.2314],[b;w;bd;d;V;t],'<=');

resetA = eye(6);

resetc = zeros(6,1);

reset = struct('A', resetA, 'c', resetc);

guard = levelSet(-V+13,[b;w;bd;d;V;t],'<=');

trans = transition();

loc(2) = location('loc2',inv,trans,m2);

% Parameters --------------------------------------------------------------

params.tFinal = 50;

params.R0 = zonotope(interval([88;0;0;0;12;0],[90;0.1;0;0;12;0]));

params.startLoc = 1;

% Reachability Settings ---------------------------------------------------

options.timeStep = 0.1;

options.taylorTerms = 5;

options.zonotopeOrder = 5;

options.alg = 'lin';

options.tensorOrder = 2;

options.lagrangeRem.simplify = 'simplify';

% settings for hybrid systems

options.guardIntersect = 'polytope';

options.enclose = {'box'};

% Reachability Analysis ---------------------------------------------------

HA = hybridAutomaton(loc);

R = reach(HA, params, options);

% Simulation --------------------------------------------------------------

%simOpt.points = 10;

%simRes = simulateRandom(loc, params, simOpt);

% Visualization -----------------------------------------------------------

dims = {[6 1],[6 2],[6 3]};

figure;

for k = 1:length(dims)

subplot(1,3,k); hold on; box on;

projDim = dims{k};

% plot reachable sets

useCORAcolors("CORA:contDynamics")

plot(R,projDim,'DisplayName','Reachable set');

%

%plot initial set

%plot(R(1).R0,projDim, 'DisplayName','Initial set');

%

% % plot simulation results

% %plot(simRes,projDim,'DisplayName','Simulations');

%

% % label plot

% xlabel(['x_{',num2str(projDim(1)),'}']);

% ylabel(['x_{',num2str(projDim(2)),'}']);

% legend()

end

[Torsten]

Contact

• Prof. Dr.-Ing. Matthias Althoff
• Dr. Niklas Kochdumper
• Dr.-Ing. Mark Wetzlinger
• Tobias Ladner, M.Sc. (primary contact)
• Research Group