Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN ODE15s

Hi everyone!

I'm writing a code to simulate a car's 2D rigid motion: I need to contemporarily solve 35 equations, 7 differential and 28 algebraic, so I'm using ODE15s.

Simpler version of this model (up to 19 equations) worked (after long fights), but now I'm stuck: after declaring where differential equations are in the M matrix, Matlab warns me that:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
> In ode15s (line 589) 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

What shall I do to sort this out?

Here's the code (pacejka and slittamenti are functions):

M(13,13)=1;
M(14,14)=1;
M(15,15)=1;
M(24,24)=1;
M(25,25)=1;
M(26,26)=1;
M(27,27)=1;

options = odeset('Mass',M,'RelTol',1e-8,'AbsTol',tolleranze);

y0_=[alfa_1_0;
    alfa_1_0;
    alfa_2_0;
    alfa_2_0;
    F_y_1_0;
    F_y_1_0;
    F_y_2_0;
    F_y_2_0;
    F_x_1_0;
    F_x_1_0;
    F_x_2_0;
    F_x_2_0;
    v0;
    r0;
    u_0;
    Z_statico_1;
    Z_statico_1;
    Z_statico_2;
    Z_statico_2;
    slip_ratio_1_0;
    slip_ratio_1_0;
    slip_ratio_2_0;
    slip_ratio_2_0;
    omega_1_0;
    omega_1_0;
    omega_2_0;
    omega_2_0;
    schiacciamento_1_0;
    schiacciamento_1_0;
    schiacciamento_2_0;
    schiacciamento_2_0;
    R_puro_rotolamento_1_0;
    R_puro_rotolamento_1_0;
    R_puro_rotolamento_2_0;
    R_puro_rotolamento_2_0];

[t,y] = ode15s(@(t,y) longitudinale_puro_nomi_veri(t,y,J_1,J_2,...
    delta_v_0,tau_1,tau_2,a_1,a_2,t_1,t_2,m,J_z,J_zx,...
    dati_pacejka_r,dati_pacejka_f,R_1,R_2,g,h,Cz_1,Cz_2,...
    Z_statico_1,Z_statico_2,...
    theta_11,theta_12,theta_21,theta_22,l,tspan),tspan,y0,options);

Here's the function where all the equations are written:

function out=longitudinale_puro_nomi_veri(t,y,J_1,J_2,...
    delta_v_0,tau_1,tau_2,a_1,a_2,t_1,t_2,m,J_z,J_zx,...
    dati_pacejka_r,dati_pacejka_f,R_1,R_2,g,h,Cz_1,Cz_2,...
    Z_statico_1,Z_statico_2,...
    theta_11,theta_12,theta_21,theta_22,l,tspan)

delta_v=0;

momento_11=0; % Torque
momento_12=0;
momento_21=0;
momento_22=0;

delta_11=tau_1*delta_v;       % Wheel angles
delta_12=tau_1*delta_v;       
delta_21=tau_2*delta_v;
delta_22=tau_2*delta_v;

% Unknowns:
alfa_11=y(1);
alfa_12=y(2);
alfa_21=y(3);
alfa_22=y(4);
F_y_11=y(5);
F_y_12=y(6);
F_y_21=y(7);
F_y_22=y(8);
F_x_11=y(9);
F_x_12=y(10);
F_x_21=y(11);
F_x_22=y(12);
v=y(13);
r=y(14);
u=y(15);
F_z_11=y(16);
F_z_12=y(17);
F_z_21=y(18);
F_z_22=y(19);
epsilon_11=y(20);
epsilon_12=y(21);
epsilon_21=y(22);
epsilon_22=y(23);
omega_11=y(24);
omega_12=y(25);
omega_21=y(26);
omega_22=y(27);
schiacciamento_11=y(28); 
schiacciamento_12=y(29); 
schiacciamento_21=y(30);
schiacciamento_22=y(31);
R_pr_11=y(32);
R_pr_12=y(33);
R_pr_21=y(34);
R_pr_22=y(35);

out=[-alfa_11+delta_11-atan((v+r*a_1)/(u-r*t_1/2)); % Congruenza: deriva alfa_11
-alfa_12+delta_12-atan((v+r*a_1)/(u+r*t_1/2));      % Congruenza: deriva alfa_12
-alfa_21+delta_21-atan((v-r*a_2)/(u-r*t_2/2));      % Congruenza: deriva alfa_21
-alfa_22+delta_22-atan((v-r*a_2)/(u+r*t_2/2));      % Congruenza: deriva alfa_22
-F_y_11+pacejka('Fy0',F_z_11,0,-alfa_11*180/pi,-epsilon_11,dati_pacejka_f); % Legame costitutivo F_y_11
-F_y_12+pacejka('Fy0',F_z_12,0,-alfa_12*180/pi,-epsilon_12,dati_pacejka_f); % Legame costitutivo F_y_12
-F_y_21+pacejka('Fy0',F_z_21,0,-alfa_21*180/pi,-epsilon_21,dati_pacejka_r); % Legame costitutivo F_y_21
-F_y_22+pacejka('Fy0',F_z_22,0,-alfa_22*180/pi,-epsilon_22,dati_pacejka_r); % Legame costitutivo F_y_22
-F_x_11+pacejka('Fx0',F_z_11,0,-alfa_11*180/pi,-epsilon_11,dati_pacejka_f); % Legame costitutivo F_x_11 
-F_x_12+pacejka('Fx0',F_z_12,0,-alfa_12*180/pi,-epsilon_12,dati_pacejka_f); % Legame costitutivo F_x_12
-F_x_21+pacejka('Fx0',F_z_21,0,-alfa_21*180/pi,-epsilon_21,dati_pacejka_r); % Legame costitutivo F_x_21 
-F_x_22+pacejka('Fx0',F_z_22,0,-alfa_22*180/pi,-epsilon_22,dati_pacejka_r); % Legame costitutivo F_x_22 
(F_y_11+F_y_12+F_y_21+F_y_22)/m-u*r;                                     % Equilibrio trasversale
((F_y_11+F_y_12)*a_1-(F_y_21+F_y_22)*a_2)/J_z;                           % Equilibrio imbardata 
(F_x_11+F_x_12+F_x_21+F_x_22)/m-(F_y_11*delta_11+F_y_12*delta_12)/m+v*r; % Equilibrio longitudinale 
-F_z_11+Z_statico_1-theta_11*(F_y_11+F_y_12+F_y_21+F_y_22)-theta_12*((F_y_11+F_y_12)*a_1-(F_y_21+F_y_22)*a_2)-(F_x_11+F_x_12+F_x_21+F_x_22-(F_y_11*delta_11+F_y_12*delta_12))*h/l;  
-F_z_12+Z_statico_1+theta_11*(F_y_11+F_y_12+F_y_21+F_y_22)+theta_12*((F_y_11+F_y_12)*a_1-(F_y_21+F_y_22)*a_2)-(F_x_11+F_x_12+F_x_21+F_x_22-(F_y_11*delta_11+F_y_12*delta_12))*h/l;                        
-F_z_21+Z_statico_2-theta_21*(F_y_11+F_y_12+F_y_21+F_y_22)-theta_22*((F_y_11+F_y_12)*a_1-(F_y_21+F_y_22)*a_2)+(F_x_11+F_x_12+F_x_21+F_x_22-(F_y_11*delta_11+F_y_12*delta_12))*h/l;                       
-F_z_22+Z_statico_2+theta_21*(F_y_11+F_y_12+F_y_21+F_y_22)+theta_22*((F_y_11+F_y_12)*a_1-(F_y_21+F_y_22)*a_2)+(F_x_11+F_x_12+F_x_21+F_x_22-(F_y_11*delta_11+F_y_12*delta_12))*h/l;                       
-epsilon_11+scorrimento_11(u,R_pr_11,omega_11);  % Congruenza: slip x_11    
-epsilon_12+scorrimento_12(u,R_pr_12,omega_12);  % Congruenza: slip x_12    
-epsilon_21+scorrimento_21(u,R_pr_21,omega_21);  % Congruenza: slip x_21
-epsilon_22+scorrimento_22(u,R_pr_22,omega_22);  % Congruenza: slip x_22
(momento_11-F_x_11*(R_1-schiacciamento_11))/J_1;  % Bilancio ruota 11
(momento_12-F_x_12*(R_1-schiacciamento_12))/J_1;  % Bilancio ruota 12
(momento_21-F_x_21*(R_2-schiacciamento_21))/J_2;  % Bilancio ruota 21
(momento_22-F_x_22*(R_2-schiacciamento_22))/J_2;  % Bilancio ruota 22
-schiacciamento_11+F_z_11/Cz_1;   % Schiacciamento 11
-schiacciamento_12+F_z_12/Cz_1;   % Schiacciamento 12
-schiacciamento_21+F_z_21/Cz_2;   % Schiacciamento 21
-schiacciamento_22+F_z_22/Cz_2;   % Schiacciamento 22
-R_pr_11+R_1-schiacciamento_11/3;  % Raggio puro rotolamento 11
-R_pr_12+R_1-schiacciamento_12/3;  % Raggio puro rotolamento 12
-R_pr_21+R_2-schiacciamento_21/3;  % Raggio puro rotolamento 21
-R_pr_22+R_2-schiacciamento_22/3]; % Raggio puro rotolamento 22

Answer 1

[Francesco Ardia]

SOLVED! Equations and initial conditions were correct from the beginning. Assigning null torque on all 4 wheels caused the ODE to generate a singular matrix, and thus the failure. Any different value, even for just one wheel, works fine!

So when solving DAEs be careful to the forcing terms!

Answer 2

[Walter Roberson]

You define matrix M but do not show it being used. You show m (lower case) being used.

If we assume that M and m are the same then look inside the routine to see

(F_y_11+F_y_12+F_y_21+F_y_22)/m-u*r

That is a matrix right division by m, similar to matrix multiplication by inv(m). That would not work if m is singular.

If you are intentionally using matrix right division then it would be clearer if you were to recode using \ instead of /

If you are not intentionally using matrix right division then ./ might be a more appropriate operator.

[Francesco Ardia]

options = odeset('Mass',M,'RelTol',1e-8,'AbsTol',tolleranze);

I defined it in the options, but hadn't posted it. m is the car's mass and is a scalar

[Walter Roberson]

Try setting MassSingular 'yes'

[Francesco Ardia]

it doesn't work, unfortunately... this is driving me mad!

[Walter Roberson]

Sorry, I do not know.

[Francesco Ardia]

I managed to solve it with a little trick. Thank you very much for your interest and time! :)

[ahmad eldeeb]

could you tell?

[Francesco Ardia]

If I put 0.000001 as one of the single torques instead of 0, the code runs

Answer 3

[Torsten]

Start with consistent initial conditions, i.e. the algebraic equations should be satisfied at t=0. If necessary, use "fsolve" to get such a consistent set of initial values.

Best wishes

Torsten.

[Francesco Ardia]

Hi Torsten: I derived the initial conditions from physical considerations, and applying fsolve to the algebraic equations only seems to confirm that they are correct, as I get this message:

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

Elapsed time is 0.048910 seconds.

What else could it be?

[Francesco Ardia]

The warning message is actually as follows:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
> In ode15s (line 589) 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

So perhaps there's something to sort within the ODE solver, but I can't interpret it...

[Torsten]

Output the vector "out" you obtain in "longitudinale_puro_nomi_veri".

Are there NaN values involved ?

Best wishes

Torsten.