ODE45 - Modeling GHR Microspic Car following Model

Hi everyone,

Can someone help with modeling the following Gazis–Herman–Rothery (GHR) model in Matlab:

I have the data for the leader car stored in mytrafficdata5 files as folllows:

I have edited the following code from a similar model (OVM). A difference here is the time lag (tau=1 sec) between the leader and follower vehicles. I was unable to run the code and wonder if someone can help to fix it as per the given model. Thanks

function basic_ghr
% y is velocity
% t is time
% Using interpolant for distance
close all;
load mytrafficdata5;
time_vector = Time_t;
init_speed = 0;
dist_vector = Distance_n;
vel_vector = Velocity_n
% dist_interp = interp1(time_vector, dist_vector,);
% vel_interp = interp1(time_vector, vel_vector,);
F_time = time_vector(end);
[T,Y] = ...
ode45( @(t,y) ...
some_system(t,y,time_vector,dist_vector),[0,F_time],init_speed);
plot(T,Y(:,1),'b-',time_vector,vel_vector,'r--')
legend('By GHR','Real Velocity')
title('Comparisam of velocity with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Velocity (m/sec)')
function dy = some_system(t,y,time_vector,dist_vector, vel_vector)
c = 125;
m = 0.2;
l = 1.6;
%%%%%%%%%%%%%%%%%
% Using interpolant for distance and velocity 
dist_interp = interp1(time_vector, dist_vector, t);
vel_interp = interp1(time_vector, vel_vector, t);
dy = c *((y)^m)/(vel_interp)^l))*(dist_interp);

[Stephan]

Attaching the data would help

[Khalilullah Mayar]

I have attached the data file now. Not quite sure if the interp1 function was the right thing to deal with delta x and delta v quantities in this model.

Answer 1

[Stephan]

This code runs, but still has a problem:

load mytrafficdata5;
basic_ghr(Distance_n,Time_t,Velocity_n)

function basic_ghr(Distance_n,Time_t,Velocity_n)
% y is velocity
% t is time
% Using interpolant for distance
close all;
time_vector = Time_t;
init_speed = 0;
dist_vector = Distance_n;
vel_vector = Velocity_n;
% dist_interp = interp1(time_vector, dist_vector,);
% vel_interp = interp1(time_vector, vel_vector,);
F_time = time_vector(end);
[T,Y] = ode45(@(t,y)some_system(t,y),[0,F_time],init_speed);
plot(T,Y(:,1),'b-',time_vector,vel_vector,'r--')
legend('By GHR','Real Velocity')
title('Comparisam of velocity with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Velocity (m/sec)')
function dy = some_system(t,y)
c = 125;
m = 0.2;
l = 1.6;
%%%%%%%%%%%%%%%%%
% Using interpolant for distance and velocity 
dist_interp = interp1(time_vector, dist_vector, t);
vel_interp = interp1(time_vector, vel_vector, t);
dy = c *(y^m)/(vel_interp)^l*(dist_interp);
end
end

The result of Y(:,1) is a vector of one zero and then all NaN - so you will have to check your equation.

[Khalilullah Mayar]

I have attached my excell worksheet solution to the probelm with relevant plots. Your code produces same result for leader car (given data vectors) but not reacting to the follower car (n+1). Do you think the interp1 function is the right choice for calculating delta x and delta v which are the the difference of distances and velocities one step apart (tau=1sec)? Plus I have an intial value of 0.2 for accelration of the follower car [a (n+1) = 0.2 ] which is not considered in the code.

Note: In excell, I couldn't get any plot for the follower car (n+1) velocity or acceleration with zero intial values. An initial value other than zero, either for the velocity (v n+1=y) or acceleration (a n+1 = dy) produced the intended results.

How could we plot the follower car accleration vs the leader one.

Thanks

load mytrafficdata5;
basic_ghr(Distance_n,Time_t,Velocity_n, Acceleration_n)
function basic_ghr(Distance_n,Time_t,Velocity_n, Acceleration_n)
% y is velocity
% t is time
% Using interpolant for distance
close all;
time_vector = Time_t;
init_speed = 4;
dist_vector = Distance_n;
vel_vector = Velocity_n;
ac_vector = Acceleration_n;
% dist_interp = interp1(time_vector, dist_vector,);
% vel_interp = interp1(time_vector, vel_vector,);
F_time = time_vector(end);
[T,Y] = ode45(@(t,y)some_system(t,y),[0,F_time],init_speed);
figure (1);
plot(T,Y(:,1),'b-',time_vector,vel_vector,'r--')
legend('By GHR','Real Velocity')
title('Comparisam of velocity with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Velocity (m/sec)')
figure (2);
plot(T,z(:,1),'b-',time_vector,ac_vector,'r--')
legend('By GHR','Real Acceleration')
title('Comparisam of Acceleration with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Acceleration (m/sec2)')
function dy = some_system(t,y)
c = 125;
m = 0.2;
l = 1.6;
%%%%%%%%%%%%%%%%%
% Using interpolant for distance and velocity 
dist_interp = interp1(time_vector, dist_vector, t);
vel_interp = interp1(time_vector, vel_vector, t);
ac_interp = interp1(time_vector, ac_vector, t);
dy = c *(y^m)/(vel_interp)^l*(dist_interp);
z=dy;
end
end

[Stephan]

I think your equation is not correct coded:

dy = c *(y^m)/(vel_interp)^l*(dist_interp);

should be:

dy = c *(y^m)/(dist_interp)^l*(vel_interp);

as far as i could find out about this model - but im not familar with this - i saw this in a bachelor thesis. Because of the structure of your equation using zero as initial value leads to a result of all zeros - i suggest to use for example:

init_speed = 0.00001;

These changes will give some kind of solution. Playing around with l and setting:

l = 1.9;

leads to the following result:

I have no idea if it is this what you expect to get - as said before you have to check if your equation is correct and if you use the correct data as input. The code is now free of errors and Matlab does his job...

[Khalilullah Mayar]

Thanks Stephan for picking the error in model equation. I have fixed that, still while the same equations are used in both excell and MATLAB, the plots are not the same. Excell gives an excellet fit while MATLAB gives quite a different result.

load mytrafficdata5;
basic_ghr(Distance_n,Time_t,Velocity_n)

function basic_ghr(Distance_n,Time_t,Velocity_n)
% y is velocity
% t is time
% Using interpolant for distance
close all;
time_vector = Time_t;
init_speed = 0.2;
dist_vector = Distance_n;
vel_vector = Velocity_n;
% dist_interp = interp1(time_vector, dist_vector,);
% vel_interp = interp1(time_vector, vel_vector,);
F_time = time_vector(end);
[T,Y] = ode45(@(t,y)some_system(t,y),[0,F_time],init_speed);
plot(T,Y(:,1),'b-',time_vector,vel_vector,'r--')
legend('By GHR','Real Velocity')
title('Comparisam of velocity with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Velocity (m/sec)')
function dy = some_system(t,y)
c = 125;
m = 0.2;
l = 1.6;
%%%%%%%%%%%%%%%%%
% Using interpolant for distance and velocity 
dist_interp = interp1(time_vector, dist_vector, t);
vel_interp = interp1(time_vector, vel_vector, t);
dy = c *((y^m)/(dist_interp)^l)*(vel_interp);
end
end

Excell file is attached

[Khalilullah Mayar]

Just wondering if putting function seperately in the following format will help:

function dydt = ghr1(t,y)
dydt= zeros(3,1)
dydt(1)= y(2);
dydt(2)= c *((y(2))^m)/(dist_interp)^l)*(vel_interp);
dydt(2)= y(3)
c = 125;
m = 0.2;
l = 1.6;
dist_interp = interp1(time_vector, dist_vector, t);
vel_interp = interp1(time_vector, vel_vector, t);
acc_interp = interp1(time_vector, acc_vector, t);

and Code:

close all;
load mytrafficdata5;
time_vector = Time_t;
dist_vector = Distance_n;
vel_vector = Velocity_n;
acc_vector= Acceleration_n
F_time = time_vector(end);
tspan=[0,F_time];
y0=[0;0;0.2];
[t,y] = ode45(@ghr1,tspan,y0);
figure (1);
plot(t,y(:,1),'b-',time_vector,dist_vector,'r--')
legend('By GHR','Real Distance')
title('Comparisam of distance with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Distance (m)')

figure (2);
plot(t,y(:,2),'b-',time_vector,vel_vector,'r--')
legend('By GHR','Real Velocity')
title('Comparisam of velocity with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Velocity (m/sec)')

figure (3);
plot(t,y(:,3),'b-',time_vector,acc_vector,'r--')
legend('By GHR','Real acceleration')
title('Comparisam of acceleration with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Acceleration (m/sec2)')

[Stephan]

You are using the wrong dataset i guess. I attached a modified .mat-file taking the delta-values for X and V from your Excel file. Those are introduced to the function as S and V. This appears to be more logic to me if i look at your equation:

load mytrafficdata5;
basic_ghr(S,Time_t,V,Distance_n,Velocity_n)
function basic_ghr(S,Time_t,V,Distance_n,Velocity_n)
% y is velocity
% t is time
% Using interpolant for distance
close all;
time_vector = Time_t;
init_speed = 0.2;
dist_vector = S;
vel_vector = V;
% dist_interp = interp1(time_vector, dist_vector,);
% vel_interp = interp1(time_vector, vel_vector,);
F_time = time_vector(end);
[T,Y] = ode45(@(t,y)some_system(t,y),[0,F_time],init_speed);
plot(T,Y(:,1),'b-',time_vector,Velocity_n,'r--')
legend('By GHR','Real Velocity')
title('Comparisam of velocity with time bewteen Real Data Set VS. GHR model')
xlabel('Time')
ylabel('Velocity (m/sec)')
legend('Location','best')
function dy = some_system(t,y)
c = 125;
m = 0.2;
l = 1.6;
%%%%%%%%%%%%%%%%%
% Using interpolant for distance and velocity 
dist_interp = interp1(time_vector, dist_vector, t);
vel_interp = interp1(time_vector, vel_vector, t);
dy = c *(y^m)/(dist_interp)^l*(vel_interp);
end
end

With this values and the initial value of 0.2 m/s it appears to be correct:

I can only repeat myself: It is not a Matlab problem, but a problem of a correct equation and the usage of the correct data...

If my help was useful, dont miss to accept this answer ;-)

[Khalilullah Mayar]

Thanks Stephan, that was helpul.

[Khalilullah Mayar]

Hi Stephan,

Any comment on this would be appreciated:

https://www.mathworks.com/matlabcentral/answers/481828-idm-car-following-microspoic-model