The analysis is aborted because the time step size converges to 0

[akaph2p]

Hi,

I have a nonlinear ODE system with 8 states. When i try running CORA I receive this error:

Warning: Check computation for order of abstraction error 
Warning: The analysis is aborted because the time step size converges to 0.
         The reachable sets until t = 1.4257e-17 are returned. 

Are there ways to fix this while keeping the ODEs nonlinear? Or is this a limitation of MATLAB itself?

[akaph2p]

clear
close all
clc
format 'long'

% Parameters
alpha = 2e-5;  
T0 = 652;    
BETA = 0.0067;   
LAMBDA = 10e-2;      
lambda1 = 0.08;
beta1 = 0.0067; 
%T_in = 650;   % C
P0 = 42e9; % Watt

mp = 2215;       
Cp = 42e6/17; %J/kg K, FLiBe    

H = 0.9104;                                                                  % Heat transfer coefficient, KW/m^2-K
A = 2*pi*0.02385/2*4.5*1777;                                                % (2*pi*r*h)*number of tubes, m^2
r = 1.1; % pressure ratio across the turbine = inlet pressure /outlet pressure
eta = 1; % Turbine isentropic efficiency
Cturb = 0.55;
ms = 3047.4; %Kg
Cs = 4.3; %KJ/Kg-K

%% Highpower ss, controlling Trx

% Parameters --------------------------------------------------------------
params.tFinal = 200;
n0 = 1;
C10 = n0*beta1/LAMBDA/lambda1;
Trx0 = 675;
eint0 = 4.6e-4/0.0001;
mdotp0 = 340;
Tp0 = 650;
Ts0 = 400;
mdots0 = 340;

params.R0 = zonotope([  [n0; C10; Trx0; eint0; mdotp0; Tp0; Ts0; mdots0], diag([0,0,0,0,0,0,0,0])   ]);
params.U = zonotope([0,0]);

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

options.timeStep = 0.01;
options.alg = 'lin-adaptive';
options.maxError = [0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01]; 

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

reactor = nonlinearSys(@modular_system_odes,8,0);

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

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

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

[akaph2p]

function dx = modular_system_odes(x,u)

% Parameters
alpha = 2e-5;  
T0 = 652;    
BETA = 0.0067;   
LAMBDA = 10e-2;      
lambda1 = 0.08;
beta1 = 0.0067; 
%T_in = 650;   % C
P0 = 42e9; % Watt

mp = 2215;       
Cp = 42e6/17; %J/kg K, FLiBe    

H = 0.9104;                                                                  % Heat transfer coefficient, KW/m^2-K
A = 2*pi*0.02385/2*4.5*1777;                                                % (2*pi*r*h)*number of tubes, m^2
r = 1.1; % pressure ratio across the turbine = inlet pressure /outlet pressure
eta = 1; % Turbine isentropic efficiency
Cturb = 0.55;
ms = 3047.4; %Kg
Cs = 4.3; %KJ/Kg-K

% Control
n = x(1);
C1 = x(2);
Trx = x(3);
error_int = x(4);
mdotp = x(5);
Tp = x(6);
Ts = x(7);
mdots = x(8);

% Reactivity control
Tref = 675;
error = Tref - Trx;
KP = 0.2;
KI = 0.0001;
rhoc = KP*(error) + KI*error_int;
rho = rhoc -alpha*(Trx-T0);

% Primary flow rate control
% mdotpref = 340;
% error2 = mdotref - mdotp;
% KP2 = 60;
% KI = 100;
% mdotp = KP2*erro2 + KI2*error_int2;

Thot = 2*Trx - Tp;
Tpin = Thot;
Tsin = (Ts)*(1-eta*(1-(1/r)));

% differential equations
% Reactor
dx(1,1) = n*(rho - BETA)/LAMBDA + lambda1*C1;
dx(2,1) = beta1*n/LAMBDA - lambda1*C1;
dx(3,1) = P0*n/(mp*Cp) - 2*mdotp*(Trx-Tp)/mp;
dx(4,1) = Tref-Trx;
dx(5,1) = 0;
% Heat exchanger
dx(6,1) = mdotp*(Tpin - Tp)/mp - H*A*(Tpin - Tsin)/(mp*Cp);
dx(7,1) = mdots*(Tsin - Ts)/ms + H*A*(Tpin - Tsin)/(ms*Cs);
dx(8,1) = 0;
% turbine
%dx(9,1) = mdots*Cturb*(Ts - Tsin);
end

[wetzlingerm]

Hi,

Thank you for contacting us! Unfortunately, the adaptive algorithm is not well-equipped to deal with analyses that have a degenerate initial set, in your case, a single point. Also, note that setting options.timeStep = 0.01 has no effect when setting options.alg = 'lin-adaptive' since the time step size is tuned automatically.

The manually tuned algorithm parameters

options.timeStep = 0.2;
options.taylorTerms = 6;
options.zonotopeOrder = 100;
options.alg = 'lin';
options.tensorOrder = 2;

yield a tight reachable set enclosure. The computation time was about 12 seconds on my local machine.

Please note that plotting of a large number of reachable sets is prohibitively slow in CORA v2024.1.5, which we have fixed in v2024.1.6.

[akaph2p]

Thank you for your prompt response. In the past, i have found it difficult to manually find the optimal parameters for the 'lin' algorithm, so Ive let the adaptive algorithm find them.

Is there a systematic method to tune the linear algorithm myself or is it a trial-and-error process?

[wetzlingerm]

There is no clear-cut way of finding suitable values for the algorithm parameters for nonlinear systems. Generally, larger zonotope orders options.zonotopeOrder, options.errorOrder, and options.intermediateOrder yield tighter enclosures, and options.tensorOrder = 3 outperforms options.tensorOrder = 2 in tightness at the cost of more computational effort. The time step size is the most important algorithm parameter; here one usually starts with large time step sizes and then reduces the value until no significant improvement is noticeable. Since this knowledge is fairly case-dependent, we have developed our adaptive algorithm to aid users by automatically finding suitable values for all algorithm parameters. Unfortunately, we cannot yet cover all cases.

I hope this helps in future investigations!

[akaph2p]

I applied the adaptive algorithm on my system with one state being a set with range 10. My resulting plots look incomplete. Seems like the reachable set is only computed for the first ~ 9 seconds, but the simulations are computed for the whole time span. Could you please help me understand if im doing something wrong?

[akaph2p]

Here is my updated code and ode function. Looks my the reachable set is going exponential.

clear
close all
clc
format 'long'

% Parameters
alpha = 2e-5;  
T0 = 652;    
BETA = 0.0067;   
LAMBDA = 10e-2;      
lambda1 = 0.08;
beta1 = 0.0067; 
P0 = 42e9; % Watt

mp = 2215;       
Cp = 42e6/17; %J/kg K, FLiBe    

H = 420000;                                                            % Heat transfer coefficient, W/m^2-K
A = 1000;                                             % (2*pi*r*h)*number of tubes, m^2
r = 1.1; % pressure ratio across the turbine = inlet pressure /outlet pressure
eta = 1; % Turbine isentropic efficiency
Cturb = 0.55;
ms = 3055.7; %Kg 
Cs = 42e6/8.5; %J/Kg-K

%% Highpower ss, controlling Trx

% Parameters --------------------------------------------------------------
params.tFinal = 200;
n0 = 1;
C10 = n0*beta1/LAMBDA/lambda1;
Trx0 = 675;
eint0 = 4.6e-4/0.001;
mdotp0 = 340;
Tp0 = 650;
Ts0 = 625;
mdots0 = 340;

params.R0 = zonotope([  [n0; C10; Trx0; eint0; mdotp0; Tp0; Ts0; mdots0], diag([0,0,0,0,0,0,0,50])   ]);
params.U = zonotope([0,0]);

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

options.alg = 'lin-adaptive';
%options.maxError = [100;100;100.001;100.001;100.001;100.001;100.001;100.001]; 
% % 
% options.timeStep = 0.2;
% options.taylorTerms = 6;
% options.zonotopeOrder = 100;
% options.alg = 'lin';
% options.tensorOrder = 2;

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

reactor = nonlinearSys(@reach_modular_system_odes,8,0);

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

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

% Simulation --------------------------------------------------------------
tic
simOpt.points = 10;
simRes = simulateRandom(reactor, 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','LineWidth',2);
sgtitle('n vs time','FontSize', 20)
ylim([0 2.5])
xlim([0 200])
set(gca,'FontSize',20)
toc
disp(['computation time of n: ',num2str(toc/60) ' mins']);

%%
tic
figure(2)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,2,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,2, 'DisplayName','Initial set');
% plot simulation results      
plotOverTime(simRes,2,'DisplayName','Simulations','lineWidth',2);
ylim([0 2.5])
xlim([0 200])
sgtitle('C1 vs time', 'Fontsize', 20')
set(gca,'FontSize',20)
disp(['computation time of C1: ',num2str(toc/60) ' mins']);

%%
tic
figure(3)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,3,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,3, 'DisplayName','Initial set');
% plot simulation results     
plotOverTime(simRes,3,'DisplayName','Simulations', 'LineWidth',2);
sgtitle('Trx vs time','Fontsize', 20)
ylim([674 676])
xlim([0 200])
set(gca,'FontSize',20)
%ylim([670 680])
disp(['computation time of Trx: ',num2str(toc/60) ' mins']);

%%
tic
figure(4)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,4,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,4, 'DisplayName','Initial set');
% plot simulation results     
plotOverTime(simRes,4,'DisplayName','Simulations', 'LineWidth', 2);
sgtitle('rhoc/Ki vs time', 'FontSize', 20)
ylim([0 1.1])
xlim([0 200])
set(gca,'FontSize',20)
disp(['computation time of rhoc/Ki: ',num2str(toc/60) ' mins']);

%%
tic
figure(5)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,5,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,5, 'DisplayName','Initial set');
% plot simulation results     
plotOverTime(simRes,5,'DisplayName','Simulations', 'LineWidth', 2);
sgtitle('mdotp vs time', 'FontSize', 20)
set(gca,'FontSize',20)
disp(['computation time of mdotp: ',num2str(toc/60) ' mins']);

%%
tic
figure(6)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,6,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,6, 'DisplayName','Initial set');
% plot simulation results     
plotOverTime(simRes,6,'DisplayName','Simulations', 'LineWidth', 2);
sgtitle('Tp vs time', 'FontSize', 20)
ylim([615 660])
xlim([0 200])
set(gca,'FontSize',20)
disp(['computation time of Tp: ',num2str(toc/60) ' mins']);

%%
tic
figure(7)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,7,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,7, 'DisplayName','Initial set');
% plot simulation results     
plotOverTime(simRes,7,'DisplayName','Simulations', 'LineWidth', 2);
sgtitle('Ts vs time', 'FontSize', 20)
ylim([550 650])
xlim([0 200])
set(gca,'FontSize',20)
disp(['computation time of Ts: ',num2str(toc/60) ' mins']);

%%
tic
figure(8)
hold on; box on;
% plot reachable sets
useCORAcolors("CORA:contDynamics")
plotOverTime(R,8,'DisplayName','Reachable set');
% plot initial set
plotOverTime(R.R0,8, 'DisplayName','Initial set');
% plot simulation results     
plotOverTime(simRes,8,'DisplayName','Simulations', 'LineWidth', 2);
sgtitle('mdots vs time', 'FontSize', 20)
set(gca,'FontSize',20)
disp(['computation time of mdots: ',num2str(toc/60) ' mins']);

function dx = reach_modular_system_odes(x,u)

% Parameters
alpha = 2e-5;  
T0 = 652;    
BETA = 0.0067;   
LAMBDA = 10e-2;      
lambda1 = 0.08;
beta1 = 0.0067; 
P0 = 42e9; % Watt

mp = 2215;       
Cp = 42e6/17; %J/kg K, FLiBe    

H = 420000;                                                            % Heat transfer coefficient, W/m^2-K
A = 1000;                                             % (2*pi*r*h)*number of tubes, m^2
r = 25/24; % pressure ratio across the turbine = inlet pressure /outlet pressure
eta = 1; % Turbine isentropic efficiency
Cturb = 0.55;
ms = 3055.7; %Kg 
Cs = 42e6/8.5; %J/Kg-K 

% Control
n = x(1);
C1 = x(2);
Trx = x(3);
error_int = x(4);
mdotp = x(5);
Tp = x(6);
Ts = x(7);
mdots = x(8);

% Reactivity control
Tref = 675;
error = Tref - Trx;
KP = 0.2;
KI = 0.001;
rhoc = KP*(error) + KI*error_int;
rho = rhoc -alpha*(Trx-T0);

Tpin = 2*Trx - Tp;
Tsin =(Ts)*(1-eta*(1-(1/r))); % 600 at ss

% differential equations
% Reactor
dx(1,1) = n*(rho - BETA)/LAMBDA + lambda1*C1;
dx(2,1) = beta1*n/LAMBDA - lambda1*C1;
dx(3,1) = P0*n/(mp*Cp) - 2*mdotp*(Trx-Tp)/mp;
dx(4,1) = Tref-Trx;
dx(5,1) = 0;
% Heat exchanger
dx(6,1) = mdotp*(Tpin - Tp)/mp - H*A*(Tpin - Tsin)/(mp*Cp);
dx(7,1) = mdots*(Tsin - Ts)/ms + H*A*(Tpin - Tsin)/(ms*Cs);
dx(8,1) = 0;
% turbine
%dx(9,1) = mdots*Cturb*(Ts - Tsin);
end

[wetzlingerm]

When I run your updated code, I get the warning below:

Warning: The analysis is aborted because the time step size converges to 0.
         The reachable sets until t = 8.555 are returned.

The time step size converges to 0 because the size of the reachable set explodes (this effect can be observed in the plot you attached). To avoid making large approximation errors, the time step size has to be chosen very small, but the damage is already done, so the algorithm exits and returns the reachable sets until that point in time.

Reachability analysis for nonlinear systems is an active area of research. Unfortunately, no algorithm to date always returns a tight reachable set enclosure for any system under any uncertainty. The most common failure is that the approximation error grows exponentially, which happened in your case. This issue can be mitigated by finding better algorithm parameter valuations, e.g., for the time step size, or by decreasing the uncertainty in the model parameters, such as the initial set. As a last resort, one may also split the initial set into smaller initial sets and compute them in parallel or sequentially, but this scales exponentially in the state dimension due to the curse of dimensionality.

[akaph2p]

I see, thank you for the explanation. I will try the mitigation techniques you mentioned.