First, the decision variables are defined. They are the real-valued matrices P _c e ℝ^nxv that represent the charging power for each EV in the set V = {1,...,v}, where v is the total amount of vehicles in the fleet considering each time step in the set N = {1,..., n} and n is the number of steps in the calculation period; and P _d e ℝ^nxsv that represents the discharging power for each EV. Additionally, the optimization has other two decision variables with binary values. The first one is the matrix A _v e ℝ^rxv whose elements (i, j) have unitary values if the j-th vehicle performs the i-th travel, which belongs to the set of travels R = {1,..., r}. Hence, if vehicle number 4 performs the first travel, then A _v (1,4) = 1. The second one is the matrix A _c e ℝ^nxv that indicates if a vehicle is charging and contains unitary values in the positions corresponding to the time steps at which a vehicle is charging.

The objective function is the minimization of the EVs operative costs, which include the expenses for purchasing energy to charge the EVs and the batteries wear cost. Additionally, the revenues generated by delivering power from EVs to the electric network are included. This objective function is indicated in Equation (1).

where p _c(i) and p _{d (i)} are the column vectors of the matrices P _c and P_d, respectively; w is the wear cost function, Ep e Rⁿ is the electricity price at each time step, and S(¿) are the column vectors of the state of charge matrix S e ℝ ^{(n+1) xv}, which is also called SOC and is calculated using the Equation (2) for the i-th vehicle at the j-th time step.

P _r e ℝ^nxv is the matrix that contains the power consumption for the travels performed by each EV at each time step, after each travel is assigned one of the EVs. σ is the self-discharging factor of the batteries and E_m, the maximum energy that can be stored in them.

The matrix P _r depends on the travels assignment. Equation (3) indicates how to calculate this variable.

C _t ∈ ℝ^nxr is the matrix that contains the power consumption in each travel, which is the output of the urban traffic microsimulation,

Furthermore, the optimization constraints, are shown in Equations (4)-(12):

P_r, P _c , and P _d must be positive real-valued matrices:

The state of charge of the batteries must be kept between the minimum allowed energy content S_m, and 1.

The total power consumption in the travels at each time step, must be equal to the total power that the EVs spend during their assigned travels each time step:

Travels must be assigned only to one vehicle:

The charging power of the EVs must be kept between 0 and its rated value defined by the manufacturer P_c,m, at the time steps when EVs are charged (indicated by ones in the A _c matrix).

The discharging power of the EVs must be kept between 0 and its rated value defined by the manufacturer P_cm, at time steps when EVs are discharged, i.e. when they are not performing a travel or charging

where U ∈ ℝ ^nxv is a binary matrix calculated with Equation (10). The rows of U correspond to the time steps and its columns to the ID number for each EV in the fleet. This matrix has ones in the time steps when a vehicle is performing a travel, and 0 otherwise.

Being A_t e ℝ^nxr a binary matrix that contains ones at time steps where travel is required, it has as many columns as required travels. This matrix is obtained from the output data of the urban traffic microsimulator and entered as the input of the optimization. Furthermore, the sum of the elements of the binary matrix that indicate the charging status and those that show the running status for a vehicle performing travel must be less than 1, as indicated in Equation (11)

Finally, the power charged or discharged from the EVs in the fleet must be lower than the power capacity of the feeder (P;) to which the charging station is connected:

Batteries wear modeling:

The function w in Equation (1) is calculated in two ways:

a) The wear is calculated with the methodology presented by Xu et al (2018) using the Rain-flow algorithm.

This algorithm identifies the peak points of the charging/discharging profile. Then, from, these peak data, the half and complete equivalent cycles from the profile are calculated, using the algorithm detailed by Blumenthal (1935), and entered as inputs of the wear function that is described in Equation (13)

where the constant C¡, takes into account the equivalent present value of the replacement costs for a defined Capital Recovery Factor, and the constants k₁, k ₂ are defined by the battery manufacturer. The variable A takes a value of 1 for complete cycles, and a value of 0,5 for half cycles.

This study proposes a novel approach to calculate the peak data from the charging/discharging profile in order to compute the Rain-flow algorithm and calculate the sub-gradient of this function for solving the optimization problem with nonlinear methods based on the gradient. The peak data from the charging/discharging profile for the i-th EV and the j-th time step can be extracted using Equation (14)

The logic of the Rain-flow (RF) can be implemented using conditional operators, to identify the halfand complete cycles from the peak data. The objective function is not differentiable at some points, the cycle junction points. Therefore, the subgradient function is defined in Equation (15) based on the approach presented by Xu et al. (2018).

where the points j ₁ and j ₂ are the cycle junction points between the point j.

However, considering a fleet of more than 1 EV, the nonlinear optimization problem described in Equations (1)-(15) becomes a mixed-integer nonlinear problem. Hence, in order to achieve results with a lower complexity and computation time, the wear calculation is reformulated.

b) The wear is calculated with the Achievable Cycle Count (ACC), based on the approaches presented by Choi and Kim (2016) and Han et al. (2014).

The wear can also be calculated using the cycle life curves of the battery, given by the manufacturer. This cycle life is determined by how deeply the battery is used, i.e. by the depth of the discharge (D = 1 - S). Hence, a relationship between D and the life cycle must be established. In the literature, this relationship is commonly modeled by a polynomial function as in the studies of Choi and Kim (2016) and Han et al. (2014), where the cycle life N _cycle is determined with the polynomial function of Equation (16).

where constants a and b are estimated from the manufacturer data. The wear cost is calculated as Choi and Kim (2016):

w _d (s) is the wear-out density function, which is calculated as Equation (18) shows:

In Equation (18), the constant Cb takes into account the equivalent present value of the replacement costs for a defined Capital Recovery Factor. Hence, the complete expression for the battery wear is given in Equation (19).

A linear approximation for polynomial expressions is proposed on the integral calculated in the Equation (19). Thus, the approximated wear is:

where p ₁ is the slope of the linear function approximating the polynomial expression of Equation (19).

Once the wear function of Equation (20) is considered, the optimization problem becomes a mixed-integer convex problem, whose solution is presented in the Results and analysis section.

The first study case considers one electric vehicle, a Renault Twizy. The travels performed by the EV in the evaluation period are predicted from the behavior of co-housing users described by Semanjski and Gautama (2016). The operating hours of this vehicle are between 5:00 and 16:00. The energy consumption in each travel is calculated with the daily travels, which are generated randomly. The random generation of daily travels uses a normal distribution for the distance of each travel and the number of travels per day, taking into account the next parameters: a daily average traveled distance of 8,4 km, and a daily average of travels performed of 3. Three days were simulated. The parameters considered for the wear function were taken from the study of Correa-Florez et al. (2018) and are presented in Table 1.

Table 1

Source: (Correa-Florez et al., 2018)

Hence, with the energy consumption data already calculated, this study case only requires the execution of the steps indicated in the Control Algorithm block in Figure 1. The optimization, including the wear cost, was solved using Matlab 2015b. Three optimization cases were taken into account: the first case without considering the wear cost, the second case using the Rain-flow algorithm to calculate the wear cost, and the third case that includes the wear cost and uses the ACC function to calculate it. The solutions of the first and third optimization problems were found with the CVX optimization tool and the Gurobi solver, as both are convex problems. Gurobi uses the Branch and Cut algorithm. The solution of the second case was found using the fmincon function of the MATLAB Optimization Toolbox with the Active-Set optimization algorithm.

The study case used to apply the algorithm of Figure 1 includes the traffic network, and the passengers demand described by Ruiz et al. (2018), which presents a study case applicable for the city of Medellfn, Colombia. The urban traffic microsimulation to calculate the required travels and power consumption in them was performed using the software SUMO. The electric vehicles considered are BRTs Type 1, as described in Ruiz et al. (2018), whose energy capacity is 324 kWh. The data of the required travels and power consumption in them were taken from Ruiz et al. (2018), taking into account only the travels performed by BRTs Type 1. The electricity price time series and the power limit for the feeder of the charging station were also taken from the aforementioned study.

The parameters considered for the evaluation of the batteries aging are presented in Table 1 with a replacement cost of batteries of USD 5 400. This wear calculation was done only with the ACC approach, using Equation (20). Time steps of 10 minutes were considered, a simulation period of 24 hours, a fleet size of 5 BRTs, and 21 travels required to cover the passenger demand. The software used to compute the optimization was Matlab 2015b with the Gurobi solver.