The shark has a new position *Yk*+1 *i* due to moving forward, and *Yk*+1 *i* is determined by the previous speed and position, which is expressed as

$$Y\_i^{k+1} = X\_i^k + V\_i^k \cdot \Delta t\_k \tag{15}$$

where Δ*tk* represents the time interval of the *k* th iteration. In addition to moving forward, sharks usually rotate along their path to look for stronger odor particles and improve their direction of movement, which is a real way of moving.

The rotating shark moves in a closed interval which is not necessarily a circle. From the point of view of optimization, sharks implement local search at each stage to find better candidate solutions. The search formula for this position is as follows,

$$Z\_i^{k+1,m} = Y\_i^{k+1} + R\mathbf{3} \cdot Y\_i^{k+1} \tag{16}$$

where *m* = (1, 2, ···, *M*) presents the number of points at each stage of the location search; *R*3 is the random number between [−1, 1]. If the shark finds a stronger scent point in the rotation, it moves toward the point and continues the search path. The location search formula is described as follows,

$$X\_{i}^{k+1} = \arg\max\left\{ OF(Y\_{i}^{k+1}), OF(Z\_{i}^{k+1,1}), \\ \cdots, OF(Z\_{i}^{k+1,M}) \right\} \tag{17}$$

As can be seen from the above formula, *Yk*+1 *i* is obtained from the linear movement and *Zk*+1,*M i* is obtained from the rotation movement. Sharks will choose the candidate solution with higher evaluation index value as shark's next location *Xk*+1 *i* .

#### *3.2. Preference Information and Composite Angle Cosine*

To improve the optimization effect, the preference phenomenon should be considered, and the data or information used to quantify the impact of preferences on evaluations is called preference information [38]. The preference vector angle is a kind of preference information, which is used to reflect the degree of preference between the solution vector and the target demand vector. The schematic diagram of the preference phenomenon classification and the preference angle is shown in Figure 2.

**Figure 2.** Schematic diagram of the preference phenomenon classification and the preference angle.

As can be seen from Figure 2, according to the preference phenomenon, preference event Ω is divided into three categories: (*ClassI*, *ClassI I*, and *ClassI I I*), and the corresponding preference angles are ( *PI*, *PI I*, and *PIII*).

If only the numerical angle cosine is used as evaluation index in optimization, the preference phenomenon will be ignored, and it is easy to lead the evaluation result unreasonable. To take account of the numerical difference and preference difference between the solution vector and the target demand vector, a new evaluation index is proposed in this paper, that is, the compound angle cosine. The compound angle cosine is the cosine of the sum of the numerical angle and the preference angle. The schematic diagram of the compound angle cosine is shown in Figure 3.

**Figure 3.** The schematic diagram of the compound angle cosine.

In Figure 3, the angle *N* represents the numerical angle; it is used to reflect the numerical calculated result between solution vector and the target demand vector. The angle *N* represents the preference angle; it is used to reflect the preference value between solution vector and the target demand vector; the angle *C* represents the composite angle, *C* = *P* + *N*. The composite angle cosine is used as the evaluation index for the optimization algorithm proposed in this paper. The specific calculation formula of the composite angle cosine is as follows.

$$\begin{aligned} \angle \mathcal{C} &= \angle P + \angle \mathcal{C} \\ \cos \mathcal{C} &= \cos P \cos N - \sin P \sin N \end{aligned} \tag{18}$$

In some optimization problems, there are several the preference events. The specific calculation formula of the composite angle is as follows,

$$
\angle P = \sum\_{ip=1}^{np} \angle P\_{ip} \tag{19}
$$

where *ip* represents the preference event index, *np* represents the preference event number, and *Pip* represents the preference angle of the *i* th preference event.

There are three preference events in the automatic train operation velocity trajectory optimization with the multi-objective comprehensive performance index *F*(*u*) ,(*F*(*u*)=(*KE*(*u*), *KJerk*(*u*), *KT*(*u*))). The specific preference categories categories circumstances in automatic train operation velocity trajectory optimization are shown in Table 1. Where, "Prefect", "Qualified", and " Unqualified" represents the valuation level of preference event according to preference phenomenon, the preference angle for "Prefect", "Qualified", and "Unqualified" are 0, *π*/24, *π*/12, respectively; *E*1, *E*2, and *E*3 represent the boundary values of valuation level for energy consumption, which are decided by train parameters, line conditions, prescribed running time, and running distance; the boundary values of valuation level for "Comfort level" and "Time error" are decided by train operation regulation.

**Table 1.** The preference categories circumstances in automatic train operation velocity trajectory optimization.

If a certain performance index of the solution vector *T* cannot reach the valuation level "Unqualified", solution vector *T* is impermissible. At present, only valuation level "Prefect" and "Qualified" are permitted for most of automatic train operation scenarios.

#### *3.3. Fusion Distance*

Mahalanobis distance is the data covariance distance defined by Mahalanobis, which can accurately calculate the covariance distance between two samples. The formula of Mahalanobis distance between the sample *X* to be examined and the basic space set *Y* is expressed as

$$\begin{aligned} \Sigma &= \operatorname{Cov}(\mathbf{X}, \mathbf{Y}) = E\left[ (X - E\left(\mathbf{X}\right)) \left( \mathbf{Y} - E\left(\mathbf{Y}\right) \right) \right] \\ &= \begin{bmatrix} \operatorname{Cov}\left(\mathbf{x}\_{1}, y\_{1}\right) & \operatorname{Cov}\left(\mathbf{x}\_{1}, y\_{2}\right) \cdots \operatorname{Cov}\left(\mathbf{x}\_{1}, y\_{j}\right) \\ \operatorname{Cov}\left(\mathbf{x}\_{2}, y\_{1}\right) & \operatorname{Cov}\left(\mathbf{x}\_{2}, y\_{2}\right) \cdots \operatorname{Cov}\left(\mathbf{x}\_{2}, y\_{j}\right) \\ \vdots & \vdots & \ddots & \vdots \\ \operatorname{Cov}\left(\mathbf{x}\_{j}, y\_{1}\right) & \operatorname{Cov}\left(\mathbf{x}\_{j}, y\_{2}\right) \cdots \operatorname{Cov}\left(\mathbf{x}\_{j}, y\_{j}\right) \end{bmatrix} \end{aligned} \tag{20}$$

where Σ is the expected matrix of the covariance matrix for the basic space set *Y*.

The fusion distance is the combination of Mahalanobis distance and Euclidean distance, taking into account the independence and relevance of the characteristic variables, which can

effectively improve the accuracy of distance calculation [39]. The specific formula for calculating fusion distance is expressed as

$$\begin{cases} \begin{aligned} d\_{\text{Mix}} &= \omega \times MD(\mathbf{X}, Y) + (1 - \omega) \times ED(\mathbf{X}, Y) \\ \mathbf{C}\_{Y} &= \begin{bmatrix} \rho\_{Y\_1 Y\_1} & \rho\_{Y\_1 Y\_2} \cdots \rho\_{Y\_1 Y\_n} \\ \rho\_{Y\_2 Y\_1} & \rho\_{Y\_2 Y\_2} \cdots \rho\_{Y\_2 Y\_n} \\ \vdots & \vdots & \ddots & \vdots \\ \rho\_{Y\_n Y\_1} & \rho\_{Y\_n Y\_2} \cdots \rho\_{Y\_n Y\_n} \end{bmatrix} \\ \omega &= \sqrt{1 - |\mathbf{C}\_Y|} \end{aligned} \tag{21}$$

where *dMix* represents the fusion distance, *MD* represents the Mahalanobis distance, *CY* represents the correlation coefficient matrix for the sample set *Y*, *n* is the size of sample set *Y*, *Yi* (*i* = 1, ··· , *n*) is the corresponding elements of sample set *Y*, and *ρ* is the correlation coefficient. The fusion distance is fused by the weight *ω* with the relevant information, and the Euclidean distance is fused by 1 − *ω* .

#### *3.4. Particle Swarm Optimization*

Particle swarm optimization (PSO) is an optimization algorithm proposed by American psychologist Kennedy and electrical engineer Eberhart in 1995. The update formula of the velocity and position of the particle in dimensional space is described as follows,

$$\begin{cases} \boldsymbol{v}\_{ip,t+1}^{d} = \boldsymbol{\omega} \times \boldsymbol{v}\_{ip,t}^{d} + c1 \times rand \times (p\_{ip,t}^{d} - \boldsymbol{\chi}\_{ip,t}^{d})\\ \boldsymbol{+} \, \boldsymbol{c2} \times rand \times (p\_t^{\mathcal{S}} - \boldsymbol{\chi}\_{ip,t}^{d})\\ \boldsymbol{\omega}\_{ip,t+1}^{d} = \boldsymbol{\omega}\_{ip,t}^{d} + \boldsymbol{v}\_{ip,t+1}^{d} \end{cases} \tag{22}$$

where *ip* ∈ [1, 2, ..., *N*] is *ip* th particle of the particle population; *N* represents the size of particle population; *d* ∈ [1, 2, ..., *D*] is the *d* th dimension of the particle; *D* represents the number of dimension; *t* ∈ [1, 2, ..., *T*] is the *t* th iteration; *T* represents the number of iteration; *c*1 and *c*2 represents the acceleration constants; *rand* is the random real number of the interval (0,1); *ω* is the weight coefficient, which is used to balance the degree of global search and local search; the position vector is represented as → *Xip* = (*xip*,1, *xip*,2, ··· , *xip*,*d*, ··· , *xip*,*D*); the velocity vector is expressed as → *Vip* = (*vip*,1, *vip*,2, ··· , *vip*,*d*, ··· , *vip*,*D*); the optimal location of the particles' individual is recorded as → *Pip*,*t* = (*p*1 *ip*,*t* , *p*2 *ip*,*t* , ··· , *pd ip*,*t* , ··· , *pD ip*,*t* ); the best position of the swarm is denoted as → *Pg*,*t* = (*p*1 *g*,*t* , *p*2 *g*,*t* , ··· , *pd g*,*t* , ··· , *pD g*,*t* ).

#### *3.5. Dual-Population Evolution Mechanism*

Due to the limitations of both the evolutionary environment and the initial population, the problem of slow evolution and even stagnant evolution will appear during late iteration [40]. In the long process of iteration, the optimal individual will dominate all the population to some extent, making it difficult to converge globally. The proposal of dual-population strategy makes the optimal individuals of two populations exchange with each other, and the long-term dominance of the optimal individual in the original population is easily lost due to the change of the population environment [41,42].

To improve the optimization performance of SSO algorithm, an improved strategy combining genetic algorithm, particle swarm algorithm, and SSO algorithm based on dual-population Evolution Mechanism is proposed in this paper. In the process of iteration, the SSO algorithm brings each individual of the shark population into the optimal position rapidly. Particle swarm algorithm has the same defect of local convergence as SSO algorithm due to its fixed foraging behavior. At the same time, the crossover and mutation of genetic algorithm can prevent the SSO algorithm based on evolutionary from stalling immediately when it falls into the local optimum, and it can cause certain disturbance and mutation, which can help the SSO algorithm to jump out of the local optimum dilemma. In addition, the composite angle cosine is used as the evaluation index and the parallel evolutionary mechanism of dual-population is adopted to prevent the population from being dominated by extreme individuals. The flowchart of improved shark smell optimization algorithm proposed in this paper is shown in Figure 4.

**Figure 4.** The flowchart of improved shark smell optimization algorithm proposed in this paper.

As can be seen form Figure 4, the dual-population strategy uses two populations (shark population *S*1 and particle swarm *S*2) to evolve at the same time, and compares the optimal individuals of the two populations, so as to break the equilibrium state in the population and to jump out of the local optimum.

In the optimization process, some generated particle may be beyond the boundary, the specific treated formula is described as follows,

$$\begin{array}{l} \neg \mathbf{x}\_{i}^{k} > \mathbf{x}\_{\text{max}} \Rightarrow \mathbf{x}\_{i}^{k} = \mathbf{x}\_{\text{max}}\\ \neg \mathbf{x}\_{i}^{k} < \mathbf{x}\_{\text{min}} \Rightarrow \mathbf{x}\_{i}^{k} = \mathbf{x}\_{\text{min}} \end{array} \tag{23}$$

where ¬... ⇒ ... represents the sign of "if ... then ..." ;*xk i* is the updated particle; *x*max is the maximum value of particle boundary; *x*min is the minimum value of particle boundary.

#### **4. Experimental Simulation**

#### *4.1. Data and Parameters for Experimental Simulation*

In this paper, the scenarios about Jinpu line No.1 and rail transit line No.12 in Dalian, China are selected as the research objects. Jinpu line No.1 is the intercity railway line that is 46.76 km long and has 11 stations in the initial stage, which extend from Jiuli light rail station to the terminal of Zhenxing road in Dalian under construction. Rail transit line No.12 is an urban rail transit line that is 40.38 km long and has 8 stations, which extend from Hekou station to terminal of Lvshun New Port. The simulation running line of scenario about Jinpu line No.1 is from the Jiuli to the 19th bureau, and the interval length between the above two station is 2.74 km. The running simulation line of scenario about rail transit line No.12 is from Lvshun New Port to Tieshan Town, and the interval length between the above two station is 2.94 km, with two long downhill and a long uphill ramps. The main parameters of the above scenarios are shown in Tables 2 and 3, ramp parameters and velocity limit for automatic train operation are shown in Figure 5.

**Table 2.** The main parameters of the scenario about Jinpu line No.1 in Dalian.

**Table 3.** The main parameters of the scenario about rail transit line No.12 in Dalian.

The calculation formula of traction characteristics is expressed as follows,

$$\begin{cases} F(\upsilon) = F\_{\text{max}} & \upsilon\_{\eta} \le \upsilon \le \upsilon\_{c} \\\ F(\upsilon) = \, ^{P}\_{\text{max}}\prime\_{\upsilon} & \upsilon\_{c} \le \upsilon \le \upsilon\_{d} \\\ F(\upsilon) = \, ^{P}\_{\text{max}} \times \upsilon\_{d}\prime\_{\upsilon} 2 & \upsilon\_{d} \le \upsilon \le \upsilon\_{\text{max}} \end{cases} \tag{24}$$

where *F*(*v*) is the instantaneous traction of vehicles, *F*max is the vehicle's maximum traction, *P*max is the maximum traction power of the vehicle, *vd* is the switching velocity of the constant power zone and the reduced power zone, *vq* is the switching velocity of the traction starting region and constant torque region, and *v*max is the maximum design speed of the vehicle.

**Figure 5.** The ramp parameters and the speed limit of experiment simulation for two ATO scenarios. (**a**) Jinpu line No.1 in Dalian scenario. (**b**) rail transit line No.12 in Dalian scenario.

Simulation optimization results must satisfy the following conditions; the train instantaneous speed cannot surpass the speed limit, the train must finish the whole process, and the parking error is less than 0.2 m. The improved shark smell optimization algorithm parameters are set as follows; the particle swarm size is 40, the weight coefficient is 0.9, the acceleration coefficients *C*1 and *C*2 are 0.5, crossover probability is 0.7, mutation probability is 0.09, selection probability is 0.5, the iteration number is 200, the shark population size is 40, the random number *R*1 = 0.4, *R*2 = 0.3, *R*3 = 0.25, *ηk* = 0.2, speed limit factor *βk* is 1.3, and the weight coefficient *αk* is 0.15. Setting parameters for the optimization algorithm is necessary to consider the convergence speed and optimization effect, such as population size and iteration number, with the increase of population size, the global search ability of the optimization algorithm will be enhanced, but the convergence speed will be reduced; similarly, with the increase of the number of iterations, the finding opportunity of optimal solution will be increase, and the more computing time and resource will be expend. The parameters characteristics and multiple experimental simulation results are taken into account for the above parameters setting. The multi-objective optimization parameters of the ATO scenario simulation of Jinpu line No.1 in Dalian are set as follows; the scheduled running time is 177*s*; *KE* ∈ [90,000, 150,000]; *KT* ∈ [0, 0.2]; *KJerk* ∈ [6, 10]; the intrinsic weight factors *ω* 1; *ω* 2 and *ω* 3 are, respectively, 0.5, 0.3, and 0.2; the target demand vector is [98,000, 6.2, 0.01]. The multi-objective optimization parameters of the ATO scenario simulation of rail transit line No.12 in Dalian are set as follows; the scheduled running time is 180 s, *KE* ∈ [80,000, 130,000], *KT* ∈ [0, 0.2], *KJerk* ∈ [5, 10], the intrinsic weight factors *ω* 1, *ω* 2, and *ω* 3 are, respectively, 0.5, 0.3, and 0.2, the target demand vector is [90,000, 5.2, 0.01].

#### *4.2. Matlab/simulink Simulation Results for Automatic Train Operation Scenarios*

According to the automatic train operation scenarios of rail transit line No.12 and Jinpu line No.1 in Dalian, China, the approximate optimal solutions are obtained by using the improved algorithm proposed in this paper, traditional improved shark smell optimization algorithm (chaotic shark smell optimization) [17] and traditional improved particle swarm optimization (dynamic multiple populations particle swarm optimization algorithm based on decomposition) [43], the above Matlab/simulink platform is chosen as a simulation platform. The specific configuration of the Matlab/simulink platform used in this paper is described as follow: the Matlab/simulink revision is Matlab GUI 2016b; the major computer configuration is "CPU Core i9-7920X @ 2.9 GHZ" and "Windows 10". The specific Matlab/simulink optimization simulation results are shown in Figures 6 and 7 and Tables 4–7, and three different algorithms are recorded as ISSO, CSSO, and dMOPSO.

**Figure 6.** The Matlab/Simulink optimization curves of different algorithms for Jinpu line No.1 in Dalian scenario. (**a**) Jinpu line No.1 in Dalian scenario. (**a**) Target velocity trajectory profiles. (**b**) Operating condition distance curves. (**c**) Iterative convergence curves of each optimization objectives. (**d**) Iterative convergence curves of unified goals.

**Table 4.** The Matlab/simulink optimization results of different algorithms for Jinpu line No.1 in Dalian scenario.

**Figure 7.** The Matlab/simulink optimization curves of different algorithms for rail transit line No.12 in Dalian scenario. (**a**) Target velocity trajectory profiles. (**b**) Operating condition distance curves. (**c**) Iterative convergence curves of each optimization objectives. (**d**) Iterative convergence curves of unified goals.

**Table 6.** The Matlab/simulink optimization results of different algorithms for rail transit line No.12 in Dalian scenario.

**Table 7.** The HILS evaluation of the results of different algorithms for Jinpu line No.1 in Dalian scenario.

As can be seen in Tables 4–7, the optimization solution obtained by the improved ISSO is superior to that of CSSO and dMOPSO, and three indexes of energy saving, punctuality, and comfort have been improved considerably. The ATO experiment simulation scenario for Jinpu line No.1 from Jiuli to the 19th bureau is a typical slope with ups and downs. It is necessary to keep the train at high speed before driving down the long down ramp. The ATO experiment simulation scenario for rail transit line No.12 from Lvshun New Port to Tieshan Town is in the hilly of Dalian Lvshunkou district, and the hilly region is the typical characteristics of Dalian. When the train is running in such a terrain, the control sequence should be concise. The train speeds up in the long down slope and slows down in the long up slope as much as possible, which can save energy and avoid turbulence. As can be seen from Figures 6a,b and 7a,b ( the target velocity trajectory and control sequence ), the improved algorithm ISSO can obtain extremely simple control sequence and make the most of long down and up slopes, so as to obtain the target velocity trajectory as smooth as possible. As can be seen from Figures 6c,d and 7c,d (the iterative convergence curves), the convergence rate of ISSO is faster than that of CSSO and dMOPSO. Even in the late iterations, ISSO has a strong ability of global convergence performance.

Compared with other comparison optimization algorithms, ISSO has several obvious advantages in matlab/simulation environment, as there is huge difference between matlab/simulation environment and actual scenario yet, the effectiveness of ISSO is necessary to further test and verify.

#### *4.3. HILS Results for Automatic Train Operation Scenarios*

Matlab/Simulink is a simulation technology that is completely separated from the real train operation environment. Therefore, some problems (such as real-time velocity sampled inaccurately, signal delay and packet loss in transmission, a certain degree of unstable in tracking control system, etc.) need to be considered in the actual control process and are cannot be truly reflected. To more effectively test the performance of the optimization algorithm in the actual train operation environment, dSPACE hardware-in-the-loop simulation (HILS) technology is adopted to write the verified optimization algorithm or control algorithm to the chip of the optimizer or controller. Compared to traditional simulation platform based on pure software environment, dSPACE HILS platforms contain the real on-board equipments, which can truly reflect the the real situation for automatic train operation. Due to the restriction of funds and experiment conditions, it is difficult to implement the corresponding actual automatic train operation experiment. Moreover, compared with real actual experiments, dSPACE HILS has the advantages of low experimental cost, implement easily, and high security protection of personal and equipment [33]. Therefore, HILS is highly valued by researchers and developers, abundant HILS-related products are researched, developed, and applied in various fields and numerous relative research results have been achieved [44,45]. The semi-physical simulation equipments in automatic train operation mainly include optimizer, controller, sensors, simulator, conditioning circuit, signal conditioning unit, and connector. The simulation topology diagram of automatic train operation HILS platform, the structure diagram for automatic train operation HILS, and the physical diagram of controller cabinet and simulation cabinet for automatic train operation HILS are shown in Figures 8 and 9.

**Figure 8.** The simulation topology diagram of automatic train operation hardware-in-the-loop simulation (HILS) platform.

**Figure 9.** The structure diagram for automatic train operation HILS.

As can be seen from Figure 8, the Controller/Optimizer and its Service equipment make up the automatic train operation HILS platform, the Monitor serves as the monitoring center for human–computer interaction, and the Alarm set and Circuit breaker are the warning and protection center. Alarm set detects the working status of each HILS device in real-time, and the Circuit breaker is used to protect equipment by breaking when a certain kind of exception is detected. As can be seen from Figure 9, the automatic train operation HILS platform contains various actual hardware, and simulation hardware for automatic train operation. The "lower-layer control loop" based on controller and "upper-layer optimization loop" based on optimizer of HILS constitute the two independent communication systems, respectively, and there is a certain connection between the two loops. "Controller" is also named as traction control unit (TCU), which could provide control commands for corresponding equipments in real-time using a proper control algorithm in real-time, enable the the urban rail vehicle to track the ideal profile curve; "Optimizer" is also named as main processor unit (MPU), which could provide the velocity ideal trajectory profile (target speed curve of automatic train operation) for "lower-layer control loop" tracking control based on 'dSPACE emulator'. Second, the "dSPACE emulator", "conditioning circuit", "signal processing unit", "sensors", "connectors", and so on are service equipments for ATO HILS: "dSPACE emulator" provides some correlative simulation environments for the automatic train operation HILS, the related models included such as accurate braking model, traction transformer model, running line model, velocity fluctuation model, etc.; "conditioning circuit" can regulate electrical waves properly for "Controller" appropriately; "signal processing unit" can regulate net signals properly for "Optimizer" appropriately; and the "sensors" and "connectors" are used to feed electrical waves and net signals back to the "Controller" and "Optimizer" in real-time. Third, the "DC power source", "Converter system", "Electric motor", "Digital rheostat box", and "Gear box" are simulation electric hardware equipments of in place of the actual: "DC power source" acts as the actual pantograph; "Converter system" transfer the electric energy from "DC power source" to "Electric motor" through a series of current-voltage conversion processes, it includes AC–DC converter, DC–AC converter, low-pass filters, etc.; "Electric motor" acts as the actual urban rail vehicle motors set; the capacity of simulation circuit loop is smaller than actual. In Figure 11, "train controller cabinet" and "train emulator cabinet" are the vital equipments for automatic train operation HILS, except controller and emulator, conditioning circuit, signal processing unit and corresponding switch groups are included.

Obviously, compared with these automatic train operation scenarios based on pure software (such as Matlab/simulink simulation), their identical scenarios based on HILS are closer to actual. Therefore, based on the same scenarios, tracking control algorithm and HILS platform, the comprehensive performance index for automatic train operation obtained by optimization algorithms could reflect the optimization performance of these algorithms effectively. To further verify the effectiveness of the algorithm, according to the automatic train operation scenarios of rail transit line No.12 and Jinpu line No.1 in Dalian, China, the approximate optimal solutions are

obtained by using ISSO, CSSO, and dMOPSO, the above ATO HILS platform is chosen as simulation platform. The specific configuration of the ATO HILS platform used in this paper is described as follow: the Matlab/simulink revision is Matlab GUI 2016b; the major computer configuration is "CPU Core i9-7920X @ 2.9 GHZ" and "Windows 10"; the core chip of "Controller" and "Optimizer" is "TMS320F28335"; the simulation software of "dSPACE emulator" is dSPACE software(revision is control desk 6.1); the communication protocol of the ATO HILS is MVB (multifunction vehicle bus); the MPC (model predictive control) algorithm is adopted as tracking control algorithm; the three optimization algorithms (ISSO, CSSO, and dMOPSO) used are written in the kernel chip of "Optimizer" in "upper-layer optimization loop" for contrasting. The specific HILS optimization results obtained by "Optimizer" in "upper-layer optimization loop" are shown in Figures 10 and 11 and Tables 8–11.

**Figure 10.** The HILS optimization curves of different algorithms for Jinpu line No.1 in Dalian scenario. (**a**) Target velocity trajectory profiles. (**b**) Operating condition distance curves. (**c**) Iterative convergence curves of each optimization objectives. (**d**) Iterative convergence curves of unified goals.

**Table 8.** The Matlab/simulink evaluate results of different algorithms for rail transit line No.12 in Dalian scenario.

**Table 9.** The HILS optimization results of different algorithms for Jinpu line No.1 in Dalian scenario.

**Figure 11.** The HILS optimization curves of different algorithms for rail transit line No.12 in Dalian scenario. (**a**) Target velocity trajectory profiles. (**b**) Operating condition distance curves. (**c**) Iterative convergence curves of each optimization objectives. (**d**) Iterative convergence curves of unified goals.

**Table 10.** The HILS optimization results of different algorithms for rail transit line No.12 in Dalian scenario.

**Table 11.** The HILS evaluate results of different algorithms for rail transit line No.12 in Dalian scenario.

In Figures 10a,b and 11a,b (which is related to the optimization effect of automatic train operation process), the power is switched on, the pantograph is raised, the dSPACE simulator is in the working state (the "dSPACE" button is pressed), the human–computer interaction signal is normal (the "Design" button is green), and the circuit breaker is normally closed, which is in a normal optimization simulation state.

According to optimization simulation results of different algorithms in Tables 8–11, compared with the traditional improved particle swarm optimization algorithm (dMOPSO) and traditional improved shark smell optimization algorithm (CSSO), the improved algorithm proposed in this paper (ISSO) has the obvious performance improvement effect, the three indexes of energy saving, punctuality and comfort of target velocity trajectory have been improved considerably; meanwhile, and the convergence evolution generations and computation time have been reduced considerably. As can be seen from Figures 10a,b and 11a,b, the ideal target velocity trajectory obtained by ISSO was the smoothest, compared with the contrasted algorithms (CSSO and dMOPSO), ISSO obtained extremely simple control sequence and made use of the most of long down and up slopes sufficiently, it enables the train to reduce unnecessary traction and braking status and to make full use of coasting status as much as possible. As can be seen from Figures 10c,d and 11c,d, compared with contrasted algorithms (CSSO and dMOPSO), the three optimization objective indexes and two unified goals obtained by ISSO have been improved to a considerable extent not only in the optimization effect but also in the computation efficiency.

Compared with the improving optimization effect of ATO target velocity trajectory optimization, improving the actual tracking control effect is more practical. The optimization effect calculation speed could estimate the performance of the ATO target velocity trajectory optimization algorithms, achievable difficulty for tracking control is also a significant evaluation index. To better verify the effectiveness of the ISSO, the "lower-layer control loop" is used to tracking control according to optimal target velocity trajectory generated by the "upper-layer control loop". The MPC (model predictive control) algorithm has some its own advantages such as good performance in tracking precision, powerful robustness, fast tracking speed, etc. and widely in ATO traction system, so it is selected as tracking control algorithm in this paper. The specific HILS tracking control results obtained by "Controller" in "lower-layer control loop" are shown in Figures 12 and 13 and Tables 12 and 13.

**Table 12.** The HILS tracking control results of different algorithms for Jinpu line No.1 in Dalian scenario.

**Table 13.** The HILS tracking control results of different algorithms for rail transit line No.12 in Dalian scenario.

In Figures 12 and 13 (which relate to the tracking control effect of automatic train operation process according to optimization result), the power is switched on, the pantograph is raised, the dSPACE simulator is in the simulation state (the "dSPACE" button is pressed), the design parameters cannot be changed (the "Design" button is white), and the circuit breaker is normally closed, which is in a normal tracking control simulation state.

**Figure 12.** The HILS tracking control curves of different algorithms for Jinpu line No.1 in Dalian scenario.

**Figure 13.** The HILS tracking control curves of different algorithms for rail transit line No.12 in Dalian scenario.

According to optimization simulation results of different algorithms from Tables 12 and 13, compared with the traditional improved particle swarm optimization algorithm (dMOPSO) and traditional improved shark smell optimization algorithm (CSSO), the corresponding tracking control curve tracking the target velocity trajectory, obtained by improved algorithm proposed in this paper (ISSO), has the obvious performance improvement effect, the its three indexes of energy-saving, punctuality, and comfort have been improved considerably. As can be seen from Figures 12 and 13, compared with the contrasted algorithms (CSSO and dMOPSO), the corresponding tracking control curve for ISSO has better tracking control effects. As can be seen from the enlarged areas of Figures 12 and 13, the velocity fluctuation of the velocity distance trajectory curves corresponding to ISSO is weaker, because the target velocity trajectory obtained by ISSO was more ideal and tracking controlled easily in the automatic train operation simulation scenarios.

To verify that the improved algorithm (ISSO) in this paper has more superiorities, such as optimization effect, calculation speed, easily achievable in tracking control, etc., two scenarios about rail transit line No.12 and Jinpu line No.1 in Dalian, China are chosen as optimized objects for ISSO and some comparison algorithms, the corresponding optimization and tracking control results of matlab/simulation and hardware-in-the-loop simulation (HILS) as shown in Tables 4–13 and Figures 5 and 6, and 10–13. Obviously, compared with the contrasted algorithms (CSSO and dMOPSO), the corresponding optimization and tracking control results obtained by ISSO are more ideal. It is indicate that ISSO is a appropriate algorithm with powerful optimization capability, so as to solving practical automatic train operation more effectively.

#### **5. Conclusions**

Automatic train operation target velocity trajectory optimization is a very complex issue that needs to take into account multiple objectives, and the ideal optimization solution is not easy to be obtained. An improved multi-objective dual-population shark smell optimization algorithm using incorporated composite angle cosine for automatic train operation velocity trajectory optimization is proposed in this paper, and the specific advantages are described as follows.

For the multi-objective optimization problem, the evaluation index of the solution is very important. Nonetheless, the linear weighted target belongs to the common multi-objective unified target, and there is a problem that subjective parameters are selected blindly. In this paper, the multi-objective angle cosine is used as evaluation index, which can effectively avoid the blind selection of subjective parameters. To make the evaluation index more reasonable, on the basis of the preference information, the composite angle cosine which takes into account both the preference difference and the numerical difference is be used as the evaluation index of the solution in this paper.

Because the updating rules of velocities and positions in the SSO and PSO make all individuals close to the extreme individuals, after a long iteration, the extreme individuals will form a certain degree of dominance over the population, which makes it difficult to converge globally. First, the dual-population evolution mechanism is used to jump out of the local optimum. Second, to suppress the local convergence of optimization algorithm, it is necessary to determine accurately whether the individual aggregation occurs in shark population. In this paper, the fusion distance can be used to take into account the relativity and independence of speed and time, which can detect whether there is the phenomenon of individual aggregation preciously, so the local convergence is better suppressed. At the same time, this paper also introduces the dual-population collaborative optimization mechanism of SSO algorithm and particle swarm algorithm to further improve the optimization performance of the algorithm.

According to the the Matlab/simulink results and HILS results about automatic train operation scenarios, compared with the conventional improved shark smell optimization algorithm and the conventional improved particle swarm optimization algorithm, the improved algorithm (ISSO) proposed in this paper improves the calculation accuracy and optimization ability of the optimization algorithm to some extent, so that more ideal target velocity trajectory can be obtained. This clearly shows that ISSO have been improved to a considerable extent for the automatic train operation target velocity trajectory optimization not only in the pure software scenarios but also in the semi-physical scenarios. As the automatic train operation HILS is close to its actual experiment, the problem that the ISSO has worse improvement optimization effect in the actual than anticipatory could be avoid to certain extent.

**Author Contributions:** The work presented here was performed in collaboration among all authors. L.W. designed, analyzed, and wrote the paper, and completed the simulation experiment. X.W. guided the full text and provided simulation conditions. Z.S. conceived idea and involved simulation experiment. S.L. involved simulation experiment and analyzed the data. All authors have read and agree to the published version of the manuscript.

**Funding:** This research was funded by The Nature Science Foundation of China (grand number 51609033 and 61773049).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

HILS hardware-in-the-loop simulation

ISSO improved shark smell optimization

CSSO chaotic shark smell optimization

dMOPSO multi-objective particle swarm optimization based on decomposition

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Method for the Optimization of Daily Activity Chains Including Electric Vehicles**

#### **Dimitrios Rizopoulos \* and Domokos Esztergár-Kiss**

Department of Transport Technology and Economics (KUKG), Faculty of Transportation Engineering and Vehicle Engineering (KJK), Budapest University of Technology and Economics (BME), 1111 Budapest, Hungary<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>Tel.: +30-697-865-6250

Received: 29 January 2020; Accepted: 14 February 2020; Published: 18 February 2020

**Abstract:** The focus of this article is to introduce a method for the optimization of daily activity chains of travelers who use Electric Vehicles (EVs) in an urban environment. An approach has been developed based on activity-based modeling and the Genetic Algorithm (GA) framework to calculate a suitable schedule of activities, taking into account the locations of activities, modes of transport, and the time of attendance to each activity. The priorities of the travelers concerning the spatial and temporal flexibility were considered, as well as the constraints that are related to the limited range of the EVs, the availability of Charging Stations (CS), and the elevation of the road network. In order to model real travel behavior, two charging scenarios were realized. In the first case, the traveler stays in the EV at the CS, and in the second case, the traveler leaves the EV to charge at the CS while conducting another activity at a nearby location. Through a series of tests on synthetic activity chain data, we proved the suitability of the method elaborated for addressing the needs of travelers and being utilized as an optimization method for a modern Intelligent Transportation System (ITS).

**Keywords:** daily activity chains; electric vehicles; optimization; charging stations; intelligent transportation systems; ITS

#### **1. Introduction**

In most urban environments, there is an ever-growing need for navigation through transportation networks. Although more and more mobility services are offered to travelers, the complexity of their use has risen alongisde the offer and demand for such services. Thus, the utilization of Intelligent Transportation Systems (ITS) and the development of travel-related services have become an immense need for the facilitation of the everyday life of citizens.

With the expected rise of the use of Electric Vehicles (EVs) all over the world [1] and the positive impacts that they are expected to bring [2], advanced EV travel planning algorithms are now essential more than ever [3].  Not only their development and availability are vital, but also the user experience and the satisfaction of specific user preferences are two crucial points when creating such services [5]. Nonetheless, those solutions can help to mitigate the risk of electrification of transportation and deal with one of the enormous obstacles in EV adoption that is called range anxiety, a term used to describe the psychological worry of EV drivers that the remaining electric energy in the vehicle is not enough to attend their destination. While the battery capacity of the EVs is growing at a satisfactory rate and is even expected to double until 2030 [6], it is proven that range anxiety is still one of the most significant concerns of EV drivers and is one of the most critical obstacles towards universal adoption and mainstream use [7]. While it has been shown that EV users that have a charger

available at home mainly cover their needs at a satisfactory rate, for users that solely depend on public charging infrastructure, charging can be a more challenging task [8].

Alongside the ability of those methods to address problems of EV drivers on an individual level, they can improve the situation for the system-level [3], too. Mitigation strategies on electrification on a regional level, as their importance as highlighted in recent articles [9], can benefit from driver assistance ITS, since the underlying methods can be independent of the spatial-context and are reproducible over several regions according to needs of travelers of those specific regions.

Apart from the work of Cuchý et al. [3], little research has been conducted to study the flexibility or willingness to use travel planning services for EV drivers. However, the authors believe that based on other previous research on other modes of transport, and mainly conventional automobiles [10], a driver assistance system customized for EV users can be of crucial importance. This observation can be further highlighted by other literature studies that have shown that travel patterns can vary further when influenced by technologic systems [11]. While traveler groups' flexibility can differ according to age groups [12] and social characteristics [13], if the benefits acquired from the planning process are perceived as worthy [14], benefits emerge for planning trips in pre-trip phase [15] or in real-time [16].

The goal of this research is to provide a realistic approach to the Daily Activity Chain Optimization (DACO-EV) problem by exploiting activity-based modeling and the Genetic Algorithm (GA) framework. While the incorporation of the EV range constraints is a significant part of our work, a particular focus is also given to mechanisms that enable the implementation of such a system in a real-world ITS system that will serve the travelers. Section 2 presents a literature review on related topics. Section 3 provides a detailed description of the DACO-EV method, its parameters, and its attributes. Section 4 discusses the implementation of the method and to provide realistic solutions to travelers. Section 5 presents the results, including use cases. Section 6 serves as a conclusion to the article and its contribution.

#### **2. Literature Review**

To the best of our knowledge, an approach to the DACO-EV method, including enough parameters to render its solution useful in realistic settings, does not exist in the scientific literature prior to the publication of this article. However, insights and valuable research directions can be extracted from articles regarding the incorporation of electromobility into transportation problems.

The activity-based analysis of transportation systems stems from the idea that the demand for transportation between locations is interconnected with the demand for participation in some activity. As indicated by a series of works in the literature [17–19], it has been extensively studied and used for the analysis of travel behavior and the planning of transportation systems. Our work on the optimization of daily activity chains falls under the umbrella of activity-based research in transportation. Recent notable articles [20,21] show that activity-based modeling has been successfully applied to cases where electromobility is incorporated into transportation systems. In the first paper by Dong et al. [20], the impact of different types of deployment of charging infrastructure is analyzed in regard to the range anxiety of EV users. An activity-based method is used to evaluate the deployments according to travel patterns. In the work by Kontou et al. [21], activity-travel patterns from a National Household Travel Survey data are used, in order to assess two schemes for centralized charging management of EVs, one scheme that regards preferences of individuals and another that considers the government perspective. The results indicate the differing nature of the interests of those two stakeholders, and confirm that the availability of charging stations at the workplace can greatly affect the charging profiles of EV travelers. The latter observation is one of great importance that can be found in other studies [8], too.

In regard to the problem of Daily Activity Chains Optimization (DACO), the basic instance of DACO-EV, it has been previously addressed in the literature, and a GA approach has been introduced [22,23]. Charypar et al. [24] have addressed a similar problem and introduced a similar solution approach. The authors have considered several simplifications, such as the computation of distances according to geometric distance and not based on the real transportation network. Another similar work was conducted by Abbaspour et al. [25], where the authors tackled the

same problem but had a specific focus on the touristic aspect of the optimization of activity chains. Although the aforementioned papers regard heuristic optimization methods that have been developed to solve the DACO problem, different modelings could also be considered. A modern approach was presented by Liao et al. [26], who modeled the problem as a graph super-network. They included space-time constraints to enable the valid selection of locations for activities in time, multi-modal transport options and modeled the parking choice. While their activity-travel scheduling algorithm computes the final solutions to the problem in minimal time, in the scheduling problem and the modeling, the priorities and the flexibilities of the travelers are not included, like in the case of DACO.

There is a lack of exact optimization approaches for the DACO or DACO-EV in the literature. Exact optimization methods, where the gap between the solution calculated and the globally optimal solution of the problem is zero, are based on mathematical models of the transportation problem. Although exact approaches tend to provide better solutions in terms of reduction of cost according to criteria, they cannot yield solutions in the amount of time that heuristics can, and usually take much more time to complete the calculations. Works can be found in the literature where routing in multi-modal networks has been addressed with exact approaches [27,28].

One motivation for the extension of our method was to provide an efficient tool dealing with the range anxiety of EV drivers in urban environments. Range anxiety is one of the main bottlenecks of EV adoption and affects both inner-city and inter-city trips. Many studies [29] exist in the literature that addressed this topic and provided useful insights regarding key factors that can evaluate its effect and help deal with it. Although technology, and its advancement, is one of the main ways to deal with range anxiety, other key factors are considered as important for the medium-term horizon, such as battery costs, coverage of the charging stations network and CO2 vehicle standards [7]. Range anxiety has also been studied in the context of daily schedules by Neubauer et al. [8]. Significant findings in their studies showed that workplace charging could play an essential role in promoting the utility of EVs for high mileage commuters. Additionally, the broadly available charging infrastructure has been shown to be important to both high and low mileage commuters. Range anxiety importance has also been shown in computational studies that discussed its relationship with the CS network. In the work by Guo et al. [30], the battery charging station location problem was addressed where range anxiety and distance convenience were taken into account. While the KSIGALNS algorithm that they developed to solve their bi-level integer programming model was proven to be very effective compared to standard and previous solutions approaches, the analysis of range anxiety as a parameter shows that it is an essential factor to the location strategy of the CS network. Other interesting computational studies exist in the context of range anxiety, such as the one by Esmaili et al. [31], where the authors examine EVs as distributed energy storage units and their potential contribution to microgrids when vehicle-to-grid service is considered. An interesting result that emerges from the solution of their Mixed-Integer Linear Programming (MILP) model is that when range anxiety, as modelled in their mathematical optimization program, is in higher levels, then the average State of Charge of the EV drivers rises and, as a result, the total cost of the microgrid is higher, too.

Range anxiety, however, is a phenomenon that is being progressively dealt with on both a personal and systemic level. It is more evident than ever that not all EVs are the same regarding experiencing range anxiety [32]. There are observations that range anxiety is minimally experienced by users of specific types of EVs in specific regions, such as in the study by Gorenflo et al. [33], in which, their analysis of the e-bike usage and battery charging data supports that range anxiety does not exist among the participants of the survey. Finally, in order to underline the importance of guidance systems, like the one presented in this paper, we further present the study of Cuchý et al. [3]. In their work, they were able to evaluate multi-destination transportation scenarios according to single-user perspective and infrastructure perspective by using AgentPolis simulation framework. Except for the obvious benefits for EV users, they showed that when users plan their mobility up-front, the total utilization of infrastructure can increase by more than 100%.

Particular focus was given to articles dealing with EV routing and the parameters that were considered when it comes to the range that vehicles can cover. Extensive work on the subject has been conducted by the team of Baum. Notably, in work conducted in 2013 [34], they extended their algorithm, named Customizable Route Planning (CRP), to calculate fast queries on graphs that are suited to Electric Vehicles (EV). Their goal was to calculate energy optimal routes, where several parameters can be included, such as the recuperation of the electric vehicle, the battery capacity constraints, and the dynamic behavior in energy consumption of the vehicle. In their next work conducted in 2016 [35], they proposed a framework with a more holistic approach to the problem, where parameters are included, such as the use and location of charging stations, turn costs, and the adjustment of speed in order to save energy. In their work conducted in 2015 [36], they introduced an approach named CHArge, which solves the EV routing problem in realistic settings. They discussed the properties of charging functions, which are used to map an initial State of Charge (SoC) and the duration of charging sessions with a resulting SoC. The works by Baum et al. [34–36], are very close to our approach in terms of methods and calculation mechanisms that were used to deal with the incorporation of electromobility. Additionally, useful information about the charging network, charging sessions, plugs, and strategies were also extracted from the work of Moghaddass et al. [37].

There is a series of articles that are connected to our work since they address the optimization of activity chains considering the use of EVs. The work of Liao et al. [38] is very close to our approach, but they do not solve the same problem. The authors provide modeling and solution to the EV shortest travel time path problem and the fixed tour EV touring problem. They consider a battery swap system for the touring problem, where two cases are discussed: the on-site station model, where each city is considered as a swapping station, and the off-site model, where the swapping station is further away from the city. Although their approach addresses a similar problem to ours, the presented work is not directly linked to a real-world application. In the works of Cuchý et al. [39,40], the authors tackled the Whole Day Mobility Planning with Electric Vehicles (WDMEV) problem and modeled it according to graph theory. The definition of the problem is similar to DACO-EV, but the authors included fewer parameters than in our problem definition. They did not include the priorities of the travelers, and they are neglecting the time windows of the desired attendance of the traveler. Furthermore, the authors did not include in the models the consumption according to Worldwide Harmonized Light Vehicle Test Procedure (WLTP), different EV models, and it is not clear if they consider the charging connector types. In the modeling, the overall weight of the vehicle and the option to have a final desired amount of energy in the battery when the tour ends were not included. Finally, the availability of the chargers was not taken into account in the calculation of the routes, but they considered waiting time at Charging Station (CS) according to randomly parameterized models. The aforementioned parameters may lead to very different solutions spaces and computational results, so the two problems and methods cannot be directly compared.

Overall, the main contribution of this article lies in the elaboration of a method that addresses the DACO-EV problem, including all the necessary parameters that allow the solution of the problem to be personalized to the specific needs of travelers, and enabling the use of the method in realistic settings. The method includes:

#### **3. Method Description**

#### *3.1. Parameters of the DACO-EV*

The DACO-EV method incorporates the use of an EV; therefore, it requires modeling of the constraints regarding the energy equilibrium that takes place in the energy storage mechanism of the car while visiting several locations. Because EVs are still in an early phase in their adoption by the public, the available charging facilities are not always optimally distributed, and the minimum distance between available charging points can reach up to several kilometers depending on the type of road network [41]; this situation renders a demanding reality for EV drivers that need to utilize the public charging network to cover their charging needs. The extra parameters that affect this energy equilibrium can be grouped into three major categories, such as vehicle design, driver behavior, and the environment [42]. In our attempt to provide a realistic and practical method for the elaboration of meaningful solutions to the EV drivers, we have included several of those parameters, while we ignored others of lesser importance.

The problem and the elaborated method are based on three sets of parameters. The first set is defined by the traveler and refers to his or her schedule and overall preferences concerning the activities. The second group of parameters is about the constraints related to electromobility when the traveler utilizes an EV. The third set consists of the parameters that do not depend on the user and refer to the existing network, the timeframe of operation of facilities and services. Table 1 refers to the parameters that are set by the traveler.

**Table 1.** The first group of parameters.

Particular emphasis should be put on the second group of parameters about the EV and its use. Since the traveler specifies the used EV model, the parameters of this second group can be considered as sub-parameters of the first group (Table 2).

**Table 2.** The second group of parameters concerning the characteristics of the Electric Vehicle (EV).

Finally, there is the third group of parameters that are needed for the problem to be fully described. Those parameters are not defined for each passenger, but instead, they are specified for the network, the modes of transport, and the types of activities available. They can also be described as the static parameters of the DACO-EV problem (Table 3).

**Table 3.** The third group of parameters that are independent of the user's input.

#### *3.2. Degrees of Freedom of the Problem and Decision Variables of the Solution*

At this point, we should make a reference to the decision variables of the DACO-EV that fully define a solution for the travelers. Those decision variables are essential for the modeling as they will form the chromosomes of the solution encoding of the Genetic Algorithm. In Table 4, the decision variables and their descriptions are provided.

All of the aforementioned variables (ones included in Table 4) are the independent variables of the problem. By deciding on the order, the locations, the time windows of attendance, the modes of transport, and the charging time, we decide on a specific way that the first tour changes and gets optimized. Of course, that information is not enough in order to fully articulate a meaningful solution for a traveler in a human-readable format. From the independent variables, we can extract a series of information, such as exact paths to be followed in between activities, battery power after the visits at each of the locations, starting and ending time of the whole tour, travel time of the tour, starting and final state of charge of the EV.

**Table 4.** Decision variables of DACO-EV.

#### *3.3. Charging Scenarios*

To fully cover the description of the DACO-EV problem, we also need to deal with the possible types of charging. As a first charging scenario, we consider the case where the EV users want to charge their vehicle at a CS, and they want to conduct one of their activities at another location. This is called the *en passant* charging scenario. Usually, it can occur when a charger and a location for an activity are close to each other, and the user can just walk a few hundred meters to the facility and back to the charging station after completing the activity. As a second charging scenario, the classic scenario of charging is included, where the traveler again deviates from the first tour and embarks on a subtour to a CS, but in this case, he or she just waits in the car until the charging is finished.

First, we provide a graphical example of the classic charging scenario of the DACO-EV (Figure 1). Given a set of five initial activities, the activity chain is optimized. For simplicity reasons, we only provide an optimized tour where nothing changes in the spatial context except for the inclusion of a charging session. The charging station is indicated by the blue circle, which is a new activity in the chain.

**Figure 1.** The classic charging scenario of the DACO-EV problem.

In Figure 2, present the *en passant* charging scenario of the DACO-EV. The charging detour happens between Activities 4 and 5, and the transfers from the CS to the location of Activity 5 are represented by the dotted lines. In this case, the EV remains at the CS, while the user processes Activity 5 using another mode of transportation (e.g., walking).

**Figure 2.** The *en passant* charging scenario of the DACO-EV problem.

#### *3.4. Optimization Method*

The method is based on a Genetic Algorithm (GA) that iteratively solves the Travelling Salesman Problem with Time-Windows (TSP-TW) for different combinations of the characteristics of the possible solutions and evaluates each of them according to several criteria. The efficiency of the proposed method has been previously tested and proved for the primary instance of the problem in the articles by Esztergár-Kiss et al. [22,23]. After running the algorithm for several iterations, a more optimal solution for the initial schedule and the preferences of the traveler is derived. In the algorithm, we make sure that after the traveler executes the tour suggested, and given the EV usage constraints, the EV has the desired FSoC when it reaches home.

By defining the properties (i.e., decision variables) of an individual (i.e., candidate solution) of the problem, we can form the populations of solutions in the solution space. Then the genetic operators of selection, mutation, and crossover are utilized to search the available space of solutions and derive the ones that are more suitable for each case of traveler. Except for the genetic operators that have been used, and according to GA literature [43], the method is characterized by a set of parameters, such as:

#### 3.4.1. The Genetic Algorithm

Although the GA framework contains more or less standard data structures, functions, and operations utilized for each application, it was considered appropriate that a more formal specification of the algorithm is presented. This algorithm describes how temporal and spatial flexibility were handled and includes several technical details in order to enable its reproduction.

#### 3.4.2. Fitness Function and Optimization Criteria

For the efficient solution of the DACO-EV, several criteria have been considered. While most of them regard travel time, whose reduction has been our primary aim, a few more criteria were considered as important to include. The method is capable of handling five criteria. They were weighted by the corresponding parameters *a*, *b*, *c*, *d*, *e* to create an efficient multi-criteria fitness function for the problem.

$$\mathcal{L}I(X) = (a \ast TT) + (b \ast TTST) + (c \ast AT) + (d \ast \mathcal{W}T) + (e \ast CT) \tag{1}$$

Fitness function *U* has been applied to every candidate solution *X* that has been calculated by the genetic algorithm. This evaluation, as indicated in the algorithm at *Step 3.4*, was applied to each potential solution, when not infeasible, and has been utilized to label and prioritize according to the potential needs of travelers. For the first criterion, *TT*, its actual value is given in regard to the time that the user has spent traveling within his or her tour. The second criterion, the Travel Time in the Subtour (TTST), regards travel time only in the charging subtour. This specific type of travel time has been calculated as an extra variable, in order to be able to be separately included in cases, where the effect of the subtour needs to be minimized in the overall solutions. Arrival Time (AT) at the final destination is also another criterion considered that can be essential for some users, as well as, the Walking Time(WT), when it occurs in the sub-tour for charging, and the en-passant scenario. Finally, Charging Time (CT) has been additionally considered as an extra criterion, and has been added to the function as a possibility for optimization. The reason for including all five criteria was to provide a robust optimization approach that is applicable to a wide set of real-world applications. While usage of the method is independent of spatial context, the potential needs of travelers vary across different cities, countries, and the available charging stations network.

For the values of the parameters *a*, *b*, *c*, *d*, *e*, different solutions are derived. The proper choice of the values is of great importance to the type of solutions that will be calculated. A sensitivity analysis of the effect of the parameters on the solutions produced is one of the future steps to take in order to assure suitable performance for each real-world case (i.e., region of application) and understand the effect of each criterion on the activity chain.

#### **4. Implementation of the Method**

For the efficient functioning of a modern application that will support the travelers in their daily commutes, the implementation of the method has been realized. Because of the complexity of the problem addressed and the vast amount of calculations needed, a series of heuristic rules were applied to the method in order to direct the search of the solution space and yield meaningful results within realistic computation times. This addition to the method is called the pre-optimization phase, which creates a smaller solution space for the algorithm to search. However, in this paper, we focus on the charging related functions of the method.

For the implementation of the method, Python programming language was utilized. Data from OpenStreetMap and Geofabrik [44] was used for establishing the road network, OpenTripPlanner (OTP) engine [45], was the primary tool for the calculation of road distances and routes, and calculations for the availability of CS were conducted based on data from e-MOBI [46] about the charger network. The chargers based on which the availability was calculated were 89. Elevation data are taken from NASA's Shuttle Radar Topography Mission (SRTM). Regarding the EV model used, the traveler is allowed to enter the vehicle that he or she is going to use, and we were able to extract the properties (battery capacity, available plug types, consumption according to WLTP, vehicle weight) of each vehicle from an electric vehicle database [47]. A simplification was made regarding the charging rate, which was considered according to the available plug types at both the charging station and the EV. A matchmaking mechanism was developed to match EVs with fast charging plugs if they are likely to be available, or with standard charging plugs, otherwise. Regarding the priorities of the travelers, we included them in four levels which are described in Table 5:

**Table 5.** Priorities of the travelers considered.

When an activity is assigned with Priority Label 1, it means that the traveler considers this activity as both spatially and temporally fixed. Priority Label 2 is assigned by the travelers to activities that are spatially flexible, which means they can be conducted at several locations, but temporally not, meaning that the algorithm is bound to calculate the visit within the desired time windows. With Priority Label 3, the traveler indicates that the activity must be conducted at the provided location, but it can happen in any time-window within the day. Oppositely to Priority Label 1, Priority Label 4, is the case where the traveler indicates that the activity is totally flexible, which means that it can be both conducted in another location (providing the same type of service) and at any time (within the overall timeframe of the user and the operating hours of venues and shops).

Before actually providing information about the implementation of the GA and its mechanisms that allow the effective solution of the DACO-EV, an overview of the work is presented in Figure 3 that depicts the method, its parameters, external tools, and data used. On the one hand, the traveler provides an initial activity chain and information about his or her preferences by setting the first, and thus the second, group of parameters. On the other hand, the third group of parameters provide external information, such as the network topology (road network and elevation), the electric vehicle

database for consumption calculations, and the alternative activities' location database for calculating the optimized activity chain by examining the full spectrum of the spatial possibilities for the activity chain. Finally, based on the charging station usage data and other related information to the CS (i.e., charging stations database), the availability of charging stations was calculated using data mining techniques. Before using the heuristic optimization model, a pre-optimization is performed, and the routing engine is used to calculate optimal routes between activity locations. As a result, the genetic algorithm calculates optimized activity chains for users.

**Figure 3.** The method elaborated.

#### *4.1. Explanation of the Pre-Optimization Phase*

The pre-optimization of each solution is a useful tool that can help make the GA approach more efficient. In the context of this work, several techniques were implemented in order to make the GA compute better results faster. The pre-optimization phase aims at reducing the initial solution space of the algorithm, which is created every time that an instance of the problem is created according to the static and dynamic parameters of the problem, from a generalized vast solution space to, what one can call, a personalized solution space.

The personalized network is achieved by deciding on some aspects of the solution before triggering the actual search for the optimal solution, thus creating a smaller solution space. It considers the activities with spatially fixed priorities to serve as gravity points for the optimization procedure. More specifically, after the traveler provides the input of the preferences, we have an initial set of choices that the traveler initially intends to follow. This cutting of the solution space happens around the locations of activities that the travelers consider as totally or spatially fixed. This creates a subset of alternative locations of activities and charging stations. Then this subset of the initial space is searched with the GA in the main optimization phase.

For example, if an activity is fixed (Priority Label 1) in the search for an optimal solution, all other temporal and spatial choices concerning this activity can be ignored. In the same way, if the activity is temporally fixed, we include all alternative facilities (locations), but we do not permit the algorithm to search outside the time windows specified by the traveler. If the activity is spatially fixed, then we cut out the parts of the solution space that has to do with other location of facilities for the activity, and we only allow the algorithm to search for different time windows. If we have the fourth case of total flexibility, then we have to search both for alternative facilities and time windows of attendance, and we cannot pre-optimize.

#### *4.2. Calculation of the Charging Time and Deciding the Charging Scenario*

Charging time is a deterministic variable and a function of the rest of the attributes. The amount of time spent at the CS is calculated with a forward lookup of the energy requirements of the EV to cover the full tour. This calculation occurs by comparing the *SSoC*, the *FSoC*, and the remaining energy in the battery of the vehicle at every location. The candidate solutions (i.e., individuals) to the problem fully describe the locations to visit and the order of attendance to the activities. From those two attributes, we can extract the exact routes between locations of activities, as well as the kilometric distance, elevation gained and lost, and charging time at the CS.

Several aspects of the detour are modeled as independent variables and can stochastically change. Thus, we included a heuristic rule in the solution process that further directs our search. This rule defines that if the CS is within a direct distance *K* from the location of the activity after the trip to the CS, then the charging scenario is set to *en passant* charging scenario. When the CS is at a direct distance greater than *K*, then the charging scenario is set to a classic charging scenario. By setting *K* to the appropriate value, the algorithm can lead to better solutions by avoiding assigning the *en passant* scenario to cases where the location of the CS is too far away from the next location of the activity. The values of *K* can vary according to what modes are included and are available in the detour when the EV is left for charging. In our implementation, we only added walking mode for the *en passant* scenario, and we have noticed that the implementation yields satisfying results for the values of 800 to 1000 m for urban environments.

In the cases where the solution gets mutated to include the detour right at the end of the tour before the traveler goes to the desired final destination, the charging scenario is always set to the classic charging scenario. In Figure 4, a combination is given, where the heuristic rule is applied, and charging station A is utilized with the *en passant* scenario, while charging station B is utilized with the classic charging scenario. For the optimization of this initial activity chain (i.e., green markers), charging stations that fall within a distance K from the location of Activity 5, can be utilized under the en passant scenario (dotted lines) since the walking distance is considered walkable. Otherwise, the algorithm sets the charging method to the classic charging scenario.

**Figure 4.** The difference between the classic and the *en passant* charging scenarios, and the role of distance K.

#### *4.3. Availability of Charging Stations*

Regarding the availability of charging stations, a method has been elaborated in order to incorporate this parameter of the problem. We used the data provided by e-Mobi [46], which is the organization responsible for e-mobility investments (especially for deploying charging stations) in Hungary, to derive a suitable metrics that could be included in the optimization of the DACO-EV problem for travelers in Budapest. The dataset provided information about the names, locations, and plug types of chargers in Budapest. Furthermore, we had a list of charging sessions for those chargers that happened between the 23rd of February 2018 and the 12th of February 2019. For this period, we applied data mining techniques to derive average statistical properties that lead to a percentage of the occupation, in similar terms, probability for the charger to be occupied for each hour of a day. If the charger was occupied for 20 min on average in a specific hour in the examined period, then it was assumed that this charger was occupied for 33.33% on average.

These statistics were included in the method for the optimization of the DACO-EV. The main component of the method that was utilized for including the availability of CS was the evaluation operator of the GA. Every time we needed to evaluate a solution that emerged from the mutation and crossover operators of the GA, we applied the evaluation operator and checked whether the CS included in the solution is likely to be occupied. We aimed to prioritize solutions in the population that include chargers that are less likely to be occupied. The prioritization occurs only by labeling the solutions, and the corresponding charger as feasible or infeasible. Other than that, the percentage to be occupied did not have any impact on the fitness function and did not change its value. In order to avoid empirical bias in how the solutions get labeled, we implemented a labeling mechanism based on stochasticity. The method produced a random percentage in the positive continuous set of [0,100], and if that random percentage was lower than the calculated percentage of the occupation, then the solution was marked infeasible. In that regard, CS, whose rate of occupation was high were more likely to be labeled as infeasible and were not included in the population of solutions.

In Figure 5, we provide example graphs for the availability of Charging Stations (CS) as calculated by the method. The percentage of occupation for four random CS in Budapest can be seen for 24 h. The percentage of hourly occupation can significantly vary between the chargers. In the case of CS 1 and CS 3, nearly 0% occupancy rate can be derived, while in the case of CS 2 and CS 4, the occupancy rate reaches almost 50% in the morning hours.

**Figure 5.** Examples of the availability of charging stations throughout the day.

As a future step of the data mining method for the calculation of the availability of charging stations, we would include a calculation of the rate of occupation also according to the day of the week. A more enormous dataset than the one available would be more suitable for this kind of extension. For example, this extension would help the method performing better on special days like weekends, where the availability may be different from weekdays.

#### *4.4. Energy Consumption Model and the Incorporation of Elevation*

The incorporation of elevation of the road segments was a significant part of our work, and it was based on several articles [48–51], with the work of Fišer [52] being a vital starting point to our work. In his thesis, he proposed an Approximated Consumption Function:

$$\log(e, m) = \begin{cases} \kappa \cdot \delta(e) + \lambda \cdot \Delta \text{elev}\_{\mathfrak{c}\mathfrak{c}} & \text{if } \Delta \text{elev}\_{\mathfrak{c}} \ge 0 \\\quad \kappa \cdot \delta(e) + a \cdot \Delta \text{elev}\_{\mathfrak{c}\mathfrak{c}} & \text{otherwise} \end{cases} \tag{2}$$

where *g* is the consumption function, δ is a function that assigns a distance to an edge, *e* and *m* are the edge traversed and the EV model correspondingly, Δ*eleve* is the difference in elevation from the start of the edge *u* to end of it, *v.* Furthermore, κ, λ, and *a* are tuning constraints through which we can control the values of function *g* in order to yield realistic results. However, to cover the needs of our more practical approach, we moved on to introduce the Practical Consumption Function, which can be described by the following equation:

$$\mathbf{g}(\mathbf{e}, m) = \begin{cases} \text{Cons}\,\mathsf{WLTP}\_{\varepsilon, \mathsf{m}} + \mathsf{coef}\, f\_{\mathsf{pull}} \* \left( m \* \mathbf{g} \* \mathsf{Adv} \mathbf{v}\_{\varepsilon} \right), \,\mathrm{if}\,\Delta \mathrm{dev}\_{\varepsilon} \ge 0 \\\quad \mathrm{Cons}\,\mathsf{WLTP}\_{\varepsilon, \mathsf{m}} - \mathsf{coef}\, f\_{\mathsf{follow}\,\mathsf{IM}} \* \left( m \* \mathbf{g} \* \mathsf{Adv} \mathbf{v}\_{\varepsilon} \right), \,\mathrm{otherwise} \end{cases} \tag{3}$$

where *ConsWLTP* is the consumption according to EV model and the *WLTP* standard, and the Electric Vehicle database [47]. The reason for using *coe*ff*uphill* and *coe*ff*downhill* is to emulate the performance of the motor of the vehicle when going uphill and the performance of the recuperation system going downwards. The energy recuperated while braking on a flat road is ignored. In comparison to other approximation functions used for considering the consumption for EV routing, our Practical Consumption Function avoids using further tuning parameters and includes the actual interchange between the potential and kinetic energy of the vehicle. Nonetheless, the results derived from our Practical Consumption Function are realistic and very close to findings reported in the literature of EV routing. The incorporation of more advanced formulas like ones presented in [53,54] our method can be considered as the next step in our research.

In order to provide the necessary background for the use of the Practical Consumption Function, the data for the distance and elevation are derived from the OpenTripPlanner API, which bases its calculations on OpenStreetMap data for the road network and SRTM data for the elevation. The aim of this model is not to introduce an EV consumption calculation model as compared to recent approaches [53,54], but to provide a mechanism for our method calculating meaningful results considering limited computation times. The model itself is based on WLTP standard for the consumption of EV models that already includes several parameters and laboratory tests and is only extended to include the effect of elevation, for the road segments traversed in the tour of the traveler. The consideration of more complex approaches are not necessary and not within the scope of this research.

#### **5. Results**

#### *5.1. Use Cases*

The experiments were run on ten instances of problems, where the solution space can significantly differ according to different preferences of the user. While one traveler may need to attend a set of fixed activities, another may be flexible concerning the locations of the activities or the temporal schedule of the activities. This means that the first traveler's initial schedule and preferences would not allow the method to search for alternatives since there is no flexibility concerning the schedule, while the second traveler's preferences would enable the method to search through a vast solution space for other options.

For all of the cases used for our tests, the activity chains contained at minimum one out-of-the-house activity and at maximum four out-of-the-house activities. For the first six instances of problems, cases 1 to 6, we considered synthetic cases of the problem, where the locations, types of activities, and flexibility were based on real activity chains of travelers, but also randomized in some of their attributes to a certain degree. For the rest of the examined cases, cases 7 to 10, we kept the set of locations and types of activities the same as in case 6, and then we included potential combinations of flexibility flags. We consider cases 1 to 6 the balanced cases, which could be input from real-world travelers, and cases 7 to 10, the extreme cases where the computational limits of the algorithm are tested.

In our use cases for the DACO-EV, the same vehicle was considered, and the *SSoC* was set to a medium-range percentage while the *FSoC* was required to be at a higher percentage than the *SSoC*, which means that a detour for charging is needed in all cases of problems displayed. For the results presented, and regarding case-specific characteristics, and the results presented below, for the provided instance of tests, for cases 1, 2, and 5 to 10 we have four out-of-the-house activities, whereas for cases 3 and 4, three out-of-the-house activities were conducted. Finally, let us state that the results provided here are calculated based on one of the criteria that are available in the method described and in Section 3.4.2, that of travel time. Thus, the parameters in the fitness function are all set to *0* except for parameter *a* that is set equal to *1.*

#### *5.2. DACO-EV Results*

The results for the evaluation of the performance of the method in addressing the DACO-EV problems are presented. The performance of the algorithm is compared in regard to the travel time in the initial activity chain of the travelers and the resulting activity chain in both cases when precomputation is applied and when not, are given. The columns of the tables are the following:

In Table 6, the results addressing the several cases of experiments are presented.

**Table 6.** Comparisons for the runs addressing the DACO-EV—effect of precomputation on solutions and computation time.

It can be stated that the developed approach provides appropriate results that can serve EV drivers in their daily commutes. On the one hand, the overall travel time gets optimized(travel time for attending activities plus travel time in the introduced subtour), while the computation times remain at acceptable levels for a driver assistance system, that can be run by the driver for a few seconds in his or her pre-trip planning phase.

The optimized overall travel time is not reduced in all cases. This is because for the tests that we have run, in which charging is required because of the SSoC and FSoC levels, an extra charging activity or charging session was added. Although activity chains are optimized to reduce travel time, in some cases (e.g., cases 1 and 2), the travel time is raised. When precomputation is applied, the algorithm still finds only a solution with more travel time than the initial one (e.g., case 1), but it can calculate solutions that have less travel time than the initial one (e.g., case 2). It can be noticed that the travel time for the subtour, which, in most cases, is a big part of the overall travel time. A positive observation is that the computation time in balanced cases was 19.51 seconds(s) on average, which is a reasonable solution time if the method can be possibly reproduced to serve real travelers in their daily lives.

In Figures 6 and 7, and Tables 7 and 8, we present an instance of case 3 in more detail to show how an activity chain of the traveler is initially set and how it can change after applying our algorithm. In Figure 6, we can see the initial activity chain, and in Figure 7, we can see how the algorithm produces an optimal activity chain. The travel time of the tour depicted in Figure 6 is 70.75 minutes (min), while the tour depicted in Figure 7 takes 55.37 min. The figures present the spatial aspects of the initial and optimized activity chain, and the tables present the temporal attributes of the two activity chains.

**Figure 6.** Initial activity chain of the traveler for case 3.

**Figure 7.** Optimized activity chain of the traveler for case 3.

**Table 7.** Part of the input of the traveler regarding the spatial and temporal aspect of activities.

**Table 8.** Attributes of the solution regarding the visits to the activities.

Regarding the initial activity chain of the traveler, it can be noticed that it includes three out-of-the-house activities, each with different levels of flexibility. Except for the first visit, that is conducted at the university, the other two activities are spatially flexible. This fact, in combination with the significant difference between the SSoC and the FSoC that forces the algorithm to add a visit to a CS, changes the spatial attributes of the initial activity chain to a great extent. In other words, the method proposes different locations to be visited by the traveler. In this specific example, the order of the visits to the activities does not change, but a charging session is added between the visits to the locations for the second and the third activities. Finally, we can notice that while the optimized activity chain includes an extra activity that corresponds to the detour for charging, the selection of alternative locations and suitable CS enables the reduction of the overall travel time.

#### *5.3. Research Implications for EV Usage and Energy Systems*

Apart from the obvious benefits that the method presents on an individual level for travelers, such as decreasing the travel time in their daily activity chains and dealing with range anxiety, there are two main takeaways regarding the use of EVs. Firstly, on an individual level, the method can be an effective supporting tool that helps minimize consumption, distributing drivers to the suitable CS according to availability and allowing them to find possible charging opportunities according to two realistic charging scenarios. The travelers that do not own chargers at their homes can greatly benefit from this method since it allows the calculation of tours and routes, that guarantee a FSoC for the drivers when reaching their home. On a system-level, the use of such a method, although it is aimed for personal usage, can provide benefits for cities. The better distribution of drivers to charging stations is a major possible benefit, as well as, the more efficient tours that are to be followed will minimize the consumed energy and the load on the charging stations network as a whole. As shown in work by Cuchý et al. [3], and their simulation efforts based on multi-destination planning tools like the one provided in this article, those methods can allow the better alignment of supply and demand and can help to make the overall transition to electromobility smoother.

#### **6. Conclusions**

A method has been elaborated for the solution DACO-EV problem, which many EV drivers face in their daily lives. Given the initial schedule of activities of the traveler, their locations, and priorities concerning each activity, our method calculates an optimal tour that guides the travelers in a better way. At the same time, travelers do not need to face range anxiety since the solution guarantees that the EV battery level remains above the desired energy level. In the modeling, two charging scenarios were included, where the traveler can charge the EV either while staying in the EV at the CS or leaving it for charging at a CS while conducting another activity. This article describes the following contributions of the method:

Our research efforts aimed to develop an optimization model that is capable of modeling the real-world problem of optimization of daily activity chains for EV users and can serve as a driver assistance solution that reduces the range anxiety of EV users. The model itself can be used for studying the utilization of CS networks in urban environments and the study of the expected behavior of EV drivers. As discussed earlier in this article, although those methods address the problem on an individual level, they can be beneficial for the system-level, the CS network infrastructure, and the road network.

As a few next steps in this research, we propose further investigation of the incorporation of the availability of charging stations. It should be studied how solution times and the quality of solutions delivered by the algorithm change when the availability is considered according to different data mining methods. A worth mentioning research direction is the inclusion of the effect of traffic on the tours that happen on rush-hours and may get delayed travel times in regard to traffic. The utilization of this algorithm as part of a simulation framework based on activity-based modeling is a long-term research goal.

**Author Contributions:** Conceptualization, D.R. and D.E.-K.; methodology, D.R. and D.E.-K.; writing—original draft, D.R. and D.E.-K.; data curation, D.R. and D.E.-K.; software, D.R.; supervision, D.E.-K. All authors have read and agreed to the published version of the manuscript.

**Funding:** Electric Travelling project has received funding from the ERA-NET COFUND Electric Mobility Europe (EMEurope), on the basis of grant contract No. ZFF24/2018-ITM\_SZERZ concluded with the Ministry for Innovation and Technology.

**Acknowledgments:** The research reported in this paper was supported by the Higher Education Excellence Program in the frame of the Artificial Intelligence research area of Budapest University of Technology and Economics (BME FIKP-MI/SC).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article*

## **New Dispatching Paradigm in Power Systems Including EV Charging Stations and Dispersed Generation: A Real Test Case**

#### **Fabio Cazzato 1, Marco Di Clerico 1, Maria Carmen Falvo 2,\* , Simone Ferrero 1 and Marco Vivian 2**

Received: 22 January 2020; Accepted: 18 February 2020; Published: 20 February 2020

**Abstract:** Electric Vehicles (EVs) are becoming one of the main answers to the decarbonization of the transport sector and Renewable Energy Sources (RES) to the decarbonization of the electricity production sector. Nevertheless, their impact on the electric grids cannot be neglected. New paradigms for the management of the grids where they are connected, which are typically distribution grids in Medium Voltage (MV) and Low Voltage (LV), are necessary. A reform of dispatching rules, including the management of distribution grids and the resources there connected, is in progress in Europe. In this paper, a new paradigm linked to the design of reform is proposed and then tested, in reference to a real distribution grid, operated by the main Italian Distribution System Operator (DSO), *e-distribuzione*. First, in reference to suitable future scenarios of spread of RES-based power plants and EVs charging stations (EVCS), using Power Flow (PF) models, a check of the operation of the distribution grid, in reference to the usual rules of management, is made. Second, a new dispatching model, involving DSO and the resources connected to its grids, is tested, using an Optimal Power Flow (OPF) algorithm. Results show that the new paradigm of dispatching can effectively be useful for preventing some operation problems of the distribution grids.

**Keywords:** Dispatching Service; Dispersed Generation; Distribution Grid; Electrical Vehicle; Optimal Power Flow; Power Flow; Power System; Simulation Models

#### **1. Introduction**

Stopping climate change requires a revolution in the way we produce and consume energy. On the production side, the most important step is to reduce the generation from carbon and other fossil fuel to near-zero, and to increase the amount of the Renewable Energy Sources (RES). On the consumption side, the revolution requires an improvement of the efficiency and in the electrification of the energy demand. In the transport sector, these two objectives can be reached thanks to the spread of Electric Vehicles (EVs), whose Charging Stations (EVCS) power is enhancing (up to 250 kW) to allow a faster recharge [1–4]. All these changes affects the operation of the power systems, at transmission grids level in High Voltage (HV) and distribution grids level in Medium Voltage (MV) and Low Voltage (LV) [5,6]. The issues are different in the two stages of the power system: most of the RES-based power plants are connected to the distribution grids (Dispersed Generation, DG) and all the EVCS are connected to the distribution grids, thus, the first impact is on the distribution level. However, at the same time, the DG, typically RES-based, also has indirect impacts on the transmission level. some problems related to the operation of the transmission grids, due to the DG, are:

On the other hand, some problems related to the operation of the distribution grid, due to DG and EVCS are:

Many of the listed issues are due to the fact that the RES-based DG and EVCS are still managed with a "fit and forget" approach: the amount of electricity available in production or needed in consumption is put into or delivered by the system, without any limits. This type of approach is no longer feasible, especially when the rate of DG and EVs become high, as is the case in Europe.

A solution of most of the listed issues could be an active management of distribution grids and of the sources, there connected, in charge to the Distributor System Operators (DSOs), similar to the management done by the TSOs with the HV grids. For all these reasons, a reform of the electrical dispatching is in progress in Europe, including Italy [7]. The common objective of the review process in all the European Countries is to test new ways for getting the necessary resources to guarantee the reliability of the whole power system, through a new form of dispatching service. It also includes taking in new entities, such as DG, MV and LV flexible end-users (such as EVCS) and their combinations, and thus involving in the dispatching functions the DSOs.

In literature, many recent works are present on the impact on the transmission and distribution grids of RES-based DG and EVCS. Some paper are dedicated only to the impacts and possible management solution for the RES-based DG on the distribution grids [8–11]. Others are related only to the influence and possible control strategies of EVCS, in distribution grids, according to a Vehicle-to-Grid approach, in reference to a smart grid environment [11–21]. Other papers are dedicated to analyze the combination of the effect of the two (DG and EVCS) and of their suitable management in the distribution grids operation and planning [22–25]. Few papers deal with the problems related both to the distribution and transmission grids, linked to the DG and EVCS [26,27].

In this context, the originality of the present work is to propose and test a new dispatching model on a real distribution system, in line with the reform that is in progress in Europe and in Italy [7], that implies a smart management of EVCSs, coordinated with DG, mainly RES-based. Analyses and simulations are done on the new model in reference to a real distribution system, in reference to future scenarios (2030), owned and operated by the most important DSO in Italy, *e-distribuzione*.

The work is so organized. Section 2 includes a short description of new dispatching models proposed by the reform in progress. The main features of the real grid object of the analysis are in Section 3. Section 4 shows the methodologies and the assumptions for the simulations. Section 5 describes and discusses the obtained results. In Section 6, conclusions are drawn.

#### **2. Dispatching Models in the European Reform**

The Italian TSO (TERNA) is moving forward in the dispatching reform guided by the Italian Authority for electricity (ARERA) in accordance with the European policies and rules [7]. In 2013, the debate concerning new dispatch models started. The investigation was through three possible models managing the dispatching sources that exploit RES-based generators, which can be summarized as follows:

Nowadays the reforms are heading toward Model 1, where the TSO can directly draw new dispatching sources from entities connected to the distribution system. The requests of the TSO do not consider the issues that may arise in the distribution grids, such as congestions, under-voltage and over-voltage phenomena and so on; thus, the DSOs are not only suffering the unavoidable spread of DG, but also an external control, which may lead to further problems. The control of the DG and of the flexible loads, by the DSO over its system, could be the best solution to avoid technical issues (Model 2 or 3). In the meantime, to guarantee a predictable exchange profile at primary ESS, reducing the global dispatching cost could be easy in this case for the DSO and could be useful for TSO (Model 3). Obviously, this can be done in a smart grid, where suitable communication [15] and control technologies are available between the DSO and the dispersed resources. This is the reason why the most important Italian DSO, *e-distribuzione*, is investigating the opportunity to apply Model 2 or 3 instead of Model 1, which is currently tested by the Italian TSO.

#### **3.** *E-Distribuzione* **Test Case**

The study has been carried out on a real distribution grid belonging to *e-distribuzione*, which is the most important DSO in Italy, and is involved in many studies on the impact of EVCS and DG on distribution grids [28–32]. The main features of the distribution grid are summarized in Figure 1.

**Figure 1.** *e-distribuzione* distribution grid.

The system is located in Southern Italy and it is subtended to a primary Electrical Sub-station (ESS) with a 25 MVA power transformer and includes four lines. The main data on the lines are summarized in Table 1.

All the MV generators are Photovoltaic (PV) systems with a total peak power of 8.76 MW. The powers and the lines to which they are connected are detailed in Table 2.

Big consumers connected directly to MV lines are detailed in Table 3 with their available power and delivery line.

For each passive customer (industrial and residential) in MV, *e-distribuzione* has provided data about the active and reactive power profile in time, every 15 min, in reference to the four seasons and the two typical days of the week (workdays and weekends). The LV loads and generators are represented, in their respective secondary ESS, as aggregate profiles of active and reactive power, because the information about the LV network are not given. The production profile of the PV systems have been calculated in reference to the irradiation values in the analyzed area through System Advisor Model (SAM) software [33,34]. In addition, *e-distribuzione* provides the technical features of the MV grid components.

All these data, useful for the simulations, are not reported here, both for lack of space and privacy policies.

**Table 3.** MV passive customers.

#### **4. Models for Simulations**

The system has been implemented in a MATLAB environment and simulated through MATPOWER Power Flow (PF) and Optimal Power Flow (OPF) algorithms [35,36]. Twelve simulations have been analyzed according to season and the day of the week was considered. The PF and OPF run every 15 min. The simulations take into account only the active powers as control variables, because it is assumed to be a dispatching model that deals only with the control of the active power of the sources involved with a different objective, according to the reform proposals in progress in Europe and in Italy.

A first PF analysis in the present scenario was performed, to check the operation of the grid in the present and without any type of dispatching action implemented. The results show that the peak demand is reached during the summer evening. Nevertheless, the grid constraints, such as voltage and current limits, are widely respected all year long and the grid turns out to be discharged and highly suitable for new sources integration.

#### *4.1. Assumptions for Simulations*

An evaluation was carried out about possible future scenarios in 2030 considering the targets imposed by the European Commission and Italian Authority [1–3]. The future scenarios have been chosen based on the experience on the field of *e-distribuzione* authors and their expertise coming out by the management and operation of the Italian distribution grid [28]. It is well known that the forecasts are dependent on many different factors, such as investments, policies and politics of the selected nation. For these reasons, three different scenarios have been deduced according to the speed of the growth, as shown in Table 4.

**Table 4.** Percentage of growth for 2030.

Since 31 GW of PV power has to be installed at distribution level to reach 51 GW, for the LV size, the installed power of domestic PV plants should at least double and in the fast growth scenario, it should triple. Thus, the variation will be considered within 100 and 200%. Concerning the PV installed at MV, the current plants are left as they are, while new PV plants are supposed to rise in the proximity of the MV customers, affecting between 50% and 90% of them. Most of the studies highlight a sudden boom in the sale of EV between 2020 and 2025, when the technology would be mature to be worldwide adopted. Unfortunately, nowadays, Italy holds one of the last positions in the European ranking of EV readiness, hence, the forecasts are discouraging; even in the best scenario the percentage goal to contain the emission is not achieved. The Italian percentages of EV share in 2030 according to *e-distribuzione* fluctuate between 10 and 15%.

The integrated EVCSs are assumed to be:

The number of every kind of EVCS is proportional to the EVs share and the number of vehicles circulating in the area, but it is obtained following different methodologies and considering a battery capacity of 80 kWh. It is worth to note that the study is focused on how the distribution grid could be affected by the increase of installed PV power, EVCSs and their management. Thus, any kind of technical information, such as kind of socket connection, details of the recharge and battery properties, is avoided, because they do not concern this dissertation.

From the DSO point of view, the EVCS could be seen as a flexible load. Indeed, the charging stations can not only withdraw energy from the grid, as is pretty much common today, but they will be able to inject power, as is shown in many pilot projects, with so-called Vehicle to Grid (V2G) technology [13–21]. In this framework, an assumption for simulations is to consider that 25% of the batteries are not accessible, because they need to be charged.

#### *4.2. OPF Objetive Functions*

In reference to these scenarios and assumptions, different OPF analyses have been performed in reference to different objective function. In particular, to minimize the losses, to get a peak shaving function and to get a voltage drop containment.

The higher the power crossing the primary ESS, the higher are the losses in the system, and consequently, the operation cost. This is the reason why the DSOs are often concerned about the minimization of the active power flow in the primary ESS. The peak shaving concept lies on the assumption that some loads or generator can be shifted in time to avoid high losses on the grid. Any storage available in the grid could be a source to fulfil this need, since it can inject and withdraw energy whenever it is needed, respecting the capacity and maximum power constraints. The batteries of the EV connected to the distribution grid can be a useful storage for this task. In particular, the peak shaving corresponds to an energy injection (or withdraw in case of negative peak) that has to be withdrawn (or injected) within the day to allow the daily cycle of an EV battery. The proposed EV battery cycle considers that the EV drivers change the range setting of the lower limit of the State of Charge (SOC) to 50%. Furthermore, the hypothetic daily trip (9% of discharge) is divided into the route from home to office and vice versa (Figure 2).

**Figure 2.** EV battery cycle of use hypnotized in the simulations.

The energy injected to the grid has to be the 82% of the energy withdrawn during the day to have a perfect charging cycle. The energy spent during the route is seen as a loss from the grid point of view, but it is in fact the energy that makes the EVs move, thus, the most useful from the driver's point of view. Two algorithms have been carried on with two different primary goals: to reduce the evening peak and to reduce the reverse power flow. These two purposes cannot be completely decoupled from each other; indeed, as an EV battery is going to be discharge in the evening, it will unavoidably recharge during the day, flatting the negative peak and vice versa. The goal of the analysis is to evaluate the power, the energy and the percentage of EVs involved guaranteeing a given peak shaving service. This evaluation of peak shaving potential is based on thresholds that limit the positive active power, which can flow through the primary ESS; once this limit is broken, the algorithm is set to exploit the charging station to balance the grid. The evening peak set point is 4.7 MW, because it allows a realistic exploitation of the EVs batteries without requiring more power than the available one. The minimum active power threshold is set to –2.7 MW.

As the voltage drop is remarkable during summer evening, the home stations are programmed to inject power to the grid when violation occurs. An important assumption is made: the stations belong to the DSO or the EVs users are inclined to accept less power availability, the available power indeed can be modified to provide local services, only if there is a will from the battery side. The percentage of the availability of the battery to be charged or discharged and the presence of the vehicles have been considered equal to 45% of the total installed power. This analysis aims to verify whether the defined percentage of dispatchable power can guarantee the containment of the voltage drop within ±6.3%.

#### **5. Results and Discussion**

#### *5.1. PF Analysis*

A first PF analysis has been made in reference to the future scenarios, for checking the performance of the grid, without any dispatching model implemented. The results are here reported in reference to the case more stressing, which is the one for the fast growth scenario in spring and summer. The impact of increased PV systems, home stations and pole stations on the grid in terms of active power at the MV bus-bar is shown, for the fast growth scenario, in Figure 3a for the spring season and in Figure 3b for the summer season, compared to the real profile of nowadays. It is evident that the presence of the EVCS has as results an increasing of the amount of active power at MV bus-bar, only during evening hours, linked to the evening use of the home stations. During the noon hours, the prevalent effect is of increasing of the negative active power at MV bus-bar (up to 5 MW), due to the higher production by new PV plants.

**Figure 3.** Active power at MV side of the primary ESS for fast growth scenario during spring (**a**) and summer (**b**) considering the impact of Home stations and pole stations.

Figure 4 shows the trend of the active power in the primary ESS adding also the super-fast EVCS, which could require up to 750 kW, during spring season (a) and summer season (b). The effect is the same described for the case in Figure 3.

**Figure 4.** Active power at MV side of the primary ESS for fast growth scenario during spring (**a**) and summer (**b**) considering the impact of Home stations, pole stations and super-fast stations.

From the simulation results, the voltage was also checked: a drop of 6.71% was tested during summer, but it was confined within 5% during the rest of the year. Nevertheless, the grid still did not experience any congestion and it is highly suitable for further installations.

#### *5.2. OPF Analysis In Case Of Evening Peak Shaving Service*

The most interesting scenario analyzed is the most stressful, that is, the fast growth one during a summer weekday. A threshold assumed to be of 4.7 MW implicates an energy injection of 760 kWh, which stands for an exploitation of 29 home stations more than the foreseen ones. The peak power indeed reaches 800 kW (Figure 5), which means a contemporary recharge of about 111 EVs. The recharges consider the energy spent for the daily trip and they occur during the rest of the day, preferably during the time of reverse power flow. The 111 vehicles are 60% of the total EVs fleet, but 79% of the available EVs.

**Figure 5.** Active power at the MV side of the primary ESS and power injected by the EV home stations for the evening peak shaving (**a**); power available by the home stations (**b**).

#### *5.3. OPF Analysis In Case Of Reverse Power Flow Shaving*

From the energy point of view, the system will need to charge 22 EVs more to increase the energy demand of 880 kWh. On the other hand, 949 kW of peak power (Figure 6) requires 21 EVCS more, with the strong assumption that they are used at the same time. In this case, the peak power dictates the right spread of EVs which is equal to 49%. Nevertheless, this percentage is not actually a needed condition for the right functioning of the system, it should need just 42 EVs plugged at the same time. The recharge implies the discharge of part of the battery during the evening; this helps to smooth out the peak power and the voltage.

**Figure 6.** Active power at MV side of the primary ESS and power injected by the EV pole station for the reverse power flow shaving (**a**); power withdrawn by the EV pole stations (**b**).

Table 5 shows the results of the peak shaving analysis: the power and energy required by the DSO vary according to the goal set. The EV drivers' awareness is calculated as the number of EV drivers that allow the grid to use their own battery over the total number of EV drivers in the area. The easier objective to implement appears to be the noon negative peak shaving, which requires only 22% of drivers' awareness. Nevertheless, this solution implies that installing 21 EVCS is too many compared to the circulating fleet (420 more EVs should be bought).

**Table 5.** Power and energy required for peak shaving services.

The noon shaving goal instead could be reach without any economic effort; it simply depends on the EVS drivers' availability, which should reach a quite high percentage.

#### *5.4. OPF Analysis In Case Of Voltage Drop Containment*

The active power through the primary ESS has not changed a lot, except for the evening peak between 8 pm and 10 pm, when the home stations inject active power to compensate the voltage violation (Figure 7). The home stations are indeed the only source that can be exploited because it is the only one available during those hours.

**Figure 7.** Profile of the minimum voltage in the whole network (**a**); power injected by the EVs home station to contain the voltage violation (**b**).

The total power requested by the grid is 273 kW, thus, 57% of the total dispatchable power. The violation is highly contained within the threshold; this scenario requires only 37 available EVs batteries for a total of 160 kWh. The energy needed is much lower than the one that the involved EVs can provide; thus, the reduction of energy in each car will be almost derisory.

#### **6. Conclusions**

A new paradigm for the management of distribution grids, with big amount of EVCS and RES-based DG, is investigated in reference to a real distribution grid, belonging to the main Italian DSO, *e-distribuzione*. First, in reference to future scenarios of spread of RES-based power plants and EVCS, some PF analyses have been made to check the performance of the distribution grid in reference to the usual operation. Second, considering a new dispatching model involving DSO and EVCS, the opportunity to provide different local and system services are tested, using OPF algorithms.

The impact of EVCS and DG power growth, in 2030, is modelled and simulated in MATLAB environment, for three future scenarios, according to the intensity of the diffusion. After running a first PF simulation, the network detected a high peak power, given by the evening loads and a negative peak power provided by the DG plants. Thus, with OPF algorithms, the flexible management of the installed EVCS for providing services to the DSO was tested. The first algorithm was intended to cut off the voltage violations that occur if the voltage limits are at ±6.3%, as local service. The OPF algorithm tends to maximize the available power at the EVCS and acts only if a violation occurs. The study proved that it is possible to successfully contain the voltage drop within 6.3%, assuming an EVs availability of 45%. Another OPF algorithm had the goal to restrict the active power flow in the primary ESS within an upper and lower limit. This is effectively done to lessen the grid losses and provides a given value of maximum input/output power to the transmission grid, as system service. The tests have shown that the availability of EVs drivers to provide this kind of services is fundamental rather than high spread of EVs. In general, all the results show that the proposed paradigm of dispatching involving DSO and their resources (EVCS and DG) can effectively improve the distribution grid performance, to prevent its operation problems (such as under-voltages or congestions) and to support the whole power system management with system services.

**Author Contributions:** Conceptualization, M.C.F., F.C., M.D.C.; methodology, M.V., M.D.C. and S.F.; software, M.V.; validation, M.V., M.D.C. and S.F.; formal analysis, M.V. and M.C.F.; investigation, M.V.; resources, M.V. and M.C.F.; data curation, M.V., M.D.C. and S.F.; writing-original draft preparation, M.V.; writing, review and editing, M.V. and M.C.F.; visualization, M.V.; supervision, M.C.F. and F.C.; project administration, M.C.F. and F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Optimal Design for a Shared Swap Charging System Considering the Electric Vehicle Battery Charging Rate**

#### **Lingshu Zhong 1 and Mingyang Pei 2,3,\***

Received: 23 February 2020; Accepted: 1 March 2020; Published: 6 March 2020

**Abstract:** Swap charging (SC) technology offers the possibility of swapping the batteries of electric vehicles (EVs), providing a perfect solution for achieving a long-distance freeway trip. Based on SC technology, a shared SC system (SSCS) concept is proposed to overcome the difficulties in optimal swap battery strategies for a large number of EVs with charging requests and to consider the variance in the battery charging rate simultaneously. To realize the optimal SSCS design, a binary integer programming model is developed to balance the tradeoff between the detour travel cost and the total battery recharge cost in the SSCS. The proposed method is verified with a numerical example of the freeway system in Guangdong Province, China, and can obtain an exact solution using off-the-shelf commercial solvers (e.g., Gurobi).

**Keywords:** shared swap charging system; electric vehicle; operational design; battery charging rate; binary integer programming

#### **1. Introduction**

Electric vehicles (EVs) are a promising technology for reducing the environmental impacts of road transport [1] and have increased rapidly in number over the past ten years [2]. However, there are several barriers to overcome for expanding the adoption of EVs. One problem with large-scale EV adoption is the limited maximum driving range [3,4] and range anxiety [5–8], which may make it difficult to complete some long-distance tours. The other problems are high battery purchase cost [9–11] and long charging time [11–13]. To solve the problems above, an increasing number of researchers have focused on deploying EV charging systems, which will significantly shape current EV coverage [10,14]. These infrastructures can generally be divided into three categories [15]: plug-in EVs (PEVs, i.e., slow chargers and fast chargers) [16], wireless charging EVs (WCEVs, i.e., inductive charging during driving) [17,18], and swap charging EVs (SCEVs) [19,20]. Table 1 shows the comparisons between these infrastructures.

charger (22 kW) provides 75 miles of driving range in 1 h of charging. The charging rate of the rapid charger (50 kW or more) can provide up to 90 miles of driving range in 30 min. Most often used in these proposed scenarios, other scenarios are omitted due to limited space. Cost abbreviations: L—low; M—medium; H—high; VH—very high. Abbreviations: EV—electric vehicle; PEV—plug-in EV; SCEV—swap charging EV; WCEV—wireless charging EV. Data source: the Renault website of https://www.renault.co.uk/electric-vehicles/, the NIO website of https://www.nio.com/nio-power, and the WIRED website of https://www.wired.com/.

Battery swap stations (BSSs) were originally implemented by the company Better Place, which went out of business in 2013 [15]. Then, five cities in China start testing BSS technology, where they serve personal vehicles, and commercial vehicles [21]. Although BSSs best replicate the experience of existing gas stations, there are still issues that prevent their wide-scale implementation. These include battery ownership, high battery purchase cost, complicated battery standardization issues, and safety issues in the swap and charge process. Since not only do an increasing number of EV consumers expect charging approaches that include short charging times (similar to refueling their current fuel vehicles) [22], but also global economic growth means more people can afford high-cost options, the SCEV mode is becoming increasingly popular [4,19,23–25]. SCEVs are good in that they have both fast and economical charging modes [4,26,27]. As shown in Figure 1, a driver can drive into a battery swap station, and a robot replaces the depleted battery with a fully charged spare [28,29]. This swap time could be very short (e.g., less than one minute based on a report from Tesla) with further automation and refinements on the vehicle [30].

**Figure 1.** Battery swap automation and refinements on the vehicle (figure source: SUN mobility and Tesla).

Swap charging (SC) can reduce the peak consumption of electricity by centralized charging [30,31] and avoid grid overloading due to mass EV charging [32] because the empty batteries that are swapped out can be charged when electricity is cheap or demand is low. Since SCEVs are considered to be a suitable EV mode, an increasing number of studies on the SC system (SCS) have emerged worldwide [4,19,25,27,33]. The Fluence Z.E. was the first electric car enabled with battery swapping technology and deployed within the Better Place network in Israel and Denmark in 2012 [4,20,27]. Then, with the advanced SC technology, fully automatic battery swapping was even faster than refueling at gas stations. NIO proposed the smallest power swap station in the world which only took up three parking spaces [2,31]. Based on these state-off-the-art battery swapping technologies' tests, some researchers have proposed an advanced concept called shared SCS (SSCS) [31]. The SSCS is an SCS that can provide heterogeneous services and requires online reservations in advance. The SSCS has a lot of differences from the regular SCS mode, and the comparisons are shown in Table 2.

The SCS and the SSCS proposed in this paper are both used for SCEVs, which separate the batteries from the vehicles and allow the SC mode. The SSCS has a few new features, as listed below:

\* The symbol • in this table denotes that the factor is considered and symbol − denotes otherwise. \*\* Battery types: single type—fully charged battery; multiple types—varying state of charge (SOC) batteries. Cost abbreviations: L—low; M—medium. Abbreviations: BSS—battery swap station; SC—swap charging; SCS—swap charging system; SSCS—shared SCS.

#### *1.1. Literature Review*

Since public power charging infrastructure plays a critical role in EV systems [7,14,34], an increasing number of researchers have begun to focus on EV routing problems under SC technology and with the battery charging dispatch model [1,4,20,24,25,27,29,30,35], which holds promise to realize long-distance EV travel [4,20]. Here, we summarize some applications and modeling attempts to develop SC in recent years, as shown in Table 3, and the findings can be briefly synthesized as follows.

**Table 3.** Comparison between existing related SC models and our model \*.

 Thesymbolsymbol BatterytypeM—multipletypesvaryingS—single type (i.e., fully charged battery). Abbreviations: DP—dynamic programming; IP—integer programming; LRP—location routing problem; MINLP—mixed-integer nonlinearprogramming;MISOCP—mixed-integersecond-orderconeprogramming;MIP—mixed-integerprogramming;VRP—vehicleroutingproblem.

\*

Although research focusing on BSS strategies has been ongoing, the results are fragmented. Currently, an integrated way of considering the VRP, battery dispatching, and battery charging efficiency (considering the battery charging rate) has not been fully investigated. To bridge these research gaps and realize the vision of the SSCS, this paper proposes an exact approach to describe the EV routing problem and BSS battery dispatching and determine the optimal SSCS design to minimize the overall system operational cost. We formulate this problem into a binary programming model so that it can deal with the various large-scale strategy issues. This model has a binary decision variable and thus quickly solves an exact solution by off-the-shelf commercial solvers (e.g., Gurobi).

#### *1.2. Contributions*

This paper focuses on SC technology and proposes a new structured SSCS to overcome the difficulties in optimal swap battery strategies for a large number of EVs with charging requests and simultaneously considers the varying battery charging rate. The contributions of this paper are mainly three-fold.

The remainder of this paper is organized as follows. Section 2 introduces the operation characteristics, notation, and concept of the proposed SSCS. Section 3 formulates the SSCS model with alternative systems. Section 4 tests the proposed model with a numerical example in China and conducts corresponding sensitivity analyses. Finally, Section 5 provides conclusions and recommends future research directions.

#### **2. Model Description**

This section introduces the operational process of the SSCS and underlying assumptions. For the convenience of readers, we list some notation frequently used in the paper in Table 4.

Consider a set of vehicle stations I{1, ... , *I*} in space. For each vehicle station *i* ∈ I, there is a BSS. These stations can also be the ODs of vehicles. Consider a set of batteries with varying SOC Q{1, ... *Q*} that a shared BSS can provide. Let *q* ∈ Q denote the battery SOC. For each station, the number of battery types can be different. Consider a set of the vehicle trip characteristic index U{1, ... , *U*} which has a series of various travel demands (i.e., origin station *i* + *u* , destination station *i* − *u* , and the initial state of the battery charge *q*0 *u*). Let *xujq* denote whether vehicle *i* heads to station *j* and replaces the depleted battery with a well-charged battery in the state of *q* ∈ Q.

To fully understand the operation process of an SSCS, Figure 2 shows an example with shared BSS stations I = {1, ... , 5} and three types of battery SOCs *q* = {1, 2, 3}. In this figure, on each link between two stations, the segment of a different number represents the travel distance between the stations. The different combinations of colors for the stations represent the battery types they provide.

**Table 4.** Notation.

**Figure 2.** Example network with different battery mode supplies in the SSC station.