a nonlinear system and its control input in simulink

hi guys the system is as follows ;

which k(a,h) is defined as follows ;

I generally implement its block - daigram in simulink as follows ; I'm not sure whether it's correct or not ; I need your help on defining the functions

Answer 1

[Sam Chak]

Hi @controlEE

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.

[controlEE]

Hi @Sam Chak , thanks , how we can plot control signal in the code by the way , besides I develped a simmulink model based on the pic , but actully there are some errors related to sympref and heaviside , I attached simulnik here on comment

[Sam Chak]

Hi @controlEE

I forgot to mention that the heaviside() function doesn't work in Simulink for an unspecified reason. Therefore, it needs to be replaced with a formulated sign() function. The syntax for this replacement is as follows:

Rh = ((sin(pi*t/h)).^2).*((sign(t - h) + 1)/2) - ((sin(pi*t/h)).^2).*((sign(t - 2*h) + 1)/2);

When it comes to plotting simple figures, I prefer to avoid loops and lengthy code. To plot the figures in the simplest manner, I generate the responses in Simulink and then send the signal using two separate 'To Workspace' blocks. You can view the block diagram in the attached Simulink file, 'firstSimu_edited.slx'.

In MATLAB, I can plot the signals and using these two simple lines:

subplot(2,1,1), plot(tout, x), grid on, 
subplot(2,1,2), plot(tout, u), grid on

Please let me know if the replacement works.

[controlEE]

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]

Hi @controlEE

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]

Hi @controlEE

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)

[controlEE]

Hi Mr.Chak , Thank you for your answer

Yes , the control signal which designed here is based on a prescribed time control theory .

No I didn't stimulate this part yet . I also need help in defining functions , or instead of functions can we use simulink blocks ?

Is there any way I can reach you ?

thanks

[Sam Chak]

Hi @controlEE

Apart from the bell-shaped function and the unrecognized function, I don't observe any other special mathematical functions that might require assistance. The function appears to be a piecewise function, with the second sub-function is a sine-squared function. These elementary functions are supported within the Simulink math blocks.

If you search for 'sig' in Wolfram MathWorld, you won't find anything useful. The term was likely informally coined by non-mathematician researchers in the past, and some control practitioners adopted it for ease of reference.

[controlEE]

Hi Mr.Chak

the sig is a discrete function SIGN which can be described as follows :

I want to implement state space for the function on the right and the control signal for the function on the left .

thank you

[Sam Chak]

Hi @controlEE

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)

[controlEE]

Hi Mr.Chak , thank you for the clarification , I used functions to implement the model and control signal : which leads to the below :

this is the complete block - diagram :

I define control function :

and also function on the right for xdot :

I assume d(t,x) which is the uncertainty of the system equals to 0.1x ,

But I have got sevel errors and aslo do not get the ideal results ; based on the article the result should be as follows ;

[Sam Chak]

Hi @controlEE

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

[controlEE]

hi Mr.Chak , thank u , It's a good example you brought here

but I do need help in my model , I mean it's way more complicated .

[Sam Chak]

Hi @controlEE,

The 1st-order system you mentioned doesn't appear to be more complex than what you have described within the MATLAB function block. When you are ready to proceed, I suggest posting the MATLAB code using the template provided in my Example. This will assist in streamlining the process and obtaining the desired outcome.

[controlEE]

Hi Mr.Chak

the block - diagram I implement , I'm not utterly sure about it , I think it would be better if we use only a matlab function for defining the state space , and u see the term d(t,x) in the xdot to , based on the article we should get the results for d(t,x) = 0.1x and also d(t,x) = 0.1exp(-t)x , so first we all integrate xdot to optain the term x , but we should first implement u as an input to the matlab function , should we use matlab function again or can we implement it by blocks ?

[Sam Chak]

Hi @controlEE

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]

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]

despite all of the challenges , 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)) , it would come to be wrong

[controlEE]

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]

Hi @controlEE,

Good effort with the starting code. Could you find and fix the missing bracket(s) in this line first (scroll up)? We will need to troubleshoot this step-by-step. I previously overlooked the time-delayed term ) Later on, we will need to switch from the ode45 solver to the dde solver.

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)

By the way, you should break down the control signal u into multiple parts so that you can troubleshoot the issue. Additionally, is computed as a constant and is not dynamically linked to the nonlinear system. Therefore, you can directly assign "Wc = 1.8269" to save computational effort.

Lastly, I addressed the expression "(sign(x)).^expn" in my previous comment. It's possible that you overlooked it. When x is greater than 0, will always produce 1, and raising 1 to the power of 'expn' will always yield 1. That's why I mentioned that the authors didn't intend "sig(x)" to be "(sign(x)).^expn".

[Sam Chak]

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]

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                                

[Sam Chak]

Hi @controlEE,

You've made good progress so far. Let's now shift our focus to getting the code to run smoothly without any error messages. I've highlighted the parts that need your attention in order to fix them.

By the way, I'm curious. Have you previously implemented any linear full state feedback control algorithms in your prior experience before diving into the realm of nonlinear control?

[controlEE]

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]

Hi @controlEE

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]

Hi mr.Chak

thank you for your help and also important tips that you mentioned here

actully I attach the paper here , yes , Kinda I do read it 2 or 3 times , but I didn't notice any usefull tips around sig function

besides , in the code I also want to plot u , I mean our control signal , How can I do that ? I tried to plot in the form same as x , but it came with an error and didn't show me anything .

[Sam Chak]

@controlEE, Thank you for sharing the article. To find all instances of characters with "sig" in the PDF, you can utilize the "SEARCH" function. Once you have located them, you can proceed with working on the code.

[controlEE]

Hi , you're welcome Mr.Chak , thank u I searched it ,

the thing that I like to mention here is that , the control signal is different in two disturbances , I'd like to plot it too , and also u mentioned that because of sig function the results we concluded here from the code aren't accurate , so how we can troubleshoot the code to get ideal results , by the way the control signal that is mentioned here in the article I should apply it to one - llimk manipulator too on page 7 .

[Sam Chak]

Could you please share the updated Divide-n-Conquer code that addresses the bullet issues #1 and #2 mentioned in my previous comment? This will enable me to troubleshoot and identify any potential errors that occurred on your side and your team members.

If you are encountering difficulties with the Divide-n-Conquer code, let's divide the simulation problem into two sub-problems for the time being. [Ignore this if you already know how to proceed]

Task #1: Plot the piecewise function .

Task #2: Plot the sigmoidal function .

This approach will allow your team member to address each task separately and ensure the accuracy of the plotted functions in part 2 (p2) and part 3 (p3) of the Divide-n-Conquer code.

[controlEE]

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);

>>

[Sam Chak]

Hi @controlEE,

Could you please include the time vector variable 't' in the code? It is necessary for the function . Additionally, you can click the 'Play' button on the MATLAB Answers forum to execute the code. This will help you identify any missing elements in the code when it throws error messages in red. In fact, it's also the method I use to troubleshoot your code.

If you are unfamiliar with the plot function, I suggest referring to the examples provided in this link. None of the examples utilize a for-loop approach.

https://www.mathworks.com/help/matlab/ref/plot.html#btzitot-4

It's important to note that this plotting code is separate from the Divide-n-Conquer code and is not directly related to it. You need to see what the piecewise function looks looks like before inserting the code in . Refer to the PDF, not the coded formula in the my Answer.

[Sam Chak]

Hi @controlEE

I noticed that you made some revisions to your comment by adding additional code. However, when I clicked the Play button to run the simulation, it seemed to complicate the whole matter.

Here is an example of how I plot a continuous-time W-shaped function over the range of . This will give you an idea of how I approach plotting in a simpler manner.

However, do remember that is NOT a continuous-time function but a piecewise function.

%% Define the W-shaped function
tStart  = -2;                              
tFinal  =  2;
numpts  = (tFinal - tStart)*1000 + 1;       % number of data points
t       = linspace(tStart, tFinal, numpts); % time vector
W       = (cos(pi/2*t)).^2;                 % W-shaped function of time

%% Plot the graph
plot(t, W), grid on, ylim([-0.5, 1.5])      % requires inputs from t and W vectors
xlabel('t'), ylabel('W(t)')
title('W-shaped function')

[controlEE]

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]

Hi @controlEE

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]

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.

[controlEE]

Hi @Sam Chak , good evening Thank you very much for clarification I saw how true Rh looks like , now for proceeding with the code , should we use dde instead of ode solver and bring delay in our computations ? I mean to have state and control signals same as the article .

[Sam Chak]

Hi @controlEE

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:

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

[controlEE]

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]

Hi @controlEE

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

[controlEE]

Hi @Sam Chak

I've voted , I really appreciate your help , Thank you very much

I think the reason that our control and state don't reach as same as atricle might be for the delay term we didn't consider in p4 of our controller , and afterwards utilize dde instead of ode .

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

d = 0.1*x -------> exact fixed time stable ------>T=4s

d=0.1exp(-t)*x ------->fixed time stable --------> T=3.5

which they both are equal or smaller than T=2h(h=2s)

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

and also thanks for the H(x).

[Sam Chak]

Hi @controlEE

You're correct that the h-second-delayed state is not considered in part 4. However, the equations in the code are 99% complete. The remaining step is to convert the ODE function into a DDE function so that you can utilize the dde23 solver. DDE requires past history so that the solver can access the solution for times before the initial integration point.

If the paper doesn't mention it, you can assume a constant delay that is the same as the initial value in the state. In other words, since , you can assume that .

To proceed, follow the steps outlined in the link below to code the DDE function. Once you've written the code, please post it, and I'll review it for any errors.

Link: DDE with Constant Delays

[controlEE]

Hi @Sam Chak

as we do have a constant delay in p4 I write it same as the help document as follows ;

tspan = [0, 6];
sol = dde23(@ddefun, lags, @history, tspan);
% Parameters
lags = 2 ; 
tau = 5/7;
den = 2*(1 - tau);              % denominator
% Code equation 
function dydt = ddefun(~ , x , z)
xlag1 = z(:,1);
dydt = abs(xlag1).^den;
end 
% Code Solution history 
function s = history(t)
s = ones(1);                    % assumed 
end 

I wonder we do have this delayed - state only in p4 , are there any other ways to calculate delay ?

[Sam Chak]

You did the 1% in DDE, but forgot to insert 99% of the code from the ODE function into the DDE function.

[controlEE]

kinda confused how to bring this to our code

[Sam Chak]

Just as shown in the DDE example. The entire differential equation shown in the image of your question above is DDE due to a single time delayed state.

There is NO separate partial ODE and partial DDE. Everything (99% + 1%) is grouped under a single DDE function.

[controlEE]

I got it , but when it comes to implement it insted of ode , we should apply some changes in the code and also variables

[controlEE]

I got it , but when it comes to implement it insted of ode , it might demand some serious changes in the code and also variables .

[Sam Chak]

Um... What kind of significant changes might be needed in the code and variables? Are these changes performed in accordance to the example shown in the "DDE with Constant Delays" article?

The derivatives and incorporate time-delayed states on the right-hand side, while represents a regular ordinary differential equation. However, all three dynamic states are grouped together within the 'ddefun()' function. Interestingly, the code remains largely unchanged for the non-time-delayed state, resembling the ODE function.

[controlEE]

alright ,So how do we define our equatios in ddefun ?

I mean first we should determine our delay constants which we only have one delay = 2

so

lags = 2 ;

function dxdt = ddefu(t,x,z)

xlag1 = z(:,1);

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

end

somehow we should bring u signal into our ddefun as the delay was defined in p4

[Sam Chak]

Hi @controlEE

Referring to my previous comments, I have consistently emphasized the importance of placing the "entire code in the ode function" within the dde function. The example has two DDEs, but you only have a single 1st-order delayed differential equation

,

where the ONLY time-delayed state denoted by 'xlag1' lies in part 4 of u:

In the previous code, part 4 was implemented within the nonlinearSystem() function as:

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

However, it should have been coded within the ddefun() function as:

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

This is the reason why I mentioned that 99% of the code remains unchanged.

Please prioritize this step, and I will assist you in identifying the necessary adjustments.

[controlEE]

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
 

[Sam Chak]

Hi @controlEE

While you did add the 99% of the code, you also unintentionally removed the crucial 1% that is necessary for solving delay differential equations using the dde23 solver.

In the example, the ddefun function requires a third parameter, Z, which represents the time-delayed state. Initially, you included this correctly in your comment, but it seems to have been accidentally omitted during the copy/paste process. when changing the function names from 'nonlinearSystem(t, x)' to 'ddefun(t, x)'. Moreover, the assignment of Z to the time-delayed state 'xlag1' has been overlooked.

Please rectify these issues and plot the response using the following syntax:

plot(sol.x, sol.y)

[controlEE]

I got it , thank you , by the way , do the errors arise only because of not defining z in ddefun ?

[Sam Chak]

Indeed, similar to running the ode45 solver, the odefun(t, x) function requires the inclusion of the time parameter t and the state variable x in the input arguments. However, when adhering to the guidelines of the dde23 solver, the ddefun(t, x, z) function requires an additional input argument, specifically the time-delayed state variable z, as the third parameter.

[controlEE]

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 

[Sam Chak]

Hi @controlEE

To troubleshoot the code, it is important to start by identifying the error message and understanding its meaning. To view the error message, simply click the run button, and the message will be displayed.

In this specific case, you have provided 'lags' as an input argument to the dde23 solver, but no variable named 'lags' has been assigned prior to this line. Do you happen to know what value should be assigned to 'lags'?

Additionally, unused variable 'x0' in the code should be removed, and the annotation "Call ode45 solver" should be revised to adhere to Good Engineering Programming Practice.

Update: Previously, you assigned the value of '2' to the variable 'lags'.

[controlEE]

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]

Hi @controlEE

Congratulations! You have successfully implemented the Prescribed-Time Control in the DDE, and the response closely matches the one depicted in the article. It seems that you prefer a more minimalist coding style, then you can utilize the following code instead.

sympref('HeavisideAtOrigin', 1);    
sol = dde23(@ddefun, 2, @history, [0, 6]);
plot(sol.x, sol.y), grid on, xlabel t, ylabel x(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 

[controlEE]

Hi @Sam Chak

thank you very much

I modified the plot term a little to get the control too ;

sympref('HeavisideAtOrigin', 1);    
sol = dde23(@ddefun, 2, @history, [0, 6]);
%plot(sol.x, sol.y), grid on, xlabel t, ylabel x(t)
plot(sol.x, sol.y, sol.x, sol.yp(1, :), 'r'), grid on, xlabel('t'), ylabel('x(t) and 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;    
    %d= 0.1*x
    dx  = x.^3 + (1 + x.^2).*u + d;
end

function s = history(t) 
    s   = 1;
end 

now I want to proceed with the second and last example of the paper .

[controlEE]

Hi @Sam Chak

I wanted to take a moment to sincerely thank you for your help with the problem . Your explanation and guidance were incredibly helpful, and I feel much more confident with the process now.

I was wondering if I could impose on your expertise again. I have another example I'd like to stimulate and learn from, but I'm not quite sure how to approach it. It's similar to the previous one ; I Attached a pdf in the comment and also ;

Would you be available to take a look and offer some pointers on how to set up the simulation for this new example?

Thanks again for your time and support!

Best regards

[Torsten]

If you have no further constraints or expressions to minimize, you can choose phi = phi(t,x1,...,xn) to be any function with phi(t,0,...,0) = 0 for all t >=0. And ode45 will integrate system (22).

So either this is the solution or you did not state your problem properly.

[Sam Chak]

Hi @controlEE

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:

1. The mathematical model of the system,
2. The system's parameters,
3. The time-delayed state-feedback control equation, u, and
4. 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.

[controlEE]

hi , thank you @Sam Chak and @Torsten

thank you for the tips

I asked a new question which provide the information needed ;

https://www.mathworks.com/matlabcentral/answers/2110591-stimulation-of-a-3-order-system-and-its-control-signal-in-simulink

I sincerely appreciate your help , thanks

[controlEE]

Hi @chak , can u help me proceed with second example .

[Sam Chak]

Hi @controlEE

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]

If you observe the number of views on your post, you will notice that the count likely includes views from both @Torsten and myself.

By addressing the three key points ☝️ mentioned earlier and ensuring that you properly incorporate them into your post, you should be able to improve the situation and attract more attention 📈 from the experts. "Marketing skills" are also taught in Engineering Schools.

[controlEE]

hi @Sam Chak

thanks for the tips , I was wondering actually the author didn't mention the control signal in the second example in details unlike first one, by the way I read the part agian, in fact it brought a control signal in theorem 2 ; I attached a pic , I mean it says by adpoting backstepping precedure they overgeneralize this algorithm for higher order systems , so my point is we have the model I mean its dynamics , it mentions them both in state space way and also in a way of defining Fn and Gn and such functions, for modelling or coding where should we start , and also the control signal I do have some confusion about it , it has some terms like zn and ... , by califing these vague things , I'm struggling to code it first and afterwrds go for a simulink model .

[Sam Chak]

Hi @controlEE

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