# Giving a rectangular pulse in ode45

6 views (last 30 days)
Kashish Pilyal on 15 Aug 2022
Edited: Kashish Pilyal on 15 Aug 2022
I am trying to use a rectangular pulse in my system through ode 45, so I decided to split my code and run it 3 times. Once at 0, then at 1 and at 0 again. I am not getting results as expected. As you can see, I have a system code as:
clear all
%%
q=0; %position
v=0; %velocity
a=0; %acceleration
x_0=[q;
v;
a]; %initial conditions
t=linspace(0,5,100);
%parameters of second vehicle
q1=0; %position
v1=0; %velocity
a1=0; %acceleration
x_1=[q1;
v1;
a1]; %initial conditions
%parameters of third vehicle
q2=0; %position
v2=0; %velocity
a2=0; %acceleration
x_2=[q2;
v2;
a2]; %initial conditions
%initial conditions of the platoon
x=[x_0; x_1; x_2];
[t,y1]=ode45(@code,t,x);
%% second step
q=1; %position
v=0; %velocity
a=0; %acceleration
x_0=[q;
v;
a]; %initial conditions
t=linspace(5,50,100);
%parameters of second vehicle
q1=0; %position
v1=0; %velocity
a1=0; %acceleration
x_1=[q1;
v1;
a1]; %initial conditions
%parameters of third vehicle
q2=0; %position
v2=0; %velocity
a2=0; %acceleration
x_2=[q2;
v2;
a2]; %initial conditions
%initial conditions of the platoon
x=[x_0; x_1; x_2];
[t,y2]=ode45(@code,t,x);
%% third step
q=0; %position
v=0; %velocity
a=0; %acceleration
x_0=[q;
v;
a]; %initial conditions
t=linspace(50,60,100);
%parameters of second vehicle
q1=y2(end,4); %position
v1=0; %velocity
a1=0; %acceleration
x_1=[q1;
v1;
a1]; %initial conditions
%parameters of third vehicle
q2=y2(end,7); %position
v2=0; %velocity
a2=0; %acceleration
x_2=[q2;
v2;
a2]; %initial conditions
%initial conditions of the platoon
x=[x_0; x_1; x_2];
[t,y3]=ode45(@code,t,x);
%% concating
time=linspace(0,20,300);
%parameters of 1st
q_1=[y1(:,1);y2(:,1);y3(:,1)];
%parameters of 2nd
q_2=[y1(:,4);y2(:,4);y3(:,4)];
%parameters of 3rd
q_3=[y1(:,7);y2(:,7);y3(:,7)];
plot(time',q_1,time',q_2,time',q_3)
The system response is not as expected from ode45. I found out that the reason it does not match is due to the vehicles depending upon the predecessor's acceleration which does not happen as I shift to ode 45 mid code. I tried giving initial condition as some acceleration but it just messed up the rectangular pulse input that I am giving. Is there any other way to feed the ode45 solver a rectangular pulse?

Torsten on 15 Aug 2022
Edited: Torsten on 15 Aug 2022
I can't evaluate the equations and initial conditions you use. This is what the solver returns - slight over-and undershooting.
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
q=0; %position
v=0; %velocity
a=0; %acceleration
x_0=[q;
v;
a]; %initial conditions
t=linspace(0,5,100);
%parameters of second vehicle
q1=0; %position
v1=0; %velocity
a1=0; %acceleration
x_1=[q1;
v1;
a1]; %initial conditions
%parameters of third vehicle
q2=0; %position
v2=0; %velocity
a2=0; %acceleration
x_2=[q2;
v2;
a2]; %initial conditions
%initial conditions of the platoon
x=[x_0; x_1; x_2];
[t1,y1]=ode45(@code,t,x,options);
%% second step
q=1; %position
v=0; %velocity
a=0; %acceleration
x_0=[q;
v;
a]; %initial conditions
t=linspace(5,50,100);
%parameters of second vehicle
q1=0; %position
v1=0; %velocity
a1=0; %acceleration
x_1=[q1;
v1;
a1]; %initial conditions
%parameters of third vehicle
q2=0; %position
v2=0; %velocity
a2=0; %acceleration
x_2=[q2;
v2;
a2]; %initial conditions
%initial conditions of the platoon
x=[x_0; x_1; x_2];
[t2,y2]=ode45(@code,t,x,options);
%% third step
q=0; %position
v=0; %velocity
a=0; %acceleration
x_0=[q;
v;
a]; %initial conditions
t=linspace(50,70,100);
%parameters of second vehicle
q1=y2(end,4); %position
v1=0; %velocity
a1=0; %acceleration
x_1=[q1;
v1;
a1]; %initial conditions
%parameters of third vehicle
q2=y2(end,7); %position
v2=0; %velocity
a2=0; %acceleration
x_2=[q2;
v2;
a2]; %initial conditions
%initial conditions of the platoon
x=[x_0; x_1; x_2];
[t3,y3]=ode45(@code,t,x,options);
%% concating
time=[t1;t2;t3];
%parameters of 1st
q_1=[y1(:,1);y2(:,1);y3(:,1)];
%parameters of 2nd
q_2=[y1(:,4);y2(:,4);y3(:,4)];
%parameters of 3rd
q_3=[y1(:,7);y2(:,7);y3(:,7)];
plot(time,[q_1,q_2,q_3]) function y=code(t,x)
%-----------legend------
% x(1),x(2),x(3) [position, velocity, acceleration] 1st
% x(4),x(5),x(6) [position, velocity, acceleration] 2nd
% x(7),x(8),x(9) [position, velocity, acceleration] 3rd
tau_i=0.1;
h_i=0.5;
k_d=0.68;
k_p=0.2;
T=tau_i/h_i;
u_i=( (1-T)*x(3) );
%x(2)=0.5*t;
% CACC equations
x_1dot=0; %q(position)dot
x_1ddot=x(3); %v(velocity)dot
x_1dddot=(-(1/tau_i)*x(3))+((1/tau_i)*u_i); %a(accelration)dot
%% second vehicle
tau_i=0.2;
h_i=0.5;
T=tau_i/h_i;
e_i1_1=x(1)-( x(4)+(h_i*x(5)) );
e_i2_1=x(2)-( x(5)+(h_i*x(6)) );
%e_i0_1=int(e_i1_1,t);
u_i=T*[k_p k_d]*[e_i1_1;e_i2_1]+ ( (1-T)*x(6) ) + ( T*x(3) );
% CACC equations
x_2dot=x(5); %q(position)dot
x_2ddot=x(6); %v(velocity)dot
x_2dddot=(-(1/tau_i)*x(6))+((1/tau_i)*u_i); %a(accelration)dot
%% third vehicle
tau_i=0.1;
h_i=0.5;
T=tau_i/h_i;
e_i1_2=x(4)-( x(7)+(h_i*x(8)) );
e_i2_2=x(5)-( x(8)+(h_i*x(9)) );
u_i=T*[k_p k_d]*[e_i1_2;e_i2_2]+ ( (1-T)*x(9) )+ ( T*x(6) );
% CACC equations
x_3dot=x(8); %q(position)dot
x_3ddot=x(9); %v(velocity)dot
x_3dddot=(-(1/tau_i)*x(9))+((1/tau_i)*u_i); %a(accelration)dot
y=[x_1dot; %position of 1st
x_1ddot; %velocity of 1st
x_1dddot; %acceleration of 1st
x_2dot; %position of 2nd
x_2ddot; %velocity of 2nd
x_2dddot; %acceleration of 2nd
x_3dot; %position of 3rd
x_3ddot; %velocity of 3rd
x_3dddot]; %acceleration of 3rd
end
Kashish Pilyal on 15 Aug 2022
Edited: Kashish Pilyal on 15 Aug 2022
Yes, it does. I think the issue here is at the time there is a step increase (at 5 secs from 0 to 1) and at the step decrease. The leader has a change in acceleration which I am not able to put in code properly as the acceleration allows the others to follow suit. The graph should be somewhat like this: No overshoot just proper settling. This is why I am confused how to add rectangular pulse.
Edit:I may have a solution that I can run the ode solver in between again when the leader goes from 0 to 1. I need to put some value of acceleration and make the solver stop when it reaches 1. Any ideas on how to do that. I checked there is something called event function but I am not sure how to implement it here as the documentation does not show it being used in the ode45 code.

Torsten on 15 Aug 2022
Edited: Torsten on 15 Aug 2022
time=linspace(0,20,300);
is wrong.
You must concatenate the t-vectors returned from ode45 in the three calls in the same way as you did it to define q_1, q_2 and q_3.
We cannot test your code since you didn't include "code.m".
Kashish Pilyal on 15 Aug 2022
The code you are asking for is:
function y=code(t,x)
%-----------legend------
% x(1),x(2),x(3) [position, velocity, acceleration] 1st
% x(4),x(5),x(6) [position, velocity, acceleration] 2nd
% x(7),x(8),x(9) [position, velocity, acceleration] 3rd
tau_i=0.1;
h_i=0.5;
k_d=0.68;
k_p=0.2;
T=tau_i/h_i;
u_i=( (1-T)*x(3) );
%x(2)=0.5*t;
% CACC equations
x_1dot=0; %q(position)dot
x_1ddot=x(3); %v(velocity)dot
x_1dddot=(-(1/tau_i)*x(3))+((1/tau_i)*u_i); %a(accelration)dot
%% second vehicle
tau_i=0.2;
h_i=0.5;
T=tau_i/h_i;
e_i1_1=x(1)-( x(4)+(h_i*x(5)) );
e_i2_1=x(2)-( x(5)+(h_i*x(6)) );
%e_i0_1=int(e_i1_1,t);
u_i=T*[k_p k_d]*[e_i1_1;e_i2_1]+ ( (1-T)*x(6) ) + ( T*x(3) );
% CACC equations
x_2dot=x(5); %q(position)dot
x_2ddot=x(6); %v(velocity)dot
x_2dddot=(-(1/tau_i)*x(6))+((1/tau_i)*u_i); %a(accelration)dot
%% third vehicle
tau_i=0.1;
h_i=0.5;
T=tau_i/h_i;
e_i1_2=x(4)-( x(7)+(h_i*x(8)) );
e_i2_2=x(5)-( x(8)+(h_i*x(9)) );
u_i=T*[k_p k_d]*[e_i1_2;e_i2_2]+ ( (1-T)*x(9) )+ ( T*x(6) );
% CACC equations
x_3dot=x(8); %q(position)dot
x_3ddot=x(9); %v(velocity)dot
x_3dddot=(-(1/tau_i)*x(9))+((1/tau_i)*u_i); %a(accelration)dot
y=[x_1dot; %position of 1st
x_1ddot; %velocity of 1st
x_1dddot; %acceleration of 1st
x_2dot; %position of 2nd
x_2ddot; %velocity of 2nd
x_2dddot; %acceleration of 2nd
x_3dot; %position of 3rd
x_3ddot; %velocity of 3rd
x_3dddot]; %acceleration of 3rd
Concating the time vector made no difference. If you run the whole code, you will see an overshoot which should not be there. Also the system settles in less than a second which is not the case here.

### Categories

Find more on Matrix Indexing in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!