search

TUMcps / CORA

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

Array indices error in hybrid automaton reach #49

Closed akaph2p closed 9 months ago

akaph2p commented 9 months ago

I created a linear hybrid automaton for my system. But when I try to run reach i receive this error

Array indices must be positive integers or logical values.

Error in reach_adaptive>aux_ebar (line 1430)
ebarredt_max = ebar.red_e(t_idx-1) + ...

Error in reach_adaptive (line 149)
        ebar = aux_ebar(t,k,e,ebar,timeStep,options.tFinal);

Error in linearSys/reach (line 83)
            [timeInt,timePoint,res] = reach_adaptive(obj, 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_linear_hybrid_automaton (line 202)
R = reach(HA, params, options);

I have attached my .m file here. This error has not happened in the past when i created other automatons.

akaph2p commented 9 months ago 
% Parameters 
T_in = 650;
alpha = 2e-5;  
T0 = 652;    
BETA = 0.0067;   
LAMBDA = 10e-2;      
lambda1 = 0.08;
beta1 = 0.0067; 
P0 = 42e9; % Watt
m = 2215;       
Cp = 42e6/17; %J/kg K, FLiBe    

%% Low power, controlling n

% xdot = Ax + Bu;
% where x = [n~, C1~, Trx~, eint1, eint2, mdotp~]
% and u = [nref~, mdotref~]
% I'm modulating rhoc~ and mdotp~
% error1 = nref~ - n~
% error2 = mdotref~ - mdotp~

% n~` = (rho0-BETA - n0 Kp1)/LAMBDA)n~ + lambda1*C1~ - alpha n0/LAMBDA Trx~ + n0/LAMBDA Ki1 eint1 + n0/LAMBDA Kp1 nref~
% C1~` = n~ beta1/LAMBDA - lambda1 C1~
% Trx~` = p0n~/mCp - 2/m mdotp0 Trx~ - 2/m (Trx0-Tin)Ki2 eint2 + 2/m (Trx0-Tin) Kp2 mdotp~ - 2/m (Trx0-Tin) mdotref~
% eint1` = -n~ + nref~
% eint2` = -mdotp~ + mdotref~
%mdotp~` = 0

% Operating point
n_0 = 0.1; % at n = 0.1, Trx = 652.5 
C1_0 = n_0*beta1/lambda1/LAMBDA; 
Trx_0 = 652.5; 
rhoc_0 = 1e-5;
rho_0 = rhoc_0 - alpha*(Trx_0-T0);
mdotp_0 = 340; 

% PI parameters
Kp1 = 0.2;
Ki1 = 0.0001;
Kp2 = 0.5;
Ki2 = 0.15;

A_lp = [(rho_0-BETA - n_0*Kp1)/LAMBDA,     lambda1,    -alpha*n_0/LAMBDA,  n_0/LAMBDA*Ki1,     0,                      0;
    beta1/lambda1,                      -lambda1,   0,                  0,                  0,                      0;
    P0/(m*Cp),                          0,          -2*mdotp_0/m,       0,                  -2/m*(Trx_0-T_in)*Ki2,  2/m*(Trx_0-T_in)*Kp2;
    -1,                                 0,          0,                  0,                  0,                      0;
    0,                                  0,          0,                  0,                  0,                      -1;
    0,                                  0,          0,                  0,                  0,                      0];

B_lp = [n_0/LAMBDA*Kp1, 0;
    0,               0;
    0,               -2/m*(Trx_0-T_in);
    1,               0;
    0,               1;
    0,               0];

reactor_lp = linearSys(A_lp,B_lp);

%% High power approach, controlling n
% Operating point
n_0 = 1; % at n = 1, Trx = 675
C1_0 = n_0*beta1/lambda1/LAMBDA; 
Trx_0 = 675; 
rhoc_0 = 4.6e-4;
rho_0 = rhoc_0 - alpha*(Trx_0-T0);
mdotp_0 = 340; 

% PI parameters
Kp1 = 0.006;
Ki1 = 0.0000001;
Kp2 = 0.5;
Ki2 = 0.15;

A_hpa = [(rho_0-BETA - n_0*Kp1)/LAMBDA      lambda1     -alpha*n_0/LAMBDA                   n_0/LAMBDA*Ki1  0                       0;
    beta1/lambda1                       -lambda1    0                   0               0                       0;
    P0/(m*Cp)                           0           -2*mdotp_0/m        0               -2/m*(Trx_0-T_in)*Ki2   2/m*(Trx_0-T_in)*Kp2;
    -1                                  0           0                   0               0                       0;
    0                                   0           0                   0               0                       -1;
    0                                   0           0                   0               0                       0];

B_hpa = [n_0/LAMBDA*Kp1 0;
    0               0;
    0               -2/m*(Trx_0-T_in);
    1               0;
    0               1;
    0               0];

reactor_hpa = linearSys(A_hpa,B_hpa);

%% High power steady state, controlling Trx

% xdot = Ax + Bu;
% where x = [n~, C1~, Trx~, eint1, eint2, mdotp~]
% and u = [Tref~, mdotref~]
% I'm modulating rhoc~ and mdotp~
% error1 = Tref~ - Trx~
% error2 = mdotref~ - mdotp~

% n~` = (rho0-BETA/LAMBDA)n~ + lambda1*C1~ + (- alpha n0 - Kp1 n0)/LAMBDA Trx~ + n0/LAMBDA Ki1 eint1 + n0/LAMBDA Kp1 Tref~
% C1~` = n~ beta1/LAMBDA - lambda1 C1~
% Trx~` = p0n~/mCp - 2/m mdotp0 Trx~ - 2/m (Trx0-Tin)Ki2 eint2 + 2/m (Trx0-Tin) Kp2 mdotp~ - 2/m (Trx0-Tin) mdotref~
% eint1` = -Trx~ + Tref~
% eint2` = -mdotp~ + mdotref~
%mdotp~` = 0

% Operating point
n_0 = 1; 
C1_0 = 0.8375; 
Trx_0 = 675; 
rhoc_0 = 4.6e-4;
rho_0 = rhoc_0 - alpha*(Trx_0-T0);
mdotp_0 = 340; 

% PI parameters
Kp1 = 0.0005;
Ki1 = 0.00001;
Kp2 = 0.5;
Ki2 = 0.15;
Tref = 0;
mdotref = 0;

A_hpss = [(rho_0-BETA)/LAMBDA    lambda1     (-alpha*n_0/LAMBDA)-n_0/LAMBDA*Kp1     n_0/LAMBDA*Ki1  0                       0;
    beta1/lambda1           -lambda1    0                   0               0                       0;
    P0/(m*Cp)               0           -2*mdotp_0/m        0               (-2/m)*(Trx_0-T_in)*Ki2   2/m*(Trx_0-T_in)*Kp2;
    0                       0           -1                  0               0                       0;
    0                       0           0                   0               0                       -1;
    0                       0           0                   0               0                       0];

B_hpss = [n_0/LAMBDA*Kp1 0;
    0               0;
    0               -2/m*(Trx_0-T_in);
    1               0;
    0               1;
    0               0];

reactor_hpss = linearSys(A_hpss,B_hpss);

%% Hybrid automaton 

% Invariant set
inv = polytope([-1,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,0; 
                0,0,0,0,0,-1;],[-1;-1;-5;-1;-1;-5]); % n~>=-1, C1>=-1, Trx>=-5, eint1 >= -1, eint2 >= -1, mdotp >= -5

% Locations/ nodes
guard1 =  conHyperplane([-1,0,0,0,0,0],0,[-1,0,0,0,0,0],0.1); % transition to 2 when n >= 0.1, n~ >= 0
reset1.A = [1 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 0; 0 0 0 0 0 1];
reset1.c = [0, 0, 0, 0, 0, 0];
%reset1.f = @(x,u)([x(1); x(2); x(3); x(4); x(5); x(6)]);
trans1 = transition(guard1,reset1,2); % --> next loc: 2
loc(1) = location('lowpower',inv,trans1,reactor_lp);  % low power

guard2 =  conHyperplane([-1,0,0,0,0,0],0,[-1,0,0,0,0,0],0); % transition to 3 when n >=1, n~>=0
reset2.A = [1 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 0; 0 0 0 0 0 1];
reset2.c = [0, 0, 0, 0, 0, 0];
%reset2.f = @(x,u)([x(1); x(2); x(3); x(4); x(5); x(6)]);
trans2 = transition(guard1,reset2,3); % --> next loc: 3
loc(2) = location('highpowerapproach',inv,trans2,reactor_hpa);  % high power approach

guard3 =  conHyperplane([-1,0,0,0,0,0],10,[-1,0,0,0,0,0],10); % transition to 1 when n >=10, n~>=10
reset3.A = [1 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 0; 0 0 0 0 0 1];
reset3.c = [0, 0, 0, 0, 0, 0];
%reset3.f = @(x,u)([x(1); x(2); x(3); x(4); x(5); x(6)]);
trans3 = transition(guard3,reset3,1); % --> next loc: 1
loc(3) = location('highpowerss',inv,trans3,reactor_hpss);  % high power ss

HA = hybridAutomaton(loc);

%% Reachability

% Parameters --------------------------------------------------------------
clc
params.tFinal = 200;
params.R0 = zonotope([[-0.1; -0.08375; -2.5; -0.3; 0; 0], diag([0,0,0,0,0,0])]); % states: n~, C~, Trx~, eint1, eint2, mdotp~. Initially, n = 0, C1 = 0, Trx = 650, eint1 = -0.04e-3/0.001, eint2 = 0, mdotp~ = 0. Operating point n0 = 0.1, C10 = 0.0375, Trx0 = 652.5
params.startLoc = 1;

% Reachability Settings ---------------------------------------------------
%options.timeStep = 0.1;
%options.taylorTerms = 4;
%options.zonotopeOrder = 20;
options.linAlg = ('adaptive');

% settings for hybrid systems
options.guardIntersect = 'zonoGirard';
options.enclose = {'box'}; 
options.error = 1;

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

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

% Simulation --------------------------------------------------------------
simOptions.startLoc = params.startLoc; 
tic
simOpt.points = 5;
simRes = simulateRandom(HA, params, simOpt);
tComp = toc;
disp(['computation time of simulation: ',num2str(tComp/60) ' mins']);

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

tic
figure(1)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,1,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,1, 'DisplayName','Initial set');
% plot simulation results      
plotOverTime(simRes,1,'DisplayName','Simulations');
sgtitle('n~ vs time')
%toc
disp(['computation time of n: ',num2str(toc/60) ' mins']);
wetzlingerm commented 9 months ago

There was indeed a bug in reach_adaptive, which is now fixed in CORA v2024.1.3.

I would also like to point out a modeling error in the reset functions: The vectors in reset1.c, reset2.c, and reset3.c should be column vectors, not row vectors. The new version also detects that error now. FYI, I have noticed that guard2 is not used in trans2.

Thank you for helping us improve our code!

akaph2p commented 9 months ago

Thank you for catching my modeling errors, and thank you for your prompt response.