Giving a rectangular pulse in ode45

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
%%

%paramters of leader vehicle
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

%paramters of leader vehicle
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

%paramters of leader vehicle
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?

Answer 1

[Torsten]

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);
%paramters of leader vehicle
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
%paramters of leader vehicle
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
%paramters of leader vehicle
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;
%% leader
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]

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.

Answer 2

[Torsten]

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]

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;

%% leader
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.