Effacer les filtres
Effacer les filtres

a nonlinear system and its control input in simulink

120 vues (au cours des 30 derniers jours)
controlEE
controlEE le 16 Mar 2024
Commenté : Sam Chak le 27 Juin 2024 à 8:33
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

Réponse acceptée

Sam Chak
Sam Chak le 14 Avr 2024
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 commentaires
controlEE
controlEE le 18 Avr 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
Sam Chak le 18 Avr 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.

Connectez-vous pour commenter.

Plus de réponses (3)

Sam Chak
Sam Chak le 16 Mar 2024
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 commentaires
controlEE
controlEE le 2 Avr 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
Sam Chak le 5 Avr 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)')

Connectez-vous pour commenter.


Sam Chak
Sam Chak le 2 Avr 2024
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 commentaires
Sam Chak
Sam Chak le 5 Avr 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:
  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
controlEE le 5 Avr 2024
Modifié(e) : controlEE le 5 Avr 2024
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 .

Connectez-vous pour commenter.


Sam Chak
Sam Chak le 5 Avr 2024
Modifié(e) : Sam Chak le 5 Avr 2024
Great job on your progress shown in the comment. Now, I'd like to share my piecewise function formula with you:
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 commentaires
controlEE
controlEE le 26 Juin 2024 à 15:36
Modifié(e) : controlEE le 26 Juin 2024 à 15:38
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 ?

Connectez-vous pour commenter.

Catégories

En savoir plus sur Simulink dans Help Center et File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by