a nonlinear system and its control input in simulink

Sign in to answer this question.

Follow Question

Answer 1

Sam Chak on 14 Apr 2024

1 vote

The control signal plotted in your comment is incorrect. Since you have completed the learning tasks and are capable of solving the DDE using the dde23 solver in MATLAB, I believe it's time to proceed with running the simulation in Simulink, as this was the original problem you asked about. Visualizing the plots of the system response

and the control signal

in Simulink is much easier.

Below is the block diagram in Simulink, which is self-explanatory. The plots of

and

seem to resemble the results depicted by the blue dashed lines in Figures 3 and 4 of Ding et al., 2023. If you find the Simulink solution provided in this response to be satisfactory, you can "unaccept" my previous answer and consider "Accepting" ✔ this answer and giving it a thumbs-up vote 👍.

Best of luck with working on the second and final example from the paper. If you encounter any technical issues while solving them, feel free to post a new question to prevent this thread from becoming too lengthy, as there are already over 50 comments.

4 Comments
Show 2 older comments Hide 2 older comments

controlEE on 18 Apr 2024

Hi @Sam Chak , Thank you very much , I ran the model in simulink and due to replacement of heaviside() with sign() , It performs correctly .

but I wonder what the signal in blue is for ?

and also based on the tip you mentioned here I modified the plot part and ;

sympref('HeavisideAtOrigin', 1);    
sol = dde23(@ddefun, 2, @history, [0, 6]);
%plot(sol.x, sol.y), grid on, xlabel t, ylabel x(t)
subplot(2,1,1), plot(tout, x), grid on,xlabel t , ylabel x(t) 
subplot(2,1,2), plot(tout, u), grid on,xlabel t , ylabel u(t) 
function [dx, u] = ddefun(t, x, z)
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    Rh  = ((sin(pi*t/h)).^2).*heaviside(t - h) - ((sin(pi*t/h)).^2).*heaviside(t - 2*h);
    Kah = Rh.*(1/1.8269).*exp(-a*(h - 2*t));
    sig = sign(x).*(abs(x)).^(2*tau - 1);
    u   = (((-a/(2*(1 - tau)) - 0.1)*x - (Kah/(2*(1 - tau))).*sig.*(abs(z).^(2*(1 - tau)))) - x.^3)./(1 + x.^2);
    d   = 0.1*exp(-t)*x;    % 0.1*x
    dx  = x.^3 + (1 + x.^2).*u + d;
end
function s = history(t) 
    s   = 1;
end 

all clear and correct both in matlab and simulink .

THANKS

Sam Chak on 18 Apr 2024

Congratulations on your success! The blue curve displayed in the Scope represents the solution response, which is delayed by 2 seconds,

.

If you find the Simulink solution provided in this response to be satisfactory, you have the option to "unaccept" my previous answer and consider "Accepting" ✔ this answer instead.

Answer 2

Sam Chak on 16 Mar 2024

1 vote

While I cannot evaluate the MATLAB code from the image, it seems that

in Equation (6) represents a definite integral. The integrand is integrated from h to

, and the resulting function

is plotted below. It exhibits a distinct bell-shaped curve, which is commonly utilized in Prescribed-Time Control algorithms.

Did you obtain the same result?

syms t

a = 0.1;

h = 2;

tau = 5/7;

Rh = (sin(pi*t/h))^2;

expr= exp(2*a*t)*Rh

expr =

Wc = int(expr, [h, 2*h])

Wc =

Wc = double(Wc)

Wc = 1.8269

W = inv(Wc);

Kah = Rh*W*exp(-a*(h - 2*t));

fplot(Kah, [h, 2*h]), grid on, xlabel('t')

title({'$K_{a, h}(t)$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

28 Comments
Show 26 older comments Hide 26 older comments

Sam Chak on 17 Mar 2024

Thank you for the update. I have plotted the function you mentioned.

If

is indeed intended to represent the mathematician's signum function,

, then when the argument x is a negative real number, it would result in complex numbers due to the fractional exponent in the expression. It's important to note that

represents the nth root of

. Therefore, I believe the authors did not intend

to be mathematically interpreted as

. It should be defined somewhere in the paper, perhaps in a numbered or unnumbered equation.

It is worth noting that in the field of control, it is not uncommon for some practitioners with a practical focus (or those lacking an extensive mathematical background) to uncritically adopt the approaches of others.

x = linspace(-1, 1, 20001);

tau = 5/7;

expn= 2*tau - 1 % exponent is a fraction 2143/5000

expn = 0.4286

y = (sign(x)).^expn;

plot(x, y, 'linewidth', 2), grid on, ylim([0, 1.2])

Warning: Imaginary parts of complex X and/or Y arguments ignored.

xlabel('x'), ylabel('y')

title({'$\left[\mathrm{sign}(x)\right]^{2 \tau - 1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

Sam Chak on 20 Mar 2024

I would suggest simulating the nonlinear system in MATLAB initially, before implementing the design model in Simulink. This approach makes it relatively easier to run the simulation on the MATLAB online platform and identify any errors that may arise from your code. You can click on the indentation icon

to copy/paste the code into the grey field.

Below is an example written as a 4-line script, with a 'function' named 'nonlinearSystem()' placed at the end of the script (from Line 7 to Line 18).

,

where the control input u is given by

tspan = [0, 5]; % simulation time % Line 1

x0 = 1; % initial value % Line 2

[t, x] = ode45(@nonlinearSystem, tspan, x0); % Line 3

plot(t, x), grid on, xlabel t, ylabel x(t) % Line 4

% Line 5

%% Function that describes the Nonlinear System % Line 6

function dx = nonlinearSystem(t, x) % Line 7

% parameters % Line 8

a = 2; % Line 9

b = 3; % Line 10

c = 5; % Line 11

% Line 12

% control input % Line 13

u = - (c*x)/b; % Line 14

% Line 15

% differential equation % Line 16

dx = a*sin(x) + b*u; % Line 17

end % Line 18

Sam Chak on 21 Mar 2024

It seems there is some confusion. My suggestion is to start by implementing the algorithm in MATLAB first, as shown in the second sample code below. If you achieve success with the MATLAB implementation, you can then copy the same code (from Line 7 to Line 16) into your Simulink model, and it should work accordingly.

However, if you encounter any issues or errors, feel free to investigate and share the error message here for further assistance. Take a moment to contemplate and determine what should be written in Line 9, which corresponds to the control input u.

tspan   = [0, 5];   % simulation time           % Line 1
x0      =  1;       % initial value             % Line 2
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   % Line 3
plot(t, x), grid on, xlabel t, ylabel x(t)      % Line 4
                                                % Line 5
%% Function that describes the Nonlinear System % Line 6
function dx = nonlinearSystem(t, x)             % Line 7
    % control input                             % Line 8
    u   = (...);                                % Line 9
                                                % Line 10
    % disturbance                               % Line 11
    d   = 0.1*x;                                % Line 12
                                                % Line 13
    % differential equation                     % Line 14
    dx  = x^3 + (1 + x^2)*u + d;                % Line 15
end                                             % Line 16

controlEE on 22 Mar 2024

Edited: Torsten on 22 Mar 2024

Hello Mr.Chak

first thank for your help

I define the control input as follows , first I determine the constants that are used in it ;

tspan   = [0, 5];   % simulation time           
x0      =  1;       % initial value             
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   
plot(t, x), grid on, xlabel t, ylabel x(t)      
%% Function that describes the Nonlinear System 
function dx = nonlinearSystem(t, x)             
% control input                            
syms t
a   = 0.1;
h   = 2;
tau = 5/7;
Rh  = (sin(pi*t/h))^2;
expr= exp(2*a*t)*Rh
Wc  = int(expr, [h, 2*h])
Wc  = double(Wc)
W   = inv(Wc);
Kah = Rh*W*exp(-a*(h - 2*t));
expn= 2*tau -1 ;
y   = (sign(x)).^expn;
u   = 1/(1 + x^2)*((-a/2*(1-tau)-0.1))*x(t)+(-Kah/2(1-tau))*y*(abs(x(t-h)))^2(1-tau)-x^3)                               
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.
% disturbance                               
d   = 0.1*x;                               
% differential equation                    
dx  = x^3 + (1 + x^2)*u + d;                
end                

but it doesn't show me any results on x .

controlEE on 22 Mar 2024

tspan   = [0, 5];   % simulation time           
x0      =  1;       % initial value             
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   
plot(t, x), grid on, xlabel t, ylabel x(t)      
                                               
%% Function that describes the Nonlinear System 
function dx = nonlinearSystem(t, x)             
    % control input                            
    syms t
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    Rh  = (sin(pi*t/h))^2;
    expr= exp(2*a*t)*Rh
    Wc  = int(expr, [h, 2*h])
    Wc  = double(Wc)
    W   = inv(Wc);
    Kah = Rh*W*exp(-a*(h - 2*t));
    expn= 2*tau -1 ;
    y   = (sign(x)).^expn;
    u   = 1/(1 + x^2)*((-a/2*(1-tau)-0.1))*x(t)+(-Kah/2(1-tau))*y*(abs(x(t-h)))^2(1-tau)-x^3)                                % Line 9
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.
                                            
    % disturbance                               
    d   = 0.1*x;                               
                                               
    % differential equation                    
    dx  = x^3 + (1 + x^2)*u + d;                
end                                   

Sam Chak on 22 Mar 2024

Edited: Sam Chak on 23 Mar 2024

OP: for instance we are working with the term x and we have a term in control x(t) , are they both same because we want to write the abs(x(t-h))

Regarding

and

, they represent the same state of the nonlinear system. However,

specifically refers to the state that is delayed by 2 seconds. It's important to note that you cannot directly use the mathematical notation 'x(t - 2)' in the code because MATLAB will interpret it as accessing the element of the vector x based on its index location (at t - 2) in the vector.

To handle the time-delayed state

, you need to implement a special assignment. For now, you can assume temporarily that

is equal to

without any delay.

% – - – - – - – - – - – - – - – - – -

Update: I have included a colored graphic to assist you in coding the four main parts of the controller individually, using just four lines of code. Once you have completed these four parts, you can combine them in the fifth line to create the complete controller.

% – - – - – - – - – - – - – - – - – -

I have managed to get the code running using a simple Proportional Control scheme, without incorporating the Prescribed Time Control law. To implement the 'non-time-delayed' Prescribed Time Control law, please modify the code by replacing 'uu' with the corresponding expression from your Prescribed Time Control law.

tspan = [0, 5]; % simulation time

x0 = 1; % initial value

[t, x] = ode45(@nonlinearSystem, tspan, x0);

plot(t, x), grid on, xlabel t, ylabel x(t)

%% Function that describes the Nonlinear System

function dx = nonlinearSystem(t, x)

% control input

% syms t % unnecessary if directly assigning Wc

a = 0.1;

h = 2;

tau = 5/7;

Rh = (sin(pi*t/h))^2;

expr= exp(2*a*t)*Rh;

% Wc = int(expr, [h, 2*h])

Wc = 1.8269; % direct assignment but Kah is not used

W = inv(Wc);

Kah = Rh*W*exp(-a*(h - 2*t)); % not used in this run

expn= 2*tau -1 ;

y = (sign(x)).^expn; % not used in this run

k = 3; % larger value makes convergence faster

uu = - k*x; % proportional control (for testing purposes)

u = (uu - x^3)/(1 + x^2); % cancel out the destablizing nonlinearity x^3

% state-dependent disturbance

d = 0.1*x;

% differential equation

dx = x^3 + (1 + x^2)*u + d;

end

controlEE on 24 Mar 2024

Moved: Sam Chak on 24 Mar 2024

Hi Mr.Chak , Thank u very much for the tips you mentioned here , I did some research upon dde and ode and ode45 too , I really appreciate that .

ok , alright the uu should be replaced by :

and I rewritted down in the code :

There are several things I wanted your help with , First u said we can assum x(t-2) as x(t) , " h =2 " , I mean as a non-delayed term , the secoond thing is that u itself has dynamics and we would have a curve for it which depeds on the time as I Sent the results of the example before , I mean it wouldn't have a certain value or being as a constant function , The third is when I'm trying to run the code to see the results actullay a window came up which is called profiler , It's about a few days that It doesn't work for me , and the fourth is that there is an error with the code which keeps coming up about parentheses , I do consider them but I don't know what the problem accurately is .

tspan   = [0, 5];   % simulation time           
x0      =  1;       % initial value             
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   
plot(t, x), grid on, xlabel t, ylabel x(t)      
                                               
%% Function that describes the Nonlinear System 
function dx = nonlinearSystem(t, x)             
    % control input                            
    % syms t                        % unnecessary if directly assigning Wc
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    Rh  = (sin(pi*t/h))^2;
    expr= exp(2*a*t)*Rh;
    % Wc  = int(expr, [h, 2*h])
    Wc  = 1.8269;                   % direct assignment but Kah is not used
    W   = inv(Wc);
    Kah = Rh*W*exp(-a*(h - 2*t));   % not used in this run
    expn= 2*tau -1 ;
    y   = (sign(x)).^expn;          % not used in this run
    
    k   = 3;                        % larger value makes convergence faster
   % uu  = - k*x;                   % proportional control (for testing purposes)
    uu =  (-a/2(1 - tau) - 0.1)*x - kah/2(1-tau)*y*abs(x).^2(1-tau) ; 
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.
    u   = (uu - x^3)/(1 + x^2);     % cancel out the destablizing nonlinearity x^3
                                            
    % state-dependent disturbance                          
    d   = 0.1*x;                               
                                               
    % differential equation                    
    dx  = x^3 + (1 + x^2)*u + d;                
end                                

controlEE on 25 Mar 2024

Edited: controlEE on 25 Mar 2024

Hi Mr.Chak , How are you ? I hope you feel fresh and good .

I troubleshooted the code based on your refers as follows :

tspan   = [0, 5];   % simulation time           
x0      =  1;       % initial value             
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   
plot(t, x), grid on, xlabel t, ylabel x(t)      
                                               
%% Function that describes the Nonlinear System 
function dx = nonlinearSystem(t, x)             
    % control input                            
    % syms t                        % unnecessary if directly assigning Wc
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    Rh  = (sin(pi*t/h))^2;
    expr= exp(2*a*t)*Rh;
    % Wc  = int(expr, [h, 2*h])
    Wc  = 1.8269;                   % direct assignment but Kah is not used
    W   = inv(Wc);
    Kah = Rh*W*exp(-a*(h - 2*t));   % not used in this run
    expn= 2*tau -1 ;
    y   = (sign(x)).^expn;          % not used in this run
    
    k   = 3;                        % larger value makes convergence faster
    % uu  = - k*x;                   % proportional control (for testing purposes)
    uu = (-a/2*(1 - tau) - 0.1)*x - Kah/2*(1-tau)*y*abs(x).^2*(1-tau);
    u   = (uu - x^3)/(1 + x^2);     % cancel out the destablizing nonlinearity x^3
                                            
    % state-dependent disturbance                          
    d   = 0.1*x;                               
    % d = 0.1 * exp(-t) * x;                                           
    % differential equation                    
    dx  = x^3 + (1 + x^2)*u + d;                
end                                  

and the state ;

and for the d = 0.1e^(-t) * x :

based on the article , states should converge to zero (our goal) before t = 5 sec , I Increased the length of time , but states didn't converge to zero yet , I'm not sure but mabe becuse of delay term ;

for the tspan = [ 0 , 30 ]

for d = 0.1e^(-t) * x ;

by the way , I studied on power electronics for my Bachelor's degree , and I did some research for fuzzzy control , but kind of I'm interested in nonlinear control too .

thank you

Sam Chak on 25 Mar 2024

Edited: Sam Chak on 26 Mar 2024

By the way, what is the title of the paper? Have you attempted to search for the paper to check if the author who wrote the code has mathematically defined the "sig" function, such as

?

Previously, you made a mistake by assuming that 'a/b*c' is equal to

. However, I have fixed that issue in the code. Nevertheless, the response is still unstable due to the incorrect computation of the 'sig' function. Once you fix that, the response should become stable, as demonstrated by the coding author in the paper.

Transitioning from linear control to nonlinear control can be quite challenging even for those studying for a PhD, if you are unfamiliar with various aspects such as:

Mathematical notations (e.g., piecewise function in Rh and the 'non-recognized' sig function)
Math operator precedence (+, −, ×, ÷, ^, parentheses),
Elementary math functions (, , , , , , , , etc.)
Differential equations (ordinary vs. delay),
Divide and Conquer strategy for coding.

Mastering these concepts is crucial for making the prescribed time control algorithm works.

tspan = [0, 5]; % simulation time

x0 = 1; % initial value

[t, x] = ode45(@nonlinearSystem, tspan, x0);

plot(t, x), grid on, xlabel t, ylabel x(t)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

%% Function that describes the Nonlinear System

function dx = nonlinearSystem(t, x)

% parameters (minor divide-and-conquer)

a = 0.1;

h = 2;

tau = 5/7;

Rh = (sin(pi*t/h))^2; % Minor issue: incorrect math function

Wc = 1.8269;

W = inv(Wc);

Kah = Rh*W*exp(-a*(h - 2*t));

expn= 2*tau -1 ;

%% DIVIDE: Four main parts of the controller

den = 2*(1 - tau); % denominator

p1 = (-a/den - 0.1)*x; % part 1

p2 = Kah/den; % part 2

p3 = (sign(x))^expn; % part 3 (MAJOR issue: incorrect math function)

p4 = abs(x)^den; % part 4 (assume no time-delay in x)

%% and CONQUER: Control signal

uu = p1 - p2*p3*p4; % dxdt = uu + d

u = (uu - x^3)/(1 + x^2); % control container

% state-dependent disturbance

d = 0.1*x;

% d = 0.1*exp(-t)*x;

% differential equation

dx = x^3 + (1 + x^2)*u + d;

end

controlEE on 2 Apr 2024

Edited: Sam Chak on 5 Apr 2024

Hi , Mr.Chak , Yes for sure , I writed this code to plot Rh

tspan   = [0, 5];   % simulation time
h       = 2;
for i = 1:length(t)
    Rh_values(i) = (sin(pi*t(i)/h))^2;  % Calculate Rh at each time step
end
plot(t, Rh_values);
xlabel('Time');
ylabel('Rh(t)');
title('Plot of Rh vs Time');

and this is the full code :

tspan   = [0, 5];   % simulation time           
x0      =  1;       % initial value             
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   
plot(t, x), grid on, xlabel t, ylabel x(t)      
%% Function that describes the Nonlinear System 
function dx = nonlinearSystem(t, x)             
    % parameters (minor divide-and-conquer)
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    Rh  = (sin(pi*t/h))^2;          % Minor issue: incorrect math function
    Wc  = 1.8269;
    W   = inv(Wc);
    Kah = Rh*W*exp(-a*(h - 2*t));
    expn= 2*tau -1 ;
    %% DIVIDE: Four main parts of the controller
    den = 2*(1 - tau);              % denominator
    p1  = (-a/den - 0.1)*x;         % part 1
    p2  = Kah/den;                  % part 2
    p3  = (sign(x))^expn;           % part 3 (MAJOR issue: incorrect math function)
    p4  = abs(x)^den;               % part 4 (assume no time-delay in x)
    %% and CONQUER: Control signal
    uu  = p1 - p2*p3*p4;            % dxdt = uu + d
    u   = (uu - x^3)/(1 + x^2);     % control container
    % state-dependent disturbance                          
    d   = 0.1*x;                               
    % d = 0.1*exp(-t)*x;                                
    % differential equation                    
    dx  = x^3 + (1 + x^2)*u + d;                
end                                  

and also :

tspan = [0, 5]; % simulation time

x0 = 1; % initial value

[t, x] = ode45(@nonlinearSystem, tspan, x0);

% Plot the state variable x

subplot(2,1,1);

plot(t, x);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

grid on;

xlabel('t');

ylabel('x(t)');

title('State Variable x(t)');

% Plot the control signal u

u = zeros(size(t));

for i = 1:length(t)

%% ------ added by Sam ------

a = 0.1;

h = 2;

tau = 5/7;

Rh = (sin(pi*t(i)/h))^2;

Wc = 1.8269;

W = 1/Wc;

Kah = Rh*W*exp(-a*(h - 2*t(i)));

%% ------ added by Sam ------

u(i) = (-0.1*x(i)/(2*(1 - 5/7)) - Kah/(2*(1 - 5/7)) * (sign(x(i)))^(2*5/7 - 1) * abs(x(i))^(2*(1 - 5/7)) - x(i)^3)/(1 + x(i)^2);

end

subplot(2,1,2);

plot(t, u);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

grid on;

xlabel('t');

ylabel('u(t)');

title('Control Signal u(t)');

function dx = nonlinearSystem(t, x)

% parameters

a = 0.1;

h = 2;

tau = 5/7;

Rh = (sin(pi*t/h))^2;

Wc = 1.8269;

W = inv(Wc);

Kah = Rh*W*exp(-a*(h - 2*t));

expn = 2*tau - 1;

% Four main parts of the controller

den = 2*(1 - tau);

p1 = (-a/den - 0.1)*x;

p2 = Kah/den;

p3 = (sign(x))^expn;

p4 = abs(x)^den;

% Control signal

uu = p1 - p2*p3*p4;

u = (uu - x^3)/(1 + x^2);

% State-dependent disturbance

d = 0.1*x;

% Differential equation

dx = x^3 + (1 + x^2)*u + d;

end

the error :

>> untitled5

Warning: Imaginary parts of complex X and/or Y arguments ignored.

> In untitled5 (line 8)

Unrecognized function or variable 'Kah'.

Error in untitled5 (line 17)

u(i) = (-0.1*x(i)/(2*(1 - 5/7)) - Kah/(2*(1 - 5/7)) * (sign(x(i)))^(2*5/7 - 1) * abs(x(i))^(2*(1 - 5/7)) - x(i)^3)/(1 + x(i)^2);

>>

controlEE on 2 Apr 2024

hi @Sam Chak , thank you for the links , Why do we need to plot Rh ? we have Kah already plotted , I think we might have problem with sign function causing inaccurate results , Do you have any idea how to plot u signal in Divided code not using a loop ?

Sam Chak on 5 Apr 2024

The error message "Unrecognized function or variable 'Kah'" indicates that the 'Kah' variable cannot be found in the for-loop script, which is responsible for computing the control signal u(i). In simpler terms, the algorithm is unable to locate the necessary information it needs to perform the calculations. To resolve this, you need to ensure that the required input is provided to the algorithm in the correct format.

I have made some modifications to the for-loop section based on your coded functions for 'Rh', so that the control signal u can be plotted. However, please note that the

function in the code is NOT a piecewise function, which causes the control signal u to oscillate within the time frame of

. During this period,

should not produce any output.

OP's question: Do you have any idea how to plot u signal in Divided code not using a loop?

If you don't wish to use the for-loop approach, then directly use the Vectorization method.

function dx = odefun(t, x)

u = - 2*x;

dx = sin(x) + u;

end

[t, x] = ode45(@odefun, [0 10], 1);

ut = - 2*x; % <-- vectorize the control algorithm

plot(t, x, '-o', t, ut, '-o'), grid on, legend('x(t)', 'u(t)')

Answer 3

Sam Chak on 2 Apr 2024

1 vote

Hi @controlEE,

If you define

solely as a continuous-time function, it would appear as displayed in Figures 1 and 2.

h = 2;

Rh = @(t) (sin(pi*t/h)).^2;

%% Plot over time range 2 < t < 4 ... (h < t < 2*h)

t = linspace(2, 4, 2001);

plot(t, Rh(t)), grid on, ylim([-0.5, 1.5])

xlabel('t'), ylabel('R_{h}'), title('Fig.1: R_{h} defined for 2 < t < 4')

%% Plot over time range from 0 to 6

t = linspace(0, 6, 6001);

plot(t, Rh(t)), grid on, ylim([-0.5, 1.5])

xlabel('t'), ylabel('R_{h}'), title('Fig.2: R_{h} plotted over 0 < t < 6')

Piecewise function

However, I have been trying to emphasize that the true

is actually a piecewise function, as defined in the PDF. This needs to be addressed before moving forward, as what you have implemented in the ODE solver may not qualify as Prescribed-Time Control.

%% Plot True Rh according to the definition in the PDF

t1 = linspace(0, 2, 2001);

t2 = linspace(2, 4, 2001);

Rh1 = zeros(1, numel(t1)); % time interval t1 vector is used

Rh2 = (sin(pi*t2/h)).^2; % time interval t2 vector is used

plot(t1, Rh1, 'color', '#265EF5'), hold on

plot(t2, Rh2, 'color', '#265EF5'), hold off, grid on, ylim([-0.5, 1.5])

xline(2, '--');

text(0.7, 1.25, 'Interval 1')

text(2.7, 1.25, 'Interval 2')

xlabel('t'), ylabel('R_{h}'), title('Fig.3: True R_{h} defined according to PDF')

By the way, the straight-line triangle graph shown is definitely not representative of

. You should be able to recognize that something is amiss, as

contains purely sinusoidal functions, resulting in a curvaceous graph.

3 Comments
Show 1 older comment Hide 1 older comment

Sam Chak on 5 Apr 2024

You're welcome! I simply provided the essential guidance for you to gain insight into how mathematical functions behave across a function's domain. This understanding is crucial for studying dynamics and control. Merely observing the math notations of equations won't facilitate learning as it cannot be taught.

If you have these 3 items checked, you'll be prepared to proceed with solving the delay differential equation:

Is the 'Rh' function correctly coded? ✔️
Is the 'sig' function correctly coded? ✔️
Check if the ode45 solver generates a result that closely resembles the one depicted in Figure 3 of the article. ✔️

controlEE on 5 Apr 2024

Edited: controlEE on 5 Apr 2024

Hi @Sam Chak

ok then , for Rh we figured out that we should utilize it as a piecewise function with the second sub-function is a sine-squared function , in the code it should be divided into two terms Rh1 and Rh2 , in other hand we can say Rh should yield to zero in t = (0,2) , which leads to some changes in the parameter part ;

%% Plot True Rh according to the definition in the PDF
h = 2 ; 
t1  = linspace(0, 2, 2001);
t2  = linspace(2, 4, 2001);
Rh1 = zeros(1, numel(t1));  % time interval t1 vector is used
Rh2 = (sin(pi*t2/h)).^2;    % time interval t2 vector is used
plot(t1, Rh1, 'color', '#265EF5'), hold on
plot(t2, Rh2, 'color', '#265EF5'), hold off, grid on, ylim([-0.5, 1.5])
xline(2, '--');
text(0.7, 1.25, 'Interval 1')
text(2.7, 1.25, 'Interval 2')
xlabel('t'), ylabel('R_{h}')

----------------------------------------

% parameters

a = 0.1;

h = 2;

tau = 5/7;

Rh = (sin(pi*t/h))^2; incorrect -----> Rh1 and Rh2 correct

Wc = 1.8269;

W = inv(Wc);

Kah = Rh*W*exp(-a*(h - 2*t));

expn= 2*tau -1 ;

-----------------------------------------

for the sign function as defined in notation (aslo pointed out by you ) :

as you found out that , and It would definitely lead to some changes in the p3

p3 = (abs(x))^expn * sign(x)

let's first see the sign itself ;

before :

x   = linspace(-1, 1, 20001);
tau = 5/7;
expn= 2*tau - 1     % exponent is a fraction 2143/5000
y   = (sign(x)).^expn;
plot(x, y, 'linewidth', 2), grid on, ylim([0, 1.2])
xlabel('x'), ylabel('y')
title({'$\left[\mathrm{sign}(x)\right]^{2 \tau - 1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

and now ;

x   = linspace(-1, 1, 20001);
tau = 5/7;
expn = 2*tau - 1     % exponent is a fraction 2143/5000
y   = sign(x).*(abs(x)).^expn ;
plot(x, y, 'linewidth', 2), grid on, ylim([0, 1.2])
xlabel('x'), ylabel('y')
title({'$\left[\mathrm{sign}(x)\right]^{2 \tau - 1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

and by editing this part our state gets closer to the article ;

tspan   = [0, 5];   % simulation time           
x0      =  1;       % initial value             
[t, x]  = ode45(@nonlinearSystem, tspan, x0);   
plot(t, x), grid on, xlabel t, ylabel x(t)      
                                               
%% Function that describes the Nonlinear System 
function dx = nonlinearSystem(t, x)             
    % parameters (minor divide-and-conquer)
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    Rh  = (sin(pi*t/h))^2;          % Minor issue: incorrect math function
    Wc  = 1.8269;
    W   = inv(Wc);
    Kah = Rh*W*exp(-a*(h - 2*t));
    expn= 2*tau -1 ;
    %% DIVIDE: Four main parts of the controller
    den = 2*(1 - tau);              % denominator
    p1  = (-a/den - 0.1)*x;         % part 1
    p2  = Kah/den;                  % part 2
    p3  = sign(x).*(abs(x))^expn;   % part 3 
    p4  = abs(x)^den;               % part 4 (assume no time-delay in x)
    
    %% and CONQUER: Control signal
    uu  = p1 - p2*p3*p4;            % dxdt = uu + d
    u   = (uu - x^3)/(1 + x^2);     % control container
                                            
    % state-dependent disturbance                          
    d   = 0.1*x;                               
    % d = 0.1*exp(-t)*x;                                
    
    % differential equation                    
    dx  = x^3 + (1 + x^2)*u + d;                
end                                  

lastly I think we should change Rh and seperate it in Rh1 and Rh2 . and check the ode45 solver result .

Answer 4

Sam Chak on 5 Apr 2024

Edited: Sam Chak on 5 Apr 2024

1 vote

Great job on your progress shown in the comment. Now, I'd like to share my piecewise function formula with you:

where

is the Heaviside Step Function.

Update: The article also discusses the piecewise function for the time interval

in Definition 1. Although the system remains stable, it is important to note that the settling time is significantly longer than the prescribed-time control performance, which is indicated as 3.4 seconds.

If you find the piecewise formula and the code provided helpful, I kindly request your consideration to vote 👍 for the solution presented in this answer. Your feedback and support are greatly appreciated.

%% Setting Unit Step function (to be used in constructing Rh)

old = sympref('HeavisideAtOrigin', 1);

%% Call ode45 solver

tspan = [0, 6]; % simulation time

x0 = 1; % initial value

[t, x] = ode45(@nonlinearSystem, tspan, x0);

%% Pass time vector t and solution x to nonlinearSystem to get Rh and u

[~, Rh, u] = nonlinearSystem(t, x);

%% Plot results

subplot(311)

plot(t, x, 'linewidth', 1.5, 'Color', [0.8500, 0.3250, 0.0980])

grid on, ylabel x(t), title('Response, x(t)')

subplot(312)

plot(t, Rh, 'linewidth', 1.5, 'Color', [0.9290, 0.6940, 0.1250])

grid on, ylabel Rh(t), title('Piecewise function, Rh')

subplot(313)

plot(t, u, 'linewidth', 1.5, 'Color', [0.4660, 0.6740, 0.1880])

grid on, ylabel u(t), title('Control action, u'), xlabel t

%% Nonlinear System

function [dx, Rh, u] = nonlinearSystem(t, x)

% parameters (minor divide-and-conquer)

a = 0.1;

h = 2;

tau = 5/7;

%% Piecewise function Rh using the piecewise formula

Rh1 = 0; % sub-function at interval 1

Rh2 = (sin(pi*t/h)).^2; % sub-function at interval 2

Rh3 = 0; % sub-function at interval 3

Rh = Rh1 + (Rh2 - Rh1).*heaviside(t - h) + (Rh3 - Rh2).*heaviside(t - 2*h);

Wc = 1.8269;

W = inv(Wc);

Kah = Rh.*W.*exp(-a*(h - 2*t));

expn= 2*tau - 1;

den = 2*(1 - tau); % denominator

%% DIVIDE: Four main parts of the controller

p1 = (-a/den - 0.1)*x; % part 1

p2 = Kah/den; % part 2

sig = sign(x).*(abs(x)).^expn; % part 3 renamed to 'sig'

p4 = abs(x).^den; % part 4

%% and CONQUER: Control signal

uu = p1 - p2.*sig.*p4; % dxdt = uu + d

u = (uu - x.^3)./(1 + x.^2); % control container

% state-dependent disturbance

d = 0.1*x;

% d = 0.1*exp(-t)*x;

% differential equation

dx = x.^3 + (1 + x.^2).*u + d;

end

45 Comments
Show 43 older comments Hide 43 older comments

controlEE on 9 Apr 2024

Moved: Sam Chak on 9 Apr 2024

Hi @Sam Chak , I edited it as bellow ; there are some problems : 1- we should replaced our nonlinear system with ddefun 2- define a history (noticeably it's defined for two delays ) 3- defining dde23 in the first

%% Setting Unit Step function (to be used in constructing Rh)
old     = sympref('HeavisideAtOrigin', 1);
%% Call ode45 solver
tspan   = [0, 6];       % simulation time           
x0      =  1;           % initial value      
lags = 2 ; 
sol = dde23(@ddefun,lags,history,tspan);
%% Pass time vector t and solution x to nonlinearSystem to get Rh and u
[~, Rh, u] = ddefun(t, x);
%% Plot results
subplot(311)
%plot(t, x,  'linewidth', 1.5, 'Color', [0.8500, 0.3250, 0.0980])  sol.x
%grid on, ylabel x(t),  title('Response, x(t)')
%subplot(312)
%plot(t, Rh, 'linewidth', 1.5, 'Color', [0.9290, 0.6940, 0.1250])
%grid on, ylabel Rh(t), title('Piecewise function, Rh')
%subplot(313)
%plot(t, u,  'linewidth', 1.5, 'Color', [0.4660, 0.6740, 0.1880]) sol.u
%grid on, ylabel u(t),  title('Control action, u'), xlabel t
% history 
function s = history(t) % history function for t <= 0
  s = ones(1);
end
%% Nonlinear System 
function [dx, Rh, u] = ddefun(t, x)
    % parameters (minor divide-and-conquer)
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    
    %% Piecewise function Rh using the piecewise formula
    Rh1 = 0;                        % sub-function at interval 1
    Rh2 = (sin(pi*t/h)).^2;         % sub-function at interval 2
    Rh3 = 0;                        % sub-function at interval 3
    Rh  = Rh1 + (Rh2 - Rh1).*heaviside(t - h) + (Rh3 - Rh2).*heaviside(t - 2*h);
    
    Wc  = 1.8269;
    W   = inv(Wc);
    Kah = Rh.*W.*exp(-a*(h - 2*t));
    expn= 2*tau - 1;
    den = 2*(1 - tau);              % denominator
    %% DIVIDE: Four main parts of the controller
    p1  = (-a/den - 0.1)*x;         % part 1
    p2  = Kah/den;                  % part 2
    sig = sign(x).*(abs(x)).^expn;  % part 3 renamed to 'sig'
    p4  = abs(xlag1).^den;              % part 4
    
    %% and CONQUER: Control signal
    uu  = p1 - p2.*sig.*p4;         % dxdt = uu + d
    u   = (uu - x.^3)./(1 + x.^2);  % control container
                                            
    % state-dependent disturbance
    d   = 0.1*x;
    % d = 0.1*exp(-t)*x;
    
    % differential equation
    dx  = x.^3 + (1 + x.^2).*u + d;
end
 

controlEE on 10 Apr 2024

Edited: controlEE on 10 Apr 2024

Hi @Sam Chak , I modified it as bellow , but steel ecountring some syntax errors :

%% Setting Unit Step function (to be used in constructing Rh)
old     = sympref('HeavisideAtOrigin', 1);
%% Call ode45 solver
tspan   = [0, 6];       % simulation time           
x0      =  1;           % initial value             
sol = dde23(@ddefun,lags,@history,tspan);
Unrecognized function or variable 'lags'.
%% Pass time vector t and solution x to nonlinearSystem to get Rh and u
[~, Rh, u] = ddefun(t, x , z);
%% Plot results
subplot(311)
plot(t, x,  'linewidth', 1.5, 'Color', [0.8500, 0.3250, 0.0980])
grid on, ylabel x(t),  title('Response, x(t)')
subplot(312)
plot(t, Rh, 'linewidth', 1.5, 'Color', [0.9290, 0.6940, 0.1250])
grid on, ylabel Rh(t), title('Piecewise function, Rh')
subplot(313)
plot(t, u,  'linewidth', 1.5, 'Color', [0.4660, 0.6740, 0.1880])
grid on, ylabel u(t),  title('Control action, u'), xlabel t
%% Nonlinear System 
function [dx, Rh, u] = ddefun(t, x,z)
    % parameters (minor divide-and-conquer)
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    lags = 2 ; 
    xlag1 = z(:,1);
    %% Piecewise function Rh using the piecewise formula
    Rh1 = 0;                        % sub-function at interval 1
    Rh2 = (sin(pi*t/h)).^2;         % sub-function at interval 2
    Rh3 = 0;                        % sub-function at interval 3
    Rh  = Rh1 + (Rh2 - Rh1).*heaviside(t - h) + (Rh3 - Rh2).*heaviside(t - 2*h);
    
    Wc  = 1.8269;
    W   = inv(Wc);
    Kah = Rh.*W.*exp(-a*(h - 2*t));
    expn= 2*tau - 1;
    den = 2*(1 - tau);              % denominator
    %% DIVIDE: Four main parts of the controller
    p1  = (-a/den - 0.1)*x;         % part 1
    p2  = Kah/den;                  % part 2
    sig = sign(x).*(abs(x)).^expn;  % part 3 renamed to 'sig'
    p4  = abs(xlag1).^den;              % part 4
    
    %% and CONQUER: Control signal
    uu  = p1 - p2.*sig.*p4;         % dxdt = uu + d
    u   = (uu - x.^3)./(1 + x.^2);  % control container
                                            
    % state-dependent disturbance
    d   = 0.1*x;
    % d = 0.1*exp(-t)*x;
    
    % differential equation
    dx  = x.^3 + (1 + x.^2).*u + d;
end
function s = history(t) 
s = ones(1);
end 

controlEE on 10 Apr 2024

Edited: controlEE on 10 Apr 2024

ok , I define lags in the dde function , but I moved it before dde23 and errors went off

yes Mr.Cahk I define lags = 2 , because h = 2 , now if I remove subplot I will get different answer , but how can I plot Rh and u here

syms t;
syms x; 
syms z;
%% Setting Unit Step function (to be used in constructing Rh)
old     = sympref('HeavisideAtOrigin', 1);
%% Call ode45 solver
tspan   = [0, 6];       % simulation time           
x0      =  1;           % initial value 
lags = 2;
sol = dde23(@ddefun,lags,@history,tspan);
%% Pass time vector t and solution x to nonlinearSystem to get Rh and u
[dx, Rh, u] = ddefun(t, x , z);
%% Plot results
t = sol.x;
x = sol.y;
subplot(311)
plot(t, x,'b')
%grid on, ylabel x(t),  title('Response, x(t)')
%subplot(312)
%plot(t, Rh, 'linewidth', 1.5, 'Color', [0.9290, 0.6940, 0.1250])
%grid on, ylabel Rh(t), title('Piecewise function, Rh')
%subplot(313)
%plot(t, u,  'linewidth', 1.5, 'Color', [0.4660, 0.6740, 0.1880])
%grid on, ylabel u(t),  title('Control action, u'), xlabel t
%% Nonlinear System 
function [dx, Rh, u] = ddefun(t,x,z)
    % parameters (minor divide-and-conquer)
    a   = 0.1;
    h   = 2;
    tau = 5/7;
    %lags = 2 ; 
    xlag1 = z(:,1);
    %% Piecewise function Rh using the piecewise formula
    Rh1 = 0;                        % sub-function at interval 1
    Rh2 = (sin(pi*t/h)).^2;         % sub-function at interval 2
    Rh3 = 0;                        % sub-function at interval 3
    Rh  = Rh1 + (Rh2 - Rh1).*heaviside(t - h) + (Rh3 - Rh2).*heaviside(t - 2*h);
    
    Wc  = 1.8269;
    W   = inv(Wc);
    Kah = Rh.*W.*exp(-a*(h - 2*t));
    expn= 2*tau - 1;
    den = 2*(1 - tau);              % denominator
    %% DIVIDE: Four main parts of the controller
    p1  = (-a/den - 0.1)*x;         % part 1
    p2  = Kah/den;                  % part 2
    sig = sign(x).*(abs(x)).^expn;  % part 3 renamed to 'sig'
    p4  = abs(xlag1).^den;              % part 4
    
    %% and CONQUER: Control signal
    uu  = p1 - p2.*sig.*p4;         % dxdt = uu + d
    u   = (uu - x.^3)./(1 + x.^2);  % control container
                                            
    % state-dependent disturbance
    %d   = 0.1*x;
    d = 0.1*exp(-t)*x;
    
    % differential equation
    dx  = x.^3 + (1 + x.^2).*u + d;
end
function s = history(t) 
s = ones(1);
end 

Sam Chak on 21 Apr 2024

The 2nd problem also involves solving a time-delayed differential equation (DDE), but this time, the system in question is a third-order system, which adds complexity to the task. Therefore, I recommend posting it as a New Question so that other experts can contribute their insights and assistance.

To ensure effective communication and engagement, it would be helpful to clearly describe the DDE problem accompanied by relevant images. These images should include:

The mathematical model of the system,
The system's parameters,
The time-delayed state-feedback control equation, u, and
All design parameters in u.

It is worth noting that many people may lose interest in reading full articles (PDFs) to search for equations and information if they are unfamiliar with the symbols and terminology. Hence, including code snippets of your initial attempt using the dde23 solver would be beneficial.

You can mention that the inspiration for this approach came from a similar solution provided in this thread (don't forget to include the URL link). People are often more inclined to point out errors or suggest improvements in existing code rather than writing the code from scratch.

Sam Chak on 24 Apr 2024

I noticed a few key points regarding your recent post that may have contributed to it not receiving sufficient attention from the experts:

Key #1:

In your post, item 3 (the control equation for the second example) and the requested MATLAB code snippet are missing. Including these essential details could help generate more interest from the experts.

Key #2:

Based on the information provided in item 2, it appears relatively straightforward to code item 1. If you revise your post to include all the recommended materials (items 1–4), it will enable experts to run your code and provide valuable feedback, such as identifying errors or suggesting improvements.

Key #3: (Very important!)

Additionally, it's important to remember to work out your code in MATLAB before transitioning to Simulink. If your code functions correctly in MATLAB, it is highly likely to work in Simulink as well. Since most users are familiar with MATLAB but not necessarily Simulink, including the keyword "Simulink" in your title may cause MATLAB experts to skip your question. Similarly, Simulink experts may lose interest if your post does not include a block diagram or an .slx file.

Sam Chak on 25 Apr 2024

Upon checking the keyword "Backstepping" on Wikipedia, I discovered that it is not a controller in itself but rather a recursive many-step control design procedure. This procedure can be quite tedious, which may explain why the authors did not provide the detailed workings in their paper. Since it involves multiple steps (in this case, three steps due to the system being third-order), it can consume a lot of space to present all the iterations in IEEE paper.

Unfortunately, my professor did not cover this specific topic. However, you have the option to implement other controllers that your professor taught you. If you are interested, I can provide you with a coding template below so you can design a state-feedback controller, via feedback linearization. Consider it a mini challenge to further develop your skills in controller design.

Also, please insert the provided code template in your new question and ensure that you revise the title to exclude the word "Simulink" if you would like to attract MATLAB experts to review your code.

q0 = -0.2;

Dq0 = -0.1;

I0 = 0.1;

tspan = [0, 20]; % set the simulation time

x0 = [q0; Dq0; I0]; % initial condition

[t, x] = ode45(@onelink, tspan, x0);

plot(t, x), grid on, xlabel('t'), title('Response')

%% One-link Manipulator

function [dx, u] = onelink(t, x)

% reference signals

q = (pi/2)*sin(t)*(1 - exp(- 0.1*t^2));

dq = exp(-0.1*t^2)*((-pi/2 + pi/2*exp(0.1*t^2))*cos(t) + 0.1*pi*t*sin(t));

Vp = 0.1*sin(50*pi*t);

% parameters

B = 1;

N = 1;

M = 1;

R = 1;

KB = 1;

L = 1;

% controller

u = 0; % to be designed by the User

% differential equations

dx(1,1) = x(2);

dx(2,1) = (N/M)*sin(q) - (N/M)*sin(x(1) + q) - (B/M)*x(2) + (1/M)*x(3);

dx(3,1) = - (R/L)*x(3) - (KB/L)*x(2) + (1/L)*u + (1/L)*Vp;

end

Sam Chak on 24 May 2024

I have reviewed the Simulink model (Fig. 1) and set the initial value to 1 to test if your designed controller can track the reference signal

. The result, as shown in Fig. 2, indicates that the system's output failed to track the reference signal

.

I understand that the Integral action in your PI controller can reduce or eliminate the steady-state error when properly tuned. However, in theory, in the absence of any unknown disturbances, such a single integrator system

does not require a PI controller. A Proportional-only controller would be sufficient to perfectly track the given reference signal

.

Can you redesign the Proportional-only controller based on the concept of driving the error

to zero? You should be working on the error dynamics

, and not the state dynamics

.

Please provide the relevant design equations so that we can review the design steps as you would publish in the Automatica Journal or IEEE Control Systems Letters. Your MATLAB code (preferred) or Simulink model should be based on the derived equations.

- - - - - - - - - - - - - - - - - - - -

Perhaps, start from here...

- - - - - - - - - - - - - - - - - - - -

Figure 1: Block diagram

Figure 2: Result

Figure 3: If

and

.

controlEE on 13 Jun 2024

Edited: Sam Chak on 14 Jun 2024

Hi @Sam Chak

I kinda get to the following code , can u troubleshoot and help me modify it based on the example 2 , I mean I defined all the parameters but when I run it and it kinda deosn't respond properly , if u have any suggestion I would like to hear and really appreciate it ,We are on this step for about a mounth .

thanks

% Parameters from the article

J = 1.625e-3; % kg·m²

m = 0.506; % kg

L0 = 0.305; % m

R0 = 0.023; % m

B0 = 16.25e-3; % N·m·s/rad

L = 25e-3; % H

R = 5; % Ω

K_tau = 0.09; % N·m/A

Ke = 0.09; % V·s/rad

G = 9.81; % m/s²

% Derived parameters

M = J / K_tau + m * L0^2 / (3 * K_tau) + 2 * m * R0^2 / (5 * K_tau);

N = m * L0 * G / (2 * K_tau) + m * L0 * G / K_tau;

B = B0 / K_tau;

% Initial conditions and parameters

tspan = [0, 10]; % Example time span

x0 = [-0.2, -0.1, 0.1]; % Initial conditions

a = -0.15; % Example value

tau = 81/83; % Example value

h = 5; % Example value

lags = [h]; % Delay values

% Solve the DDE

sol = dde23(@(t, x, xtau) manipulator_dynamics(t, x, xtau, a, tau, h), lags, x0, tspan);

% Plot results

figure;

subplot(3,1,1);

plot(sol.x, sol.y(1,:));

xlabel('Time (s)');

ylabel('x1 (error in q)');

title('State x1 vs Time');

subplot(3,1,2);

plot(sol.x, sol.y(2,:));

xlabel('Time (s)');

ylabel('x2 (error in dq)');

title('State x2 vs Time');

subplot(3,1,3);

plot(sol.x, sol.y(3,:));

xlabel('Time (s)');

ylabel('x3 (error in I)');

title('State x3 vs Time');

% R_h(t) function

function R_h = Rh(t, h)

if t >= 0 && t <= h

R_h = 0;

elseif t > h && t <= 2*h

R_h = (t - h)^5 * (2*h - t)^5;

else

R_h = 0;

end

% K_n(t) function

function K_n = Kn(t, a, h)

if t <= 0

K_n = 0;

else

K_n = a * Rh(t, h);

end

% rho(t, x) function

function rho_val = rho(t, x)

rho_val = 0.1 * (1 + norm(x));

end

% Control input function

function u = control_law(t, x, xtau, a, tau, h)

% Retrieve parameters

M = 1.625; % kg·m²/N·m

N = 0.5; % kg·m/s²/N·m

B = 0.01; % N·m·s/rad/N·m

R = 5; % Ω

L = 25e-3; % H

Ke = 0.09; % V·s/rad

% Desired trajectory

q_d = (pi/2) * sin(t) * (1 - exp(-0.1 * t^2));

dq_d = (pi/2) * cos(t) * (1 - exp(-0.1 * t^2)) + pi * sin(t) * exp(-0.1 * t^2) * t;

ddq_d = (pi/2) * (-sin(t) * (1 - exp(-0.1 * t^2)) - 0.1 * t^2 * sin(t) * exp(-0.1 * t^2)) + pi * cos(t) * exp(-0.1 * t^2) * t + pi * sin(t) * exp(-0.1 * t^2);

% Errors

e1 = x(1) - q_d;

e2 = x(2) - dq_d;

e3 = x(3) - xtau(3);

% Control input

Kn_t = Kn(t, a, h);

rho_val = rho(t, x);

z3 = e3; % Assuming z3 is equal to e3, you may need to define it based on the article

u = ((-a / (2 * (1 - tau))) * z3 ...

+ ddq_d ...

+ (N / M) * sin(q_d) - (N / M) * sin(x(1) + q_d) ...

- rho_val * sign(z3) ...

- (Kn_t / (2 * (1 - tau))) * sign(z3 * abs(z3 - h)^(2 * (1 - tau)))^(2*tau-1))*1/L;

end

% DDE function

function dxdt = manipulator_dynamics(t, x, xtau, a, tau, h)

% Retrieve parameters

M = 1.625e-3; % kg·m²/N·m

N = 0.506; % kg·m/s²/N·m

B = 0.01; % N·m·s/rad/N·m

R = 5; % Ω

L = 25e-3; % H

Ke = 0.09; % V·s/rad

% Control input

u = control_law(t, x, xtau, a, tau, h);

% Dynamics

dxdt = [x(2);

(N/M) * sin(xtau(1)) - (N/M) * sin(x(1) + xtau(1)) - (B/M) * x(2) + (1/M) * x(3);

-(R/L) * x(3) + (Ke/L) * x(2) + (1/L) * u];

end

Sam Chak on 26 Jun 2024

Exmp2.m

At the moment, your new post of the question (Example 2) has not yet appeared. I sincerely hope you will be able to share it soon, as this request has been postponed several times. Long thread with too many posts, lengthy codes and graphical figures makes loading the page relatively slow and uneasy to navigate.

I have reviewed the provided PDF. However, I was unable to locate your personalized step-by-step control design for the specific Example 2. The PDF in Section IV-A only outlines the general procedure for designing the proposed controller for an nth-order nonlinear system. This poses a technical challenge for me to thoroughly evaluate your code.

That said, one obvious issue I have observed is that you have implemented the controller using Eq. (26) without first designing Eq. (23). Since Example 2 is a 3rd-order system, Eq. (26) requires you to design

using Eq. (23). However, you have assigned a value of zero to this parameter 'x3_d', without any accompanying annotations in the code to explain the rationale behind this decision.

function dxdt = manipulator_dynamics(t, x, Z, K_a_h)
    x1_d = 0;
    x2_d = 0;
    x3_d = 0;
    ...
    z1 = x(1) - x1_d;
    z2 = x(2) - x2_d;
    z3 = x(3) - x3_d;
    z3lag = Z(3,1) - x3_d;
    % Control input
    u = L * ((-a / (2 * (1 - tau))) * z3 + x3_d ...
        + R / L * x(3) + K_B / L * x(2) ...
        - 5 * sign(z3) ...
        - (K_a_h(t) / (2 * (1 - tau))) * (sig(z3 * (abs(z3lag)^(2 * (1 - tau)))))^(2 * tau - 1));
    ...
end

controlEE on 26 Jun 2024

Edited: controlEE on 26 Jun 2024

Hi @Sam Chak , you're right , Now by assigning values for x2d and x3d , I can see that our controller works , here I canculate values for these two based on the Eq.(23) , and F , G based on the Eq.(30) ;

and I came to these results , could you please check them whether I calculated them correctly or not ?

and how can I assign x2d and x3d in my code , I think firstly we need to define k1,k2 and ....

Besides that , here for manipulator I have 2 small question questions ;

1- If we conside Vp/L as the uncertainty , it equals with a constant value , by asumming x = 0 , it isn't equal with zero , which was a condition mentioned below of the Eq.(22) .

2- here in x2(dot) , we can see there is a sine term which is a function of x2 , so is it still in the Strict-Feedback form ?

Sam Chak on 27 Jun 2024

https://www.mathworks.com/matlabcentral/answers/2132131-design-of-strong-prescribed-time-spt-controller-for-manipulator-example-2

Could you move your "Comment" to your New Question?

Answer 5

Elijah on 19 Feb 2025

0 votes

what the sigma am i looking at

2 Comments
Show None Hide None

Sam Chak on 19 Feb 2025

Hi @Elijah

It is not sigma, but rather an Absurd Sigmoidal function. The term 'Absurd' in this context is a compound word derived from 'Absolute' and 'Surd'. Essentially, it is defined as

.

However, this notation is not commonly used in the field of mathematics and may cause confusion for some individuals, as is evident in your case. Typically, control theorists who are familiar with Super Twisting Control, Terminal Sliding Mode Control, and Prof. Han's Nonlinear ADRC (Active Disturbance Rejection Control) will recognize this notation.

Some control theorists, such as Prof. Leonid Fridman, use the notation

instead.

x = linspace(-1, 1, 2001);

y = (abs(x).^(1/4)).*sign(x);

plot(x, y), grid on, xlabel('x'), ylabel('y')

●

controlEE on 10 Nov 2025

Edited: controlEE on 10 Nov 2025

Hi @Sam Chak, I hope you are feeling sound and good. I have designed a prescribed time disturbance observer based on the idea that was used to design a prescribed time controller. The problem I have is that the disturbance observer does not work in a prescribed time, I change the settling time, but it still works in a fixed settling time that is around 1.5 seconds. I don't know where the problem comes from.