TUMcps / CORA

Toolbox for Reachability Analysis
GNU General Public License v3.0
108 stars 35 forks source link

Non-linear hybrid automaton - error in abstractionError_adaptive function #45

Closed akaph2p closed 11 months ago

akaph2p commented 12 months ago

Im trying to run reach on a nonlinear hybrid automaton, but i receive a matrix dimension error in the abstractionError_adaptive function.

My code:

% Dynamics 
pke = nonlinearSys(@pke,4,0);

% Invariant set
inv = polytope([-1,0,0,0; 0,-1,0,0; 0,0,-1,0; 0,0,0,-1],[0;0;0;-1]); % n>=0, C1>=0, Trx>=0, rhoc>=-1

%%
% Locations/ nodes
guard1 = conHyperplane([0,0,-1,0],675); % transition to 2 when Trx >= 675
reset1.f = @(x,u)([x(1); x(2); x(3); 0.66 + 1.2*(675-x(3))]);
trans1 = transition(guard1,reset1,2); % --> next loc: 2
loc(1) = location('startup',inv,trans1,pke); % startup 

guard2 =  conHyperplane([0,0,-1,0],810); % transition to 3 when Trx >= 1.2*675
reset2.f = @(x,u)([x(1); x(2); x(3); -1]);
trans2 = transition(guard2,reset2,3); % --> next loc: 3
loc(2) = location('steadystate',inv,trans2,pke);  % steady state

guard3 =  conHyperplane([0,0,1,0],650); % transition to 1 when Trx <= 650
reset3.f = @(x,u)([x(1); x(2); x(3); 0.7]);
trans3 = transition(guard3,reset3,1); % --> next loc: 1
loc(3) = location('scram',inv,trans3,pke);  % SCRAM

HA = hybridAutomaton(loc);

%% Reachability
clc
% Parameters --------------------------------------------------------------

params.tFinal = 400;
params.R0 = zonotope([[0; 0.0012; 650; 0.7], diag([0,0,0,0])]);
params.startLoc = 1;
parama.U = zonotope(0,0);
% Reachability Settings ---------------------------------------------------

options.timeStep = 1;
options.taylorTerms = 10;
options.zonotopeOrder = 50;
options.alg = 'lin-adaptive';
options.tensorOrder = 2;
options.lagrangeRem.simplify = 'simplify';

% settings for hybrid systems
options.guardIntersect = ['polytope'];
options.enclose = {'box'}; 

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

tic
R = reach(HA, params, options);
tComp = toc;
disp(['computation time of reachable set: ',num2str(tComp)]);

The error:

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise
multiplication, use '.*'.

Error in abstractionError_adaptive (line 101)
            err(i,1) = 0.5 * dz' * options.hessianConst{i} * dz;

Error in nonlinearSys/linReach_adaptive>aux_initStepTensorOrder (line 419)
[~,~,L0_2,options] = abstractionError_adaptive(obj,options,Rdelta);

Error in nonlinearSys/linReach_adaptive (line 98)
    options = aux_initStepTensorOrder(obj,options,Rstart);

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

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

Error in location/reach (line 56)
R = reach(loc.contDynamics,params,options_,spec);

Error in hybridAutomaton/reach (line 136)
            [Rtemp,Rjump,res] = reach(HA.location(locID),R0,tStart,options);

Error in reactor_HA (line 53)
R = reach(HA, params, options);

Related documentation
akaph2p commented 12 months ago

my nonlinear system is:

% Inputs:
%    x - [n, C1, Trx, rhoc]
% 
%
% Outputs:
%    dx - [n',  C1', Trx', rhoc']

% Parameters
alpha = 0.001;  
T0 = 0;    
BETA = 0.01;   
l = 1;      
lambda1 = 1;
beta1 = 0.0248; 
Tc = 650;   
P0 = 250;

mdotp = 5;  
m = 1;       
Cp = 1;    
q0 = 0;

% Control
n = x(1);
C1 = x(2);
Trx = x(3);
rhoc = x(4);
% differential equations
dx(1,1) = n*(rhoc - alpha*Trx - alpha*T0 - BETA)/l + lambda1*C1 + q0;
dx(2,1) = beta1*n/l - lambda1*C1;
dx(3,1) = -2*mdotp*Trx/m + 2*mdotp*Tc/m + P0*n/(m*Cp);
dx(4,1) = 0;
wetzlingerm commented 12 months ago

Thank you for reporting this issue. Please define your dynamics-file (pke.m) as function dx = pke(x,u) if you haven't already. Then, instantiate the nonlinear system using pke = nonlinearSys(@pke,4,1); with one input. CORA always expects at least one input so that internal computations don't have to differ between systems with and without input. We acknowledge that this is confusing since your system is obviously not influenced by external inputs -- a change will be made in the next update so that your code executes with your current setting as well.

akaph2p commented 11 months ago

That fixed the issue, thank you