Predictive Maintenance of Bus Fleet by Intelligent Smart Electronic Board Implementing Artiﬁcial Intelligence

: This paper is focused on the design and development of a smart and compact electronic control unit (ECU) for the monitoring of a bus ﬂeet. The ECU system is able to extract all vehicle data by the on-board diagnostics-(ODB)-II and SAE J1939 standards. The integrated system Internet of Things (IoT) system, is interconnected in the cloud by an artiﬁcial intelligence engine implementing multilayer perceptron artiﬁcial neural network (MLP-ANN) and is able to predict maintenance of each vehicle by classifying the driver behavior. The key performance indicator (KPI) of the driver behavior has been estimated by data mining k-means algorithm. The MLP-ANN model has been tested by means of a dataset found in literature by allowing the correct choice of the calculus parameters. A low means square error (MSE) of the order of 10 − 3 is checked thus proving the correct use of MLP-ANN. Based on the analysis of the results, are deﬁned methodologies of key performance indicators (KPIs), correlating driver behavior with the engine stress deﬁning the bus maintenance plan criteria. All the results are joined into a cloud platform showing ﬂeet e ﬃ ciency dashboards. The proposed topic has been developed within the framework of an industry research project collaborating with a company managing bus ﬂeet.

Introduction
The controller area network (CAN) bus and on-board diagnostics (ODB) communication interfaces II are standards typically used to extract information from a vehicle [1][2][3], to control the vehicle conditions, and to deduce any anomalies by accessing the electronic control unit (ECU). The CAN system is characterized by a relatively low cost per node when compared with other information systems in the automotive bus systems [4]. The CAN standard can be integrated with mobile information systems [5] and data mining algorithms such as artificial neural networks (ANN) [6]. In [7], some researchers have provided important indications about the procedures being performed for the diagnosis and prognosis of the vehicle starting from the analysis of the ECU data, thus suggesting adopting vehicle data for maintenance plan. In [8], the integrability of the sensors in a complete diagnostic network has been demonstrated, thus suggesting ideating a smart and compact unit adaptable to different types of vehicles by using specific connectors. In this scenario, ODB II connector has been adopted in [9,10] to acquire data for the fuel consumption trend, which could therefore indirectly provide indications on the analysis of drivers behavior and in a certain way is able to predict the vehicle wear (the fuel consumption is a parameter that is a function of the acceleration and of the vehicle speed). The analysis, mainly derived from the speed and consumption indications of a bus, can also be carried out by analyzing the data acquired by global positioning systems (GPS) [11][12][13][14][15][16][17][18][19], IoT 2020, 1 181 able to track in real time the vehicle route and to estimate vehicle speed. GPS technology can be integrated with global system for mobile (GSM) [11,12], providing additional support for data transmission, for anti-theft systems [13], and for signaling emergency situations or engine failures [14]. GPS tracking is strategic also for logistics [16]. Machine learning (ML) algorithms are good candidates to support data analysis, especially for engine failures prediction [20]. On the other side, data mining could provide decision-making in fleet management, estimating economical maintenance by using k-means algorithm [21]. Other unsupervised algorithms have been adopted for predictive maintenance in the automotive sector [22]. Moreover, bus surveillance systems by camera could provide further information about violators and traffic [23], thus supporting drivers. Recently, the use of Internet of Things (IoT) and microcontroller technologies enabled automatic systems predicting fleet health and maintenance [24]. A full architecture able to implement all the facilities fleet management, including the use of external sensors and microcontrollers, has been analyzed in [25]. Different frameworks have been proposed for decision-making processes in bus fleet maintenance [26]. Concerning driver behavior, clustering techniques are suitable to classify characteristics such as safe and ecological driving [27]. Driver behavior features can also be extracted by deep learning (DL) approaches analyzing road complexity [28]. Following the state-of-the-art, are formulated the main specifications of an industry project, concerning the development of a software and hardware platform oriented on predictive maintenance and driver behavior estimation. Specifically, the research project concerns the design and development of an engineered system for the acquisition of bus fleet data and for the management of their maintenance, using predictive analysis. The project also provides the monitoring of each individual bus vehicle by means of the GPS system and the integration of surveillance system. Specifically, the main specifications of project include: • the use of data acquisition interfaces (electronic interfaces) implementing on-board diagnostics II (OBD-II) communication standard diagnostics II (OBD-II) communication standard (the data are extracted from the control units by means of scheduled procedures); • the use of a central database (MySQL) for the collection of data, which are processed by data mining and artificial intelligence algorithms (e.g., clustering, artificial neural networks), supporting the formulation of the fleet maintenance plan, based on wear prediction deriving from the analysis of vehicle data such as revolutions per minute (RPM), accelerations (throttle position), stops, refueling, fuel consumed, inconsistencies between loaded values and actually consumed volumes, etc.; • the creation of dashboards indicating the wear levels of the single bus, and the predictive scheduling of the maintenance plan based on the outputs of the artificial intelligence algorithms; • a camera module: two night vision cameras are mounted on each vehicle (front and rear view) for video streaming and for recording. The cameras, in addition to the linking of the GPS signal, can be also used to verify the driving style of the driver, and to analyze the status of the itineraries that could influence the vehicle wear for a long time period (intense traffic conditions, roads with potholes, etc.); • a GPS monitoring module allowing the tracking of all the movements and activities of each vehicle and to monitor any inefficiencies (for example excessive consumption due to an inappropriate driving style, risky driving styles due to speed limits not respected, etc.). GPS data can be processed by the data mining engine for the definition of the driver's reliability and efficiency indices, for the mapping of the activities of each individual vehicle, and for the support in predicting maintenance procedures. Figure 1 illustrates the electronic architecture related to the data on acquisition on board: the ECU is connected to the ODB II port transferring data to a raspberry Pi board; the boards trough an internet key, transmit data and images to the cloud. Data and images are then processed by a server where an artificial intelligence engine is implemented. The artificial intelligence algorithms provide as output the driver key performance indicators (KPIs) and predictive maintenance procedures (see Figure 2). IoT 2019, 2 FOR PEER REVIEW 3 Figure 1. Architecture of the data acquisition on board system. The software design (artificial intelligence engine) integrates the following modules: • Multilayer perceptron (MLP) artificial neural network (ANN) model providing prediction about vehicle wear; • k-means algorithm able to provide driver clusters indicating the correct and inappropriate behaviors.
The outputs of both the algorithms are combined to update the predictive maintenance procedure.  The software design (artificial intelligence engine) integrates the following modules:

GPS
• Multilayer perceptron (MLP) artificial neural network (ANN) model providing prediction about vehicle wear; • k-means algorithm able to provide driver clusters indicating the correct and inappropriate behaviors.
The outputs of both the algorithms are combined to update the predictive maintenance procedure. The software design (artificial intelligence engine) integrates the following modules:

GPS
• Multilayer perceptron (MLP) artificial neural network (ANN) model providing prediction about vehicle wear; • k-means algorithm able to provide driver clusters indicating the correct and inappropriate behaviors.
The outputs of both the algorithms are combined to update the predictive maintenance procedure. The paper is structured by the following discussions: • Bus protocol description;

Methodological Approaches and Experiments
The ODB II communication standard is applied to the BUS Iveco Crossway. The used standards for bus communication are SAE J1962 (HW connector) and SAE J1939 (PGN) standards, useful for the design and development of the IOT system recovering the parameters from the ECU. The OBD port has been designed to communicate with different transmission protocols, and therefore with different ECU models, which over time have been replaced with CAN. The use of OBD-II technology allows direct access to the data of the engine control unit (ECU), by means of the SAE J1939 standard.
The type of frame used for the SAE J1939 standard is the extended one, which provides the 29-bit identification (ID) field given by the sum of two sub-fields: The maximum payload size is 8 bytes. The protocol is shown in Through the SAE J1939 standard, the messages are identified as "parameter group" or PG (group of parameters) corresponding to a set of quantities belonging to the same topic or subsystem. To each PG is assigned a "parameter group number", or PGN ("Parameter Group Number"), which identifies the same dataset of information.
Two types of PGN are possible: • PG global PGN: it identifies a group of PG parameters, which are sent to all devices or broadcast.
Here the Protocol Data Unit (PDU) format, PDU specific, data page and extended data page are used for the identification of the corresponding PG. Global PGNs occur when the PDU format value is greater than or equal to 240. In fact, the PDU specific corresponds to the group extension. The format of the PDU used for this data is the second.

•
Specific PGN: are parameter group PG transmitted to particular devices (peer-to-peer). Here the PDU format, data page and the extended page are used for the identification of the corresponding PG. As for the PDU Format, it assumes a value less than 240 and the specific PDU is set to zero. The format of the PDU used for this data is therefore the first.
More details about SAE J1939 protocol are provided in Appendix A.
Predictive maintenance algorithm has been implemented by multilayer perceptron neural networks (MLP) implemented by Konstanz Miner (KNIME) open source tool based on the use of graphical user interfaces (GUIs) as blocks enable data processing [29,30]. The MLP is a feed-forward artificial neural network (ANN) model that maps sets of input vehicle data by providing output wear prediction. The MLP network is constituted by multiple nodes linked in different layers. Each node behaves as a calculus cell named neuron able to process data by means of a properly defined activation function. The MLP approach implements the backpropagation training computing the gradient of the loss function with respect to the weights of the network for a single input-output sample.
The approach used for the predictive maintenance is the definition of the training and testing dataset following the scheme of Figure 3: a first data partition is used for the model training and the last sample is adopted for the model testing. The detected data samples are stored into a MySQL table where each record contains the attributes to process.
IoT 2019, 2 FOR PEER REVIEW 5 PG is assigned a "parameter group number", or PGN ("Parameter Group Number"), which identifies the same dataset of information. Two types of PGN are possible:  PG global PGN: it identifies a group of PG parameters, which are sent to all devices or broadcast. Here the Protocol Data Unit (PDU) format, PDU specific, data page and extended data page are used for the identification of the corresponding PG. Global PGNs occur when the PDU format value is greater than or equal to 240. In fact, the PDU specific corresponds to the group extension. The format of the PDU used for this data is the second.  Specific PGN: are parameter group PG transmitted to particular devices (peer-to-peer). Here the PDU format, data page and the extended page are used for the identification of the corresponding PG. As for the PDU Format, it assumes a value less than 240 and the specific PDU is set to zero. The format of the PDU used for this data is therefore the first.
More details about SAE J1939 protocol are provided in Appendix A.
Predictive maintenance algorithm has been implemented by multilayer perceptron neural networks (MLP) implemented by Konstanz Miner (KNIME) open source tool based on the use of graphical user interfaces (GUIs) as blocks enable data processing [29,30]. The MLP is a feed-forward artificial neural network (ANN) model that maps sets of input vehicle data by providing output wear prediction. The MLP network is constituted by multiple nodes linked in different layers. Each node behaves as a calculus cell named neuron able to process data by means of a properly defined activation function. The MLP approach implements the backpropagation training computing the gradient of the loss function with respect to the weights of the network for a single input-output sample.
The approach used for the predictive maintenance is the definition of the training and testing dataset following the scheme of Figure 3: a first data partition is used for the model training and the last sample is adopted for the model testing. The detected data samples are stored into a MySQL table where each record contains the attributes to process. The dataset partition used for the MLP-ANN calculus is: 80% of training and 20% of testing. This proportion provides the best model for the analyzed dataset (low calculus error value).  The dataset partition used for the MLP-ANN calculus is: 80% of training and 20% of testing. This proportion provides the best model for the analyzed dataset (low calculus error value).
The MLP model has been tested using the dataset found in [31], where the following attributes were identified:  IoT 2019, 2 FOR PEER REVIEW 6 The MLP model has been tested using the dataset found in [31], where the following attributes were identified:  GPS time (time acquired by the GPS module);  Device time (internal clock time of the device);  Longitude (longitude of the GPS coordinate);  Latitude (latitude of the GPS coordinate);  GPS speed (measured in meters/second representing the speed of the vehicle);  Horizontal dilution of precision (horizontal error on the GPS position);  Altitude (altitude acquired by the GPS module);  Bearing (the horizontal angle between the direction and the north);  G(x) (angular velocity (degree per second) on the X axis acquired with gyroscope);  G(y) (angular velocity (degree per second) on the Y axis acquired with gyroscope);  G(z) (angular velocity (degree per second) on the Z axis acquired with gyroscope);  G calibrated (gyro calibration error);  Engine coolant temperature (temperature in °C of coolant engine liquid);  Engine RPM (angular speed of the motor shaft expressed in rpm);  Intake air temperature (temperature of the air entering the combustion chamber expressed in °C);  Engine load % (percentage of the maximum power supplied by the engine);  Mass air flow rate (flow rate of the air flow entering the combustion chamber of the engine expressed in g/s);  Throttle position manifold % (percentage of the accelerator position pressed). Figure 4 shows the statistical plots of the dataset adopted for the performance check of the MLP algorithm.    The MLP network has been implemented by the KNIME workflow of Figure 5 structured as follows: IoT 2020, 1

-
A data source "CSV Reader" block loading the bus data into a local repository (data extracted from the MySQL database); -A data pre-processing filtering attributes to select for the data processing ("Column Filter" block); -A data pre-processing block normalizing all numerical data of the filtered dataset (predictive model markup language (PMML) normalizer "Normalizer (PMML)"); -A data pre-processing partitioning data for the training and testing processing; -For the training of the MLP is adopted the efficient RProp algorithm [32,33] ("RProp MLP Learner" block constituting the training dataflow); - The "MultiLayerPerceptronPredictor" block model the MLP neural network merging the training workflow with the testing one; - The numeric score provides the mean squared error (MSE) defined as: where y i is the measured value andӯ is the predicted one; the "Excel Writer (XLS)" block writes the scoring results in an excel file.
IoT 2019, 2 FOR PEER REVIEW 7 The MLP network has been implemented by the KNIME workflow of Figure 5 structured as follows: -A data source "CSV Reader" block loading the bus data into a local repository (data extracted from the MySQL database); -A data pre-processing filtering attributes to select for the data processing ("Column Filter" block); -A data pre-processing block normalizing all numerical data of the filtered dataset (predictive model markup language (PMML) normalizer "Normalizer (PMML)"); -A data pre-processing partitioning data for the training and testing processing; -For the training of the MLP is adopted the efficient RProp algorithm [32,33] ("RProp MLP Learner" block constituting the training dataflow); -The "MultiLayerPerceptronPredictor" block model the MLP neural network merging the training workflow with the testing one; -The numeric score provides the mean squared error (MSE) defined as: where yi is the measured value and ӯ is the predicted one; -the "Excel Writer (XLS)" block writes the scoring results in an excel file. For clustering results, indicating driver behavior has been applied to the k-means algorithm [34,35] using RapidMiner tool (see the related workflow in Figure 6). For the analysis of correlations between the variables, the correlation matrix algorithm of RapidMiner tool has been adopted [36] (see the related workflow in Figure 7).  For clustering results, indicating driver behavior has been applied to the k-means algorithm [34,35] using RapidMiner tool (see the related workflow in Figure 6). For the analysis of correlations between the variables, the correlation matrix algorithm of RapidMiner tool has been adopted [36] (see the related workflow in Figure 7).
IoT 2019, 2 FOR PEER REVIEW 7 The MLP network has been implemented by the KNIME workflow of Figure 5 structured as follows: -A data source "CSV Reader" block loading the bus data into a local repository (data extracted from the MySQL database); -A data pre-processing filtering attributes to select for the data processing ("Column Filter" block); -A data pre-processing block normalizing all numerical data of the filtered dataset (predictive model markup language (PMML) normalizer "Normalizer (PMML)"); -A data pre-processing partitioning data for the training and testing processing; -For the training of the MLP is adopted the efficient RProp algorithm [32,33] ("RProp MLP Learner" block constituting the training dataflow); -The "MultiLayerPerceptronPredictor" block model the MLP neural network merging the training workflow with the testing one; -The numeric score provides the mean squared error (MSE) defined as: where yi is the measured value and ӯ is the predicted one; -the "Excel Writer (XLS)" block writes the scoring results in an excel file. For clustering results, indicating driver behavior has been applied to the k-means algorithm [34,35] using RapidMiner tool (see the related workflow in Figure 6). For the analysis of correlations between the variables, the correlation matrix algorithm of RapidMiner tool has been adopted [36] (see the related workflow in Figure 7).

Testing and Results
The system represented in Figure 1 and Figure 2 has been implemented by performing the following preliminary tests:  Verification of the correct functioning of the accelerometer and the GPS module (this check also implies the correct connection with the Raspberry input pins);  Checking of the auto-starting operation with .desktop files;  Firmware testing;  Verification of solution validity using the 4G WiFi router;  Verification of server data receipt. Figure 8a shows the testing circuit system assembling the components of the architecture of Figure 1, additionally Figure 8b illustrates the photo concerning server linking. The Raspberry is powered by the cigarette lighter socket, thanks to the car adapter, which has two USB 5V sockets at the output. The preliminary tests for the accelerometer firmware are performed with the engine off, while the preliminary tests of the Bluetooth OBD reader are executed with the engine running. When the Raspberry board receiver is turned on, it connects the testing laptop with the WiFi router. The Raspberry board is remotely controlled through the remote desktop control application. The OBDrelated test of the reading script proved the detection of the vehicle data. The data flow is enabled through the Bluetooth OBD, by timing and synchronizing the Raspberry acquisition every ten seconds. The command r = requests.post (URL_ODB_BUS, json = payload) provides the following server connection response checking the correct data flow into the json package file:  In Figure 9 is illustrated a preliminary test of acceleration data acquisition.

Testing and Results
The system represented in Figures 1 and 2 has been implemented by performing the following preliminary tests: Verification of server data receipt. Figure 8a shows the testing circuit system assembling the components of the architecture of Figure 1, additionally Figure 8b illustrates the photo concerning server linking. The Raspberry is powered by the cigarette lighter socket, thanks to the car adapter, which has two USB 5V sockets at the output. The preliminary tests for the accelerometer firmware are performed with the engine off, while the preliminary tests of the Bluetooth OBD reader are executed with the engine running. When the Raspberry board receiver is turned on, it connects the testing laptop with the WiFi router. The Raspberry board is remotely controlled through the remote desktop control application. The OBD-related test of the reading script proved the detection of the vehicle data. The data flow is enabled through the Bluetooth OBD, by timing and synchronizing the Raspberry acquisition every ten seconds. The command r = requests.post (URL_ODB_BUS, json = payload) provides the following server connection response checking the correct data flow into the json package file:

Testing and Results
The system represented in Figure 1 and Figure 2 has been implemented by performing the following preliminary tests:  Verification of the correct functioning of the accelerometer and the GPS module (this check also implies the correct connection with the Raspberry input pins);  Checking of the auto-starting operation with .desktop files;  Firmware testing;  Verification of solution validity using the 4G WiFi router;  Verification of server data receipt. Figure 8a shows the testing circuit system assembling the components of the architecture of Figure 1, additionally Figure 8b illustrates the photo concerning server linking. The Raspberry is powered by the cigarette lighter socket, thanks to the car adapter, which has two USB 5V sockets at the output. The preliminary tests for the accelerometer firmware are performed with the engine off, while the preliminary tests of the Bluetooth OBD reader are executed with the engine running. When the Raspberry board receiver is turned on, it connects the testing laptop with the WiFi router. The Raspberry board is remotely controlled through the remote desktop control application. The OBDrelated test of the reading script proved the detection of the vehicle data. The data flow is enabled through the Bluetooth OBD, by timing and synchronizing the Raspberry acquisition every ten seconds. The command r = requests.post (URL_ODB_BUS, json = payload) provides the following server connection response checking the correct data flow into the json package file:  In Figure 9 is illustrated a preliminary test of acceleration data acquisition. In Figure 9 is illustrated a preliminary test of acceleration data acquisition.  The MLP model has been checked by obtaining the parameters listed in Table 2, indicating the number of the hidden layers, the neuron number for the hidden layers and the MSE: by considering the testing dataset, very low MSE values are obtained, in the order of 10 −2 . The MSE results delineate a good error trend versus the variation of the parameters, such as the number of hidden layers and number of neurons for the hidden layers thus proving the correct choice of the algorithm used for the prediction.  An example of application of the KNIME MLP network is illustrated in Figure 7, where it is possible to observe that the predicted engine power (engine load) is higher than the measured engine load values thus predicting an accentuate engine wear. Moreover, as expected, Figure 10 illustrates a close correlation between engine RPM values and engine power. The MLP model has been checked by obtaining the parameters listed in Table 2, indicating the number of the hidden layers, the neuron number for the hidden layers and the MSE: by considering the testing dataset, very low MSE values are obtained, in the order of 10 −2 . The MSE results delineate a good error trend versus the variation of the parameters, such as the number of hidden layers and number of neurons for the hidden layers thus proving the correct choice of the algorithm used for the prediction. An example of application of the KNIME MLP network is illustrated in Figure 7, where it is possible to observe that the predicted engine power (engine load) is higher than the measured engine load values thus predicting an accentuate engine wear. Moreover, as expected, Figure 10 illustrates a close correlation between engine RPM values and engine power. In order to provide information about driver behavior, we executed the k-means algorithm fixing as K = three the number of clusters (three main driver behavior). Figure 11 illustrates the clusters by grouping the GPS speed and engine RPM parameters: the cluster indicated by the orange color (cluster 0) is representative of drivers tending to travel at low velocities and by accelerating slowly (prudent driving behavior). The drivers of the cluster indicated by the green color (cluster 1) travel with low speed but forcing the engine (high average RPM engine values), thus denoting an inefficient driving style, which could accelerate the engine wear. The cluster represented by the blue color (cluster 2) denotes drivers that mainly contribute to the vehicle wear. Figure 11. k-means analysis: clusters grouped by GPS speed and engine RPM parameters, and linear regression trend. Figure 12 represents the same clusters and the linear regression trends, denoting that the cluster 0 is characterized by a low percentage of the accelerator pressed (throttle position), and a low engine power (engine load). The cluster 1 do not press excessively on the accelerator but the engine is forced, denoting that the gears of the vehicle are not often changed. The cluster 2 indicate the engine stress due to a strong pression of the accelerator. In order to provide information about driver behavior, we executed the k-means algorithm fixing as K = three the number of clusters (three main driver behavior). Figure 11 illustrates the clusters by grouping the GPS speed and engine RPM parameters: the cluster indicated by the orange color (cluster 0) is representative of drivers tending to travel at low velocities and by accelerating slowly (prudent driving behavior). The drivers of the cluster indicated by the green color (cluster 1) travel with low speed but forcing the engine (high average RPM engine values), thus denoting an inefficient driving style, which could accelerate the engine wear. The cluster represented by the blue color (cluster 2) denotes drivers that mainly contribute to the vehicle wear. In order to provide information about driver behavior, we executed the k-means algorithm fixing as K = three the number of clusters (three main driver behavior). Figure 11 illustrates the clusters by grouping the GPS speed and engine RPM parameters: the cluster indicated by the orange color (cluster 0) is representative of drivers tending to travel at low velocities and by accelerating slowly (prudent driving behavior). The drivers of the cluster indicated by the green color (cluster 1) travel with low speed but forcing the engine (high average RPM engine values), thus denoting an inefficient driving style, which could accelerate the engine wear. The cluster represented by the blue color (cluster 2) denotes drivers that mainly contribute to the vehicle wear. Figure 11. k-means analysis: clusters grouped by GPS speed and engine RPM parameters, and linear regression trend. Figure 12 represents the same clusters and the linear regression trends, denoting that the cluster 0 is characterized by a low percentage of the accelerator pressed (throttle position), and a low engine power (engine load). The cluster 1 do not press excessively on the accelerator but the engine is forced, denoting that the gears of the vehicle are not often changed. The cluster 2 indicate the engine stress due to a strong pression of the accelerator.  Figure 11. k-means analysis: clusters grouped by GPS speed and engine RPM parameters, and linear regression trend. Figure 12 represents the same clusters and the linear regression trends, denoting that the cluster 0 is characterized by a low percentage of the accelerator pressed (throttle position), and a low engine power (engine load). The cluster 1 do not press excessively on the accelerator but the engine is forced, denoting that the gears of the vehicle are not often changed. The cluster 2 indicate the engine stress due to a strong pression of the accelerator. The plot of Figure 13 confirms the correlation between engine load and engine RPM as deduced by the MLP analysis: a similar trend is observed in Figure 8.

Discussion
Starting with the MLP-ANN prediction of the engine stress, it is possible to re-plan the maintenance schedule of each vehicle. The standard predictive maintenance plan could change based The plot of Figure 13 confirms the correlation between engine load and engine RPM as deduced by the MLP analysis: a similar trend is observed in Figure 8. The plot of Figure 13 confirms the correlation between engine load and engine RPM as deduced by the MLP analysis: a similar trend is observed in Figure 8.

Discussion
Starting with the MLP-ANN prediction of the engine stress, it is possible to re-plan the maintenance schedule of each vehicle. The standard predictive maintenance plan could change based

Discussion
Starting with the MLP-ANN prediction of the engine stress, it is possible to re-plan the maintenance schedule of each vehicle. The standard predictive maintenance plan could change based on the MLP-ANN prediction of the engine stress, which is a function of the driver behavior: the planned period to perform bus maintenance can be anticipated by predicting a high engine wear. Figure 14 illustrates a theoretical plot merging this concept. on the MLP-ANN prediction of the engine stress, which is a function of the driver behavior: the planned period to perform bus maintenance can be anticipated by predicting a high engine wear. Figure 14 illustrates a theoretical plot merging this concept. The driver style and behavior are represented in three cluster typologies by analyzing the more significant parameters such as the GPS speed, the engine RPM, the engine load, and the throttle position. Each cluster can be associated with a score (average, low and high) as KPI concerning the driver velocity, the engine stress and the driver caution (see Table 3). The driver velocity and the engine stress are also representative of the fuel consumption: The low, average and high scores are denoted by red, orange and green color, respectively.
The same scoring of Table 3 is deduced mainly by the regression line slopes of Figure 7, Figure  8, and Figure 9. We observe that the predicted results and the KPI can be normalized to the unit, thus achieving an estimation scale (the results can be expressed in percentage). A full scenario is provided by the correlation matrix analysis, providing possible correlations between the most significant variables. Figure 15 reports the correlation matrix calculus enhancing the high correlation between throttle position and engine load and between throttle position and engine RPM, additionally a moderate correlation is observed between engine RPM and GPS speed, indicating the correct gear use of the drivers. The driver style and behavior are represented in three cluster typologies by analyzing the more significant parameters such as the GPS speed, the engine RPM, the engine load, and the throttle position. Each cluster can be associated with a score (average, low and high) as KPI concerning the driver velocity, the engine stress and the driver caution (see Table 3). The driver velocity and the engine stress are also representative of the fuel consumption: The low, average and high scores are denoted by red, orange and green color, respectively.
The same scoring of Table 3 is deduced mainly by the regression line slopes of Figure 7, Figure 8, and Figure 9. We observe that the predicted results and the KPI can be normalized to the unit, thus achieving an estimation scale (the results can be expressed in percentage). A full scenario is provided by the correlation matrix analysis, providing possible correlations between the most significant variables. Figure 15 reports the correlation matrix calculus enhancing the high correlation between throttle position and engine load and between throttle position and engine RPM, additionally a moderate correlation is observed between engine RPM and GPS speed, indicating the correct gear use of the drivers. The data are sampled every second. All data are grouped for a basic daily analysis. The k-means and MLP-ANN algorithms are also able to process all the collected data for monthly and yearly estimation, thus providing criteria for predictive maintenance.
The graphical dashboards can be automatized by inserting in the main workflow the delay blocks, timing the reading [37], and by simply substituting the reading CSV blocks with Python-based script, enabling the automated reading from the MySQL database trough web services [38]. All the  The data are sampled every second. All data are grouped for a basic daily analysis. The k-means and MLP-ANN algorithms are also able to process all the collected data for monthly and yearly estimation, thus providing criteria for predictive maintenance.
The graphical dashboards can be automatized by inserting in the main workflow the delay blocks, timing the reading [37], and by simply substituting the reading CSV blocks with Python-based script, enabling the automated reading from the MySQL database trough web services [38]. All the results are collected into database system linked to a cloud platform with dashboards.
The platform provides online monitoring of the KPI, vehicle health status, fuel consumption efficiency, and driver efficiency. Figure 16 illustrates the implemented dashboard. The thresholds expressed in percentage defining the low, high and average KPI are low: 0% ÷ 40%, average: 41% ÷ 60%, and high: 61% ÷ 100%. The data are sampled every second. All data are grouped for a basic daily analysis. The k-means and MLP-ANN algorithms are also able to process all the collected data for monthly and yearly estimation, thus providing criteria for predictive maintenance.
The graphical dashboards can be automatized by inserting in the main workflow the delay blocks, timing the reading [37], and by simply substituting the reading CSV blocks with Python-based script, enabling the automated reading from the MySQL database trough web services [38]. All the results are collected into database system linked to a cloud platform with dashboards.
The platform provides online monitoring of the KPI, vehicle health status, fuel consumption efficiency, and driver efficiency. Figure 16 illustrates the implemented dashboard. The thresholds expressed in percentage defining the low, high and average KPI are low: 0% ÷ 40%, average: 41% ÷ 60%, and high: 61% ÷ 100%. In Appendix A and in Appendix B are reported more details about SAE J1939 protocols and IVECO vehicle parameters, respectively.
Recent studies oriented the research in maintenance procedures by considering the programming approach [39], by classifying the state of the vehicles [40], or by applying artificial  In Appendix A and in Appendix B are reported more details about SAE J1939 protocols and IVECO vehicle parameters, respectively.
Recent studies oriented the research in maintenance procedures by considering the programming approach [39], by classifying the state of the vehicles [40], or by applying artificial intelligence algorithms for the predictive maintenance of the only engine part [41]. The proposed research is oriented on a new concept of predictive maintenance merging procedures, artificial intelligence prediction and KPI driver efficiencies, thus providing a methodology that takes into account multiple weighted factors potentially influencing vehicle maintenance.

Conclusions
The proposed work shows how it is possible to combine IoT devices detecting bus status with the data mining algorithms simultaneously estimating engine status prediction and driver behavior. The compact and implemented electronic architecture can be applied for each vehicle characterized by ODB-II and SAE J1939 standards. Data of vehicles are transmitted in the cloud to a data mining engine performing driver KPI by defining a score using k-means clustering analysis, and by predicting engine stress through MLP-ANN algorithms. The proposed data mining models have been tested mainly with a stable dataset by providing a low MSE error, thus confirming the model accuracy. The output of the data mining algorithms allowed the establishment of criteria for the predictive maintenance, thus anticipating the maintenance in cases of predicted engine stress due also to incorrect driver behavior. The bus fleet efficiency has been estimated by considering the engine stress prediction and the driver KPI. The efficiency parameters are stored into a database system and remotely visualized by dashboards. The perspectives of the proposed research are mainly oriented on the automatic management of the maintenance of a large number of vehicles, and on the possibility to choose dynamically the drivers according to the KPI evaluation. The followed scientific approach is able to combine the predictive maintenance procedure updated by wear prediction with the driver efficiency, balancing the assignment of vehicles and drivers. The adopted self-learning MLP-ANN network is stable and can be improved if a large number of vehicles and drivers will be assigned according to the project perspectives. The proposed electronic components are adaptable to different types of vehicles. The limitations of the on-board IoT solution are mainly due to few connections to dedicate to other sensors or measuring devices. In particular, the Raspberry board has an I2C port (to which the 3D accelerometer is connected), a serial universal asynchronous receiver-transmitter (UART) port (to which the GPS module is connected) and a serial peripheral interface (SPI) port. In order to overcome this limit, have been added the four USB ports, two of which are used for the OBD device and for the internet key. To have other connection points, it is necessary to insert other boards such as Raspberry and Arduino, which also have analog inputs.

Conflicts of Interest:
There are no conflicts of Interest.

Appendix A
As has already been said, with the SAE J1939 standard, data frames have an extended structure that has an extended 29-bit identification ID field. The latter is essential since it provides information on the type of PGN, on the priority of the same message, on the address of the device to which it is sending, and the intended recipient. In particular, we have the following structure: indicates the format of the data frame since its value varies the structure of the same frame. • PF < 239: the data frame refers to specific devices with their 8-bit address. If it is encoded as (0xFF) then it is transmitted to all connected devices, i.e., broadcast. • PF ≥ 240: the data frame corresponds to a broadcast message.
• PDU specific (PS) (8 bit): this field is specific to the data frame to be transmitted, as it varies in meaning according to the value encoded in the PF. • PF < 239: the PDU specific field corresponds to the address of the device to which you want to send the data frame. • PF ≥ 240: the PDU specific field becomes "Group Extension" to form the PGN of the transmitted PG. • Source address (SD) (8 bit): it is the address of the device that is transmitting the CAN data frame.
The following example shows how CAN data frames are communicated and implemented with the SAE J1939 standard. Suppose we want to know the thermal parameters of the engine of the vehicle in question, which can be a bus or a road tractor. The group of parameters to be referred to corresponds to the PGN (0x00FEEE) entered in document J1939-71 relating to the same standard.
First of all, it is appropriate to request this PG by sending the request data frame to which the PGN is associated with broadcast. The structure of the ID is as follows: The other fields of the data frame: • DLC: set to a value of three. • Data filed: the first three bytes correspond to the PGN you want to request, which in the case in question is (0x00FEEE).
Consequently, the device used, or the heavy vehicle engine ECU, sends a response to the request made with the related group of parameters (0x00FEEE). The structure of the ID is as follows:

Appendix B
Tables A4 and A5 explain the variables of the IVECO vehicle.

Variable Description
Total_Fuel_Cons By running a count for each service delivery, it is possible to monitor the total fuel consumption. Fuel_Rate It takes into account the driving style of the driver and any losses due to idling. Inst_Fuel_Eco High consumption indicates an incorrect driving style or the occurrence of a fault. Pos_Valve It allows to monitor the status of the engine fuel system. Fuel_Level It is closely linked to the autonomy of travel or to a possible loss of fuel. Engine_Hours It is useful for monitoring ordinary maintenance actions based on a time scale. Total_RPM It is useful for monitoring ordinary maintenance actions based on the RPM scale. Total_Distance It verifies the ordinary maintenance actions based on kilometric scale. Speed It is closely related to the driver's driving style. RPM A high value means more waste of fuel and an incorrect driving style. Percent_Torque It takes into account the performance of the engine with a view to safety and failures. Percent_Engine It decreases due to engine deterioration. Temp_Cool Index of correct engine operation for road safety and guide efficiency. Temp_fuel It is useful for fire risk estimation. Temp_Oil It is adopted to control the correct functioning of the motion transmission parts. Volt_Batt It is used to check possible outages following a continuous engine shutdown.

SPN
Index of any anomalies that may require roadside assistance.


extfieldsofstudy:
Computer Science
Mathematics

text:
timeXplain -- A Framework for Explaining the Predictions of Time Series Classifiers

Modern time series classifiers display impressive predictive capabilities, yet their decision-making processes mostly remain black boxes to the user. At the same time, model-agnostic explainers, such as the recently proposed SHAP, promise to make the predictions of machine learning models interpretable, provided there are well-designed domain mappings. We bring both worlds together in our timeXplain framework, extending the reach of explainable artificial intelligence to time series classification and value prediction. We present novel domain mappings for the time and the frequency domain as well as series statistics and analyze their explicative power as well as their limits. We employ timeXplain in a large-scale experimental comparison of several state-of-the-art time series classifiers and discover similarities between seemingly distinct classification concepts such as residual neural networks and elastic ensembles.

Introduction
Artificial intelligence technology has become ubiquitous in industry and our everyday life in recent years. From dynamic pricing and recommender systems to autonomous vehicles and predictive policing, machine learning (ML) models continuously expand their influence on society. As they are entrusted with more and more responsibility, we should wonder whether we can actually trust them. ML systems must be safe, exhibit predictable and expected behavior, and must not discriminate [13]. Whether a system satisfies such fuzzy, yet indispensable criteria cannot be determined by computing a single "trustworthiness" metric. Instead, researchers are pushing towards technologies that explain the decisionmaking process of ML models. Through understanding the reasoning behind the predictions of complicated estimators, engineers are enabled to both spot flaws and build trust in their systems. The movement towards explainable artificial intelligence has been further fueled by a "right to explanation" recently proposed by the European Union [8]. The Communications of the ACM dedicated their  Fig. 1: Explanations of the predictions of two classifiers on time series #46 from the UCR GunPoint test set generated with timeXplain. The explanations show which temporal slices contribute towards (red) or against (blue) the correct classification as a "gun draw". Confirming previous literature, the shapelet transform (ST) classifier mostly focuses on the absence of a dip where the actor lowered their arm [12], while the dictionary-based WEASEL mainly focuses on the draw [17].
January 2020 issue to interpretable machine learning. In the cover article, Du, Liu, and Hu categorized the existing approaches by two major criteria: are the techniques intrinsic or post-hoc and are they global or local [5]. Intrinsic methods are specialized ML models that have interpretability built into them, see [2]. Post-hoc approaches treat the model as a black box. Also, the derived explanations can either be global, aiming at understanding the whole model, or they can be local, focusing on the prediction with respect to a specific input. All the categories have in common that one's ability to comprehend the decision-making process highly depends on the model's complexity. Understanding a linear regression might be feasible [15], but grasping the whole inner workings of a neural network with millions of weights is considered to be impossible today [13]. Time series are sequences of values whose fluctuations are associated with the passing of time. All sensor data is given this way and the pervasive use of monitoring devices in digital health [19], predictive maintenance [10], or cybersecurity [24] has led to an enormous increase of available time series data. Automated time series estimators (TSE) are applied to process this data. This subsumes all ML models that assign any kind of number to a time series and includes classification [3,7,21] as well as other value prediction tasks. With the rise of machine learning applications on time series, explanations of their behavior are much needed. We introduce timeXplain, a post-hoc local explanation framework for time series estimators.
To generate a local explanation, model-agnostic approaches often use perturbation in that they query the ML model in the vicinity of the input whose analysis they want to explain. We call this input the specimen. By slightly altering its features, one can gain insight into which of those features are relevant to the model, yielding an explanation as the result. To flesh out such an approach, one needs to define mapping functions detailing how to perturb the specimen. Ribeiro et al. pioneered work in this direction with LIME, introducing mappings for images, text, and tabular data [20]. Multiple extensions have subsequently been proposed, recently unified by Lundberg and Lee with their SHAP tool, in particular its model-agnostic explanation generator Kernel SHAP [16]. In the time series domain, this explanation approach has not realized its full potential yet. Guillemé et al. [9] gave some first mappings for time series and conducted a study suggesting that the derived explanations are understandable to human users. However, a collection of mappings covering all information available to time series estimators and a thorough discussion of the strengths and weaknesses of the applied methods are still mostly lacking. We introduce the timeXplain framework by defining mappings that build on the time domain, the frequency domain, and the statistics of a time series. Exemplary explanations generated with timeXplain are shown in Figure 1. They confirm previous literature on which fragments are relevant to shapelet transform (ST) [12] and WEASEL [17] classifiers. We further discuss the implicit assumptions made by our mappings and how they affect the generated explanations. Finally, we provide a reference implementation of our framework and employ it in a large-scale automated comparison of state-of-the-art time series classifiers. This comparison reveals surprising similarities between seemingly distinct classification concepts.
After reviewing the basics of SHAP in Section 2, the domain mappings of timeXplain are defined in Section 3. Section 4 shows how to make them robust in practice. Section 5 conducts an experimental comparison between several time series classifiers using timeXplain. We conclude the work in Section 6.

Post-hoc Local Explanations with SHAP
For SHAP [16], the black-box ML model to be interpreted is a function f : I → R from the input space I. For a classification task, we employ one such function f c for each class c, outputting the probability that the input belongs to that class. Intuitively, SHAP "disables" portions of the specimen x ∈ I, called fragments, yielding a perturbed specimen z ∈ I and then computes f (z). Doing this multiple times with different fragments explores the behavior of the prediction function f in the vicinity of the specimen. Using the gained insight, an interpretable linear model is built that approximates the original model near x. This is illustrated in Figure 2. The interpretable linear model then allows us to estimate the impact that each fragment of x has on the prediction f (x). One thus needs a way to divide x into some d fragments that can be disabled individually. This is formalized as the space of simplified inputs I = {0, 1} d . Each 1 or 0 in z ∈ I respectively stands for an activated or deactivated fragment of the specimen x. The actual deactivation of fragments is done by a mapping function h x : I → I; z → z. We require that evaluating h x on the all-ones vector 1 ∈ I recovers the unperturbed specimen h x (1) = x. Definition 1. An explanation model with respect to a mapping function h x is a linear regression model g : The set of all possible explanation models is denoted by G.
The impact vector φ acts as the post-hoc local explanation with respect to x. For any index k ∈ {1, . . . , d }, the so-called SHAP value φ k quantifies the impact of the k-th fragment on the prediction of the model f . To generate meaningful explanations, g should be faithful to the original model [20]. That is, it approximates f in the vicinity of x, observing g(z ) ≈ f (h x (z )) if z ∈ I disables only a few fragments. This is made formal by defining an inverse distance π 1 : I → R + 0 . Intuitively, the further away the simplified specimen z is from the all-ones vector 1, the smaller π 1 (z ) is. For an appropriate choice of π 1 , we refer the reader to [16]. The optimal explanation model is then Kernel SHAP learns the impact vector by training the linear regression model on randomly sampled simplified perturbed specimens z to approximate g * .

The timeXplain Mappings
While Kernel SHAP can generate model-agnostic explanations, it crucially depends on well-designed mappings for the specific application domain; a task that is left to the user. The fragments of the specimen x defined by h x need to be intuitively understandable for the impact values to be meaningful. Conversely, unsuitable mappings not adapted to the model or the data could produce misleading explanations, counteracting the whole idea of interpretability. This makes the construction of good mappings an important challenge.
We now describe the timeXplain framework and propose expressive mappings for time series data. It is well-known that time series can be seen not only as a sequence of values over time but also as a superposition of different frequencies, with the Fourier transform connecting these dual views [18]. Our mappings take both views as well as the statistics of the specimen into account.
Throughout this paper, we define the input space I = R d to be the set of all time series of length d. The value of a series z ∈ I at time t ∈ {1, . . . , d} is denoted by z t .

Time Slice Mapping
We first introduce the time slice mapping. The main idea is to partition the time series x into d connected subsequences, called slices, whose lengths differ by at To disable a slice of x, that is, a fragment of the specimen, one cannot just remove it or fill it with null values as most models can cope neither with series of varying lengths nor with null values. Instead, any slice is replaced with the corresponding slice of a so-called replacement time series r ∈ I. Six such replacement series will be presented below, each one giving rise to another variant of the time slice mapping. An exemplary perturbation (using the global mean replacement) is shown in Figure 3.

Definition 2.
Let r ∈ I be a replacement time series. The time slice mapping function h x : I → I using r for replacement yields the perturbed time series that is defined, for all t ∈ {1, . . . , d} and z ∈ I = {0, 1} d , as The time slice mapping makes the following two assumptions implicitly.
1. Feature space assumption: The ML model to be explained bases its decisions solely on the occurrence of patterns in the time domain. 2. Temporal coherence assumption: Neighboring points, that is, points in the same slice, have a similar impact on the predictions of the model.
Many models throughout the field of time series estimation seem to fulfill these assumptions. However, to the best of our knowledge, the effects of their violation on impact computation have not been discussed in the previous literature. We examine those effects in Section 4.2.
Synthetic Void Information. Finding a good replacement series r is nontrivial. We first propose to use synthetic void information, where replacing a slice of x with that of r removes all local structure so that the model cannot utilize any information of the "disabled" slice. Replacing a slice this way should also not insert any new structure so as not to distort the explanation. This requires r to carry as little information as possible. We propose five replacement options with void information. In view of their application in learning, they are defined with respect to a reference set S of time series, for example test data. A possible replacement is the zero series with r (1) t = 0 for all t. This, however, may lead to large jumps at the edges of the disabled slices if the neighboring values are very different from 0, which can be mitigated by setting r t to some average instead. We define the local mean by averaging over the t-th value in all series in S. Averaging over all values in all series gives the global mean.
Local mean: r For the next two replacement series, we employ Gaussian noise whose parameters are adapted to the series in S, again leading to a local and global variant.
The replacements are illustrated in Figure 4.
Authentic Opposing Information. Naturally, none of these five replacements fully succeed in removing all information in the disabled slices. Also, if the ML model has not seen synthetic void information in training, this approach might lead to distorted explanations. To tackle these issues, one can also use authentic opposing information as a sixth option for the replacement series r. Instead of avoiding any structure, we intentionally insert information different from that of the specimen. We accomplish this via sample replacement meaning that we draw the replacement series r (6) uniformly from the reference set S. This mapping has previously been suggested in [9]. The impact vector now shows how the model discriminates between x and the particular choice for r (6) , but this does not necessarily include all slices responsible for the model's prediction of x. We present a robust modification towards the latter in Section 4.1. Observe that authentic opposing information implicitly assumes reference set homogeneity, alleging that all series in S are somewhat similar, in order to ensure that multiple runs yield similar impacts. For a detailed discussion of this assumption see Section 4.2.

Frequency Band Mapping
Some time series estimators may not only work on the temporal structure of a time series but also examine its frequency spectrum. This might even be done by an explicit conversion of the series from time domain to the frequency domain via the Fourier transform. While the time slice mapping allows us to compute interpretable impacts of temporal slices of a time series, it fails to capture any frequency bands that may be considered impactful by some ML model. To explain the decisions of those models as well, we need to introduce a mapping that considers the frequency domain. For this frequency band mapping, we split the spectrum of the specimen into d frequency bands with quadratically increasing bandwidths. We use a quadratic scaling to provide greater resolution at the information-rich lower frequencies, while not loosing too much resolution at higher frequencies (as opposed to, say, a logarithmic scaling). Each band k ∈ {1, . . . , d } now corresponds to one fragment of the specimen. Its activation is governed by the entry z k of the simplified input. This mapping makes two assumptions, analogous to those for the time slice mapping. Firstly, the model solely considers frequency information in its decisionmaking (feature space assumption), and secondly, neighboring frequencies inside the same band have a similar impact on the predictions of the model (band coherence assumption). Both assumptions are discussed further in Section 4.2.
Synthetic Void Information. Proceeding like we did with the time slice mapping, we start with the synthetic void information approach. We propose the filter mapping that disables fragments by cutting the corresponding frequency bands from the spectrum of x using a so-called bandstop filter [18]. Any type of bandstop filter gives another variant of the filter mapping. For a fixed filter type, when the mapping function is queried with some simplified input z ∈ I , it cuts all the respective frequency bands k of the specimen x for which z k = 0. An illustration of the resulting perturbation is given in Figure 5.
The purpose of the filter is to disable a specified frequency band without introducing artifacts into the resulting spectrum (frequency domain) or the time series (time domain). We compare two types of filters from digital signal processing: elliptic and FIRLS filters. Their frequency responses are depicted in Figure 6. The elliptic filter [18] is an infinite impulse response (IIR) filter that is fast in execution and can be designed to feature a clear frequency response with sharp cutoffs. However, IIR filters are notorious for their instability [18]. The straightforward choice of an elliptic filter with the goal to disable the fragments might thus not be ideal as a general purpose solution for the design of a suitable frequency band mapping function. It may have its advantages, however, when the goal is to disable very narrow frequency bands to investigate their influence on the ML model. In general though, we recommend to employ a finite impulse response (FIR), in particular the FIRLS filters designed with least-squares optimization [18]. While a FIRLS filter may not demonstrate the same cutoff sharpness as the elliptic filter, its reliability makes it a robust choice for the automated generation of explanations for unknown specimens.
Authentic Opposing Information. As the ML model might be unable to handle zeroed-out frequency bands, we apply the approach of authentic opposing information in the frequency domain as well. We propose the patch mapping. It replaces a frequency band with that of a randomly picked patch time series r ∈ S. Once again, this approach assumes reference set homogeneity.
Let rDFT : I → C 1+ d /2 be the discrete Fourier transform with real input, see [23]. In the resulting series, the fluctuations of values are associated with the change of frequency. Coupled with the inverse transform (rDFT) -1 : C 1+ d /2 → I, we can manipulate time series in the frequency domain. Let the function κ assign to each frequency bin ω > 0 the corresponding band number κ(ω) = k.

Statistics Mapping
A time series estimator may base its decision on the mean and standard deviation of the input series. We employ the statistics mapping to detect this. It has always d = 2 fragments. If z 1 = 0, then the output series h x (z ) is a re-normalization of the specimen x such that its mean is now equal to some replacement mean; same for the standard deviation if z 2 = 0. The replacement values can be computed with respect to a reference set as shown in Section 3.1 (synthetic void information), or by using the mean and standard deviation of a randomly drawn time series (authentic opposing information). The transformation defined this way applies globally to the whole series. It thus only requires the feature space assumption, with the authentic opposing information approach additionally requiring reference set homogeneity in terms of mean and variance.

Robustness of timeXplain
Applying the above mappings without taking the application context into account may result in misleading interpretations. We highlight the most common pitfalls and how to avoid them when working with timeXplain, in particular regarding the assumptions that the various mappings require.

Multi-Run Averaging and Environment Classes
Most of the mappings above depend on a reference set. When randomly sampling from S for authentic opposing information, the resulting series may be uncharacteristic or, even worse, similar to x itself. The derived impacts are then not meaningful. We mitigate this issue by taking n distinct picks into account instead of relying on a single sample. Asking the model to discern x from several different sources of opposing information leads to the detection of more fragments that are relevant to the model's prediction in general. Still, sampling from or averaging over S introduces bias from patterns that occur frequently in S. This often happens in classification if there is a single class dominating S. At least for classification, we combat this problem by treating the classes separately. A classifier provides an individual probability function f c for each class. To derive the corresponding impact vector φ c , we average over intermediate impacts φ c,c with respect to a so-called environment class c .
To make these ideas precise, we need some additional notation. We model the process of computing the impact vectors of the explanation model g as an operator Φ(f, h x ) that takes the original ML model f and a mapping function h x as input. In turn, this mapping function is constructed from some reference set S or a (randomly drawn) time series r ∈ S, denoted by H x (S) or H x (r), respectively. For a classification task, let C denote the collection of all classes, and let S c contain the series in S with true class c ∈ C.
We benefit from the fact that most classifiers compute their class-specific models f c simultaneously. Thus, for some fixed environment class c , we can also compute the vectors φ c,c simultaneously for all c, speeding up the runtime.

Avoiding Artifacts from False Assumptions
The mappings all make some assumptions on the specimen or model. We explore what happens if they are violated in order to help the user accurately assess the merits of the explanations generated by timeXplain.
Feature Space Assumption. The feature space assumption states that the ML model bases its decisions on the properties of the specimen in the respective feature space. For us, that is the time domain, the frequency domain, or the statistics of a time series. What happens if one assumes the model to work in the time domain, and therefore applies the time slice mapping, while in actuality the model works in the frequency domain? When disabling a fragment in the time domain, for example removing some spike, the removal naturally affects the frequency domain as well, presumably by lowering the magnitude of certain frequencies. If these frequencies are impactful to the model, this affects the confidence of the model's prediction. However, timeXplain only witnesses the indirect effect of removing the spike and thus marks it to be impactful.
The end user who is inspecting the impact distribution needs to be aware that the marked fragments might not actually be the important features, but instead only artifacts. Conversely, the timeXplain user may need additional knowledge of the domain of the model to choose an appropriate mapping.
Coherence Assumption. Assuming coherence alleges that neighboring values have a similar impact on the model's prediction, both in the time and frequency domain. Intuitively, this is reasonable as a model that bases its decision on isolated points, ignoring their neighborhood, would be substantially overfitting.

Time Time
Specimen Slice edges Impactful slices Fig. 7: The effect of using too few time slices. The left plot assumes only shortreaching coherence, the two rising slopes have a high impact. The right plot assumes far-reaching coherence, resulting in a misleading impact distribution.
Contrarily, coherence spanning a longer and longer range of points or frequencies becomes increasingly implausible. Assuming such a far-reaching coherence may result in misleading explanations, if for example several meaningful peaks of a time series are assumed to be one single fragment as illustrated in Figure 7. Moreover, as the resolution of the mapping function implicitly defines how far-reaching the coherence is, the derived explanation cannot serve as aposteriori justification of the initial assumption. Those who view the resulting impact distribution have to take into account that the truly important portions of the specimen might only be small parts lying inside an impactful fragment. To mitigate these issues, we recommend to use a moderate number of time slices or frequency bands incorporating a small but non-empty neighborhood.
Simplified Input Dimension Trade-off. The choice of the dimension d of the simplified input space is also a source of artifacts in another regard. Above, we have argued that too small of a value results in low resolution and potentially misleading impacts. Intuitively, the higher we choose d , the more fine-grained the derived explanation becomes. Furthermore, the fidelity to the original model should increase. This intuition, however, is not entirely true. A high-dimensional simplified input space also makes the explanation susceptible to noise in practice.
Even for a moderately sized d , it is prohibitive to brute-force all 2 d perturbations of a specimen for a single explanation, let alone with multiple specimens and mapping functions. We thus resort to sampling. However, if one were to choose the dimension d too high and then sample only a few perturbations, the resulting explanation model g would be noisy and overfitted. Additionally, Alvarez-Melis and Jaakola observed fluctuations of SHAP caused by slight changes to the specimen [1]. They concluded that perturbation-based approaches are prone to unstable behavior. These results incentivize to choose d small enough such that a suitable number of samples for the linear regression can be collected.
In summary, the simplified input dimension being either too low or too high may both reduce the quality of the explanation. In this work, we choose d based on preliminary experiments. Finding general bounds on d is left as future work.
Reference Set Homogeneity. Experiments by Guillemé et al. [9] suggest that in the time slice mapping, sample replacements produce explanations at least as, but often more, faithful to the original ML model than constant replacements. This behavior likely extends to our proposed authentic opposing information approach in the frequency domain as well as to the statistics mapping. So, is there a justification to use synthetic void information instead?
The reference set S is likely heterogeneous to some degree, with the series exhibiting differences ranging from a time shift to completely distinct behavior. This may have an adversarial effect on the robustness of the impact computation over multiple tries, even when multi-run averaging is applied, see Section 4.1. When dealing with a very heterogeneous reference set or when otherwise observing unstable impacts, synthetic void information is a good fall-back approach. It applies the same replacement strategy over multiple tries and thus increases both robustness and reproducibility of the derived explanation.

Experimental Classifier Comparison
We employ timeXplain to compare a collection of state-of-the-art time series classifiers. We thereby show how local explanations can be applied to discover similarities in the decision-making processes of distinct classifiers.

Classifiers
We include one classifier from each paradigm outlined by Bagnall et al. [3].
1. Whole series techniques compare time series pairwise with a distance measure that operates on time series as a whole in the time domain. For our experiments, we choose the elastic ensemble EE of 1-nearest-neighbor classifiers using various distance measures, including dynamic time warping [3]. 2. Interval techniques analyze only training-selected, phase-independent intervals of a series, ignoring noisy regions and focusing on discriminatory ones. As a representative, we choose the time series forest classifier TSF [3]. 3. Shapelet techniques classify series based on the occurrence of class-specific, phase-independent subseries found during training, named shapelets. We choose shapelet transform ST due to its superior performance over other shapelet methods [3] and couple it with random forest for classification. 4. Dictionary-based techniques not only take into account the presence or absence of subseries, but also their repetition. Histograms are built from these frequency counts and then fed into a classifier. We choose the word extraction for time series classification WEASEL and combine it with a logistic regressor, as proposed by Schäfer and Leser [22].
Additionally, we include random interval spectral ensemble RISE, a classifier recently proposed by Lines et al. [11] that extracts features from the frequency domain and utilizes them to learn a time series forest. We also investigate whether the behavior of general purpose classifiers shares any resemblance with that of classifiers specific to time series. The rotation forest RotF is known to lead the field of general purpose classifiers applied to time series data [3]. Finally, we add a support vector machine with a linear kernel SVM/Lin and the residual network ResNet as the current best representative of deep learning on time series [7].

Experimental Setup
The classifiers were implemented in Python using the sktime [14], pyts [6], and tslearn [25] packages. The experiments were run on a server with four Intel Xeon Gold 5118 CPUs at 2.30 GHz, 48 physical cores in total, and 62 GB RAM. The data was taken from the public UCR time series classification archive [4]. We excluded datasets with missing values or varying time series lengths, as some models cannot cope with such series. Datasets with more than four classes were also excluded to keep the multi-run averaging costs low. From the remaining 67 datasets, 33 were picked at random. One instance of each classifier was trained on each dataset. Due to its high runtime, we limited the EE classifier to 10 randomly picked datasets. From each dataset's test set, five specimens were selected for the comparison. They were drawn at random in a way that specimen 1 is of class 1, specimen 2 of class 2, and so on, wrapping around once all classes were covered. We applied timeXplain to explain the behavior of the class-specific model for the true class of the specimen. For each specimen and each classifier, we computed two impact vectors, one with the time slice mapping with sample replacement and one with the frequency band patch mapping. The test sets served as reference sets. We employed multi-run averaging over n = 10 runs with environment classes. The simplified input dimension d was set to min{ d /5 , 30}, based on preliminary experiments. In each pass, the classifier was probed with perturbations based on 1000 random simplified inputs per specimen. Note that we chose authentic opposing information mappings since the impacts computed by such mappings are more faithful to the model given that the reference sets are homogeneous (see Section 4.2), which seems to be mostly true for the UCR test sets.  We used the following method to quantify the similarities between explanations from two different classifiers on the same specimen. Given impact vectors φ and ψ, we derived the set of correlation points defined as {(φ k , ψ k )} k∈{1,...,d } . The linear correlation of the points was evaluated via the Pearson correlation coefficient. This coefficient then served as the similarity score of the explanations. The process is illustrated in Figure 8. To compare the classifiers, we computed the pairwise correlations separately in the time domain and the frequency domain. For all 33 datasets and five specimens per dataset, the derived explanations were compared with the above method. The respective median correlations in the two domains are shown in Figure 9.

Results
We first focus on the median correlations obtained from the time slice mapping in the left part of Figure 9. Note that all medians are positive since negative correlations were outnumbered by positive ones. While most median correlations stay below 0.40, we can observe multiple distinctly pronounced median correlations between classifiers.
With a median correlation of 0.56, RISE is notably similar to TSF. This is quite surprising as TSF is a time-based interval technique while RISE operates exclusively in the frequency domain. TSF also shares a median correlation of 0.56 with EE. Even though EE works on whole series while TSF only looks at subseries, one can see why the two concepts may be related. In a similar way, the median correlation of 0.49 between EE and ST seems reasonable. More interesting is the median correlation of 0.5 between EE and ResNet, considering that the latter may learn a highly complex function, while the former only applies simple distance measures. The general purpose classifiers RotF and SVM/Lin do not share any notable median correlation with any other classifier, which seems reasonable as they interpret time series only as unstructured collections of points. WEASEL not sharing notable median correlations hints towards dictionary-based techniques being a truly unique classification concept.
The median correlations computed in the frequency domain are shown on the right side of Figure 9. These similarities are substantially higher overall compared to the time domain. We have conducted extensive manual spot checks to expound this observation. From examining several time series spectra and impacts, we conjecture that for the majority of UCR datasets, most of the information of a time series that is useful to classifiers is contained in its lower frequencies.
Note that all classifier pairs with pronounced median correlations in the time domain have above average median correlations in the frequency domain. This may be due to artifacts induced by violations of the feature space assumption, see Section 4.2. Furthermore, we cannot identify a noise floor in the frequency domain and, as such, cannot determine a group of pronounced median correlations. Still, the high median correlation of 0.82 respectively 0.73 between EE and RISE respectively RotF are surprising considering they do not appear in the time domain and neither EE nor RotF have explicit insight into the frequency domain.

Conclusion
The increasing adoption of machine learning systems in our everyday life raises the need to interpret them. We devised timeXplain, a model-agnostic local explanation framework for time series estimators. We introduced several domain mappings building on the SHAP framework. These mappings were embedded in a system of categories, depending on whether they use time slices, frequency bands, or statistics, and whether they employ synthetic void information or authentic opposing one. This categorization makes it easy to extend and adapt the framework. We also discussed common assumptions underlying the mappings and how they affect the robustness of timeXplain.
Finally, we presented a new use case for explainable artificial intelligence by employing timeXplain in an experimental comparison of the behavior of several time series classifiers. Our analysis found unexpected similarities between models which are very different in design and focus.
Future work should examine these similarities further and may also discover hidden parallels between machine learning systems in other areas of time series estimation. Environment classes may be refined, extending the concept to non-classifier models and analyzing the semantics of intermediate impacts, especially negative ones. Most importantly, to advance towards zero prior knowledge explanations, there is need for a technique that detects and eliminates artifact impacts induced by violations of the feature space assumption.