Using ode4 solver in a sliding mode controller

Hello,

I did a xls model in simulink and a matlab script for a SMC controller and I was expecting the same resutl for both, but is not exactly the same. After looking for both programs, I realized that I'm using a ode4 solver in the simulink xls and seems for me that this is the main issue (I saved the time vector from simuink and use it in the .mat, the results where the same)

I'm quite lost with ode's in matlab, and more if I have to solve this control algorithm... Please How can I use ode4 in my matlab script? I have searched in the forum and in matlab help and I have found this link, but I dont get the point on how to use the function:

https://ch.mathworks.com/matlabcentral/fileexchange/59044-ode4-gives-more-accurate-results-than-ode45-ode23-ode23s

I see that ode uses this fucntion to work where I have to use my starting time, time step, final time and initial conditions. I'm lostwith the F function part. I dont know how can I integrate this in my code.

yout = ODE4(F,t0,h,tfinal,y0)

attached bot scripts, Thanks in advance!

[Sam Chak]

Hi @Mikel

Correct me if I'm wrong. In SMDOB_X2_observer.m, the disturbance is disabled , and the input is changed to so that I can compare with the standard response of a stable mass-spring-damper-like system:

The solution for . But the response from your code looks different from the analytical solution.

y = @(t) exp(-t).*(1 + t);
fplot(y, [0 10]), grid on
xlabel('t'), ylabel('x_{1}')

When I run your integration code in SMDOB_X2_observer.m, the supposedly stable mass-spring-damper system becomes unstable.

[Mikel]

Hi,

I dont understand your answer quite well, sorry. The perturbance is always active until the U control law takes control of it and neutralizes the disturbance so the state variable estimations x1s and x2s converge into the system real varaibles x1 and x2. My problem is that as in simulink using the ode4 solver the results are correct, in my matlab file I integrate the space-state equations here:

%calculate x1 and x2
x1(ind1+1,1)= deltaT*x2ant+x1ant;
x2(ind1+1,1)=(f(ind1,1)+u(ind1,1))*deltaT+x2ant;

%estimates
x1s(ind1+1,1)=x1s(ind1,1)+deltaT*v(ind1,1);
x2s(ind1+1,1)=x2s(ind1,1)+deltaT*veqs(ind1,1);

All I want to do is to implement the ode4 solver in my .mat so I can see that this was exactly the problem as I guess, but I dont know how to use ode4 in the for loop. It should be something like that but everything I do is wrong. Maybe I'm completely wrong, I hope you can understand what I'm trying to explain here.

% as the space state system is defines by: 
%   x1'=x2                  x1(0)=x10
%   x2'= u + f(x1,x2,t)     x2(0)=x20   f in this case known: f=sin(2*t)

%where x1 = x and x1'= x2

%I have a mess here, I dont know how can I define my system using these
%@functions and the somehow use it into the loop. 
F1 = @(t,x1) x1;
F2 = @(t,x2) f+u;

x1(ind1+1,1) = ode4(F1,t(ind1,1),deltaT,t(ind1+1,1),x1(ind1,1));
x2(ind1+1,1) = ode4(F2,t(ind1,1),deltaT,t(ind1+1,1),x2(ind1,1));

[Sam Chak]

@Mikel, Thanks for your reply.

Have you tested the system using the ode45() solver?

I understand that the sinusoidal disturbance is always active. I purposely set to simulate the response of a Double Integrator under the influence of a full-state feedback input .

This allows me evaluate whether integration code in the m-file works as expected or not.

Using your integration code, the states appears to grow and oscillates after 40 seconds (significant to human eyes), and they explode after 60 seconds. However, this test system should be stable.

Answer 1

[Sam Chak]

Hi @Mikel

In MATLAB, I think these are the expexted results ,where the nonlinear state estimators and are able to rapidly track the true states and , respectively, within 0.2 sec.

I didn't change your design code for the high-gain sliding mode observer part, but I used the ode45() solver instead of your integration code.

% solving the ode
tspan  = linspace(0, 2, 20001);
x0     = [1 -2 0 0];
[t, x] = ode45(@smc_obsv, tspan, x0);

% plotting the solution
plot(t, x), grid on,
legend({'$x_{1}$', '$x_{2}$', '$\hat{x}_{1}$', '$\hat{x}_{2}$'}, 'Interpreter', 'latex', 'fontsize', 20)
xlabel({'$t$'}, 'Interpreter', 'latex')
ylabel({'$\mathbf{x}$'}, 'Interpreter', 'latex')

% describing the ode
function xdot = smc_obsv(t, x)
    xdot    = zeros(4, 1);	
    
    % parameters
    f       = sin(2*t);                     % disturbance
    
    c       = 1.5;                          %
    r       = 2;                            % criterion: r > max(f)
    s       = c*x(1) + x(4);                % sliding surface
    
    u       = - r*tanh(s/0.001) - c*x(4);   % control input
    
    r1      = 10;                           % high gain
    z       = x(3) - x(1);                  % estimation error
    v1      = - r1*tanh(z/0.001);           % rate of xhat_1
    LPF     = 1/0.01;                       % adapatation rate
    v2      = (v1 - x(4))*LPF;              % rate of xhat_2
    
    xdot(1) = x(2);
    xdot(2) = f + u;
    xdot(3) = v1;
    xdot(4) = v2;
end

[Mikel]

Thanks for the answer, this is what I want to do using ode4 (as I need to be a fixed step solver). I just have one question as I'm learning SMC and my understanding of concepts are not enough good. Here the theory book I'm following (yuri shtessel's sliding mode control and observation) says that in chapter 1.6 that the injection term is chosen as rho * sign(z)

 v1      = - r1*tanh(z/0.001);
 

why are you changing that to a tanh function? I mean, is the same but do you have any particular reason for that?

[Sam Chak]

You are welcome, @Mikel.

The signum, is a non-smooth discontinuous function with discontinuity at , and ode45() is not built to handle stiff equations. I usually approximate the discontinuous with a hyperbolic tangent function by setting because it is continuous. Moreover, there are no physical systems that behave like a pure . Most behave like saturation functions, with a very steep slope . Thus, it is very reasonable use the alternative continuous function .

Another continuous candidate is a sigmoidal function , though it has a smooth discontinuous derivative.

Test:

In the following simulation, though it appears that the doesn't cause any issue, the solver actually becomes inefficient when simulating for a longer time ( sec). It happens because the sliding variable .

% solving the ode
tspan  = linspace(0, 2, 20001);
x0     = [1 -2 0 0];
[t, x] = ode45(@smc_obsv, tspan, x0);

% plotting the solution
plot(t, x), grid on,
legend({'$x_{1}$', '$x_{2}$', '$\hat{x}_{1}$', '$\hat{x}_{2}$'}, 'Interpreter', 'latex', 'fontsize', 20)
xlabel({'$t$'}, 'Interpreter', 'latex')
ylabel({'$\mathbf{x}$'}, 'Interpreter', 'latex')

% describing the ode
function xdot = smc_obsv(t, x)
    xdot    = zeros(4, 1);	
    
    % parameters
    f       = sin(2*t);                     % disturbance
    
    c       = 1.5;                          %
    r       = 2;                            % criterion: r > max(f)
    s       = c*x(1) + x(4);                % sliding surface
    
    u       = - r*sign(s) - c*x(4);         % control input
    
    r1      = 10;                           % high gain
    z       = x(3) - x(1);                  % estimation error
    v1      = - r1*sign(z);                 % rate of xhat_1
    LPF     = 1/0.01;                       % adapatation rate
    v2      = (v1 - x(4))*LPF;              % rate of xhat_2
    
    xdot(1) = x(2);
    xdot(2) = f + u;
    xdot(3) = v1;
    xdot(4) = v2;
end

[Mikel]

thank you for the information! Yeah, my guess wass as well to use the tanh function. No I'm having some issues using your answer with the ode4 function as my system is constrained to use a fixed step of 1e-4.

[Sam Chak]

If you are using the signum function, in the code, the system will become very stiff as .