Case study specification (Student version)

Case study specification (Student version)#

Written by Marc Budinger, INSA Toulouse, France

We will now develop tools allowing:

generate mission profiles corresponding to a new need for an urban transport line
optimize certain mission profiles choices to minimize consumption

In order to have flexibility in the definition of the case studies, an object-oriented programming approach will be used.

Electric vehicles#

Different types of urban electric vehicles exist (tram, trolley bus, electric bus), we can define a general vehicle class wich can represent the multiple possible supports.

class vehicle:
    def __init__(self, MaxSpeed, MaxAcc,NumberPass,Weight,Crr,Cd,FrontArea):
        self.MaxSpeed = MaxSpeed # [m/s] maximal speed
        self.MaxAcc = MaxAcc # [m/s²] maximal acceleration
        self.NumberPass = NumberPass # [-] passenger capacity (3p/m²)
        self.Weight = Weight # [kg] effective weight (with passenger 62kg/p)
        self.Crr = Crr # [-] rolling resistance coefficient
        self.Cd = Cd # [-] drag coefficient
        self.FrontArea = FrontArea # [m²] frontal area
        self.Fmax = Weight*MaxAcc # [N] max traction force
        self.Pmax = self.Fmax*self.MaxSpeed # [W] max corner power

Here you will find some useful vehicle characteristics for the rest of the study.

Attribute	Tram	Trolleybus	Bus	Car (Tesla 3)
Maximum operational speed	60 km/h	60 km/h	50 km/h	100 km/h (1st gear)
Acceleration and braking	1.2 m/s²	1.2 m/s²	1.2 m/s²	4.6 m/s²
Passenger capacity @ 3p/m²	220	138	95	4
Vehicle effective weight @ 62 kg/p	57049 kg	27656 kg	19500 kg	1800 kg
Rolling resistance coefficient	0.006	0.015	0.015	0.015
Aerodynamic drag coefficient	0.6	0.6	0.6	0.23
Frontal area	8.5 m²	8.5 m²	8.5 m²	2.22 m²

Exercice: Define 4 instances of the class vehicle corresponding to the vehicles above.

Generation and mission profile simulation#

The objective now is to simulate the dynamic evolution of the main mechanical quantities (position, speed, acceleration, tractive effort, driving power) on typical sections of urban routes.

When a vehicle travels from point A (Starting point) to point B (Stopping point ) it generally passes through four stages namely:

Accelerating mode
Constant speed mode
Coasting/freewheel mode
Braking mode

The control points for switching from one mode to another are defined in the following code using 3 parameters \(k_1, k_2, k_3\) as shown in the figure below.

The following code implements a control logic comprising the 4 phases defined previously and a numerical integration of the differential equation corresponding to the longitudinal dynamics of the vehicle.

Question: What evolution should be added to take into account variations in altitude?

import numpy as np
from scipy import signal
from scipy.optimize import fmin_slsqp
from scipy.optimize import differential_evolution
import matplotlib.pyplot as plt
from scipy.integrate import odeint

class SimulSection:
    def __init__(self,Vehicle,Distance, MeanSpeed,BrakeRatioMax,dt):
        # Parameters definition
        self.Vehicle = Vehicle # vehicle parameters
        self.Distance = Distance # [m] distance to travel
        self.MeanSpeed = MeanSpeed # [m/s] Mean Speed < MaxSpeep
        self.TravelTime = self.Distance / self.MeanSpeed  # [s] Travel time
        self.BrakeRatioMax = BrakeRatioMax # [-], inferieur à 1, ratio de puissance de freinage / Puisssance d'acceleration
        self.rho = 1.25 # [kg/m3] air density
        self.g = 9.81 # [m/s²] gravity
        self.dt = dt # Time step for numerical solver
        
        # Tests sur les lois de mouvement
        # https://fr.wikipedia.org/wiki/Loi_de_mouvement
        
        if (MeanSpeed>=Vehicle.MaxSpeed):
           print("Mean Speed is to high : Mean Speed > Vehicle Max Speed") 
        
        # profil triangulaire : minimise l'acceleration
        #     amax = 4⋅xf/T2 ;    vmax = 2⋅xf/T ; amax = vmax / (T/2) ; vmean = vmax / 2
        #     amax = 4*vmean^2 / xf
        
        # Calcul de l'acceleration possible a pleine vitesse
        AmaxFullSpeed=(self.Vehicle.Fmax - 1/2*self.rho*self.Vehicle.Cd*self.Vehicle.FrontArea*self.Vehicle.MaxSpeed**2 - self.Vehicle.Weight*self.g*self.Vehicle.Crr)/self.Vehicle.Weight
        if (AmaxFullSpeed < 4*self.MeanSpeed**2 / self.TravelTime):
            print("Vehicle Max Acceleration (depending of max force) is too low or travel distance too small or mean speed too small")
            
            
    # dynamic model for acceleration, coasting and braking
    def model(self, y,F):
        # state
        x, dxdt = y
        # System: acceleration calculation
        dxdt2 = (F-self.Vehicle.Crr*self.Vehicle.Weight*self.g*np.sign(dxdt)
                 -1/2*self.rho*self.Vehicle.Cd*self.Vehicle.FrontArea*dxdt**2**np.sign(dxdt))/self.Vehicle.Weight
        dydt = [dxdt, dxdt2]
        return dydt 

  
    # https://perso.crans.org/besson/publis/notebooks/Runge-Kutta_methods_for_ODE_integration_in_Python.html
    # solver numerique
        
    def solver(self, x):
        k1, k2, k3 = x
        dt=self.dt # [s] pas de temps pour l'integration
        t=0
        self.tsection = [0]
        self.xsection = [0]
        self.vsection = [0]
        self.asection = [0] 
        self.psection = [0]
        self.Fsection = [0]
        self.NRJsection = [0]
        NRJ=0
        
        y0 = np.array([0, 0]) 
        y= y0
        
      
        while (y[1]>=0):
            # Traction/Braking force
            if (y[0]<k1*k2*k3*self.Distance):
                F=self.Vehicle.Fmax
                dydt = self.model(y, F)
            elif (y[0]<k3*k2*self.Distance):
                F=self.Vehicle.Crr*self.Vehicle.Weight*self.g+1/2*self.rho*self.Vehicle.Cd*self.Vehicle.FrontArea*y[1]**2
                dydt = [y[1], 0]
            elif (y[0]<k3*self.Distance):
                F=0
                dydt = self.model(y, F)
            else:
                F=-self.BrakeRatioMax*self.Vehicle.Fmax
                dydt = self.model(y, F)
                
            self.Fsection = self.Fsection + [F]    
                
            # Euler integration de la postion y[0] et de la vitesse y[1]
            y = y + dt * np.array(dydt)

            t = t + dt
            NRJ = NRJ+y[1]*F*dt

           
            self.tsection = self.tsection + [t]
            self.xsection = self.xsection + [y[0]]
            self.vsection = self.vsection + [y[1]]
            self.psection = self.psection + [y[1]*F]
            self.asection = self.asection + [dydt[1]] 
            self.NRJsection = self.NRJsection + [NRJ]
            
    
    def plot(self):
        fig, axs = plt.subplots(3,2)
        axs[0,0].plot(self.tsection,self.xsection,'b-',label='Simulation')
        axs[0,0].plot(self.tsection,self.Distance*np.ones(len(self.tsection)),'g-',label='Specification')
        axs[0,0].set_ylabel("Position (m)")
        axs[0,0].legend()
        axs[0,0].grid()
        axs[1,0].plot(self.tsection,self.vsection,'b-', label='Simulation')
        axs[1,0].plot(self.tsection,self.MeanSpeed*np.ones(len(self.tsection)),'g-',label='Specification')
        axs[1,0].plot(self.tsection,np.mean(self.vsection)*np.ones(len(self.tsection)),'r--',label='Mean')
        axs[1,0].set_ylabel("Speed (m/s)")
        #axs[1].sharex(axs[0])
        axs[1,0].grid()    
        axs[1,0].set_xlabel('Time (s)')    
        axs[0,1].plot(self.tsection,self.Fsection,'b-')
        axs[0,1].set_ylabel("Force (N)")
        axs[0,1].grid()
        
        axs[2,0].plot(self.tsection,self.asection,'b-')
        axs[2,0].set_ylabel("Acceleration (m/s²)")
        axs[2,0].grid()
        axs[2,0].set_xlabel('Time (s)')    
        
        axs[1,1].plot(self.tsection,np.array(self.psection)*1e-3,'b-')
        axs[1,1].set_ylabel("Power (kW)")
        axs[1,1].grid()
        axs[1,1].set_xlabel('Time (s)')   
        
        axs[2,1].plot(self.tsection,np.array(self.NRJsection)*1e-3,'b-')
        axs[2,1].set_ylabel("Energy (kJ)")
        axs[2,1].grid()
        axs[2,1].set_xlabel('Time (s)') 
        fig.tight_layout()

    def ConsumptionPerPax(self,x):
        self.solver(x)
        print("Consumption per passenger: %.2f kJ/(Pax.km)"%(self.NRJsection[-1]/self.Vehicle.NumberPass/self.Distance))
        return (self.NRJsection[-1]/self.Vehicle.NumberPass/self.Distance)
    
    def MaxEnergyPoint(self,x):
        self.solver(x)
        print("Max energy discharge: %.0f kJ"%(max(self.NRJsection)*1e-3))
        return max(self.NRJsection)*1e-3
    

Example of use of the previous code:

# Definition of the section under study
Trajet_Tram = SimulSection(Vehicle=Tram, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=1.0,dt=1)
#Trajet_Car = SimulSection(Vehicle=Car, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=1.0,dt=1)
print('Duration of traject for given mean speed: ' + str(670/(26*1e3/3600)) +'s')

# Profile integration, plotting and consummption per passenger evaluation
X=[0.036125972011676784, 1, 0.9638740279883232]
Trajet_Tram.solver(X)
Trajet_Tram.plot()
Trajet_Tram.ConsumptionPerPax(X)
Trajet_Tram.MaxEnergyPoint(X)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 2
      1 # Definition of the section under study
----> 2 Trajet_Tram = SimulSection(Vehicle=Tram, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=1.0,dt=1)
      3 #Trajet_Car = SimulSection(Vehicle=Car, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=1.0,dt=1)
      4 print('Duration of traject for given mean speed: ' + str(670/(26*1e3/3600)) +'s')

NameError: name 'Tram' is not defined

Two extrema cases are :

a triangular profile (acceleration/deceleration, without constant speed & freewheel)
a trapezoidal profile (without freewheel)

Question 1: Give the relationship between acceleration, average speed and distance traveled for triangular profile if we assume the acceleration and braking phases are perfectly linear and of the same slope (i.e. resistive force negligible): does it match a warning message of the SimulSection class?

Question 2: Give the relationship between k1 (k2=1), average speed and distance traveled if we assume the acceleration and braking phases are perfectly linear and of the same slope.

# Calculus of lower acceleration
a_min = 4*(26*1e3/3600)**2/670
print('Min acceleration is: ' + str(a_min) +'m/s²')
Tram_low_acc = vehicle(MaxSpeed = 60*1e3/3600, MaxAcc=a_min,NumberPass=220,Weight=57049,Crr=0.006,Cd=0.6,FrontArea=8.5)

# Triangular profile test
X=[1, 1, 0.6] #if resistive force not negligieble in front of propulsive force not k1>0.5
Trajet_Tram_low_acc = SimulSection(Vehicle=Tram_low_acc, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=1.0,dt=1)
Trajet_Tram_low_acc.solver(X)
Trajet_Tram_low_acc.plot()
Trajet_Tram_low_acc.ConsumptionPerPax(X)
Trajet_Tram_low_acc.MaxEnergyPoint(X)

Min acceleration is: 0.3114059332964806m/s²
Vehicle Max Acceleration (depending of max force) is too low or travel distance too small or mean speed too small

Consumption per passenger: 18.74 kJ/(Pax.km)
Max energy discharge: 7420 kJ

7419.990488324914

# Parameters calculation
# on pose x = sqrt(k1)
from scipy.optimize import root

a = Trajet_Tram.MeanSpeed/(2*Trajet_Tram.Vehicle.MaxAcc*Trajet_Tram.Distance)**0.5

f = lambda x: a+(x*(1-x))**1.5-(x*(1-x))**0.5

X=np.linspace(0,0.5,100)

plt.plot(X,f(X),'b-')
plt.grid()
plt.show()

sol=root(f, 0.1)

k1=float(sol.x)
print(k1)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 5
      1 # Parameters calculation
      2 # on pose x = sqrt(k1)
      3 from scipy.optimize import root
----> 5 a = Trajet_Tram.MeanSpeed/(2*Trajet_Tram.Vehicle.MaxAcc*Trajet_Tram.Distance)**0.5
      7 f = lambda x: a+(x*(1-x))**1.5-(x*(1-x))**0.5
      9 X=np.linspace(0,0.5,100)

NameError: name 'Trajet_Tram' is not defined

k3=(1-k1)/k1

X=[k1,1,k3]
print(X)
Trajet_Tram.solver(X)
Trajet_Tram.plot()
Trajet_Tram.ConsumptionPerPax(X)
Trajet_Tram.MaxEnergyPoint(X)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 k3=(1-k1)/k1
      3 X=[k1,1,k3]
      4 print(X)

NameError: name 'k1' is not defined

Questions: Is this type of profile the one that minimizes the energy consumed ? What rules do you recommend to minimize the energy consumed?

Optimizing power consumption using speed profile#

In a more general case and contrary to the previous hypothesis, the braking power is not necessarily identical to that possible during acceleration. The following article shows that it is possible to optimize the consumption of a tram:

Tian, Z., Zhao, N., Hillmansen, S., Roberts, C., Dowens, T., & Kerr, C. (2019). SmartDrive: Traction energy optimization and applications in rail systems. IEEE Transactions on Intelligent Transportation Systems, 20(7), 2764-2773.

The previous class is now extended to allow optimization of variables \(k_1, k_2, k_3\) in order to minimize energy consumption (even with braking power lower than acceleration power).

class OptimSection(SimulSection):
    
    # Gradient optimization
    def objectifConso(self, x):
        self.solver(x)
        NRJmin = self.xsection[-1]*self.Vehicle.Weight*self.g*self.Vehicle.Crr
        return (self.NRJsection[-1]/NRJmin)   
    
    def contraintes(self, x):
        self.solver(x)
        TimeTrajet = self.Distance / self.MeanSpeed  
        PowerAccBrakeRatio = abs(min(self.psection))/max(self.psection)
        print('Brake power/Acceleration power ratio is: '+ str(PowerAccBrakeRatio))
        # contraintes en distance, vitesse moyenne (ou duree de trajet, vitesse au freinage), ratio puissance freinage / acc 
        return [(self.xsection[-1]-self.Distance)/self.Distance, (TimeTrajet-self.tsection[-1])/TimeTrajet, 
                self.BrakeRatioMax-PowerAccBrakeRatio ]    
      
    def optimizeConso(self,x0):
        Xbound= [(0.0, 0.3), (0.1, 1), (0.5,.99)]
        Xopt=fmin_slsqp(func=self.objectifConso, x0=x0, f_ieqcons = self.contraintes, bounds=Xbound, epsilon=1e-2)      
        return Xopt
    

    # Differential evolution optimization
    def objectifGConso(self, x):
        self.solver(x)
        NRJmin = self.xsection[-1]*self.Vehicle.Weight*self.g*self.Vehicle.Crr
        pen=0
        TimeTrajet = self.Distance / self.MeanSpeed  
        PowerAccBrakeRatio = abs(min(self.psection))/max(self.psection)
        VecC=[(self.xsection[-1]-self.Distance)/self.Distance, (TimeTrajet-self.tsection[-1])/TimeTrajet, self.BrakeRatioMax-PowerAccBrakeRatio ]      
        for C in VecC:
            if (C<0):
                pen=pen-1e2*C
        return (self.NRJsection[-1]/NRJmin+pen)
    
    def optimizeGConso(self,x0):
        Xbound= [(0.0, 0.3), (0.1, 1), (0.5,0.99)]
        res=differential_evolution(func=self.objectifGConso, x0=x0, bounds=Xbound)
        print(res)
        return res.x
    
    # Differential evolution optimization
    def objectifGEmax(self, x):
        self.solver(x)
        NRJmin = self.xsection[-1]*self.Vehicle.Weight*self.g*self.Vehicle.Crr
        pen=0
        TimeTrajet = self.Distance / self.MeanSpeed  
        PowerAccBrakeRatio = abs(min(self.psection))/max(self.psection)
        VecC=[(self.xsection[-1]-self.Distance)/self.Distance, (TimeTrajet-self.tsection[-1])/TimeTrajet, self.BrakeRatioMax-PowerAccBrakeRatio ]      
        for C in VecC:
            if (C<0):
                pen=pen-1e2*C
        return (max(self.NRJsection)/NRJmin+pen)
    
    def optimizeGEmax(self,x0):
        Xbound= [(0.0, 0.3), (0.1, 1), (0.5,0.99)]
        res=differential_evolution(func=self.objectifGEmax, x0=x0, bounds=Xbound)
        print(res)
        return res.x
    

The usage of the class is quite simular that the previous one, except an optimizewhich enable to optimize the \(k_1, k_2, k_3\) set of profile parameter:

# Definition of the section under study
Trajet = OptimSection(Vehicle=Tram, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=0.6,dt=0.25)

# Initial variables vector
X=[0.036125972011676784, 1, 0.9638740279883232]

# Optimization of the consumption
Xopt1=Trajet.optimizeGConso(X) # Global optimization with differential evolution algorithm
#Xopt1=Trajet.optimizeConso(X) # Opitmization with gradient algorithm

# print and plot results
print("Optimal vector :", Xopt1)
print("Constraints vector :", Trajet.contraintes(Xopt1))

Trajet.solver(Xopt1)
Trajet.plot()
Trajet.ConsumptionPerPax(Xopt1)
Trajet.MaxEnergyPoint(Xopt1)

# Optimization of the max energy point
print('\n######################################\n')
Xopt2=Trajet.optimizeGEmax(X) # Global optimization with differential evolution algorithm

# print and plot results
print("Optimal vector :", Xopt2)
print("Constraints vector :", Trajet.contraintes(Xopt2))

Trajet.solver(Xopt2)
Trajet.plot()
Trajet.ConsumptionPerPax(Xopt2)
Trajet.MaxEnergyPoint(Xopt2)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 2
      1 # Definition of the section under study
----> 2 Trajet = OptimSection(Vehicle=Tram, Distance=670, MeanSpeed=26*1e3/3600,BrakeRatioMax=0.6,dt=0.25)
      4 # Initial variables vector
      5 X=[0.036125972011676784, 1, 0.9638740279883232]

NameError: name 'Tram' is not defined

Question: Compare this optimized profile to a profile with only 3 segments: maximum acceleration, constant speed, maximum deceleration. Evaluate effect of specifications (distance, braking power/acceleration power ratio).

Different vehicules comparison#

Exercice: For the same travel (distance, main speed) compare the consumption of different vehicles (tramway, tram, bus, car). The function ConsumptionPerPax return and print the consumption per km and per passenger. Compare the results with this publication.

print("Calculation of energy consumption of different vehicles: ")
print("")
print("Tramway :")

print("----")

print("Trolleybus :")

print("----")

print("Bus :")

print("----")

print("Car :")

Calculation of energy consumption of different vehicles: 

Tramway :
----
Trolleybus :
----
Bus :
----
Car :

Homework#

Read the following paper:

Jaafar, A., Sareni, B., Roboam, X., & Thiounn-Guermeur, M. (2010, September). Sizing of a hybrid locomotive based on accumulators and ultracapacitors. In 2010 IEEE Vehicle Power and Propulsion Conference (pp. 1-6). IEEE. [pdf]