*Article* **An Ant Colony Algorithm for Improving Energy E**ffi**ciency of Road Vehicles**

#### **Alberto V. Donati 1,\*, Jette Krause 1, Christian Thiel <sup>1</sup> , Ben White <sup>2</sup> and Nikolas Hill <sup>3</sup>**


Received: 14 May 2020; Accepted: 30 May 2020; Published: 3 June 2020

**Abstract:** The number and interdependency of vehicle CO2 reduction technologies, which can be employed to reduce greenhouse emissions for regulatory compliance in the European Union and other countries, has increasingly grown in the recent years. This paper proposes a method to optimally combine these technologies on cars or other road vehicles to improve their energy efficiency. The methodological difficulty is in the fact that these technologies have incompatibilities between them. Moreover, two conflicting objective functions are considered and have to be optimized to obtain Pareto optimal solutions: the CO2 reduction versus costs. For this NP-complete combinatorial problem, a method based on a metaheuristic with Ant Colony Optimization (ACO) combined with a Local Search (LS) algorithm is proposed and generalized as the Technology Packaging Problem (TPP). It consists in finding, from a given set of technologies (each with a specific cost and CO2 reduction potential), among all their possible combinations, the Pareto front composed by those configurations having the minimal total costs and maximum total CO2 reduction. We compare the performance of the proposed method with a Genetic Algorithm (GA) showing the improvements achieved. Thanks to the increased computational efficiency, this technique has been deployed to solve thousands of optimization instances generated by the availability of these technologies by year, type of powertrain, segment, drive cycle, cost type and scenario (i.e., more or less optimistic technology cost for projected data) and inclusion of off-cycle technologies. The total combinations of all these parameters give rise to thousands of distinct instances to be solved and optimized. Computational tests are also presented to show the effectiveness of this new approach. The outputs have been used as basis to assess the costs of complying with different levels of new vehicle CO2 standards, from the perspective of different manufacturer types as well as vehicle users in Europe.

**Keywords:** CO2 reduction; multi-objective combinatorial optimization; meta-heuristics; ant colony optimization

#### **1. Introduction**

In recent years, the automotive industry has faced increasing pressure to reduce CO2 emissions to meet regulatory targets set in both the EU and other legislations. Several new CO2 reduction and energy efficiency improvement technologies have been developed to respond to these new regulatory developments. This technology proliferation increases the complexity of finding optimal combinations to achieve substantial CO2 reductions in a cost-efficient way. The problem consists of finding feasible configurations, represented by points in the plane total CO2 reduction versus price, and identifying the limiting subset of those, which form the Pareto front, to yield the optimal configurations for the objectives defined. The optimization problem is computationally difficult because of the size of the search space. Because of incompatibilities between certain technologies, the optimization problem cannot be solved by simple sorting according to technology efficiency (CO2 reduction per cost) and constructing configurations starting by the most efficient technology, adding technologies based on their efficiency, from the most efficient to the least.

In the following, we discuss similar problems found in the literature and various approaches to tackle them. The Knapsack Problem (KP) and its multi-dimensional variant (MKP) are problems aimed to find an optimal subset of objects to be placed in one or more knapsacks (or limited capacity boxes/containers), maximizing the total profit or utility. These problems belong to the class of NP-complete problems and have been extensively explored as a fundamental area of discrete optimization and Operations Research (OR). Whilst the standard KP seeks to maximize the utility under a single constraint, other variations or generalizations can deal with the assignment of resources, such as in freight logistics planning, capital budgeting, allocation of tasks on processors, stock cutting problems or in distributed computing as presented in [1,2]. A well-studied extension of the MKP is the Multi-Objective Multi-dimensional Knapsack Problem (MOMKP), when taken with two objective functions and no constraints on the number or capacities of the knapsacks. A detailed survey of techniques to solve the Multi-Objective Knapsack Problem (MOKP) and MOMKP is presented in [3]; it discusses exact methods, approximation algorithms and heuristics for the MOKP, including Simulated Annealing, Tabu Search and other heuristic methods, subdivided into evolutionary algorithms and those based on local search (LS). Other methods, such as Fuzzy Logic and Kalman Filtering have been successfully applied in detecting failures of integrated sensors/systems on vehicles where high-dimensional data is present [4]. State-of-the-art genetic algorithms have been extensively applied to the multi-objective optimization like in [5,6], to the multi-criterion optimal design [7].

The Technology Packaging Problem (TPP) is introduced here since the problem differs considerably from a MOMKP, in that, if a knapsack is assimilated to a package or configuration:


The difference described in point (3) can be particularly noticeable if one would try to approach the problem by fixing one objective and trying to optimize the other. More specifically, this could be done by partitioning the cost (the knapsack capacity) from a minimum to a maximum cost, say in 100 slices, and solve the same number of optimization problems each with the maximum cost of the slice considered. However this would result in a loss of solutions, since the optimization would end with a pareto front with exactly the same number of solutions as slices and, further, a potential loss of efficiency in having to re-compute from scratch the solution for each slice. This would exclude for example the possibility of incrementally adding objects/technologies to already found optimal configurations. The reduction of the 2-objective problem to a sliced single objective problem has some drawbacks which will be discussed in Section 4. On the other hand, we have attempted the reduction of the current problem to other similar multi objective and constrained problems, but none brought to a satisfactory formulation without a significant loss of details or changes in the original problem. Hence, the introduction of the Technology Packaging Problem (TPP). To solve this problem with two (or more) conflicting objectives, we have devised an Ant Colony Optimization (ACO) combined with a Local Search heuristic (ACO + LS).

Bio-inspired algorithms have gained consensus as being some of the most efficient and effective to deal with complex problems and constraints in various areas of Operational Research. Examples include evolutionary algorithms, swarm intelligence, distributed systems, neural networks and similar heuristic/meta-heuristics methods for NP-complete or other computationally hard problems. In particular, the choice of ACO is motivated by the recent success of swarm intelligence approaches to solve complex and dynamical problems with a system of interacting agents, which combines exploration and learning (ACO) with simple local rules (LS). These approaches give rise to a collective and distributed intelligence through a form of system memory whereby global information is encoded locally, as presented extensively in [8,9]. ACO heuristics have proven to be particularly effective when combined with appropriate local search techniques, as evidenced in [10] and [11]. They are commonly used in operational research problems such as the classic Traveling Salesman Problem (TSP), quadratic assignment and job-shop scheduling. ACO extends to applications in problems of a dynamic nature due to the algorithm's adaptability and re-scheduling capabilities [12], for example the Vehicle Routing Problem (VRP) with variable traffic conditions [13] or with unexpected events [14]. In [15] ACO is used for the Multi-Stage Flow Shop Scheduling Problem and applied for the scheduling of real factories. The introduction of multiple colonies and different pheromones update strategies has been proved effective for multi-objective problems. An example is the application to the time and space assembly line balancing problem, presented in [16], where eight different ACO architectures are introduced and compared.

One of the first applications of ant colony inspired algorithms to the MKP and the Subset Problem (SP)—a special case of the knapsack problem—was made by [17] with an Ant System, a variation of ACO incorporating two distinct pheromone updating rules: local and global. This was followed specifically for the MKP by [18,19] (where a general ACO is introduced for multi-objective optimization) and later by [20], reporting very successful results in comparison to other evolutionary algorithms. The MKP, SP and TPP are quite different from TSP-like problems, as the ordering of objects or tasks to complete is not important: for the former the solutions are combinations, yet for the latter they are permutations (ordered sequences).

Some further considerations are therefore needed in the representation of the problem and in particular regarding the pheromone encoding. Two main approaches are present in the literature: node-based, deposited on the node or edge-based, deposited on the edges connecting nodes/objects. In this study we adopted the former, as presented in [18].

The TPP is formulated in Section 2, and the ACO + LS algorithm is discussed in Section 3. A short description of the construction of so-called cost curves for vehicle CO2 emission reduction is also given. In Section 4, the ACO + LS results are compared with those obtained previously with a Genetic Algorithm (GA). The real world optimization problem instances, discussed in Section 5, consist of a set of specific technologies which are available and a set of initialization parameters; these include the type of vehicle, its size and powertrain, the year considered and the drive cycle over which the technology is evaluated, which together act to determine the presence or absence of a technology. Finally conclusions are drawn in Section 6.

The contribution of this study is in the following:


be used for each manufacturer type's fleet composition and the combined CO2 emissions of its circulating vehicles.

With the above, once ran over the entire EU vehicle fleet typology has been used as an input for future CO2 targets, so one of its most relevant contributions has been to the support of the EU policy making to check the feasibility of the 2025–2030 CO2 reduction targets. The cases presented in Section 5 are just a part of a number of simulations/optimizations ran by the authors. Most of these, and their rationale, have been omitted since they are beyond the scope of this study. Thanks to its efficiency, the system has been used extensively for scenario analysis, for example to find configurations with certain technologies which are always present (such as where these technologies have already been embedded in the vehicles) or for the analysis in which different future scenarios for the price projection or CO2 reduction potentials were to be used. Overall, the authors have run about 60,000 problem instances in the course of these analyses.

#### **2. The Optimization Problem**

The TPP with two-objectives is defined as follows. Given:


The problem is then to find, among all the possible packages/configurations, those which are Pareto optimal with respect to two objective functions chosen. Since the technologies considered here are CO*<sup>2</sup>* reduction technologies, each with a cost *ci* and CO*<sup>2</sup>* reduction potential *ri* (which represents the CO2 emission reduction percentage value that can be attained with the technology in place), the objective functions *F* and *G* can be defined as the total cost of the package, *Ck* and its total CO*<sup>2</sup>* reduction *Rk* as follow:

$$C\_k = \sum\_{t\_i \in T\_k} c\_i \tag{1}$$

$$R\_k = 1 - \prod\_{t\_i \in T\_k} (1 - r\_i) \tag{2}$$

Equation (2) is obtained by considering the cumulative, interacting benefits of the technologies. For example, if technology 1 has *r*<sup>1</sup> = 0.1 (i.e., 10%) and technology 2 has *r*<sup>2</sup> = 0.05 (i.e., 5%), the total reduction of the combined technologies is obtained by applying technology 2 only after we have applied technology 1. Since technology 1 will provide a CO2 emission of (1−*r*1), applying the second will result in a total CO2 emission of (1−*r*1) × (1−*r*2). The total reduction of applying the two technologies is therefore *R* = 1− (1−*r*1) × (1−*r*2) (hence the form of Equation (2)), which given *r*<sup>1</sup> and *r*<sup>2</sup> is 0.145 or 14.5%.

For this specific application, the optimization is aimed at the minimization of Equation (1) and the maximization of Equation (2). Then the Pareto optimality condition between any two configurations *Tn* and *Tm* (independently of their size) is that: the configuration *Tn* is said to be dominating the configuration *Tm* if one of the two following conditions holds:

if:

$$\mathbb{C}\_{\text{ll}} < \mathbb{C}\_{\text{ll}} \text{ and } R\_{\text{ll}} \ge R\_{\text{ll}} \text{ or if } \mathbb{C}\_{\text{ll}} = \mathbb{C}\_{\text{ll}} \text{ and } R\_{\text{ll}} > R\_{\text{ll}} \tag{3}$$

This means that either *Tn* has lower total cost and same or higher total CO2 reduction, or it has the same cost but higher total CO2 reduction. The configuration *Tn* is said to be dominated by *Tm* if one of the following two conditions holds:

if:

$$R\_{\text{ll}} > C\_{\text{ll}} \text{ and } R\_{\text{ll}} \le R\_{\text{ll}} \text{ or if } C\_{\text{ll}} = C\_{\text{ll}} \text{ and } R\_{\text{ll}} < R\_{\text{ll}} \tag{4}$$

If a configuration dominates another one, this cannot belong to the set of the Pareto optimal front. On the other hand, *Tn* is Pareto optimal with *Tm* if one of the two following conditions hold:

if:

*Cn* > *Cm* and *Rn* > *Rm* or if *Cn* < *Cm* and *Rn* < *Rm* (5)

Geometrically, these conditions are equivalent to consider in the (*R,C*) plane, the relative position of the two configurations: *Tm* = (*Rm, Cm*) in which to center Cartesian axes, and *Tn* = (*Rn, Cn*). If the latter falls in the lower-right quadrant, including the vertical and horizontal semi-axes, then *Tn* is dominating (condition (3)) *Tm*, if it falls in the upper-left quadrant, including semi-axes, it is dominated (condition (4)), while in all other cases they are pareto optimal (condition (5)), so unless other configurations are found to dominate one or both, they are both included in the Pareto optimal set. The set of all Pareto optimal configurations is also referred to as Pareto front.

The size of a configuration is the number of technologies it contains. An exhaustive exploration of solution space, in case of *N* distinct technologies, would imply the evaluation of (1) and (2) for all possible configurations. The number of all possible configurations of size *n*, is the combinations of *<sup>N</sup>* objects taken *<sup>n</sup>* at a time: *<sup>N</sup> n* = *<sup>N</sup>*! *n*! (*N*−*n*)! . The total search space size, or the total number of all possible configurations (without considering incompatibilities), is then given by sum of the above expression over all possible configurations' sizes:

$$\mathbf{S} = \sum\_{n=1}^{N} \binom{N}{n} = 2^N - 1 \tag{6}$$

where the second equality derives from the Binomial theorem. Even considering incompatibilities (*I*) between technologies, which tend to diminish the number of the feasible combinations, from a rough estimation of the search space with (6), is where *N* is usually between 50 and 80 and *I* is in the order of 2 or 3, is it clear that an exhaustive search is computationally impracticable.

#### **3. The ACO Algorithm**

A simple brute force approach would imply testing all possible combinations of the compatible technologies, but this is computationally infeasible, given Equation (6). Some simple heuristics were tested. Technologies were sorted by their cost efficiency (also referred to as 'technology efficiency'), defined as the CO2 reduction divided by the cost, and to add technologies to the configuration one by one as far as they were compatible with those already present. In the construction process, when an incompatible technology is encountered, it is simply skipped, and the next is examined. This procedure is then repeated, starting each time with a different technology and adding technologies until no further technology could be added.

It was soon noticed that a procedure such as this omits some configurations, since it always starts from a pre-ordering of technologies, based on cost efficiency. Including some randomness greatly improved the number of configurations found. However, a purely random search coupled with a greedy ordering showed to use significant computational power without an efficient search strategy in

the solution space. It was observed that most of the time was spent randomly exploring, rather than reinforcing good solutions and discarding certain areas.

The Ant Colony Optimization (ACO) algorithm was therefore devised. This involves the creation of an underlying graph where artificial ants can propagate and lay pheromones, a mechanism which allows the local encoding of global information, as discussed in [8–11]. The representation of the problem is the following: each node of the graph represents a technology, and two nodes are connected by an edge if, and only if, they are compatible. If an edge is present, the ant can step from a node to the next and this walk translates into a building process, where at each step, a technology is added to the configuration completed so far. Before adding a new technology, the ant must check if it is compatible with all the technologies "visited" beforehand. To this aim, the ant needs to keep an updated incompatibility list, which contains all technologies incompatible with those previously visited. To each ant step corresponds the addition of a technology to form a configuration with one more object. As the new configuration is found, it is added to the ant's set of configurations found so far, which is then composed by a collection of configurations of increasing size. The ant's walk ends when there are no more nodes to visit, that is, when no further technology can be added to the last configuration found. At this point, all the configurations found are evaluated against those stored in the best solution found so far. If they are suboptimal, some Local Search (LS) might be attempted (generally, if the configuration size is between 6 and 12). Pheromones are updated by a uniform evaporation and a deposition, for each configuration, on the edges visited by the ant, proportionally to the quality of the configuration. The quality criterion of a configuration is its overall cost efficiency, i.e., the total CO2 reduction divided by the total cost. Otherwise, if the configuration is dominating or Pareto with respect to one or more stored in the best solution, pheromones are boosted (usually by a factor 1000) and proportionally to the quality of the configuration found. The configuration(s) of the best solution, dominated by the newly found configurations, if any, are identified and removed, and the better configuration (s) found, if any, are added to the best solution. In this way, the best solution keeps improving and stores the best configurations found so far by any ant.

#### *3.1. Technical Description*

All data pertaining to technologies and costs are fed into the system. The specific instance of the optimization problem is determined by choosing one value for each of the initialization parameters, i.e., year, powertrain, vehicle segment, drive cycle, cost scenario, cost type, and whether to include off-cycle technologies (see the list of parameter values in Section 5 for more detail). These values determine which technologies are present, their cost, their effective CO2 reduction, and their incompatibility list. With these elements, since each technology constitutes a node, the underlying graph can be constructed at this stage.

As noted before, we denote the cost efficiency of a technology its CO2 reduction per unit cost, briefly called technology efficiency (or node efficiency), as compared to the total CO2 reduction *R*, divided by the total cost *C*, called the 'configuration efficiency'. Certainly, using the configuration efficiency (*CE* = *R*/*C*) instead of a Pareto selection criteria on the two objectives would have resulted in optimized configurations as well, and namely on the *CE* itself. However, we wanted to explicitly require an improvement in each of the two objectives separately, as in criteria in Equation (5), since the aim was to have a set of configurations to choose from with separate criteria or objectives. For example, finding which configurations can reach a certain CO2 reduction level and the expected impact of the ambition of certain CO2 targets on actual costs.

The ACO algorithm proceeds as follows:

#### *3.2. Solution Construction*

(1) An ant *a* is created and placed in the first node *ti*, the first technology added to the configuration: the first node is chosen randomly in 10% of the cases and probabilistically otherwise, based on its efficiency. The ant's list of incompatible technologies is updated with the list of incompatible

technologies of this node. The node *ti* to-do-list is set "done" so it will not be added again and add the configuration to the solution the ant is constructing.

	- (a) Greedy (usually set to 20% of cases): go to node *tj* with the maximum *p* (*tj*) ∝ τ (*ti*, *tj*) × *e* (*tj*) where τ (*ti*, *tj*) is the pheromone on the edge (*ti*, *tj*), and *e* (*tj*) is the efficiency of *tj*; this is favoring exploitation and following the colony.
	- (b) Probabilistic (usually set to ~80% of cases): from node *ti* chose next node *tj* with a probability *p* (*tj*) as given above; this criteria provides a mix of exploitation and exploration.
	- (c) Random (usually set to 1% of cases): the ant goes to the next node *tj* randomly; this represents pure exploration/trial.

Then to determine *tj*, a second random number is drawn in the case of criteria b. or c., while if criterion a. is selected, *tj* is already determined, by argmax (*p* (*tj*)).

	- (a) Evaporation: all the pheromones on the graph undergo a uniform evaporation, i.e., all are decreased by a constant percentage (usually 20%) by multiplying by a constant ρ less than one; pheromones are never allowed to drop below a minimum level of 0.1 to avoid entering into stagnation,
	- (b) Deposition: first determine if a configuration is a new optimal or not (either dominating or pareto). If it is not, LS can be applied, as will be discussed shortly. Pheromones on the edges traversed by the ant are incremented proportionally to the configuration quality, as follows:
		- i. if not an optimal configuration: increase by ε = SF × *R*/*C* (*R* and *C* given by Equations (1) and (2)), and with a SF is a scale factor, depending on the problem size and costs involved.
		- ii. if a better or Pareto configuration is found, increase by ε = SF × B × *R*/*C*, where B is the boosting factor, in the order of 103.

Note that since each configuration is a combination of technologies (e.g., the order of the technologies does not matter), a deposition will be also be performed on all edges between all technologies contained. For example, if we consider the configuration {*t*3,*t*23,*t*54,*t*65} then the pheromones on edges (*t*3,*t*23), (*t*3,*t*54), (*t*3,*t*65), (*t*23,*t*54), (*t*23,*t*65), (*t*54,*t*65) and their reverse edges undergo i or ii above.


The schema of the ACO + LS algorithm is reported in pseudo code in Algorithm 1, which also calls the procedures of ant propagation, LS and pheromones update, reported in Algorithm from 2 to 4 below. The values of the ant parameters (e.g., the proportions of greedy, probabilistic and random steps, as well as the proportion of initial random step, ρ, *SF* and B) given above were chosen in the algorithm tuning phase, where several tests were conducted to obtain the highest quality solutions in the shortest time.


The ant algorithm, ant propagation and pheromone update mechanisms presented in this study are typical of those formulated in the ACO-related literature. It is worth noting that the pheromones are node based given that no ordering of the technologies in a configuration is needed. The ant choice criteria have been integrated here with the problem specific objectives, which drive the immediate reward. Also Local Search algorithm is bespoke, having been developed and integrated for this specific problem. The novelty is then how the ACO and LS have been applied to the specific problem.

Regarding the Local Search, three different methods of removing the technologies were attempted: removing randomly, removing the least efficient technologies and removing those with the largest

incompatibility list. After a few tests, substituting technologies only on the basis of efficiency appeared to be the best criteria. Once the technologies have been removed, the same number of technologies (compatible with the reduced configuration) are randomly chosen and inserted. The modified configuration is then compared to the original and, if it is better, exchanged with the previous one in the ant's solution.

**Algorithm 2**. Ant propagation


**Algorithm 3**. Local Search

```
1 • Pick the ns nodes with lowest cost efficiency and remove them from Tk:
2 for r=1,..., ns, Tr = Tk-{tr1, ... , tr} (typically nLS=1 or 2)
3 • Update the incompatibility list λr for the reduced configuration Tr
4 • repeat ns times:
5 -
         Among the compatible nodes with λr pick one randomly and insert it in Tr
6 -
         Update λr
```
**Algorithm 4**. Pheromones update

```
1 • Evaporation: ∀ i,j, with i-
                                  j: τij=max(ρ∗τij,0.1)
2 where ρ is the persistency constant (0<ρ<1, typically ρ=0.8)
3
4 • Deposition: for i,j such that ti and tj∈Tk: τij=τij + ε
5 where ε is the deposition level (ε>0)
```
The Local Search is typically run for intermediate configurations sizes, usually of 6 to 12 technologies, since this is a greatly populated area in the solution space. It is performed only if the configuration is not an optimal one. It proceeds by removing one or more technologies from some chosen configuration and re-calculating the incompatibility list without the removed technologies; finally, some alternative technologies, compatible with the reduced configuration, are added. The number of technologies to switch can be chosen by the user, but as it heavily affects the efficiency, it was typically set to 2.

Finally, we note that many tests have been conducted in order to identify reasonable values for the various algorithm parameters over a wide set of problem instances which guarantee convergence to high quality solutions in timeframes on the order of 10 min. Detailed discussion of these tests is beyond the scope of this paper.

#### **4. Results and Benchmarking**

The ACO + LS algorithm presented in the previous section, was tested for a wide range of optimization scenarios and employed to calculate optimal solutions, which are the inputs for the calculation of the vehicle CO*<sup>2</sup>* emission reduction cost curves, as introduced in [21]. For cars and vans this included optimizations for all powertrains and vehicle segments, with and without off-cycle technologies, for 2025 and 2030 under the Worldwide harmonized Light vehicles Test Procedure (WLTP) test cycle (i.e., 224 optimization runs), and a number of other cases, as presented in Section 5 and in [22,23]. Moreover, the tool was employed for a number of additional computations. Each optimization was run for 10 min, yielding high quality solutions for the given algorithm parameters. By performing these runs in parallel, it was possible to complete overnight (14 h run-time) even a larger set of runs, up to 2000, taking advantage of 32 processors-threads of the Intel® Xeon® CPU E5-2687W at 3.10 GHz.

In an earlier approach, a heuristic method was applied for assessing vehicle CO2 reduction technologies and costs as presented in [24]. In order to solve the problem, a cloud of data points (i.e., tuples of CO2 emission reduction and cost values) was created from all possible combinations of different CO2 reduction technologies. This method was able to generate results for problems up to about 30 technologies. However, it proved to be inefficient for larger and constrained problems, so to handle these, a more efficient heuristic was developed, at first utilising a procedure where the configurations are determined iteratively by using a set of constraints based on incompatible technologies, with a Genetic Algorithm (GA) to create and select optimized solutions.

Genetic Algorithms were proved to effectively solve a number of optimization problems and also in particular the MKP, as reported in [25]. In recent years, such algorithms have also been developed and made available in MATLAB® as a toolbox [26].

Since the current problem has two objectives, a possible approach with a GA is to divide the cost range into a certain number of 'slices' (usually a few hundred, from a minimum to maximum pre-determined cost range). In each cost slice, the algorithm seeks configurations with the maximum CO2 reduction value, taking into account the compatibility constraints. The problem is thus reduced to a single objective optimization in each cost slice. In a preliminary study, it appeared crucial to establish an appropriate number of these slices since it heavily influences the computational time. The authors reported, after running some tests on an Intel Core i5-3340M CPU at 2.70 GHz, that a value of 200 slices yielded a reasonable computation time of 4-5 min per instance, as compared to 100 slices (about 2 min), 400 slices (about 10 min) and 1000 slices (about 30 min). This choice was also made since further increasing the number of slices, was not adding a sufficient improvement to justify the increased computational times.

The GA approach demonstrated to be able to handle problems of up to about 50 technologies. The computational work was carried out in 2015 using MATLAB® R2013a, using the function *ga* which is contained within the MATLAB® Global Optimization Toolbox to launch the optimization [27]. This function, if not specified by the user, automatically adjusts the levels of mutation and cross-over for the constrained optimization problem, as the total cost of a configuration must lie in the cost slice. The slices are created by computing, for each case, the minimum and maximum cost. The former is the minimum cost among the technologies considered in the instance problem, while the latter is the sum of the cost of all the technologies. The schema of the algorithm is given in Algorithm 5.

Although this approach is quite efficient, as the number of CO2 reduction technologies start to increase, it becomes computationally more and more inefficient and time consuming. We recall that the search space grows very quickly, as 2*N*, where *N* is the number of technologies.

**Algorithm 5**. Genetic Algorithm


In the following section, ACO + LS outcomes have been compared against the outcomes achieved with the GA approach. The GA runs are the benchmark runs, i.e., they were not repeated for each problem instance, but rather each is the best solution found over all the runs made for that instance; for some instances different slice resolutions were ran, with different times to improve solution quality. We have reported and compared with the best solution found for a particular instance.

Figure 1 shows an example of the CO2 reduction/cost-tuples resulting from both approaches, run for a typical case. The coloured (or grey) dots represent ACO + LS solutions, where the colours (or grey shade) indicate the size of the configuration. The black squares represent the GA benchmark solutions.

**Figure 1.** Comparison of CO2 reduction and cost for optimal configurations found by GA and ACO + LS. Black squares represent GA configurations while colored (or grey) dots represent ACO + LS configurations. The color (or gray shade) represent the configuration size).

For this case, we note that there are a total of 76 Pareto optimal GA configurations (points), against the 287 optimal configurations found with ACO + LS. Moreover the GA configurations are almost always dominated by the configurations found by ACO + LS. The high number of solutions found by the ACO + LS makes the Pareto front more densely populated, which is desirable when searching differentiation or alternative configurations.

We note also a side effects of the GA slicing: the best solution found in a slice, can be dominated by solutions in other slices. In other words, dividing the 2nd objective by M slices, does not guarantee that the Pareto front will have M points. In this case, pareto-optimal configurations are not found in some slices (as their existence can be seen by the ACO + LS), therefore the GA algorithm converges to sub-optimal configurations in a number of cases.

The two noticeable discontinuities are due to two high-cost technologies, usually characterized also by a high CO2 reduction; these, once added to a lower size optimal configuration, result in sudden increase of the additional cost and CO2 reduction, thus appearing as a discontinuity of the pareto front.

For a more systematic comparison, ten problems were chosen, as listed in Table 1, with different initial parameters for year, powertrain, vehicle segment and test cycle, where GA benchmark solutions were available. For Problems 2, 8 and 9, the GA was run over 1000 iterations (1 h runs).


**Table 1.** Reference table of the benchmark problems, with typical cost scenario and total cost type. See Nomenclature for abbreviations or Section 5 for extended descriptions of the parameter values.

For each problem, ten repeated ACO + LS runs were made with a runtime of 10 min and 1 h, respectively. Then ACO + LS solutions (sets of all configurations found) were compared to GA benchmark solutions in the following way.

Accordingly, to the optimality criteria defined in Equation (3), a configuration of a solution is dominated if there is at least one configuration of the second solution that is dominating it. Therefore, to have an appropriate comparison of two Pareto optimal fronts, also called best solutions found or simply solutions, *S1* and *S2*, each configuration *Tk* of the solution *S1* has to be compared to the set of configurations of the solution *S2*. An evaluation is made to check if *Tk* is dominating, dominated, additional Pareto for *S2* or equal (exactly the same cost and CO2 reduction, with 2 precision digits for costs and 8 precision digits for CO2 reduction - these values were set considering the precision of the GA benchmark configurations at hand). If it is dominating, the number of dominated configurations of *S2* is estimated.

Then the reverse is done, to evaluate, one by one, how many configurations of *S2* are dominating, dominated, additional Pareto, or equal with respect to *S1*. Since ACO + LS runs were repeated 10 times for each problem, for each run, the number of dominating, additional Pareto, equal and dominated configurations were computed for each run. Finally, their averages were estimated, along with their 95% confidence interval over the repeated runs. The runs of 10 min and 1 h are reported in Tables 2 and 3, respectively.

*Energies* **2020**, *13*, 2850

**Table 2.** Comparison of ACO + LS 10-min runs with GA Benchmark runs. Results of 10-minute ACO + LS runs, each repeated ten times for each problem, versus GA benchmark solutions. In the first part of the table, GA configurations are compared to ACO + LS, while in the second, vice-versa. The means are computed over the ten runs, and values in italics indicate the 95% confidence interval (CI) of the estimate of the mean. Boldface is used to indicate when GA has better or additional solutions than those found by ACO + LS.


*Energies* **2020**, *13*, 2850



The tables contain the comparison of the GA benchmark solutions *S*GA versus ACO + LS solutions *S*ACO+LS (first five columns), and then the reverse (last five columns). In particular, the size of the benchmark solutions is reported (in the tables, the column "configs. found", column 1), which is measure of the search efficiency of the respective algorithm: the greater the number of configurations, the greater is their diversity and thus the possibility of combing technologies/objects in different ways. It is also reported, for each problem, the mean number of configurations of *S*GA that are better than *S*ACO<sup>+</sup>LS, and the mean number of configurations of SACO<sup>+</sup>LS that are dominated by the *S*GA. Then, in the following columns 4 and 5, the mean number of additional Pareto configurations, and the mean number of equal configurations are reported (equal configurations is the only quantity which is invariant when comparing *S1* to *S2* or vice-versa). In the second part of the tables, columns 6–10, the reverse comparison is made, for the configurations of *S*ACO+LS versus *S*GA.

In the tables, it can be seen that the GA benchmark solutions contain additional optimal configurations in just a few cases to the optimal configurations found by ACO + LS (in the tables, these cases are marked in bold).

A value of 0.10 in column 2, means that, over 10 repeated runs, one GA configuration was better than some ACO + LS configuration (s); the value of 0.10 in column 3, tells that only one ACO + LS configuration was dominated by that GA configuration. In some other cases, a few additional Pareto configurations were found by GA. Finally a few other configurations are the same exact configurations as found by ACO + LS (the columns 5 and 10 are indeed the same). Most relevantly, in most cases ACO + LS configurations dominate the GA configurations (columns 7 and 8) or are additional Pareto (column 9) to those. The number of optimal configurations found by ACO+LS and GA is reported in columns 1 and 6 respectively, and the different order of magnitude of the two is noticeable. As computation time is increased from Table 2 to Table 3, the ACO + LS algorithm finds a slightly higher number of optimal configurations and convergence is increased over the ten repeated runs, as shown by the consistent decrease of the 95% confidence interval. Also the number of the ACO + LS dominated configurations shows a clear tendency to diminish further, as well as the GA additional Pareto solutions to ACO + LS.

It is possible to run the ACO + LS algorithm for many hours in an attempt to find even better solutions, or repeat the same runs to see how consistent is the convergence to the final solution. On an Intel Xenon at 3.10 GHz, it was found that approximately 10 min is suitable to find an optimized solution using ACO + LS. Given the high number of configurations found, a substantial improvement with respect to the GA is achieved. Depending on the number of available technologies of the problem instance, the ACO + LS algorithm typically completes between 30 and 60 million iterations (ants) in 10 min.

Besides the 10 problems discussed above, an indirect comparison was made between the fits found by the new approach for several thousand input datasets and cost scenarios. For each, the optimal configurations are fitted with a Levenberg-Marquardt non-linear regression algorithm using the *nlsLM* function of the minpack.lm library in R [28]. This provides a so-called cost function, the analytical form of the Pareto front, as documented and used in [22,23], which are then used for various other applications. This extensive exercise enhanced the authors' confidence in the efficiency and robustness of the ACO + LS algorithm.

Finally, as a further evaluation of the ACO + LS method, we examined the cloud of points in the CO2 reduction/cost plane pertaining to the configurations the ants visited during their search. By inspecting the cloud, the ant colony's distance to the Pareto front can be evaluated. If the colony stays close enough to the optimal solution, it has a higher likelihood of finding improved configurations and avoiding stagnation in local minima. The search process is driven by the parameters which determine the ants' behavior, particularly exploitation versus exploration settings (the relative share between greedy, probabilistic and random choices during propagation). The following Figure 2 shows this cloud for Problem 1. Each of the small points relates to a configuration found by the ants during their search. The points are a random extract of 5000 points over the last part of the run; only this small

number of points is extracted for visualization purposes, since the total number of points are several millions, for a typical run. Similar plots were also made for initial or intermediate intervals of the run, to monitor the activities of ants in the colony.

**Figure 2.** Randomly sampled intermediate configurations visited by the ants (small blue dots) along with optimal configurations found by ACO (colored or grey shades-bigger dots of the Pareto front).

The most significant ant parameters were varied and studied. We found that for runs of 10 min or greater, there were no significant differences in the colony behavior between the first and last fifth of iterations once the algorithm is properly tuned. This means that 10 min are sufficient for the ant colony to stabilize and indicates that it is likely that even less time would be required to find optimal solutions for the problems considered. In Figure 2 the optimal configurations (the Pareto front) found by ACO + LS are plotted indicating the number of technologies present in the particular configuration. Ideally, the cloud generated by ACO + LS should stay as close as possible to the optimal solution (which also changes during the optimization process) such that small variations in the ants' paths (explorations) can make them find new optimal configurations.

We note that configurations of intermediate size, i.e., 6 to 12 technologies, constitute the most populated area. Figure 2 visually shows why it was decided to run LS only on points in this region: to enhance exploitation in this highly populated area of the search space. Mathematically, this is due to the summation term in (6), where combinations are at their maximum when *n* = *N*/2.

In summary, the ACO + LS approach has been shown to be able to handle the computational complexity of the Technology Packaging Problem quite efficiently, in particular with the following benefits:


In particular, as a consequence of (3) above, the so-called technology pathways can be identified, i.e., certain points or regions of the Pareto front can be associated with different configurations that reach them in the most cost-efficient way from sub-optimal configurations, just by adding specific technologies. In other words, it is possible to gather concrete indications as to what are the most suitable technologies to add to pre-existing packages to achieve given CO2 reduction targets at a minimum cost, so that the transitions from these packages to new optimized ones can be done with ease and continuity.

We finally note that due to the novelty of the problem, we could not compare with other algorithms in the literature. This study aimed to present a method that could solve this problem most efficiently and with high quality solutions. This method could constitute a benchmark for future algorithms to solve this problem, and certainly provides clear indications that it is exploring extensively and efficiently the search space. The same methodology was later applied with minor changes to provide input for discussing CO2 reduction technologies for heavy duty vehicles, which demonstrates the flexibility of the method.

#### **5. Applications**

We applied the method to identify optimal technology packages for reducing CO2 emissions from light duty vehicles (LDVs, i.e., passenger cars and vans) [22] and later also for the heavy duty vehicles (HDVs) [23]. The approach presented in this study was employed within the analytic work supporting the impact assessments for post-2020 LDV and HDV CO2 standards in the EU.

A dataset of more than 80 LDV CO2 emission reduction technologies, including their reduction potentials, costs and mutual compatibilities was provided in [29]. Similarly, an extensive study was carried out for HDVs with data in [30]. For each problem, once the optimal configurations are found with the ACO + LS, a so-called cost curve can be found by fitting, and it describes the Pareto front of costs versus CO2 reductions in a continuous analytical form.

A number of values for the initialization parameters exist, which define a specific instance of the optimization problem. For the LDVs these are namely (a list of abbreviations is also available in Nomenclature Section):

	- Spark Ignition (i.e., gasoline) Internal Combustion Engine and Hybrid (SI ICE + HEV)
	- Compression Ignition (i.e., diesel) Internal Combustion Engine and Hybrid (CI ICE + HEV)
	- Spark Ignition Plug-in Hybrid (SI PHEV)
	- Compression Ignition Plug-in Hybrid (CI PHEV)
	- Spark Ignition Range-Extended Electric Vehicle (SI REEV)
	- Compression Ignition Range-Extended Electric Vehicle (CI REEV)
	- Battery Electric Vehicle (BEV)
	- Fuel Cell Electric Vehicle (FCEV)
	- Direct cost, i.e., manufacturing costs
	- Total costs, i.e., direct costs plus indirect costs, where indirect costs represent additional costs related to a technology such as R&D, warranties, marketing, etc.

The possible combinations of all these parameters give rise to a total combination of 6048 different instances. For the sake of most of applications, we set the cost type equal to total cost, which results in a total of 3024 possible optimization problem instances to be run for each scenario. This is quite a computational effort: if ran in sequence they correspond to about 21-day computation time for 10 min runs, or 7 days if only interested in WLTP cycle.

It is to be noted that the set of available CO2 reduction technologies are updated infrequently, depending on the availability of new data on technologies and their costs as well as on policy needs. Therefore, the addition of new technologies doesn't constitute any computational issue, as well as the variation of their costs, having been taken into consideration via 3 technology cost scenarios in the construction of the Pareto front for modelling the Cost Curves.

#### *Cost Curve Construction*

On the basis of the optimal configurations found by the ACO optimization described in the previous sections, vehicle emission reduction cost curves can be constructed to interpolate analytically these data points. The construction of cost curves consists in the fitting of the Pareto front configurations/points. This fitting procedure was a considerable part of the work, and it is worth noticing that all the Pareto fronts were fitted with an analytical form deriving by one single generalized functional form, with at most 4 fitting parameters. The cost curves are an essential input for applications such as the evaluation of different CO2 reduction scenarios, e.g., for calculating the costs associated with a certain desired CO2 standard for vehicles, identifying cost-minimizing distributions of CO2 reduction efforts across different vehicle types and segments.

A number of post-processing steps had to be applied to the raw output from the optimization procedure, to take into account some necessary adjustments. The post-processing steps, the cost curve fitting procedure and the resulting cost curves are documented in [22] for LDVs, and in [23] for HDVs.

We have conducted extensive computational campaigns to address a number of variations and hypotheses for scenario analysis, considering variations of the original problem instances where for example a certain technology is desired to be always present, or in another case, where a certain technology could be considered to be less effective (in term of CO2 reduction) than expected or some other more costly, to derive the optimal configurations and related cost curves. These results, obtained in relative short time scales, were made possible thanks to the properties of fast convergence and efficiency of the ACO + LS methodology.

#### **6. Conclusions**

We have presented a method based on a metaheuristic with Ant Colony Optimization combined with a Local Search algorithm to efficiently solve the Technology Packaging Problem, of finding feasible configurations of CO2 reduction technologies, maximizing the CO2 reduction and minimizing costs.

The results presented show that this methodology provides solutions of increased quality and size versus other approaches, e.g., genetic algorithms, thus offering the possibility of finding improved configurations in a highly efficient and easily adaptable manner. This has enhanced the capability to explore optimal strategies under increasingly numerous and complex technology choices against various CO2 reduction targets. The method can easily accommodate more technology options or perform analyses of different technology pathways. As the runtime is greatly improved, it has been employed to run multiple scenario analyses, including e.g., variations in technology impact, costs, and test cycles.

This method has been employed for supporting the impact assessments for post-2020 CO2 standards for cars, vans and trucks in the EU, considering diverse road vehicle fleet compositions. Different scenarios for the post-2020 EU CO2 emission standards have been evaluated by constructing vehicle emission reduction cost curves as described in the present paper, and, building on them, identifying an optimal distribution of efforts among different powertrains and segments which compose the fleet, calculating additional manufacturing costs as well as fuel savings, and computing total additional costs or savings from emission reduction. Further to applications in policy support, this methodology could also be applied by the automotive industry itself.

**Author Contributions:** Conceptualization: A.V.D., J.K. and C.T.; methodology: A.V.D., B.W. and N.H.; software: A.V.D., B.W. (for GA); validation: J.K. and C.T.; formal analysis: A.V.D.; investigation: N.H.; data curation: N.H. and B.W.; writing—original draft preparation: A.V.D.; writing—review and editing: all authors; visualization: A.V.D.; supervision: J.K. and C.T.; project administration: C.T.; funding acquisition: N.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported through the Horizon 2020 direct funding of the Joint Research Centre (JRC) of the European Commission, for the JRC authors. The work of the other authors was funded by the European Commission under the framework contract CLIMA.C.2/FRA/2012/0006. The views expressed are purely those of the authors and may not in any circumstances be regarded as stating an official position of the European Commission.

**Acknowledgments:** The authors wish to thank Yannis Drossinos (JRC), for the valuable comments and suggestions.

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

#### **Nomenclature: List of Abbreviations**


#### **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*
