# Feature Reduction Method Comparison Towards Explainability and Efficiency in
Cybersecurity Intrusion Detection Systems††thanks: Research in large conducted
thanks to funding from National Security Agency grant H98230-22-1-0017 and
National Science Foundation award 1719498.

Adam Lehavi  Viterbi School of Engineering
University of Southern California
Los Angeles, CA 90089
<EMAIL_ADDRESS>Seongtae Kim  Department of Mathematics and Statistics
North Carolina Agricultural and Technical State University
Greensboro, NC 27411
<EMAIL_ADDRESS>
###### Abstract

In the realm of cybersecurity, intrusion detection systems (IDS) detect and
prevent attacks based on collected computer and network data. In recent
research, IDS models have been constructed using machine learning (ML) and
deep learning (DL) methods such as Random Forest (RF) and deep neural networks
(DNN). Feature selection (FS) can be used to construct faster, more
interpretable, and more accurate models. We look at three different FS
techniques; RF information gain (RF-IG), correlation feature selection using
the Bat Algorithm (CFS-BA), and CFS using the Aquila Optimizer (CFS-AO). Our
results show CFS-BA to be the most efficient of the FS methods, building in
55% of the time of the best RF-IG model while achieving 99.99% of its
accuracy. This reinforces prior contributions attesting to CFS-BA’s accuracy
while building upon the relationship between subset size, CFS score, and RF-IG
score in final results.

##  I Introduction

Cybersecurity is a growing field with increasing necessity and prominence. In
2022 thus far, there have been 1.364 billion identified malware programs[1].
The average cost of a data breach is $4.24 million, with an 80% cost
difference where security AI was deployed through either Machine Learning (ML)
or Deep Learning (DL)[2]. Threat detection, or finding malicious activity, is
one of the largest components of the cybersecurity field. An intrusion
detection system (IDS) is a model that detects these threats through various
means [3]. A network-based IDS (NIDS) monitors network connections to look for
malicious traffic [4]. As networks are typically the source of more damaging
and costly problems, spanning company and organization data, building an NIDS
has importance in current research, and will be the focus of this paper.

Intrusion detection systems are typically divided into two separate classes
based upon their approach and the type of attacks they aim to cover.
Signature-based IDS relies on known attacks; it generates patterns from past
or given data and then uses those set patterns to sift through current data,
such as an antivirus package[4]. Anomaly-based IDS aims to find dynamic
patterns to group data and search for deviations [4]. Hybrid approaches
combine both methods. While signature-based systems have guarantees on
detecting common malicious activity, they perform poorly on attacks that
deviate from those patterns. In addition, they often take large amounts of
labeled data beforehand, and need constant updating. As such, anomaly-based
detection can achieve better results for large amounts of data with unknown
correlations.

NIDS are fundamentally built upon data from the network. In research, data to
test particular models and methods comes from datasets of real or simulated
network data. The most common datasets include NSL-KDD, KDD-Cup’99, UNSW-NB15,
CIC-IDS2017, Kyoto, and CSE-CIC-IDS2018 [5]. Of these sets, CSE-CIC-IDS2018 is
the most current, meaning it is the most applicable. We therefore selected
this dataset to be able to add a point of comparison with needed metrics [6].
The dataset is useful for practice, because it contains a variety of attacks
including zero-day attacks, or ones that typically happen when a network is
initially set up and open to user activity [4]. Since it is not as common as
CIC-IDS2017 and lacks many results, it is a good target dataset[5].

In building these intrusion detection systems, we will look to use ML and DL
methodology. ML concerns itself with statistical methods that construct or
evaluate patterns from known functions and behavior [7]. Within all possible
ML methods, we need to look towards classification methods that specify if a
user is attacking and how, as well as one that can work regardless of our
class distribution, as we fundamentally need to explore attacks of all types,
and CSE-CIC-IDS2018 contains very skewed data. This leads us to preferring
K-nearest neighbors and Random Forest (RF) [7]. K-nearest neighbors is a non-
parametric learning algorithm classifying new points to the majority class of
the $K$ closest points. Because K-nearest neighbors is difficult to extract
information from, and the primary motivation of this research is to compare
methodologies based on a variety of factors including relative importance of
predictors, we will not use it. We will use RF, however, for a few reasons.
For data distributed with many extreme points, using a method like RF that
does not rely upon weights for given inputs but rather specific boundaries may
give it an edge in performance. Because RF has a relatively good balance of
low variance to low bias, as well as having much better results than a
standard decision tree, all while still being relatively easily interpretable,
we will use it [7].

Deep Neural Networks (DNN) are an offshoot of DL with the goal of modeling
complex nonlinear relationships using chains of neurons and activation
functions. With regards to our work, they are a strong tool to choose and
incorporate because of their fast training time given the scope of our data
and consistently strong performance when compared to nearly all ML methodology
[8].

In exploring our results, we will display the information in a manner that is
as transparent and digestible as possible. Similar to many papers, we include
the accuracy, precision, F1 score, and false alarm rate (FAR). Most other
metrics can be determined as a derivative of these, and they provide a
foundation by which to compare results [8, 9]. While [5] achieves results for
our dataset using DNN architecture, we achieve even better results and thus
choose to include it in conjunction with RF.

Feature selection, which aims to feed a model a subset of the original
features to improve, is the central focus of the paper. These improvements can
include time to build a model as well as performance results such as accuracy,
precision, F1 score, and FAR. There are three types of feature ranking methods
commonly used. Filters use predefined criteria or patterns in the data,
wrappers build many models to compare usefulness of given features, and
embedded methods rely upon training a model that then defines features for
another model [10]. For the sake of being able to compare more modern
optimizer filter methods, we will use both other filter methods and embedded
methods as benchmarks, since wrapper methods both have better accuracy and
much larger build times. A standard method for embedded feature selection is
RF information gain (IG). Because of its large amount of time to compute, it
is useful but limited in applicability [11]. The filter methods to be compared
are the correlation based feature selection Bat Algorithm (CFS-BA), and Aquila
Optimizer (CFS-AO).

Correlation based feature selection is built upon the idea that an
outperforming subset of features can be determined purely through correlation
between variables. The score derived is one that aims to maximize correlation
between predictors and final classes while minimizing correlation between
pairs of predictors [6]. This methodology performs much faster than any
embedded method, and can give more intuitive reasoning into the decision
making process, leading us to choose it. The only issue is that we cannot
explore every possible subset, and as such need to use an optimizer method to
find the closest subset to the global maxima score.

Bat Algorithm (BA) was proposed by Yang in 2010 as a metaheuristic
optimization algorithm to explore a space mimicking the behavior of microbats
using echolocation [12]. Most metaheuristic algorithms have two phases;
exploration attempts to discover the widest range of the feature space, and
exploitation aims to discover the best local solution within a given region.
Unlike most modern metaheuristic algorithms, BA gradually switches between
exploration and exploitation using a self-tuning hyperparameter [13]. CFS-BA
has performed well for CIC-IDS2017, but does not have CSE-CIC-IDS2018 results
and lacks results for DNN and RF models [6].

The Aquila Optimizer was introduced as a faster and potentially more efficient
metaheuristic algorithm than prior methodology. While slower to converge onto
a maxima, it has been shown to outperform BA when used in feature selection on
cybersecurity benchmark datasets [14]. Because of its promising nature, the
algorithm will be used and compared to BA with CFS scoring. AO has performed
well for CIC-IDS2017, but used a different fitness than CFS and does not have
CSE-CIC-IDS2018 results [14].

The paper is organized as follows: our motivation for the particular fields,
tools, and algorithms used are in section I, particular mathematical and
technical details regarding the implementation of our research are in section
II, information specific to our dataset and its distribution and processing
are in section III, experimental results are in section IV, interpretation of
results and limitations of the study are in section V, and a summary of the
work’s main achievements as well as future work are in section VI.

##  II Methodology

Our methodology will discuss the specifics of how we will build and tune RF
and DNN models for feature selection comparison.

###  II-A Classification Algorithms

Classification algorithms work by feeding a set of inputs through a model to
output a final, categorical variable. The methods we will explore for this are
RF and DNN.

####  II-A1 Random Forest

The basic building block for RF is the decision tree. Visualization is in fig.
1. The decision tree is a model that classifies data points by putting them
through a series of binary decision boundaries and then assigning them to the
same class as the majority of the training points within a final bucket. The
equation used to choose decision splits in this paper, the Gini Index [7], is
defined as eq. 1. Entropy can also be used, defined as eq. 2 [7].

$\displaystyle G$ $\displaystyle=\sum_{k=1}^{K}\hat{p}_{mk}(1-\hat{p}_{mk})$
(1) $\displaystyle D$
$\displaystyle=-\sum_{k=1}^{K}\hat{p}_{mk}\log\hat{p}_{mk}$ (2)

Both indexes are modeled to be a measure of total variance across all classes,
with a smaller value indicating a better decision tree. Here, $K$ represents
the number of classes, and $\hat{p}_{mk}$ represents the proportion of
training observations in the $m$th region from the $k$th class. The indexes
takes a smaller value for $\hat{p}_{mk}$ closer to zero or one, meaning it
acts as a measure of _purity_ for a given node.

Figure 1: Visualizations of an example RF courtesy [15].

Splits are made to minimize eq. 1. For a standard decision tree, every split
can compare all possible predictors; however, when used in conjunction with
bootstrapping, this increases variance and decreases accuracy [7]. To work
around this, RF forces each split to only compare one of $m=\sqrt{p}$ out of
the total $p$ predictors.

In Random Forest, all inputs are used to train a set of bootstrapped decision
trees. This means that some number of trees, set in this case to be 100, are
constructed from sets with the same amount of points as the original data set,
made by randomly selecting points from the base set with replacement. Each of
these decision trees will have decision splits added but never removed, and
reach a depth $d$. We tested the $d$ range of [2,5,10,20,40,100,200] building
RF models for each and calculating the OOB score, and found $d=20$ ideal for
our performance. Out-of-bag works by plugging the training set back into the
finalized model and recording the fraction of the training set that ends up
correctly labeled. This simulates predicted accuracy to a certain degree, and
thus is a substantial enough metric for our use.

Furthermore, we can use eq. 2 to determine the information gained at each
split by each variable. Working with the average value of all of these splits
over each predictor gives us our information gain (IG) scores, which we can
use to determine feature subsets.

####  II-A2 Deep Neural Networks

Neural Networks, and DNN more specifically, consist of an input layer
$X_{1},\ldots,X_{n}$, describing all predictors in numerical terms, multiple
hidden layers, consisting of fully connected nodes passed to the next layer
through nonlinear activation functions, and an output layer. In the case of
DNN for the sake of classifying, the process is as follows:

A hidden layer $n$ with $D$ nodes with prior layer $m$ having $C$ nodes can be
represented as

$A_{d}^{(n)}=h_{d}^{(n)}(X)=g(w_{d_{0}}^{(n)}+\sum\nolimits_{j=1}^{p}w_{d_{j}}^{(n)}A_{c}^{(m)})$

for some activation value $A_{d}^{(n)}$ for node $d$ defined through
activation function $h$ expanded for some nonlinear function $g$ using weights
$w$ for all $p$ predictors. For the first layer, $A_{c}^{(m)}$ can be
substituted with $X_{j}$. This is then repeated until layer $n^{\prime}$ with
$D^{\prime}$ nodes and piped to $K$ different linear models

$Z_{k}=\beta_{k_{0}}+\sum\nolimits_{d^{\prime}=1}^{D^{\prime}}\beta_{k_{d^{\prime}}}h_{d^{\prime}}^{(n^{\prime})}(X)=\beta_{k_{0}}+\sum\nolimits_{d^{\prime}=1}^{D^{\prime}}\beta_{k_{d^{\prime}}}A_{d^{\prime}}^{(n^{\prime})}.$

This value is then taken and put through the _softmax_ activation function eq.
3 for multiclass classification and the _sigmoid_ function eq. 4 for binary
classification [7]. Afterward it is classified to the $f_{k}(X)$ with the
highest probability.

$\displaystyle f_{k}(X)=Pr(Y=k|X)$
$\displaystyle=\frac{e^{Zk}}{\sum\nolimits_{l=0}^{K}e^{Zl}}$ (3)
$\displaystyle f_{k}(X)=Pr(Y=k|X)$
$\displaystyle=\frac{e^{Z}}{1+e}=\frac{1}{1+e^{-Z}}$ (4)

All other hidden layer activation functions are chosen to be ReLU (rectified
linear unit), defined as $A(z)=(z)_{+}$ that returns $0$ for negative inputs
and $z$ otherwise, as it achieves far better performance than most other
alternatives [7, 16].

Figure 2: Deep Neural Network Visualization for 2 inputs and 1 output

The number of neurons and hidden layers, while recommended to be large, may
perform best when minimized given the sparsity of smaller classes [16]. As
such, we compared two methodologies. The first assumes that we have enough
information to warrant a complex model. This model contains three hidden
layers, all of size 100. We will refer to it as the 100-100-100 model to
reinforce this. The second model will be built on the assumption that we have
a lack of information for the size of our data. This encourages a structure
that funnels information down to the final layer, and so we will use a hidden
layer of size 50 first and then of size 25. This will be the 50-25 model.
After these hidden layer passes for each respective model, the classification
will be done.

This training is done in batches, as the results and weights are updated for
each run. Common batch sizes are 32, 64, 128, and 256. The batch size
typically makes minimal difference and is limited by a resource’s RAM, and so
we choose to use batch size 32 for the sake of maximizing result
reproducability.

Overall, the 100-100-100 model underperformed the 50-25 model in build time by
large margins and accuracy by notable margins for all of the categorical
models. This led us to compare RF and DNN using the 50-25 structure.

###  II-B Correlation Feature Selection

CFS is a filter method that analyzes the fitness of a given subset of features
based solely on the intercorrelation between features. This is done through
Spearman’s correlation, which works for numerical values calculating
colinearity. To make the formula work for our classification, we used one-hot
encoding to get individualized scores.

The score for a given subset, $s$, can be calculated as eq. 5 [17].

$M_{s}=\frac{k\bar{r_{cf}}}{\sqrt{k+k(k-1)\bar{r_{ff}}}}$ (5)

$M_{s}$ is the score, valuing a higher number as a better subset, $k$ is the
number of predictors for a given subset, $\bar{r_{cf}}$ is the average
absolute intercorrelation between predictors and final features, and
$\bar{r_{ff}}$ is the average absolute intercorrelation between predictors. We
look at the average _absolute_ intercorrelation, so that the formula rewards
inverse and direct relations the same.

###  II-C Bat Algorithm

Our main goal is to use a metaheuristic optimizer to explore a feature space
for a given best point. In our case, the main metaheuristic optimizer is BA,
the feature space is all possible subsets of the 70 preprocessed predictors,
and the best point is the one with the maximum value from eq. 5.

We define some number of bats and instantiate each one according to eq. 6.
Uniform distribution between values $a$ and $b$ is denoted as $U(a,b)$.
$x_{i}^{0}$ represents the initial position for bat $i$ out of a total of
$n=100$ bats, with an $x$ value of $[0,0.5]$ at position $i$ meaning a
rejection of the predictor in location $i$ and a value of $[0.5,1]$
representing acceptance. $v_{i}^{0}$ represents the initial velocity,
$f_{i}^{0}$ represents the initial frequency, $A_{i}^{0}$ represents the
initial loudness, and $r_{i}^{0}$ represents some initial random value.

$\displaystyle x_{i}^{0}$ $\displaystyle=\begin{bmatrix}x_{1}\\\ x_{2}\\\
\cdots\\\ x_{k}\end{bmatrix}=\begin{bmatrix}U(0,1)\\\ U(0,1)\\\ \cdots\\\
U(0,1)\end{bmatrix}$ (6) $\displaystyle v_{i}^{0}$
$\displaystyle=\begin{bmatrix}v_{1}\\\ v_{2}\\\ \cdots\\\
v_{k}\end{bmatrix}=\begin{bmatrix}U(-1,1)\\\ U(-1,1)\\\ \cdots\\\
U(-1,1)\end{bmatrix}$ $\displaystyle f_{i}^{0}$
$\displaystyle=U(0,0.1),\>\>A_{i}^{0}=U(1,2),\>\>r_{i}^{0}=U(0,1)$

We then define the best point $x_{b}$ to be the position that correlates to
the highest value for eq. 5. Then for each epoch $t$, with $t_{max}=1000$, we
first update $v_{i}^{t}$ according to eq. 7, eq. 8, and eq. 9.

$\displaystyle f_{i}^{t}$ $\displaystyle=U(0,0.1)$ (7) $\displaystyle
v_{i}^{t}$ $\displaystyle=v_{i}^{t-1}+(x_{i}^{t-1}-x_{b})f_{i}^{t}$ (8)
$\displaystyle v_{i}^{t}$
$\displaystyle=\max\left(\min\left(v_{i}^{t},1\right),-1\right).$ (9)

Then we update $x_{i}^{t}$ according to eq. 10 or eq. 11. $A_{i}^{t}$
represents the loudness.

$\displaystyle x_{i}^{t}$
$\displaystyle=\max\left(\min\left(x_{i}^{t-1}+v_{i}^{t},1\right),-1\right)$
(10) $\displaystyle x_{i}^{t}$
$\displaystyle=\max\left(\min\left(x_{b}+\begin{bmatrix}U(-0.01,0.01)\\\
U(-0.01,0.01)\\\ \cdots\\\
U(-0.01,0.01)\end{bmatrix}A_{i}^{t},1\right),-1\right)$ (11)

After this, we update our frequency and loudness using eq. 12 and eq. 13 where
$\alpha$ and $\gamma$ are constants, both of which can be set to 0.95.

$\displaystyle A_{i}^{t+1}$ $\displaystyle=\alpha A_{i}^{t}$ (12)
$\displaystyle r_{i}^{t+1}$ $\displaystyle=r_{i}^{0}(1-e^{-\gamma^{t}})$ (13)

After determining all of these values, we can run the Bat Algorithm. The main
idea of this algorithm is to transition from exploration to exploitation
gradually through self-tuning behavior. We can now create algorithm 1, based
largely on [12, 6, 13].

Algorithm 1 Bat Algorithm

Input: Correlation matrix produced from training dataset
Output: Selected best feature subset $x_{b}$

1: Initialize $n$ bats using eq. 6

2: $x_{b}\leftarrow 0,M_{b}\leftarrow 0$

3: for $i\in[1,n]$ do

4: Calculate $M_{i}$ for subset $s_{i}$ determined by bat $x_{i}$ using eq. 5

5: if $M_{i}>M_{b}$ then

6: $x_{b}\leftarrow x_{i},M_{b}\leftarrow M_{i}$

7: end if

8: end for

9: $t\leftarrow 1$

10: while $t\leq t_{max}$ do

11: for $i\in[1,n]$ do

12: Update $f_{i}$ using eq. 7

13: Update $v_{i}$ using eq. 8 and then eq. 9

14: if $r_{i}\geq U(0,1)$ then

15: Update $x_{i}$ using eq. 10

16: else

17: Update $x_{i}$ using eq. 11

18: end if

19: Calculate $M_{i}$ using eq. 5

20: if $A_{i}>U(0,1)$ and $M_{i}>M_{b}$ then

21: $x_{b}\leftarrow x_{i},M_{b}\leftarrow M_{i}$

22: Update $r_{i}$ using eq. 13

23: Update $A_{i}$ using eq. 12

24: end if

25: end for

26: $t\leftarrow t+1$

27: end while

28: return $x_{b}$

###  II-D Assessment Metrics

To determine the fitness and results we need to clarify the metrics our models
will use.

####  II-D1 Binary Classification

When we get our results for both RF and DNN, we get them in the form of a
confusion matrix. This matrix describes the amount or ratio of what the data’s
class actually is compared to what it is described as. In the case of binary
classification describing if a data point is benign or malicious, there are
universal metrics. We will be describing an attack as positive and a benign
point as negative. As taken from [3], we will denote True Positive (TP), False
Negative (FN), False Positive (FP), and True Negative (TN). Standard used
metrics are defined below as

$\displaystyle Accuracy$ $\displaystyle=\frac{TP+TN}{TP+TN+FP+FN}$ (14)
$\displaystyle Precision$ $\displaystyle=\frac{TP}{TP+FP}$ (15) $\displaystyle
Recall$ $\displaystyle=\frac{TP}{TP+FN}$ (16) $\displaystyle FAR$
$\displaystyle=\frac{FP}{FP+TN}$ (17) $\displaystyle F1$
$\displaystyle=2\left({Recall}^{-1}+{Precision}^{-1}\right)^{-1}$ (18)

with FAR denoting the False Alarm Rate.

The choice of these metrics is mainly a result of suggested and common metrics
from [18] and [8].

####  II-D2 Multiclass

For multiclass classification, we want to use the same metrics, but there is
no longer a single target class. As such, all metrics need to treat every
class as a target and average in some way. Accuracy stays consistent as the
ratio of correctly labeled points to all points. For the other metrics, we
consider using either micro, macro, or weighted values.

Micro averages compute the values for the entire table, being blind to the
individual classes. Because of this, the micro equations for F1 and Precision
are equal and equate to Accuracy, defined as eq. 19.

$Accuracy=\frac{Correctly\>Labeled\>Points}{All\>Points}$ (19)

Macro averages compute the values on a class by class basis and then average
all the independent values. Weighted averages work similar to macro averages,
but instead weigh the score from each class by the proportion of that class to
the total. This means that weighted averages end up giving us slightly less
information as they mimic overall accuracy. Furthermore, our goal with using
CFS-BA is to focus on each class with equal emphasis since we use one-hot
encoding. Therefore, macro scores will tend to give us results more reflective
of how we performed compared to our expectations.

##  III Data Description

The dataset used was created by the Communications Security Establishment
(CSE) and the Canadian Institute for Cybersecurity (CIC) in 2018 to simulate
network data for the sake of aiding in the creation and testing of NIDS
construction [19, 20]. It contains simulated data for various attacks over 10
days, with each day having a different distribution. The full distribution on
a log scale is described in fig. 3. Although there are 6 types of attacks,
they are labeled as and can be expanded to 13 types of unique attack
signatures. There are 83 numerical inputs, both continuous and discrete
inputs, and each day has at least 79 inputs.

Figure 3: Log Scaled Frequency of Classes

###  III-A Data Preprocessing

To get the data into a more suitable format, it underwent preprocessing. The
first thing done to the data was to remove {Dst.IP, Flow.ID, Src.IP,
Src.Port}, features that are present for one of the days but not others. Our
model will be constructed for all of the days, and thus the large amount of NA
values that the inclusion of these would generate would cause issues.

Of the remaining 79 predictors, 1 had NA values and 2 had $\infty$ values.
There are only 59,721 rows with NA and 36,039 rows with $\infty$ out of 16.17
million. As such, it is better to remove the points instead of the predictors.
Additionally, {Timestamp} can be removed. It is unique to the simulated
environment and the inclusion of it would lead to a model built for past data,
and not projecting towards future data. While it may seem like {Dst.Port} is
similar to this, a model can be based on banning or including certain
connections, so it is included. This leaves us with 78 predictors, of which 8
more can be removed as they have uniform values for all days.

Data is then normalized using MinMax normalization in the range of [0,1], as
data is not normally distributed enough for Gaussian normalization, and is
typically skewed about an extreme or completely non-negative, not befitting
[-1,1] MinMax normalization. The full distribution for each variable can be
found through [21]. Most of the predictors have highly skewed and nonlinear
data, with a spike toward some extreme.

For our model analysis, we used a test-train split of 50/50 because our model
is large enough such that a smaller portion of it can be used to train and the
weaker classes are small enough such that we need both test and train to be a
significant portion of the data to get accurate results.

##  IV Results

Time results were achieved through Python ran in Google Colab Pro+. CFS-BA,
CFS-AO, and RF-IG were generated using both [22] and [23]. The RF results were
generated using scikit-learn [22]. The DNN results were generated using
TensorFlow [24].

The computer by which the hyperparameter results were achieved and compared in
terms of run time has a Nvidia GTX 1650 GPU, 32 GB of RAM, and has a Intel
Core i7 processor. Finalized building and running of models was done through
Google Colab Pro +, which utilized up to 53 GB of RAM for the run sessions.

###  IV-A Correlation

Correlation structure is calculated according to the methodology in section
II-B. After calculating our correlation, we stored it in a matrix and then
added a dividing line indicating where the features end and the classes begin.
The matrix is heat-encoded and shown in fig. 4. The black lines crossing the
figure represent the divide between predictors and classes, and the diagonal
line through the figure represents each variable’s intercorrelation with
itself, which is always 1. The main issue is that most of the regions that
have strong feature to class intercorrelation also have strong feature
intercorrelation. To maximize the result of $M_{s}$ from eq. 5, we want to
find features that have the lightest average coloration beneath the diagonal
line in the top left rectangle and the darkest coloration in the bottom left
rectangle. Specific predictor names are omitted for the sake of readability
and interpretability. The original $M_{s}$ score is 0.119 for multiclass and
0.248 for binary. The original information gain sum (IG) score is 1 for both
sets by definition.

Figure 4: Correlation Heatmap for Full Feature Set

After running CFS-BA, CFS-AO, and RF-IG for our tuned hyperparameters, we are
left with features that give us $M_{s}$ scores and IG scores shown in table I
and fig. 5. The four features common to all methods are {Fwd.Pkt.Len.Mean,
Flow.Pkts., Init.Fwd.Win.Byts, Fwd.Seg.Size.Min}. Note that for table I, table
III, and table II, cat. denotes categorical, bi. denotes binary, and time
denotes build time.

TABLE I: Reduced Feature Subset Scores

Methodology | K | CFS | IG
---|---|---|---
Cat. CFS-BA | 26 | 0.147 | 0.630
Cat. CFS-AO | 14 | 0.159 | 0.350
Cat. RF-IG | 14 | 0.136 | 0.678
Bi. CFS-BA | 8 | 0.323 | 0.279
Bi. CFS-AO | 5 | 0.410 | 0.232
Bi. RF-IG | 5 | 0.287 | 0.430

[labelOnlyA=5, labelOnlyB=12,
labelOnlyC=5,labelOnlyAB=5,labelOnlyAC=0,labelOnlyBC=5,labelABC=4,labelA=RF-
IG,labelB=CFS-BA,labelC=CFS-AO,showframe=false]

Figure 5: Venn Diagram of Categorical Feature Subset Overlap

###  IV-B Random Forest

Once we determined our finalized hyperparameters, we ran the models and
generated the performance results shown in table II. The performance results
are calculated as discussed in section II-D and the time results compare the
build time, describing the time taken to actually build the model. In terms of
categorical results, CFS-BA and RF-IG both outperform the full feature set,
with RF-IG being strongest in the most metrics and CFS-BA being significantly
faster. For binary classification, all models not using the full set are worse
than building a model for categorical classification and using it solely for
binary decision making.

TABLE II: Random Forest Performance Results

Methodology | Time (s) | Accuracy | Precision | FAR | F1
---|---|---|---|---|---
Cat. Full | 6954 | 98.41% |  90.45% | 22.78% | 79.64%
Cat. CFS-BA | 5727 |  98.44% | 89.5% | 21.26% | 81.02%
Cat. CFS-AO | 4146 | 98.40% | 88.60% | 21.62% | 80.21%
Cat. RF-IG | 10359 |  98.44% | 89.3% |  19.79% |  82.09%
Bi. Full | 5392 |  98.96% |  99.16% |  2.86% |  98.12%
Bi. CFS-BA | 2161 | 96.01% | 91.49% | 4.68% | 93.26%
Bi. CFS-AO | 1792 | 95.64% | 90.69% | 4.90% | 92.69%
Bi. RF-IG | 5493 | 97.19% | 98.14% | 8.08% | 94.69%

###  IV-C Deep Neural Network

We run our DNN based upon the methodology declared in section II-A2. Results
for our performance are shown in table III. With regards to performance, we
see that CFS-BA outperforms the full feature set in every binary metric. In
categorical classification, CFS-BA gets the highest accuracy and the lowest
FAR but is beat by the full set in F1 and the RF-IG model in precision.

TABLE III: DNN Performance Results

Methodology | Time (s) | Accuracy | Precision | FAR | F1
---|---|---|---|---|---
Cat. Full | 7689 | 98.01% | 79.38% | 29.66% | 70.67%
Cat. CFS-BA | 8413 | 98.29% | 84.52% | 29.44% | 71.63%
Cat. CFS-AO | 8125 | 98.11% | 81.03% | 30.08% | 71.21%
Cat. RF-IG | 13376 |  98.35% |  85.41% |  24.83% |  76.92%
Bi. Full | 9032 |  98.95% |  99.30% |  3.01% |  98.10%
Bi. CFS-BA | 6136 | 93.77% | 95.80% | 17.89% | 87.13%
Bi. CFS-AO | 6103 | 93.65% | 95.32% | 17.99% | 86.92%
Bi. RF-IG | 9382 | 97.06% | 98.04% | 8.43% | 94.43%

Looking at the confusion matrices of RF and DNN, we see some patterns that may
explain discrepancies. In both model types, DoS and BruteForce attacks were
commonly mislabeled as one another. Web attacks, being the smallest and most
difficult to catch, were actually more easily recognized than Infiltration
attacks in RF models. DNN models nearly entirely missed Web attacks, and
struggled just as much with Infiltration attacks. All other attacks performed
relatively similarly for both model types. RF models were more inclined to
mislabel to fewer types of labels, meaning confusion matrices with less 0’s
than their DNN alternatives.

##  V Discussion

Our models are able to achieve a maximum accuracy of 98.44% on the CSE-CIC-
IDS2018 dataset, using RF-IG with 14 of the 70 features for categorical
classification. Both RF and DNN values show CFS-BA and RF-IG to outperform the
full feature set. CFS-AO outperforms the full set for DNN, but not for RF.
Using F1-score as a holistic metric, RF-IG raises 2.45% and 6.25% for RF and
DNN respectively, while CFS-BA raises 1.38% and 0.96%. Using accuracy to
determine applicability, RF-IG raises 0.34% and 0.03%, and CFS-BA raises 0.28%
and 0.03%.

The information from this paper suggests subset size or IG score may be more
important than CFS. Contrasting CFS-BA to CFS-AO, the two benefits CFS-BA may
have are a larger subset size and IG score. Of these two, IG score seems to be
more likely to directly correlate to performance. The full feature set has the
largest number of features and yet underperforms. CFS-BA has more features
than RF-IG and yet generates comparable results for categorical classification
and poorer results for binary classification. CFS-BA has an IG score 92.9% of
RF-IG’s for categorical classification and performs much better than CFS-AO,
having an IG score 51.6% of RF-IG’s.

Despite all of this, CFS-BA maintains advantages to RF-IG in the trade off
between build time and performance. It builds in 62.9% of the time for DNN and
55.3% of the time for RF. Towards explainability, CFS-BA seems to be a strong
method for funneling down a set of useful features. While not as direct as RF-
IG, the models share many features, and so it can be used as a rough RF-IG
estimate.

In the realm of explainability, RF-IG provides individual scores of importance
while CFS only scores subsets. This means that RF-IG is a superior and more
accurate method when making recommendations for a particular vulnerability of
a site, while CFS may be more applicable towards reinforcing aspects of a
site, such as showing where the heaviest loads may be distributed if an attack
is slow to be caught. However, as shown in our time results, if someone wishes
to update individualized features in a rapid manner, CFS-BA can be a
relatively comparable alternative to RF-IG.

Fundamentally, CFS using Spearman’s correlation tests monotonicity, which may
be why it does not directly indicate performance. We know our data is highly
skewed, and it may have too many non-monotonic relationships. This could be
because factors have monotonic relationships with one another within a certain
classification, but the blending of all points generates noise that distorts
calculations. A rationale to this would be the performance of binary models
being weaker than categorical models. To overcome this, future models could
look at binary $M_{s}$ for each individual class, and take the union of all
feature subsets. Another piece of data reinforcing this theory is that the
100-100-100 model was worse than the 50-25 model. This suggests either that
the data is t0o sparse, or that the data appears sparse because the features
selected are not an accurate reflection of the most beneficial features.

##  VI Conclusion

In short, this paper added results for DNN and RF performance for CSE-CIC-
IDS2018 using the full feature set, CFS-BA, CFS-AO, and RF-IG for the same
model hyperparameters. The use of CFS-BA and RF-IG indisputably benefited the
model. It removed 63%+ of features and significant improved model performance.

There are many directions the work can go in to increase its applicability in
the field. CNN structure could be integrated as a replacement for DNN to
attempt to account for the highly skewed distribution in the data as well as
clear the feature obfuscation present. Shapley game theory scores could also
be used to further improve the explainability of the models and provide more
insight into why RF-IG excels where CFS may not. Lastly, scores other than CFS
could be applied for filter comparison, or methods of correlation such as
Pearson’s or individualized Spearman’s may be a better application of CFS.

## References

  * [1] Malware statistics & trends report | AV-TEST. [Online]. Available: https://www.av-test.org/en/statistics/malware/
  * [2] Cost of a data breach report 2021. [Online]. Available: https://www.ibm.com/security/data-breach
  * [3] Y. Xin _et al._ , “Machine learning and deep learning methods for cybersecurity,” vol. 6, pp. 35 365–35 381, conference Name: IEEE Access.
  * [4] I. H. Sarker, A. S. M. Kayes, S. Badsha, H. Alqahtani, P. Watters, and A. Ng, “Cybersecurity data science: an overview from machine learning perspective,” vol. 7, no. 1, p. 41. [Online]. Available: https://doi.org/10.1186/s40537-020-00318-5
  * [5] R. I. Farhan, A. T. Maolood, and N. F. Hassan, “Performance analysis of flow-based attacks detection on CSE-CIC-IDS2018 dataset using deep learning,” vol. 20, no. 3, pp. 1413–1418, number: 3. [Online]. Available: https://ijeecs.iaescore.com/index.php/IJEECS/article/view/22288
  * [6] Y. Zhou, G. Cheng, S. Jiang, and M. Dai, “Building an efficient intrusion detection system based on feature selection and ensemble classifier,” vol. 174, p. 107247. [Online]. Available: http://arxiv.org/abs/1904.01352
  * [7] G. James, D. Witten, T. Hastie, and R. Tibshirani, _An Introduction to Statistical Learning with Applications in R_ , 2nd ed.
  * [8] J. L. Leevy and T. M. Khoshgoftaar, “A survey and analysis of intrusion detection models based on CSE-CIC-IDS2018 big data,” vol. 7, no. 1, p. 104\. [Online]. Available: https://doi.org/10.1186/s40537-020-00382-x
  * [9] S. Mahdavifar and A. A. Ghorbani, “Application of deep learning to cybersecurity: A survey,” vol. 347, pp. 149–176. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0925231219302954
  * [10] Akashdeep, I. Manzoor, and N. Kumar, “A feature reduced intrusion detection system using ANN classifier,” vol. 88, pp. 249–257. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0957417417304748
  * [11] X. Li, W. Chen, Q. Zhang, and L. Wu, “Building auto-encoder intrusion detection system based on random forest feature selection,” vol. 95, p. 101851\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167404820301231
  * [12] X.-S. Yang and A. Gandomi, “Bat algorithm: A novel approach for global engineering optimization,” vol. 29.
  * [13] X.-S. Yang and X. He, “Bat algorithm: literature review and applications,” vol. 5, no. 3, pp. 141–149. [Online]. Available: https://doi.org/10.1504/IJBIC.2013.055093
  * [14] A. Fatani, A. Dahou, M. A. A. Al-qaness, S. Lu, and M. Abd Elaziz, “Advanced feature extraction and selection approach using deep learning and aquila optimizer for IoT intrusion detection system,” vol. 22, no. 1, p. 140, number: 1 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/1424-8220/22/1/140
  * [15] L. Blakely, M. Reno, and R. Broderick, “Evaluation and comparison of machine learning techniques for rapid qsts simulations,” 07 2018.
  * [16] Deep-learning: investigating deep neural networks hyper-parameters and comparison of performance to shallow methods for modeling bioactivity data | journal of cheminformatics | full text. [Online]. Available: https://jcheminf.biomedcentral.com/articles/10.1186/s13321-017-0226-y
  * [17] M. A. Hall _et al._ , “Correlation-based feature selection for machine learning,” 1999.
  * [18] Z. Ahmad, A. Shahid Khan, C. Wai Shiang, J. Abdullah, and F. Ahmad, “Network intrusion detection system: A systematic study of machine learning and deep learning approaches,” vol. 32, no. 1, p. e4150, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/ett.4150. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/ett.4150
  * [19] IDS 2018 | datasets | research | canadian institute for cybersecurity | UNB. [Online]. Available: https://www.unb.ca/cic/datasets/ids-2018.html
  * [20] A realistic cyber defense dataset (CSE-CIC-IDS2018) - registry of open data on AWS. [Online]. Available: https://registry.opendata.aws/cse-cic-ids2018/
  * [21] IDS 2018 intrusion CSVs (CSE-CIC-IDS2018). [Online]. Available: https://www.kaggle.com/datasets/solarmainframe/ids-intrusion-csv
  * [22] F. Pedregosa _et al._ , “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * [23] N. V. Thieu, “A collection of the state-of-the-art meta-heuristics algorithms in python: Mealpy,” 2020. [Online]. Available: https://doi.org/10.5281/zenodo.3711948
  * [24] M. Abadi _et al._ , “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/


# Adaptive Safe Reinforcement Learning-Enabled Optimization of Battery Fast-
Charging Protocols

Myisha A. Chowdhury†, Saif S.S. Al-Wahaibi†, and Qiugang Lu†

†Department of Chemical Engineering, Texas Tech University, Lubbock, TX 79409,
USA
Corresponding author: Qiugang Lu; Email<EMAIL_ADDRESS>
###### Abstract

Optimizing charging protocols is critical for reducing battery charging time
and decelerating battery degradation in applications such as electric
vehicles. Recently, reinforcement learning (RL) methods have been adopted for
such purposes. However, RL-based methods may not ensure system (safety)
constraints, which can cause irreversible damages to batteries and reduce
their lifetime. To this end, this work proposes an adaptive and safe RL
framework to optimize fast charging strategies while respecting safety
constraints with a high probability. In our method, any unsafe action that the
RL agent decides will be projected into a safety region by solving a
constrained optimization problem. The safety region is constructed using
adaptive Gaussian process (GP) models, consisting of static and dynamic GPs,
that learn from online experience to adaptively account for any changes in
battery dynamics. Simulation results show that our method can charge the
batteries rapidly with constraint satisfaction under varying operating
conditions.

Keywords: Battery fast-charging; Safe reinforcement learning; Action
projection; Gaussian process; Dynamic environment

##  1 Introduction

The design of fast-charging strategies for Lithium-ion batteries is of
significance to reduce the charging time, alleviate the mileage anxiety for
electric vehicles (EV) [1], and optimize EV fleet economics [2]. Fast-charging
of batteries with overly large currents may accelerate the growth of lithium
dendrite and electrolyte decomposition, leading to degraded battery lifetime
[3]. Moreover, the excessive heat generated from aggressive charging may bring
significant safety hazards, including fires and explosions [4]. Thus,
reasonable fast-charging protocols shall be deliberately designed to reduce
the charging time while accounting for the safety and degradation of batteries
[5].

Designing fast-charging strategies for batteries has been widely studied in
the literature [6, 7]. Existing approaches can be broadly categorized into
three groups: ad-hoc, model-based, and data-driven methods [8, 9]. For ad-hoc
methods, examples include the classical constant-current constant-voltage
(CCCV) charging, pulse charging, and other variants [10, 11]. These methods
are simple in implementation but overly conservative for optimality [12]. For
model-based methods, the design of fast-charging protocols is often formulated
as a constrained optimization with battery electrochemical models included in
the constraints [13]. However, such methods are restricted by the high
complexity of solving the battery models [14]. For data-driven methods,
battery models are no longer needed, and the optimal charging profile is
acquired by machine learning (ML) [15, 16] and deep reinforcement learning
(DRL) [17]. Unlike other ML methods, e.g., those with Bayesian optimization
[15], DRL-based methods can adaptively learn the optimal charging protocol as
battery ages and parameter drifts, thus receiving wide attention.

Reinforcement learning (RL) is a sequential decision-making method that
iteratively interacts with the environment by observing the state, generating
an action, and deploying the action to the environment [18]. The decision-
making of the RL agent improves by learning from the environment reward to
eventually optimize a prescribed objective. DRL-based fast-charging
optimization is firstly studied in [19, 20], where the deep deterministic
policy gradient (DDPG) method is employed for designing charging protocols
that adapt to environmental changes. In [21], Q-learning is combined with
recurrent neural network (RNN) to design charging profiles to minimize the
cost of the participation of battery storage in the power grid. The authors in
[5] study the integration of knowledge-based multiphysics constraints with DRL
to balance between the charging speed and battery thermal
constraints/degradation. Other advances include DRL-based optimal charging for
extending battery lifetime [22], battery pack fast-charging with balance
awareness [23], model-based DRL for fast-charging design [17], and RL-based
multi-stage constant current charging [24].

Despite the rapid development of DRL-based optimization of fast-charging
strategies, the majority of reported methods cannot strictly ensure constraint
satisfaction throughout the charging process due to the usage of soft penalty
for constraint violations in the reward [5]. This can lead to degraded battery
life or safety risks. Given this knowledge gap, the proposed work will develop
safe DRL-based battery fast-charging strategies to rapidly charge the battery
while meeting the (safety) constraints with a high probability. Specifically,
this work proposes to use probabilistic Gaussian process (GP) models online to
approximate the constraint functions. With the trained and updated GP models,
for any decided RL action (charging current) from the agent, the constraint
values are predicted along with an uncertainty bound. Such information will be
employed to construct a safety region, which is a set of safe actions that can
satisfy all constraints. The agent can then project its actions into the
safety region to ensure safety prior to deploying the action. The projection
step is performed by solving a constrained optimization problem.

In this work, the classical twin-delayed DDPG (TD3) algorithm will be used as
backbone RL algorithm due to its deterministic nature for fast convergence and
ability to mitigate the overestimation of Q values incurred with the DDPG
algorithm [25]. The main contributions of this work are as follows:

(1) GP-based models are presented to adaptively learn the constraint function.
Safety regions of the action can then be constructed, and any unsafe action
will be projected into this region before deployment to meet constraints.

(2) Our work also proposes an adaptive safe RL-based optimization of battery
fast-charging protocols. The optimized protocol can minimize the charging time
while respecting the constraints with a high probability even under varying
ambient and internal conditions. Note that how to ensure hard constraints of
the system is a common issue for RL-based methods (not only for battery fast
charging). Thus, our method is of significance for applications beyond
batteries where safety constraints are critical and also the system dynamics
vary over operating conditions.

The remainder of this paper is organized as follows. Section 2 presents the
fundamentals of RL, including a brief introduction to the TD3 algorithm.
Section 3 demonstrates the proposed adaptive safe RL method, and details are
given about the GP surrogate modeling and action projection. The specific
adaptive safe RL-based battery fast-charging optimization is given in Section
4. The effectiveness of the proposed method is validated on the PyBaMM
simulator in Section 5. The conclusions are given in Section 6.

##  2 Preliminaries

This section will give a brief overview of RL and the TD3 algorithm that will
be implemented in this work. A detailed description of the RL algorithms can
be found in [18].

###  2.1 Reinforcement learning

A general RL framework consists of an agent and an environment. The agent is a
continuously learning decision-maker, while the environment refers to all the
system elements that the agent interacts with. The decision-making process of
RL can be formalized as a Markov decision process (MDP),
$\mathcal{M}:(\mathcal{S},\mathcal{A},\mathcal{P},\gamma,r)$, where
$\mathcal{S}$ and $\mathcal{A}$ are the state and action space, respectively,
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\left[0,1\right]$
is the state transition probability, $\gamma\in[0,1]$ is the discount factor,
and $r:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the reward
function. In a MDP setting, at each time step $t$, the agent observes the
state vector $s_{t}\in\mathcal{S}$, generates an action $a_{t}\in\mathcal{A}$
and deploy it to the environment, which results in one-step environment
evolution to $s_{t+1}\in\mathcal{S}$ based on the transition probability
$\mathcal{P}(s_{t+1}|s_{t},a_{t})$. At any time step $t$, the policy
$\pi(a_{t}|s_{t}):s_{t}\rightarrow\mathcal{P}(a_{t})$ gives the probability of
selecting action $a_{t}$ given the environment state $s_{t}$. Additionally,
each state is evaluated by the state-value function,
$V^{\pi}:\mathcal{S}\rightarrow\mathbb{R}$, which assigns a value to each
state by calculating the expected total cumulative reward, starting from that
state and following policy $\pi$ henceforth

$V^{\pi}(s_{t})=\mathbb{E}_{\pi}\left[\sum\nolimits_{k=0}^{\infty}\gamma^{k}r(s_{t+k},a_{t+k})\right].$
(1)

where $\gamma$ is the discount factor. Moreover, the state-action value
function $Q^{\pi}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$, also
known as Q-function, evaluates each state-action pair by calculating the
expected total cumulative reward starting from a state and the selected action
by policy $\pi$ based on that state:

$Q^{\pi}(s_{t},a_{t})=\mathbb{E}_{\pi}\left[\sum\nolimits_{k=0}^{\infty}\gamma^{k}r(s_{t+k},a_{t+k})|s_{t},a_{t}\right].$
(2)

One can re-write (2) recursively with the Bellman equation

$\displaystyle Q^{\pi}(s_{t},a_{t})=$ $\displaystyle\mathbb{E}_{s_{t+1}\sim
p(\cdot|s_{t},a_{t})}[r(s_{t},a_{t})+\gamma\mathbb{E}_{a_{t+1}\sim\pi(\cdot|s_{t+1})}\left[Q^{\pi}(s_{t+1},a_{t+1})\right]],$
(3)

where $p(\cdot|s_{t},a_{t})$ is the transition probability between successive
states of the environment, provided that the current state and action are
$s_{t}$ and $a_{t}$, respectively. The objective of the MDP is to find the
optimal policy $\pi^{\ast}$ that can yield the maximal Q-value
$Q^{\ast}:=\max_{\pi}~{}Q^{\pi}(s_{t},a_{t})$, that is,

$\pi^{*}=\arg\max_{\pi}~{}Q^{\pi}(s_{t},a_{t}).$ (4)

###  2.2 TD3 algorithm

Similar to the DDPG algorithm [26], the TD3 method is also a deterministic
actor-critic method for continuous state-action space [25]. Traditional DDPG
method suffers from the overestimation bias of the Q value due to the
repetitive maximization operation over the target Q networks (critics) in
calculating the temporal difference (TD) target. Such bias leads to inaccurate
estimations of the Q value. To address this issue, the TD3 adopts two target
critic networks, and the smaller one is used to compute the TD target [25]

$y_{t}=r(s_{t},a_{t})+\gamma\min_{i=1,2}Q_{\theta^{\prime}_{i}}(s_{t+1},\mu_{\phi^{\prime}}(s_{t+1})),$
(5)

where $\theta^{\prime}$ and $\phi^{\prime}$ are the parameters of the target
critic network $Q_{\theta^{\prime}_{i}}$ and target actor network
$\mu_{\phi^{\prime}}$, and $r(s_{t},a_{t})$ is the reward. The TD3 also has
two critic networks $Q_{\theta_{i}}(s_{t},a_{t})$, whose parameters
$\theta_{i},i=1,2$, are updated via minimizing the TD error by gradient
descent

$\nabla_{\theta_{i}}\frac{1}{|B|}\sum_{(s_{t},a_{t},r_{t},s_{t+1})\in
B}\left[y_{t}-Q_{\theta_{i}}(s_{t},a_{t})\right]^{2},$ (6)

where the batch $B$ contains transition tuples $(s_{t},a_{t},r_{t},s_{t+1})$
sampled from the reply buffer. The actor-network $\mu_{\phi}(s_{t})$
approximates the policy mapping from state $s_{t}$ to action $a_{t}$, with
parameters $\phi$ whose update is via gradient ascent

$\nabla_{\phi}\frac{1}{|B|}\sum_{s_{t}\in
B}Q_{\theta_{1}}(s_{t},\mu_{\phi}(s_{t})).$ (7)

The update of target network parameters follows that of the Polyak averaging
($0<\rho<1$, e.g., $\rho=0.99$):

$\theta_{i}^{\prime}\leftarrow\rho\theta_{i}^{\prime}+(1-\rho)\theta_{i},~{}~{}\phi^{\prime}\leftarrow\rho\phi^{\prime}+(1-\rho)\phi.$
(8)

##  3 Action projection-based safe RL

###  3.1 Gaussian process (GP) model

This section first introduces the basics of GP models, and they will be
employed in our safe RL for building surrogate models of the constraint
functions. As a non-parametric method [27], GP assumes that the function to be
learned $q(x)$ is sampled from a Gaussian distribution with the mean function
$\mathbb{E}\left[q(x)\right]=m(x)$ (often assumed to be zero) and covariance
function $cov[q(x),q(x^{\prime})]=k(x,x^{\prime})$, where $k(\cdot,\cdot)$ is
the kernel function. The covariance function defines the output correlation of
the function $q(\cdot)$ at inputs $x$ and $x^{\prime}$. One commonly used
kernel function due to its smooth and differentiable form is the radial-basis
function [27],

$k(x,x^{\prime})=\sigma_{f}^{2}\exp\left[-\|x-x^{\prime}\|^{2}/{2l^{2}}\right],$
(9)

where $\sigma_{f}^{2}$ and $l^{2}$ control the vertical span of the curve and
the correlation drop speed as distance increases, respectively. Consider $n$
sampled input-output pairs $\\{x_{1:n},y_{1:n}\\}$, where $y_{i}=q(x_{i})$.
With GP, for a query point $x^{\ast}$, the function value
$y^{\ast}=q(x^{\ast})$ is inferred by conditional multivariate Gaussian
distribution $\mathcal{N}(\cdot,\cdot)$ [28]:

$\displaystyle
p(y^{\ast}|x^{\ast},x_{1:n},y_{1:n})=\mathcal{N}(y^{\ast}|\mu_{n},\sigma_{n}),$
(10)
$\displaystyle\mu_{n}:=\mu(x^{\ast},x_{1:n},y_{1:n})=K(x^{\ast},x_{1:n})K(x_{1:n},x_{1:n})^{-1}y_{1:n},$
$\displaystyle\sigma_{n}:=\sigma(x^{\ast},x_{1:n})=K(x^{\ast},x^{\ast})-K(x^{\ast},x_{1:n})K(x_{1:n},x_{1:n})^{-1}K(x_{1:n},x^{\ast}),$

where $\mu_{n}$ and $\sigma_{n}$ are the posterior mean and variance, and $K$
is the covariance matrix constructed from the kernel function $k(\cdot,\cdot)$
with entry $K_{i,j}=k(x_{i},x_{j})$.

###  3.2 Action projection-based safe RL

Although RL has been increasingly employed for sequential decision-making
problems with black-box environments, its lack of ability to ensure
constraints hinders the usage of RL in real-world safety-critical
applications. To this end, this work proposes an action projection method (see
Fig. 1) where an additional safety layer is appended after the actor network.
The action projection below will be solved by the safety layer:

$\displaystyle\min_{\tilde{a}_{t}\in\mathcal{A}}~{}~{}||a_{t}-\tilde{a}_{t}||^{2}_{2},~{}~{}s.t.,\quad
z_{t+1}=q(z_{t},h_{t},\tilde{a}_{t})\in\mathcal{Z},$ (11)

where $a_{t}$ is the raw agent action, $\mathcal{A}$ is the action space, and
$\tilde{a}_{t}$ is the (projected) closest point to $a_{t}$ in the safety
region (green area in Fig. 1 (a)), defined by next-step safety variable
$z_{t+1}$ and feasible set $\mathcal{Z}$. In general, constraint dynamics
$z_{t+1}=q(z_{t},h_{t},\tilde{a}_{t})$ is black-box, e.g., it can the
(unknown) relation between next-step battery temperature and the charging
current to be applied ($h_{t}$ are some known variables at $t$). Thus, it is
proposed to use GP models to build data-driven surrogates of $q(\cdot)$ from
historical data. Then the prediction (posterior mean and variance) of
$z_{t+1}$ from the GP model (10) can be utilized to replace the constraint in
(11). Fig. 1 (b) gives an example where the posterior mean (blue curve) and
posterior variance (blue shaded area) of $z_{t+1}$, inferred by the trained GP
model, are used to construct the safety region (red bars) where $\mathcal{Z}$
is a threshold (red dashed) on the upper uncertainty bound of the posterior
prediction.

Figure 1: Schematics of the proposed action projection-based safe RL.

##  4 Safe RL for fast charging optimization

This section details the proposed adaptive safe RL method for optimizing
battery fast-charging protocols while respecting system constraints with a
high probability. For demonstration, the TD3 algorithm is used as the backbone
RL method, which can be easily extended to other RL algorithms, such as DDPG
[26], SAC [29], etc. Note that the proposed safe and adaptive RL method falls
into the online off-policy category (similar to TD3, DDPG, etc.) since the
agent needs to interact online with the environment to gather new data and
also access history data (collected under a different policy) for the training
[30].

###  4.1 Fast-charging optimization formulation

This work aims to develop an RL-based battery fast-charging protocol while
respecting system constraints. The fast-charging optimization problem can be
formulated as [17]

$\displaystyle\max_{\underline{I}\leq I_{t}\leq\bar{I}}\quad-t_{f},$
$\displaystyle s.t.\quad T_{t}\leq\bar{T},\quad V_{t}\leq\bar{V},$
$\displaystyle\quad\quad SOC(t_{0})=SOC_{0},\quad SOC(t_{f})=SOC_{ref},$ (12)

where $t_{f}$ and $I_{t}$ is the charging duration and current, respectively.
Additionally, $T_{t}$ and $V_{t}$ are the temperature and voltage of the
battery at time $t$, respectively, $SOC_{0}$ is the initial state-of-charge
(SOC) at time $t_{0}$, and $SOC_{ref}$ is the SOC at time $t_{f}$ when
considering the charging cycle is complete. Finally, $\bar{T}$ and $\bar{V}$
are the maximum allowable operating temperature and voltage, respectively.

###  4.2 Static safe RL for fast charging optimization

This section presents our safe RL framework for battery fast-charging
optimization. The RL environment is the battery, and the action $a_{t}$ is the
charging current $I_{t}$ to be decided. The safety projection for problem (12)
can be formulated as

$\displaystyle\min_{\tilde{a}_{t}\in\mathcal{A}}~{}~{}||a_{t}-\tilde{a}_{t}||^{2}_{2},$
(13) $\displaystyle s.t.\quad$ $\displaystyle
T_{t+1}=f(T_{t},\tilde{a}_{t-1},\tilde{a}_{t})\leq\bar{T},$ (14)
$\displaystyle V_{t+1}=g(V_{t},\tilde{a}_{t-1},\tilde{a}_{t})\leq\bar{V},$
(15)

where black-box functions $f(\cdot)$ and $g(\cdot)$ predict next-step
temperature and voltage respective for the given current $\tilde{a}_{t}$.
Static GP models $\hat{f}_{s}(\cdot)$ and $\hat{g}_{s}(\cdot)$ are trained
with historical data to approximate above constraint functions. The data of
the first $M$ RL training episodes ($n$ samples) are selected as the
historical data: $\\{x^{f}_{1:n-1},T_{2:n}\\}$ where
$x^{f}_{i}=[T_{i},a_{i-1},a_{i}]^{\top}$, and $\\{x^{g}_{1:n-1},V_{2:n}\\}$,
where $x^{g}_{i}=[V_{i},a_{i-1},a_{i}]^{\top}$. Note that the training data
for static GP is collected prior to the deployment of action projection. The
L-BFGS algorithm is used to train these two GP models. Starting from episode
$M+1$, for any query RL action $\tilde{a}_{t}$ at time $t$, the predicted
posteriors by trained static GPs $\hat{f}_{s}(\cdot)$, $\hat{g}_{s}(\cdot)$,
are $p(T_{t+1}|\tilde{a}_{t})\sim\mathcal{N}(\mu_{t}^{f},\sigma_{t}^{f})$ for
temperature and
$p(V_{t+1}|\tilde{a}_{t})\sim\mathcal{N}(\mu_{t}^{g},\sigma_{t}^{g})$ for
voltage. As in (10), here $\mu_{t}^{f}=\mu(x_{t}^{f},x_{1:n-1}^{f},T_{2:n})$
and $\sigma_{t}^{f}=\sigma(x_{t}^{f},x_{1:n-1}^{f})$. Also,
$\mu_{t}^{g}=\mu(x_{t}^{g},x_{1:n-1}^{g},V_{2:n})$ and
$\sigma_{t}^{g}=\sigma(x_{t}^{g},x_{1:n-1}^{g})$, where
$x_{t}^{f}=[T_{t},\tilde{a}_{t-1},\tilde{a}_{t}]$ contains present temperature
$T_{t}$ and the previous projected action $\tilde{a}_{t-1}$, and also the
query $\tilde{a}_{t}$ to be solved. Similarly,
$x_{t}^{g}=[V_{t},\tilde{a}_{t-1},\tilde{a}_{t}]$.

With the predicted posteriors, one can approximate (14)-(15) by using the
upper uncertainty bound (UUB) for our surrogate models. Then the problem
(13)-(15) is cast as:

$\displaystyle\min_{\tilde{a}_{t}\in\mathcal{A}}~{}~{}~{}||a_{t}-\tilde{a}_{t}||^{2}_{2},$
(16) $\displaystyle s.t.\quad$
$\displaystyle\mu_{t}^{f}+\kappa\sigma_{t}^{f}\leq\bar{T},~{}\mu_{t}^{g}+\kappa\sigma_{t}^{g}\leq\bar{V},$
(17)

where $\kappa$ is a hyper-parameter. Fig. 1 (b) shows the $\pm 3$ standard
deviation uncertainty bound (blue shaded) constructed from (17) where the
upper bound is used to construct the feasible region (red bars). For this
method, once GP models are trained, they will be fixed for future RL episodes,
thereby known as static GP-based safe RL, with pseudo-code in Algorithm 1.

Algorithm 1: Static GP-based safe RL (based on TD3 algorithm)
---
1: |  Input: initialize Q-function parameters $\phi_{1}$, $\phi_{2}$, and policy parameter $\theta$, empty
| replay buffer $\mathcal{D}$ and training data storage $\mathcal{X}$;
2: | Initialize target networks $\phi_{1}^{\prime}\leftarrow\phi_{1}$, $\phi_{2}^{\prime}\leftarrow\phi_{2}$, $\theta^{\prime}\leftarrow\theta$;
3: | for episode $m=1$ to $M$, repeat
4: |  $t\leftarrow 1$;
5: |  Observe state $s_{t}$ and select $a_{t}$ randomly in $\mathcal{A}$;
6: |  Execute $a_{t}$ in the environment;
7: |  Observe the next state $s_{t+1}$ and compute the reward $r_{t}$;
8: |  Store $\left(z_{t},a_{t-1},a_{t}\right)$ and $(z_{t+1})$ in $\mathcal{X}$ to train GP models, where $z_{t}$ is the
|  variable of interest (e.g., $T_{t}$ and $V_{t}$ in (14) and (15)) in the
constraints;
9: |  Store the transition $\left(s_{t},a_{t},r_{t},s_{t+1}\right)$ into the replay buffer $\mathcal{D}$;
10: |  Update RL agent parameters:
11: |  Randomly select a batch $B=\left\\{\left(s_{t},a_{t},r_{t},s_{t+1}\right)\right\\}$ from $\mathcal{D}$ with
|  the number of transitions as $|B|$;
12: |  Update the critic network parameters:
13: |  Compute the target for the Q-functions
|
$y_{t}=r_{t}+\gamma(\min_{i=1,2}Q_{\phi_{i}^{\prime}}(s_{t+1},a_{t+1}(s_{t+1}))$;
14: |  Update the Q-function using cost function w.r.t. $\phi_{i}$:
|
$\nabla_{\phi_{i}}\frac{1}{|B|}\sum\limits_{\left(s_{t},{a}_{t},r_{t},s_{t+1}\right)\in
B}\left(Q_{\phi_{i}}(s_{t},{a}_{t})-y_{t}\right)^{2}$, $i=1,2$;
15: |  Update the actor and target networks:
16: |  if $t\quad\%\quad d=0$, where $d$ is the delayed update period
17: |  $\nabla_{\theta}\frac{1}{|B|}\sum\limits_{s_{t}\in B}Q_{\phi_{1}}(s_{t},\mu_{\theta}(s_{t}))$,
18: |  $\phi_{i}^{\prime}\leftarrow\rho\phi_{i}^{\prime}+(1-\rho)\phi_{i}$ , $i=1,2$,
19: |  $\theta^{\prime}\leftarrow\rho\theta^{\prime}+(1-\rho)\theta$;
20: |  $t\leftarrow t+1$, until $s_{t+1}$ is terminal;
21: | end for
22: | Train the GP models $\hat{f}_{s}(\cdot)$, $\hat{g}_{s}(\cdot)$ using $\mathcal{X}$ with L-BFGS;
23: | for episode $m=M+1$ to END, repeat
24: |  $t\leftarrow 1$
25: |  Observe state $s_{t}$ and use actor network $\mu_{\theta}(\cdot)$ for action $a_{t}$;
26: |  Predict the posterior mean and variance of next-step variable
|  with the trained static GPs $\hat{f}_{s}(\cdot)$, $\hat{g}_{s}(\cdot)$;
27: |  Solve the optimization problem in (16)-(17) to find $\tilde{a}_{t}$;
28: |  Execute $\tilde{a}_{t}$, observe $s_{t+1}$, and compute $r_{t}$;
29: |  Store the transition $\left(s_{t},\tilde{a}_{t},r_{t},s_{t+1}\right)$ into $\mathcal{D}$;
30: |  Update the actor and critic following steps 10-19;
31: |  $t\leftarrow t+1$, until $s_{t+1}$ is terminal or $m$ reaches END;
32: | end for

###  4.3 Adaptive safe RL for varying environments

The previously proposed RL algorithm with a GP-based safety layer assumes a
static environment. However, the GP models trained with offline data cannot
account for the changes in environments and system dynamics, and hence may not
adaptively predict future values in such situations. Specific to batteries,
changes in the operating conditions such as ambient temperature or battery
aging alter the internal dynamics of the battery [31]. Thus, it is essential
to develop adaptive GP models that capture battery dynamics in real-time to
ensure the effectiveness of action projection in the case of changes in
ambient conditions or battery aging.

Our method will employ the previous static GP models trained with the data
from the first $M$ RL episodes as baseline models to capture the overall
trends of variables to be predicted. Then the differences between true
variable values and posterior means $\mu_{t}^{f}$ and $\mu_{t}^{g}$ from the
static GPs $\hat{f}_{s}(\cdot)$ and $\hat{g}_{s}(\cdot)$, known as residuals,
are

$\displaystyle\Delta
T_{t+1}:=T_{t+1}-\mu_{t}^{f}=f_{\Delta}(T_{t},\tilde{a}_{t-1},\tilde{a}_{t}),$
(18) $\displaystyle\Delta
V_{t+1}:=V_{t+1}-\mu_{t}^{g}=g_{\Delta}(V_{t},\tilde{a}_{t-1},\tilde{a}_{t}),$
(19)

with black-box functions $f_{\Delta}(\cdot)$ and $g_{\Delta}(\cdot)$. It is
proposed to use two GP models, known as dynamic GPs, $\hat{f}_{\Delta}(\cdot)$
and $\hat{g}_{\Delta}(\cdot)$, to capture such residuals. Such dynamic GPs
will be updated at every timestep throughout an RL episode and re-trained from
scratch at the beginning of a new episode. Consider the present episode $k>M$,
and the current timestep $t$. The collected data up to $t$ for the present
episode are $\\{x_{1:t-1}^{f},\Delta T_{2:t}\\}$ and $\\{x_{1:t-1}^{g},\Delta
V_{2:t}\\}$, with $x_{i}^{f}=[T_{i},\tilde{a}_{i-1},\tilde{a}_{i}]^{\top}$ and
$x_{i}^{g}=[V_{i},\tilde{a}_{i-1},\tilde{a}_{i}]^{\top}$. The two dynamic GPs
are then trained based on these data. For time $t$, the posterior predictions
are

$p(\Delta
T_{t+1}|\tilde{a}_{t})\sim\mathcal{N}(\bar{\mu}_{t}^{f},\bar{\sigma}_{t}^{f}),~{}p(\Delta
V_{t+1}|\tilde{a}_{t})\sim\mathcal{N}(\bar{\mu}_{t}^{g},\bar{\sigma}_{t}^{g}),$
(20)

where $\bar{\mu}_{t}^{f}=\mu(x_{t}^{f},x_{1:t-1}^{f},\Delta T_{2:t})$ and
$\bar{\sigma}_{t}^{f}=\sigma(x_{t}^{f},x_{1:t-1}^{f})$. $\bar{\mu}_{t}^{g}$
and $\bar{\sigma}_{t}^{g}$ are defined analogously and are omitted here. With
the posterior distributions for $\Delta T_{t+1}$ and $\Delta V_{t+1}$, the
overall posterior predictions from the adaptive GPs, consisting static and
dynamic GPs, are shown to be

$p(T_{t+1}|\tilde{a}_{t})\sim\mathcal{N}(\hat{\mu}_{t}^{f},\hat{\sigma}_{t}^{f}),~{}p(V_{t+1}|\tilde{a}_{t})\sim\mathcal{N}(\hat{\mu}_{t}^{g},\hat{\sigma}_{t}^{g}),$
(21)

where $\hat{\mu}_{t}^{f}=\mu_{t}^{f}+\bar{\mu}_{t}^{f}$ and
$\hat{\mu}_{t}^{g}=\mu_{t}^{g}+\bar{\mu}_{t}^{g}$. The definitions of
posterior variance $\hat{\sigma}_{t}^{f}$ and $\hat{\sigma}_{t}^{g}$ can be
flexible, either the sum or other combinations of those from static or dynamic
GPs. Our work chooses that of the static GP as the posterior variance of the
adaptive GP, i.e., $\hat{\sigma}_{t}^{f}={\sigma}_{t}^{f}$ and
$\hat{\sigma}_{t}^{g}={\sigma}_{t}^{g}$, since typically the data quality for
identifying dynamic GP is low due to lack of excitation and non-uniform
scattering of applied currents $\tilde{a}_{1:t-1}$ that often cause
$\hat{\sigma}_{t}^{f}$ and $\hat{\sigma}_{t}^{f}$ to be excessively large.

With adaptive GP surrogates of constraint functions, the original problem
(13)-(15) can be approximated as

$\displaystyle\min_{\tilde{a}_{t}\in\mathcal{A}}~{}~{}||a_{t}-\tilde{a}_{t}||^{2}_{2},$
(22) $\displaystyle s.t.\quad$
$\displaystyle(\mu_{t}^{f}+\bar{\mu}_{t}^{f})+\kappa\sigma_{t}^{f}\leq\bar{T},$
(23)
$\displaystyle(\mu_{t}^{g}+\bar{\mu}_{t}^{g})+\kappa\sigma_{t}^{g}\leq\bar{V}.$
(24)

As before, the quantity of training data for dynamic GP is limited since only
the data from current episode is used. Additionally, non-uniform data acquired
during some episodes, particularly if they are sparse in certain regions, can
lead to large posterior uncertainties in those regions. This is because the GP
has less information to make predictions in those areas. In contrast, the
training data for static GP are more extensive and diverse than the dynamic
GP. Therefore, the static GP often has a more informative confidence interval.

Algorithm 2: Adaptive GP-based safe RL (based on TD3 algorithm)
---
1: | Follow (1)-(22) in Algorithm 1 to train static GPs $\hat{f}_{s}(\cdot)$, $\hat{g}_{s}(\cdot)$;
2: | for episode $m=M+1$ to END, repeat
3: |  t $\leftarrow$ 1;
4: |  Empty the data storage $\mathcal{R}$ for training dynamic GPs;
5: |  for $t=1$ to ending step of charging for episode $m$, repeat
6: |  Observe state $s_{t}$ and use actor network $\mu_{\theta}(\cdot)$ for action $a_{t}$;
7: |  if $t\leq n$, where $n$ is the timestep starting action projection:
8: |  Execute $a_{t}$, observe the next state $s_{t+1}$, and reward $r_{t}$;
9: |  else:
10: |  Train the dynamic GP models with $\mathcal{R}$;
11: |  Get adaptive GPs from static and dynamic GPs;
12: |  Predict $z_{t+1}$ using (21);
13: |  Solve the optimization in (22)-(24), to find $\tilde{a}_{t}$;
14: |  Execute $\tilde{a}_{t}$, observe $s_{t+1}$, and compute $r_{t}$;
15: |  Predict the variable of interest $\tilde{z}_{t+1}$ with static GP;
16: |  Calculate the residual, $\Delta_{t+1}=z_{t+1}-\tilde{z}_{t+1}$;
17: |  Store $[z_{t+1},a_{t-1},a_{t}]$ and $[\Delta_{t+1}]$ into $\mathcal{R}$;
18: |  Store the transition $\left(s_{t},\tilde{a}_{t},r_{t},s_{t+1}\right)$ into $\mathcal{D}$;
19: |  Update the networks using 10-19 in Algorithm 1;
20: |  $t\leftarrow t+1$, until $s_{t+1}$ is terminal;
21: |  end for
22: | end for

##  5 Simulation Results and Discussions

This section uses two case studies, one with fixed and the other with varying
operating conditions, to validate the effectiveness of the proposed methods.
For these studies, the PyBaMM [32], an open-source Python package, is used to
perform high-fidelity simulation of battery dynamics.

###  5.1 Safe RL-based fast charging with fixed conditions

This section evaluates the proposed action projection-based safe TD3 approach
for optimizing battery fast-charging protocol. It is compared with the
traditional TD3 algorithm without the added safety layer and the classical
CCCV method. Note that here the safe TD3 based on static GP in Section 4.2 is
referred to as “safe TD3” and the one based on adaptive GP in Section 4.3 as
“adaptive safe TD3”.

This study aims to optimize the charging policy to charge the battery from
$10\%$ to $80\%$ state-of-charge (SOC) in the shortest possible time while
ensuring that the temperature and voltage remain within operational
constraints. The maximum allowable temperature $\bar{T}$ and voltage $\bar{V}$
for this case are $45^{\degree}C$ and $4.3V$, respectively. Additionally, the
ambient temperature is kept constant at $25^{\degree}C$ throughout the
experiment. The charging current $a_{t}\in[0.05C,4.5C]$, where $C$ is the
charging rate. The sampling time is 10 seconds. For our experiment, the
environment state consists of SOC, voltage, temperature, and the current
applied at the previous time step:
$s_{t}=[SOC_{t},V_{t},T_{t},a_{t-1}]^{\top}$. At each timestep in an episode,
the actor network is given the state $s_{t}$ as input and outputs a charging
current $a_{t}$. The actor network has two hidden layers, each consisting of
128 nodes. The ReLU is used as the activation function and ADAM as the
optimizer. The critic network receives the state $s_{t}$ and action $a_{t}$ as
inputs and generates Q value to assess the policy. The critic network has a
similar architecture to the actor network. The hyper-parameters of the
proposed adaptive safe TD3 algorithm are fine-tuned by the trial-and-error
method and listed in Table 3.

Table 3: List of hyperparameters used by the adaptive safe TD3. Parameters | Values (Fixed condition) | Values (Varying condition)
---|---|---
Learning rate $\alpha$ for actors | 0.0005 | 0.00005
Learning rate $\beta$ for critics | 0.005 | 0.0005
Batch size $B$ | 64 | 64
Discount factor $\gamma$ in the return | 0.99 | 0.99
Polyak averaging coefficient $\rho$ | 0.006 | 0.006
Initial noise variance $\sigma^{2}$ | 0.3 | 0.3
Noise decay factor | 0.025 | 0.025
Length scale (RBF) | 1.0 | 1.0
Noise level (white kernel) | $10^{-5}$ | $10^{-5}$

The reward function is designed to penalize the charging time while respecting
the temperature and voltage constraints [17], defined as

$r_{t}=-1+\lambda_{1}[V_{t}-\bar{V}]^{+}+\lambda_{2}[T_{t}-\bar{T}]^{+},$ (25)

where the first term is the penalty on the charging time (steps). The second
and third terms penalize the voltage and temperature constraint violations.
The clip function $[\cdot]^{+}$ returns zero if constraints are satisfied, or
the amount of violations otherwise. The weights are selected as
$\lambda_{1}=15$, $\lambda_{2}=20$.

First, the accuracy of the static and adaptive GP models is tested in
predicting the next-step temperature and voltage under a static environment,
i.e., the ambient conditions are kept fixed throughout the experiments. In
this work, both the static and adaptive GP models use an RBF kernel along with
a white noise kernel. The RBF kernel learns the global trend between the input
and output variables (described in Section 4.2 and 4.3), whereas white noise
kernels capture the local fluctuations. The parameter values of the kernel
functions are provided in Table 3. Static GPs use the data from the first
$M=5$ RL episodes as the training data, and the trained static GP models are
used thereafter for the entire experiment. For adaptive GPs, one trains the
dynamic GPs from scratch at the beginning of each RL episode and keep updating
them using the latest data in that episode. Thus, the dynamic GPs can reflect
the underlying battery dynamics in real-time. Fig. 2 demonstrates the
performance of both static and adaptive GP models in predicting next-step
temperature and voltage of the battery under different currents to be
deployed. The ambient temperature is fixed as $25\degree C$. Fig. 2 shows that
both static and adaptive GP models can predict the voltage and temperature
accurately, and the true values of these variables fall within the uncertainty
bound of the posterior means. Under fixed ambient conditions, the static GP is
able to capture the dynamics of the battery well. For adaptive GP models with
fixed conditions, the residuals, $\Delta T_{t}$ and $\Delta V_{t}$, are almost
zero, and so are the predictions of the dynamic GP models. One critical
observation is that the adaptive GPs can equally perform well as the static
GPs when there are no changes in the ambient conditions.

Figure 2: The predicted next-step temperature (a) and voltage (b) by static
(red) and adaptive (blue) GP models against the true temperature and voltage
(dashed line) under different charging currents. The GP models are trained and
tested under a fixed ambient temperature $25^{\degree}C$. Solid lines:
posterior mean; Shaded areas: $\pm$3 standard deviations.

Fig. 3 shows the training results of traditional TD3 (green), safe TD3 (red),
and adaptive safe TD3 (blue) in optimizing battery fast-charging protocols.
All three methods use the immediate reward $r_{t}$ defined in (25), and the
episodic reward can be computed as $\sum_{t=1}^{K}r_{t}$, where $K$ is the
length of each episode. The episodic rewards in Fig. 3 (a) show that the
training of traditional, safe, and adaptive safe TD3 converges to similar
solutions. However, safe TD3 and adaptive safe TD3 exhibit a slower
convergence rate compared to the standard TD3 method. This slowdown can be
attributed to a potential discrepancy between the actions used for training
the critic and actor networks within the safe TD3 algorithms (Algorithm 1).
Specifically, line 30 and line 14 of Algorithm 1 employs projected actions for
updating the critic network, whereas the actor network receives gradients
based on the pre-projected actions (line 17). This mismatch could lead to the
generation of noisy Q-values for the pre-projected actions [33]. An attempt to
address this issue is implemented by incorporating the pre-projected action
$a_{t+1}$ based on the next state $s_{t+1}$ (Algorithm 1, line 13) into the
target value $y_{t}$ used for critic network updates. While partially
effective, this approach might introduce a slight delay in convergence, as
observed in Fig. 3 (a). Fig. 3 (b) shows the consumed charging steps of every
training episode from each method to charge the battery from 10%-80% SOC.
Overall, the optimization makes progress as the charging time of each episode
eventually decreases over episodes. Specifically, the ultimate charging time
optimized by the TD3, safe TD3, and adaptive safe TD3 is shown to be 16.83,
18.83, and 22 minutes, respectively. That is, protocols optimized by safe and
adaptive safe TD3 present a slower charging speed than that by traditional
TD3. Fig. 3 (c)-(d) presents the maximum temperature and voltage reached
during each episode, and the red dashed line indicates the allowable maximum
temperature and voltage. It is seen that the protocols delivered by TD3 can
lead to significant temperature and voltage violations in almost all episodes.
In contrast, the protocols by safe and adaptive safe TD3 can well respect the
constraints throughout the entire training. The only violations occur at the
beginning when the GP models are being trained, and action projections have
not been activated at that moment.

Figure 3: Training performance of traditional (green), safe (red), and
adaptive safe (blue) TD3 in optimizing battery fast-charging protocols. (a)
Cumulative rewards; (b) Charging time; (c) Maximum temperature; and (d)
Maximum voltage, of each training episode. Magenta dashed: the allowed upper
bounds of temperature and voltage.

Fig. 4 (a) illustrates the optimal charging profiles returned by traditional
(green), safe (red), and adaptive safe (blue) TD3, as well as classical CCCV
(malibu). For all protocols, a higher current is initially applied to the
battery and then gradually decreases. The red and blue dashed lines are the
raw currents from the actor networks, while the red and blue solid lines
represent the projected currents by safe and adaptive safe TD3, respectively.
Fig. 4 (a) shows that for the safe TD3 method (red), the raw charging current
stays close to the upper bound for approximately 100 steps. Such a large
current causes a steep rise in temperature and voltage, shown in red lines in
Fig. 4 (c)-(d). As a result, the temperature of the battery approaches the
temperature constraint (purple dashed) within 20 charging steps, and
continuing such a large charging current onward may cause a violation. Thus,
safe TD3 starts to project the current into the safety region after 20
charging steps to ensure no constraint violation. For adaptive safe TD3
(blue), the final current profile is more conservative than the safe TD3 and
only requires action projection during the last 30 episodes. Fig. 4 (b) shows
that the charging protocols based on safe and adaptive safe TD3 require a
slightly longer charging time than traditional TD3, however, the short
charging time from traditional TD3 is at the expense of significant constraint
violations (see Fig. 4 (c)-(d)). Finally, despite satisfying the operating
conditions, the CCCV method takes much longer to charge the battery than all
three presented RL methods.

Figure 4: The (a) charging current, (b) SOC, (c) temperature, and (d) voltage
profiles, obtained by deploying the optimized protocols from the traditional
(green), safe (red), adaptive safe (blue) TD3, and from classical CCCV
(malibu). Magenta dashed lines in (c) and (d): allowed upper bounds of the
temperature and voltage. Solid lines in (a): safe current profiles after
projection.

###  5.2 Safe RL-based fast charging with varying conditions

This section tests the performance of safe and adaptive safe TD3 for battery
fast-charging optimization under varying operating conditions. To this end,
the ambient temperature is gradually increased from $10^{\degree}C$ to
$36^{\degree}C$ with an increment of $0.145\degree C$ per episode starting
from the 100-th episode. The diffusion-limited SEI growth model is also
included to introduce battery degradation and aging [31]. The maximum
allowable temperature $\bar{T}$ and voltage $\bar{V}$ for this case are
$45^{\degree}C$ and $4.4V$, respectively, and the sampling time is 15 seconds.
The architecture of the actor and critic networks is similar to the previous
subsection.

Figure 5: The predicted next-step temperature (a) and voltage (b) by static
(red) and adaptive (blue) GP models against the true temperature and voltage
(dashed line) under different charging currents. The GP models are trained at
$10^{\degree}C$ and tested at $36^{\degree}C$ of the ambient temperature.
Solid lines: posterior mean; Shaded areas: $\pm$3 standard deviations.

This section first validates the effectiveness of static and adaptive GPs in
modeling battery dynamics under varying conditions. For static GPs, the
previously trained static GPs serve as the baseline (trained at $10\degree
C$), and for each episode, dynamic GPs, $\hat{f}_{\Delta}(\cdot)$ and
$\hat{g}_{\Delta}(\cdot)$, are constantly updated with the online data to
capture the real-time residuals. Such dynamic GPs are expected to capture the
changes in battery dynamics due to variations in ambient temperatures and
battery aging. Fig. 5 reveals that static GPs cannot capture changes in
battery dynamics and fail to give accurate predictions once operating
conditions change. That is, using a simple fixed constrained action space (can
be defined by a fixed battery model) may not guarantee system constraints once
the battery dynamics change. In contrast, adaptive GPs can well capture
battery dynamics under varying conditions and battery aging.

Figure 6: Training performance of safe (red) and adaptive safe (blue) TD3 for
optimizing battery fast-charging protocols. (a) Cumulative rewards; (b)
Charging time; (c) Maximum temperature; and (4) Maximum voltage, of each
training episode. Magenta dashed: the allowed upper bounds of temperature and
voltage. Vertical black dashed: the episode at which the battery operating
conditions start to change.

Fig. 6 shows the training curves for the safe (red) and adaptive safe (blue)
TD3 as the ambient conditions change and the battery degrades. The episodic
reward curves in Fig. 6 (a) for both methods show similar behaviors. The
environment change occurs at the 100-th episode due to the gradual increase in
the ambient temperature from $10\degree C$ to $36\degree C$ and the addition
of the SEI growth model. After that, the rewards for both methods start to
decrease since the charging current has to be less aggressive to ensure
constraints under an increased ambient temperature, thereby the charging
duration (see Fig. 6 (b)) is significantly enlarged that leads to lower
rewards. After 250 episodes, both methods experience a sudden drop in episodic
rewards. This may be due to the high ambient temperature that causes
significant changes in system dynamics and rapid battery degradation. The safe
TD3 cannot adjust to the new conditions and violates the temperature
constraint, as in Fig. 6 (c). As a result, the temperature constraint
violations lead to a large drop in episodic reward for safe TD3. In contrast,
the adaptive safe TD3 can learn the latest battery dynamics online and can
ensure constraints even when the environment drifts. However, this is at the
price of reducing charging current and thus increasing the charging time (Fig.
6 (b)), which also causes the rapid drop in episodic reward.

The optimal profiles of charging currents designed by safe TD3, adaptive safe
TD3, and classical CCCV are given in Fig. 7 (a). Both RL methods activate the
action projection from about the 10-th step. However, as the ambient
temperature increases and the battery ages, the static GP in safe TD3 cannot
properly construct the new safety region, resulting in inadequate action
reduction. This causes the safe TD3 to violate the temperature constraint, as
shown by the red lines in Fig. 7 (c). In contrast, the adaptive safe TD3
learns the new system dynamics online and can accurately construct updated
safety regions and properly project the action to respect both temperature and
voltage constraints, as shown in Fig. 7 (c)-(d). However, the obtained
projected actions are much more conservative than those from safe TD3,
resulting in a longer charging duration, see Fig. 7 (b). Finally, for the CCCV
charging protocol, the levels of constant current and constant voltage are
carefully tuned to sufficiently meet the constraints. The designed charging
profile is shown in malibu color of Fig. 7. One can see that with such a
design strategy, the resultant protocol cannot reach the target of 80% SOC
regardless of the charging duration (Fig. 7 (b)). In other words, the CCCV
protocol is unable to tackle such harsh ambient conditions with high
temperature and battery degradation.

Figure 7: The (a) charging current; (b) SOC; (c) temperature; and (d) voltage
profiles, obtained based on the optimized protocols from safe TD3 (red),
adaptive safe TD3 (blue), and from classical CCCV (malibu). The red dashed
lines in (c) and (d) indicate the allowed upper bounds of temperature and
voltage. The solid lines in (a) are the safe current profiles after
projection.

Figure 8: (a) Total computational time; and (b) Episodic computational time
for each individual component of the adaptive safe TD3 method: RL, adaptive
GP, and constrained optimization. The inner plot in (a) shows the
computational time during the initial 200 episodes.

Finally, the computational efficiency is essential for enabling the
scalability, practicality, and sustainability of the proposed RL method. To
this end, we determine the computation time for each element of the adaptive
safe TD3 method: RL computation (action selection by the agent, execution of
the action, and the evolution of the environment to next state), adaptive GP
modeling, and constrained optimization. Fig. 8 shows the total and episodic
computation cost by each component. Specifically, Fig. 8 (a) reveals that the
RL computation dominates the total cost, accounting for roughly 65.35% of the
total computation time. This is followed by the adaptive GP modeling (26.91%)
and constrained optimization (7.45%). Fig. 8 (b) further delves into the
computation cost by each component at every episode. When the battery is new
and the ambient temperature is low (e.g., the first 200 episodes as in the
inner plot in Fig. 8 (a)), the computation cost due to GP modeling (9.39%) and
constrained optimization (3.07%) is minimal since the favorable ambient
condition allows for faster charging (few times of GP modeling and raw action
as a good initial guess for the optimization). From Fig. 8 (b), as the
training progresses, the increase in ambient temperature and battery
degradation cause a harsh condition for battery fast-charging. The charging
protocol has to be conservative to avoid constraint violations, which leads to
increased computational cost for the GP modeling (more times of GP modeling)
and constrained optimization (raw action is no longer close to the projected
action, see Fig. 7).

Remark: For the proposed scheme, the RL agent will be retrained to match the
latest changes in adaptive GP models by continuous online interactions with
the environment. As in Fig. 6, in the presence of operating condition changes,
the training of RL always goes through a transient stage during which the
adaptive GP models are updated to capture latest conditions while the RL agent
learns to adapt to new system (battery) dynamics and adaptive GP models.

##  6 Conclusion

This work proposed an adaptive safe RL approach (with TD3 as the exemplary RL
algorithm) to optimize fast-charging protocols for batteries. Our approach
introduces a safety layer after the actor network to project any unsafe
actions into a safety region before being deployed to the environment. To
construct such safety regions, adaptive GP models are developed as surrogates
of the original black-box safety constraints. The proposed adaptive GP
consists of a static GP, which captures the overall trend of battery
variables, and a dynamic GP trained online with current episode data, which
learns about the changes in battery dynamics caused by ambient conditions or
battery aging. Extensive simulation studies were conducted based on the PyBaMM
battery simulator to validate the proposed methods. Simulation results show
that the adaptive GP models are effective in capturing battery dynamics when
the environment drifts. Moreover, our proposed adaptive safe RL-based fast
charging protocol can respect constraints with a high probability regardless
of fixed or varying operating conditions. In contrast, safe RL without
adaptivity can only ensure constraints under fixed conditions. Traditional
TD3-based protocols cannot strictly ensure constraint satisfaction, while the
CCCV-based methods are overly conservative. The proposed adaptive safe RL
approach is also of theoretical and practical significance for safe RL methods
to respect hard constraints in the presence of varying system dynamics.
Finally, although this work mainly considers the temperature and voltage
constraints of batteries, the proposed safe RL framework with adaptive GP-
based surrogate modeling can be well applied to ensure satisfactions of other
system constraints. Future work includes testing the proposed safe and
adaptive RL methods for battery fast-charging on actual battery test beds to
validate their effectiveness and robustness under real operating conditions.

##  7 Acknowledgments

M. A. Chowdhury and Saif S. S. Al-Wahaibi acknowledge the support of
Distinguished Graduate Student Assistantships (DGSA) from Texas Tech
University. Q. Lu acknowledges the new faculty startup funds from Texas Tech
University. The authors acknowledge the support from the National Science
Foundation under Grant No. 2340194.

## References

  * [1] Liu K, Li K, Peng Q, Zhang C.  A brief review on key technologies in the battery management system of electric vehicles.  _Frontiers of Mechanical Engineering_. 2019;14:47–64.
  * [2] Korkas CD, Baldi S, Yuan S, Kosmatopoulos EB.  An adaptive learning-based approach for nearly optimal dynamic charging of electric vehicle fleets.  _IEEE Transactions on Intelligent Transportation Systems_. 2017;19(7):2066–2075.
  * [3] Jana A, Woo SI, Vikrant K, García RE.  Electrochemomechanics of lithium dendrite growth.  _Energy & Environmental Science_. 2019;12(12):3595–3607.
  * [4] Keyser M, Pesaran A, Li Q, Santhanagopalan S, Smith K, Wood E, Ahmed S, Bloom I, Dufek E, Shirk M, et al.  Enabling fast charging–Battery thermal considerations.  _Journal of Power Sources_. 2017;367:228–236.
  * [5] Wei Z, Quan Z, Wu J, Li Y, Pou J, Zhong H.  Deep deterministic policy gradient-DRL enabled multiphysics-constrained fast charging of lithium-ion battery.  _IEEE Transactions on Industrial Electronics_. 2021;69(3):2588–2598.
  * [6] Severson KA, Attia PM, Jin N, Perkins N, Jiang B, Yang Z, Chen MH, Aykol M, Herring PK, Fraggedakis D, et al.  Data-driven prediction of battery cycle life before capacity degradation.  _Nature Energy_. 2019;4(5):383–391.
  * [7] Ahmed S, Bloom I, Jansen AN, Tanim T, Dufek EJ, Pesaran A, Burnham A, Carlson RB, Dias F, Hardy K, et al.  Enabling fast charging–A battery technology gap assessment.  _Journal of Power Sources_. 2017;367:250–262.
  * [8] Thakur AK, Sathyamurthy R, Velraj R, Saidur R, Pandey A, Ma Z, Singh P, Hazra SK, Sharshir SW, Prabakaran R, et al.  A state-of-the art review on advancing battery thermal management systems for fast-charging.  _Applied Thermal Engineering_. 2023;226:120303.
  * [9] Bose B, Garg A, Panigrahi B, Kim J.  Study on Li-ion battery fast charging strategies: Review, challenges and proposed charging framework.  _Journal of Energy Storage_. 2022;55:105507.
  * [10] Notten PH, het Veld JO, Van Beek J.  Boostcharging Li-ion batteries: A challenging new charging concept.  _Journal of Power Sources_. 2005;145(1):89–94.
  * [11] Purushothaman B, Landau U.  Rapid charging of lithium-ion batteries using pulsed currents: A theoretical analysis.  _Journal of The Electrochemical Society_. 2006;153(3):A533.
  * [12] Jiang B, Berliner MD, Lai K, Asinger PA, Zhao H, Herring PK, Bazant MZ, Braatz RD.  Fast charging design for Lithium-ion batteries via Bayesian optimization.  _Applied Energy_. 2022;307:118244.
  * [13] Klein R, Chaturvedi NA, Christensen J, Ahmed J, Findeisen R, Kojic A.  Optimal charging strategies in lithium-ion battery.  In: _Proceedings of the 2011 american Control Conference_. IEEE. 2011; pp. 382–387.
  * [14] Zou C, Hu X, Wei Z, Wik T, Egardt B.  Electrochemical estimation and control for lithium-ion battery health-aware fast charging.  _IEEE Transactions on Industrial Electronics_. 2017;65(8):6635–6645.
  * [15] Attia PM, Grover A, Jin N, Severson KA, Markov TM, Liao YH, Chen MH, Cheong B, Perkins N, Yang Z, et al.  Closed-loop optimization of fast-charging protocols for batteries with machine learning.  _Nature_. 2020;578(7795):397–402.
  * [16] Pozzi A, Moura S, Toti D.  Deep Learning-Based Predictive Control for the Optimal Charging of a Lithium-Ion Battery with Electrochemical Dynamics.  In: _2022 IEEE Conference on Control Technology and Applications (CCTA)_. IEEE. 2022; pp. 785–790.
  * [17] Hao Y, Lu Q, Wang X, Jiang B.  Adaptive Model-Based Reinforcement Learning for Fast Charging Optimization of Lithium-Ion Batteries.  _IEEE Transactions on Industrial Informatics_. 2023;.
  * [18] Sutton RS, Barto AG.  _Reinforcement Learning: An Introduction_.  MIT press. 2018.
  * [19] Park S, Pozzi A, Whitmeyer M, Perez H, Joe WT, Raimondo DM, Moura S.  Reinforcement learning-based fast charging control strategy for li-ion batteries.  In: _2020 IEEE Conference on Control Technology and Applications (CCTA)_. IEEE. 2020; pp. 100–107.
  * [20] Park S, Pozzi A, Whitmeyer M, Perez H, Kandel A, Kim G, Choi Y, Joe WT, Raimondo DM, Moura S.  A deep reinforcement learning framework for fast charging of li-ion batteries.  _IEEE Transactions on Transportation Electrification_. 2022;8(2):2770–2784.
  * [21] Chang F, Chen T, Su W, Alsafasfeh Q.  Control of battery charging based on reinforcement learning and long short-term memory networks.  _Computers & Electrical Engineering_. 2020;85:106670.
  * [22] Kim M, Lim J, Ham KS, Kim T.  Optimal charging method for effective Li-ion battery life extension based on reinforcement learning.  In: _Proceedings of the 38th ACM/SIGAPP Symposium on Applied Computing_. 2023; pp. 1659–1661.
  * [23] Yang Y, He J, Chen C, Wei J.  Balancing Awareness Fast Charging Control for Lithium-Ion Battery Pack Using Deep Reinforcement Learning.  _IEEE Transactions on Industrial Electronics_. 2023;.
  * [24] Elouazzani H, Elhassani I, Barka N, Masrour T.  Smart Adaptive Multi Stage Constant Current Fast Charging for Lithium Ions Batteries Based on Deep Reinforcement Learning.  _Available at SSRN 4218801_ ;.
  * [25] Fujimoto S, Hoof H, Meger D.  Addressing function approximation error in actor-critic methods.  In: _International Conference on Machine Learning_. PMLR. 2018; pp. 1587–1596.
  * [26] Lillicrap TP, Hunt JJ, Pritzel A, Heess N, Erez T, Tassa Y, Silver D, Wierstra D.  Continuous control with deep reinforcement learning.  _arXiv preprint arXiv:150902971_. 2015;.
  * [27] Wang J.  An intuitive tutorial to Gaussian processes regression.  _arXiv preprint arXiv:200910862_. 2020;.
  * [28] Richardson RR, Osborne MA, Howey DA.  Gaussian process regression for forecasting battery state of health.  _Journal of Power Sources_. 2017;357:209–219.
  * [29] Haarnoja T, Zhou A, Hartikainen K, Tucker G, Ha S, Tan J, Kumar V, Zhu H, Gupta A, Abbeel P, et al.  Soft actor-critic algorithms and applications.  _arXiv preprint arXiv:181205905_. 2018;.
  * [30] Levine S, Kumar A, Tucker G, Fu J.  Offline reinforcement learning: Tutorial, review, and perspectives on open problems.  _arXiv preprint arXiv:200501643_. 2020;.
  * [31] Xiong R, Pan Y, Shen W, Li H, Sun F.  Lithium-ion battery aging mechanisms and diagnosis method for automotive applications: Recent advances and perspectives.  _Renewable and Sustainable Energy Reviews_. 2020;131:110048.
  * [32] Sulzer V, Marquis SG, Timms R, Robinson M, Chapman SJ.  Python battery mathematical modelling (PyBaMM).  _Journal of Open Research Software_. 2021;9(1).
  * [33] Kasaura K, Miura S, Kozuno T, Yonetani R, Hoshino K, Hosoe Y.  Benchmarking actor-critic deep reinforcement learning algorithms for robotics control with action constraints.  _IEEE Robotics and Automation Letters_. 2023;8:4449–4456.