Why am I recieving the following error "Warning: Failure at t=1.705994e+00. Unable to meet integration tolerances without reducing the step size below the smallest value all"

I am attempting to simulate a system of equivalent first order differential equations for three second order differential equations where I have one input function "u". I suspect that my system is becoming discontinuous as my controller is of the sliding mode category however, I would like to confirm my suspicions here and see if I could obtain some advice and help.

% Hybrid Controller
u = @(y) tanh(hybrid(y,I1,I2,Itheta1,Itheta2,L1,L2,m1,m2,l1,l2,c1,c2,M1,M2,M3,g,K1,ua1,ua2,y_com));
[t, state_values] = ode45(@(t,y)dpic(t,y,g,M1,M2,M3,Itheta1,Itheta2,c1,c2,u(y)), tspan, IC);

There is more code that represents this nonlinear model's function however I have ommited it for brevity. This is the very last line of the function which represents the dynamics of my system and also where the control "u" is applied.

Y_dot = A*Y + B*u + L; % First order nonlinear system matrix model

Below here is where I suspect my control issues are coming from. I have attempted to change the "sign" functions to "tanh" functions though not much has changed and the integration still failed. Any advice? Thank you for your time and assistance.

%% ========= Hybrid Controller  =========
function utot = hybrid(Y,I1,I2,Itheta1,Itheta2,L1,L2,m1,m2,l1,l2,c1,c2,M1,M2,M3,g,K1,ua1,ua2,Y_com)
theta1 = Y(2);
theta2 = Y(3);
% phi    = 0.3; % Boundary layer thickness
% (larger phi = less chattering, larger steady-state error) or 
% (smaller phi = more chattering, smaller steady-state error)
if cos(theta1) < 0.8
    utot = step1(Y,I1,m1,l1,ua1,K1,g);
elseif (cos(theta1) >= 0.8) && (cos(theta2) <= 0.94)
    utot = step2(Y,Itheta1,Itheta2,L1,I2,m2,l2,M1,M2,M3,c1,c2,ua2,g,Y_com);
elseif (cos(theta1) >= 0.8) && (cos(theta2) >= 0.94)
    utot = step3(Y,Itheta1,Itheta2,M1,M2,M3,c1,c2,g,Y_com);
end
end
%% ========= Step1 Controller  =========
function u1 = step1(Y,I1,m1,l1,ua1,K1,g)
% Swing up lower pendulum
% Actual system states
x          = Y(1);
theta1     = Y(2);
x_dot      = Y(4);
theta1_dot = Y(5);
phi        = 0.0001;
x1 = [x x_dot]';                                         % Cart data column vector
uz = -K1*x1;                                             % Control input to bound cart to center of rail
E1 = (1/2)*I1*(theta1_dot)^2 + m1*g*l1*(cos(theta1)-1);  % Energy of lower pendulum
utheta = ua1*sign(E1*theta1_dot*cos(theta1));            % Control input to swing up lower pend utheta=u

if (sign(utheta) == sign(uz)) || (abs(utheta) > abs(uz)) % sign(control input to swing pend) == sign(control in
    u1 = ua1*sign(E1*theta1_dot*cos(theta1)) + uz;
else, (sign(utheta) ~= sign(uz)) && (abs(utheta) <= abs(uz));
    u1 = ua1*sign(E1*theta1_dot*cos(theta1));
end
end

[Torsten]

It seems there are discontinuities in your system of differential equations. All numerical ODE integrators dislike them and often cannot continue integrating once such a discontinuity is reached.

Maybe this happens at t = 1.705, but it's only a guess.

[Akekaphop Kesavadhana]

How could I get around that? Are there any rememdies? I have looked into events for ode45, but I don't understand how I would implement it with such a complex controller. Or maybe do I need not to use if statements for my controller but use such a closed form equivalent of said controller?

[Steven Lord]

When you try to solve your system of ODEs up to t = 1.7 and plot the solution, what does it look like? Do any of the components look like they're growing to an extremely large magnitude (heading towards either +Inf or -Inf?)

Without seeing the rest of your code I'm not sure how much help we will be able to give you beyond just general advice.

[Akekaphop Kesavadhana]

This is the graph that I am observing when the two pendula start from 90 degrees counter-clockwise from the vertical 0 degrees. It does indeed look like the velocities are going to inf.

\

[Akekaphop Kesavadhana]

In the above simulation I also set the control input to 0 and I still have discontinuities.

[Sam Chak]

Hi @Akekaphop Kesavadhana,

From the plot below, it is clear that the Step1 Controller is a discontinuous function with a jump discontinuity at , which is caused by the signum function in utheta. The shape of the Step1 Controller is influenced by the first condition of the If–Else conditional statement (saturation + straight line).

Since the Step2 and Step3 controllers are not provided, only the Step1 controller can be evaluated. While the hyperbolic tangent function is continuous, applying tanh(·) to the Step1 Controller output signal does not render it continuous. Furthermore, this mathematical operation also clamps the signal to the range .

X   = linspace(-pi/2, pi/2, 36001);
Y   = ones(5, numel(X)).*X;
Y   = Y';
I1  = 1;
m1  = 1;
l1  = 1;
ua1 = 1;
K1  = 1;
g   = 9.81;

u1  = zeros(1, numel(X));
for j = 1:numel(X)
    u1(j)   = step1(Y(j,:).', I1, m1, l1, ua1, K1, g);
end
u   = tanh(u1);     % hybrid() function lacks Step2 and Step3

tl  = tiledlayout(3, 3, 'TileSpacing', 'Compact');

nexttile([2, 1])
plot(u1, u), grid on, ylim("padded") 
xlabel('u1 Amplitude'), ylabel('u Amplitude')
title('tanh effect')

nexttile([2, 2])
plot(rad2deg(X), u), grid on, ylim("padded") 
xticks(-90:30:90)
xlabel('\theta / degree'), ylabel('u Amplitude')
title('Hyperbolically-tangentized Controller, u')

nexttile(8, [1, 2])
plot(rad2deg(X), u1), grid on, ylim("padded") 
xticks(-90:30:90)
xlabel('\theta / degree'), ylabel('u1 Amplitude')
title('Step1 Controller, u1')

%% ========= Step1 Controller  =========
function u1 = step1(Y, I1, m1, l1, ua1, K1, g)
    % Swing up lower pendulum
    % Actual system states
    x           = Y(1);
    theta1      = Y(2);
    x_dot       = Y(4);
    theta1_dot  = Y(5);
    phi         = 0.0001;
    x1          = [x x_dot]';                                           % Cart data column vector
    uz          = -K1*x1;                                               % Control input to bound cart to center of rail
    E1          = (1/2)*I1*(theta1_dot)^2 + m1*g*l1*(cos(theta1) - 1);  % Energy of lower pendulum
    utheta      = ua1*sign(E1*theta1_dot*cos(theta1));                  % Control input to swing up lower pend utheta=u
    
    if   (sign(utheta) == sign(uz(1))) || (abs(utheta) >  abs(uz(1)))   % sign(control input to swing pend) == sign(control in
        u1 = utheta + uz(1);
    else (sign(utheta) ~= sign(uz(1))) && (abs(utheta) <= abs(uz(1)));
        u1 = utheta;
    end
end

[Sam Chak]

Hi @Akekaphop Kesavadhana

In the absence of a control signal, as time t approaches 1 s, the upper pendulum angle rapidly converges to approximately 1.5708 rad (90°). It seems intuitive that cos(90°) = 0, and the (minus sign) division by zero is likely the cause of the angular velocity dropping sharply and heading toward negative infinity.

Additionally, due to the tanh clamping of the control signal in your original simulation, the control torque may be insufficient to regulate the double-pendulum system, leading to instability and termination of simulation at around s.

Answer 1

[Walter Roberson]

When you call a function from ode45(), it is strictly necessary that the called function is continuous to its second derivative.

If you violate this condition and you are lucky, you will get an error message about failure to reach integration target.

If you violate this condition and you are not lucky, then ode45() will silently give you wrong answers.

if cos(theta1) < 0.8
    utot = step1(Y,I1,m1,l1,ua1,K1,g);
elseif (cos(theta1) >= 0.8) && (cos(theta2) <= 0.94)

Those look unlikely to be continuous to the second derivative.

If your discontinuities happen to be predictable based strictly on time, then loop your ode45() calls using tspans that have upper and lower bounds that match the time discontinuities.

If your discontinuities are not easily predictable based on time, then you must use event functions to detect the discontinuities and signal termination. Loop continuing the tspan from where integration left off.