# Vehicle-to-Grid Fleet Service Provision
considering Nonlinear Battery Behaviors

Joshua Jaworski, Ningkun Zheng, , Matthias Preindl, , Bolun Xu J. Jaworski, N.
Zheng, M. Preindl and B. Xu are with Columbia University, NY, USA. Contact:
{jjj2146, nz2343, mp3501<EMAIL_ADDRESS>
###### Abstract

The surging adoption of electric vehicles (EV) calls for accurate and
efficient approaches to coordinate with the power grid operation. By being
responsive to distribution grid limits and time-varying electricity prices, EV
charging stations can minimize their charging costs while aiding grid
operation simultaneously. In this study, we investigate the economic benefit
of vehicle-to-grid (V2G) using real-time price data from New York State and a
real-world charging network dataset. We incorporate nonlinear battery models
and price uncertainty into the V2G management design to provide a realistic
estimation of cost savings from different V2G options. The proposed control
method is computationally tractable when scaling up to real-world
applications. We show that our proposed algorithm leads to an average of 35%
charging cost savings compared to uncontrolled charging when considering
unidirectional charging, and bi-directional V2G enables additional 18$\%$ cost
savings compared to unidirectional smart charging. Our result also shows the
importance of using more accurate nonlinear battery models in V2G controllers
and evaluating the cost of price uncertainties over V2G.

###### Index Terms:

Energy storage, stochastic optimal control, electric vehicle charging,
vehicle-to-grid.

##  I Introduction

The International Energy Agency’s (IEA) roadmap to achieve net zero greenhouse
gas (GHG) emissions by 2050 calls for an increase in the share of renewable
energy in total global power generation and global transport sector
electrification from 29% and 2% in 2020 to 90% and 45% in 2050, respectively
[1, 2]. However, increasing energy supply intermittency and electric vehicles
(EV) charging demand can put significant stress on the grid without adequate
management and control, such as peak charging demand during periods of low
wind and solar power generation [3]. Smart charging integrates external data,
such as distribution grid constraints or time-varying electricity prices, into
unidirectional (V1G) or bidirectional (V2G) power transfer management between
the grid and the EV charging station (EVCS) [4]. Smart charging and V2G
management have emerged as a key strategy to accelerate transportation
electrification to support an increasingly renewable-powered grid operation,
minimizing EV owners’ charging cost, and leading new business models and job
opportunities [3, 5, 6, 7].

As the cost of V2G-compatible chargers continues to decline [8], software
development becomes pivotal to efficiently aggregate EVs and optimally control
their V2G responses while meeting the designated charging targets. While
plenty of works have conducted techno-economic analyses (TEA) of V2G [9, 10,
11, 12], few have considered complicating factors that practical V2G
implementations must address. We group these factors into three categories.
The first is battery model nonlinearities, in which the battery voltage,
current, efficiency, and degradation depend on the state of charge (SoC).
Controlling EV batteries accurately according to their nonlinear
characteristics is crucial to strike a balance between ensuring battery
security and economic benefits [13]. The second is grid uncertainty, that the
distribution grid load and electricity prices are time-varying and uncertain
[7]. Uncertainties are often neglected in TEA primarily due to computation
difficulties, but practical V2G implementations must consider uncertainties in
price-response applications. The last is computational scalability, the V2G
management software must manage tens to hundreds of EVs without consuming
monstrous computing power. As we will later be shown in our results, while the
aggregate benefit of V2G is pivotal for future grids, economic saving for each
individual EV is not significant to justify investment in specialized
computing hardware.

This paper presents a computation-efficient V2G management framework and a
realistic case study integrating the aforementioned complicating factors in
practical V2G implementations. Our contributions include:

  * •

We propose a computation-efficient and scalable V2G management controller
which optimizes V2G charging using accurate nonlinear battery models under
stochastic electricity prices. Combining a stochastic dynamic programming
algorithm with a least-laxity first (LLF) scheduling algorithm [14], our
proposed V2G framework minimizes charging costs for EVCS to meet charging
targets and distribution grid limits.

  * •

Using real-world electricity price and EV charging behavior data, our paper
provides a first-of-its-kind case study to demonstrate cost savings an EVCS
can realistically achieve in various V2G settings.

  * •

Our case study compares uncontrolled charging, V1G, and V2G with and without
nonlinear storage models and price uncertainties. The results quantify the
impact of various charging and model options and guide EVCS planning and
technology developments.

The remainder of this paper is organized as follows. Section II presents the
literature review. Section III describes the system model and formulates the
EV charging cost minimization problem. Section IV presents the solution
algorithm to the formulated charging problem. Section V includes simulation
results and discussion. Section VI concludes the paper.

##  II Literature Review

Previous literature has proposed multiple heuristic, optimization, or
learning-based approaches to conduct smart charging mostly considering linear
battery models with constant power rating and efficiency [15, 16]. For
example, Liu et al. [17] formulate an EVCS controller as a bi-level program
and use a genetic algorithm to minimize charging costs under a time-of-use
(ToU) tariff. Similarly, Long et al. [9] use an ordinal optimization approach
to minimize EVCS operating cost under a ToU scheme but add aggregated EV
demand, maintenance costs, V2G capability, hydrogen storage, and renewable
energy generation to the formulation. Additionally, Cao et al. [10] propose a
custom actor-critic algorithm to minimize charging costs and peak charging
load for a V2G-enabled EVCS, which results in a 24% energy cost savings when
compared to uncontrolled charging.

While most prior V2G literature included benchmarks to demonstrate the
effectiveness of the proposed algorithm [17, 9, 10], they often assume
constant EV battery parameters (power ratings, efficiency) and do not include
a penalty term to minimize battery cycling. Lab experiments and real-world
data have shown that Li-ion battery power ratings and efficiency strongly
depend on SoC, especially in nickel-cobalt-based batteries, which are the most
common choice for EVs [18, 19, 20]. In EV smart charging or V2G applications,
which aim to provide high charging power with low-cost power conversion
hardware. Battery power rating and efficiency are sensitive to the SoC, and
the charge or discharge power must be carefully controlled to ensure battery
thermal security and reduce degradation rates [21]. A common protocol for EV
charging management is the CC–CV (constant current–constant voltage) method
[22, 23, 24], that the battery charges with constant current until reaching a
high SoC level and then gradually reduces the current to maintain a constant
charging voltage to prevent over-voltage damages. Modeling battery
characteristics such as CC–CV protocols in V2G management is critical to
maximize cost savings and ensure battery security, but it requires
representing battery power and efficiencies as functions of SoC instead of
using constant values, which introduces significant computation complexities
[25, 26, 27, 28, 13].

Some EVCS smart charging algorithms partially accounted for nonlinear
charging/discharging characteristics [11, 12, 29, 30], but few were able to
model all nonlinear factors in a computation efficient approach. Starting from
an EVCS profit maximization problem formulated as mixed-integer linear
programming (MILP), Mouli et al. [11] accounts for SoC-dependent power ratings
in which the maximum power drops linearly after 80% SoC. Ebrahimi et al. [12]
models battery degradation as dependent on both SoC and depth of discharge
(DoD), and Schwenk et al. [29] incorporated both nonlinear efficiency and
nonlinear degradation terms. Lee et al. [30] implemented an adaptive
scheduling algorithm that formulates custom objectives and constraints as a
convex program and computes an optimal charging schedule in real-time while
considering feeder limits and battery tail capacity reclamation at high SoCs
with a data-driven approach. However, the reviewed control solutions do not
model all power ratings, efficiencies, and cycling penalties/degradation as
nonlinear nor do they provide a method to incorporate different behavior
curves. This property will be pivotal for the adaptability of EVCS control
algorithms to fast-charging applications and manufacturer-customized EV
battery management systems (BMS).

Besides battery models, price uncertainty is another complicating factor that
may impact the EVCS cost estimations but was rarely studied in V2G due to
computation difficulties. Most literature on smart charging assumes perfect
price forecasts or predetermined ToU tariffs [17, 9, 10, 11, 12, 29, 30]. As
EV capacity surges, future V2G projects will most likely arbitrage real-time
electricity prices that are highly volatile and uncertain, and EVCS must
consider price uncertainties. Frendo et al. and Ahmad et al. [31, 32] use day-
ahead prices to forecast the next day’s LMPs, incorporating EVCS controller
using MILP formulation. Zhang et al. [33] develop the charging control deep
deterministic policy gradient, which models EV charging as a Markov decision
process (MDP) and optimizes user satisfaction and charging costs with the
output of a long short-term memory (LSTM) network that approximates sequential
energy price dynamics. Results from this research demonstrated the
significance of modeling price uncertainties, but the computing approach is
not scalable to address nonlinear storage models.

As summarized in Table I, the reviewed literature proposes EVCS algorithms
that partially account for nonlinear EV battery behavior and price
uncertainty. We close this research gap by developing a real-time EVCS V2G
control algorithm based on analytical nonlinear stochastic dynamic programming
(SDP) and a least-laxity first (LLF) scheduling approach that adds nonlinear
EV battery behavior and price uncertainty to the EVCS control formulation
while minimizing operating costs and complying with EV charging and battery
dynamics such as CC-CV charging profiles, facility power limits, and users’
charging targets. While our proposed method accounts for system non-linearity
and uncertainty, we demonstrate that it is scalable and computationally
tractable in practical scenarios.

TABLE I: Consideration of price uncertainty and battery parameters nonlinear behavior in smart charging control formulations |  Price Uncertainty |  Nonlinear Power Rating |  Nonlinear Efficiency |  Nonlinear Cycling Penalty / Degradation Cost
---|---|---|---|---
[17, 9, 10] | ✗ | ✗ | ✗ | ✗
[11, 30] | ✗ | ✓ | ✗ | ✗
[12] | ✗ | ✗ | ✗ | ✓
[29] | ✗ | ✗ | ✓ | ✓
[31, 32, 33] | ✓ | ✗ | ✗ | ✗
Proposed | ✓ | ✓ | ✓ | ✓

##  III System Model and Formulation

We take the perspective of a public EVCS operator whose electricity cost
settles using time-varying wholesale electricity real-time prices. The
objective of the EVCS is to minimize the electricity cost under either smart
charging (V1G) or V2G operation. The EVCS has enough chargers, all with the
same specifications so that no rejection-of-service event occurs during the
simulation time frame and complies with $L$, the maximum power rating of the
charging station, at all time steps.

###  III-A EV Charging Sessions

We consider a total of $K$ EVs accessing the EVCS during the considered
period, $\mathcal{K}=\\{1,..,K\\}$ is the set of EVs. EV $k\in\mathcal{K}$
arrives on time step $A_{k}$ with a starting SoC of $S_{k}$, and departs on
time $D_{k}$ with a charging target SoC $F_{k}$, with $D_{k}>A_{k}$ and
$F_{k}>S_{k}$. Thus we use $\mathcal{T}_{k}=\\{A_{k},\dotsc,D_{k}\\}$ to
denote the time frame of the current charging session of EV $k$. At every EV
arrival, the controller updates the tuple ($A_{k}$, $D_{k}$, $F_{k}$), which
is used as input for the proposed solution algorithm. Additionally, the
controller has access to every EV SoC at all times. We also assume the EVCS
does not have information or can predict the arrival of EVs, but each EV upon
arrival will inform the EVCS its departure time and charging target.

###  III-B Battery Nonlinear Behavior

All EVs are modeled as having the same battery capacity and SoC-dependent
charge power rating ($B_{k}(e_{t-1,k})$), discharge power rating
($P_{k}(e_{t-1,k})$), single-trip efficiency ($\eta_{k}(e_{t-1,k})$) and
discharge cost penalty ($c_{k}(e_{t-1,k})$) curves. As shown in Figure 1,
efficiency, and cycling penalty are modeled as quadratic functions of SoC
[18]. The power rating curves resemble a Tesla Model S fast charging curve
[34] with custom-defined CC-CV behavior at low and high SoC. The controller
has access to a low-resolution version of the battery nonlinear parameter
curves to use as an approximation during control policy computations. This
resolution gap simulates the challenge of EV battery behavior approximation in
online smart charging control.

Figure 1: Assumed EV battery SoC-dependent parameter curves for all EVs

###  III-C Price Prediction

The EVCS operator can access the system’s day-ahead price (DAP) information
and a DAP-based real-time price (RTP) prediction tool. For this work, RTP is
modeled as a 1st-order Markov process with 12 nodes per time step trained with
historical RTP-DAP bias price data as in setting DB-Dep in [35]. The proposed
solution algorithm uses the resulting prediction to provide a control decision
that accounts for future price uncertainty.

###  III-D Formulation

We start with formulating the EVCS V2G fleet management problem including
nonlinear battery models, charging station power limits, and a causality
control policy constraint with respect to time-varying electricity price
uncertainties.

The objective of the charging station is to minimize the cost of electricity
to charge up each EV including a discharge penalty to avoid frequent cycles
that accelerate battery degradation
$\displaystyle\min_{b_{t,k},p_{t,k}}\sum_{k\in\mathcal{K}}\sum_{t\in\mathcal{T}_{k}}\lambda_{t}(b_{t,k}-p_{t,k})+c_{k}(e_{t-1,k})p_{t,k}$
(1a) in which $b_{t,k}$ is the charging power of EV $k$ during time step $t$,
while $p_{t,k}$ is the discharge power. $\lambda_{t}$ is the time-varying
price of electricity. The second term reduces battery degradation by
introducing a discharge penalty cost $c_{k}(e_{t-1,k})$ as a function of SoC.
While the EVCS operator does not assume EV battery degradation cost, the
controller incorporates an EV battery discharge penalty to avoid excessive
cycling. Note that $b_{t,k}$ and $p_{t,k}$ have been normalized by the chosen
simulation time step and have units of energy.

Each EV is subject to the following power and energy constraints
($t\in\mathcal{T}_{k}$)

$\displaystyle 0\leq b_{t,k}\leq B_{k}(e_{t-1,k})$ (1b) $\displaystyle 0\leq
p_{t,k}\leq P_{k}(e_{t-1,k})$ (1c) $\displaystyle
e_{t,k}-e_{t-1,k}=-p_{t,k}/\eta_{k}(e_{t-1,k})+b_{t,k}\cdot\eta_{k}(e_{t-1,k})$
(1d) $\displaystyle\underline{E}_{k}\leq e_{t,k}\leq\overline{E}_{k}$ (1e)
$\displaystyle e_{A_{k}}=S_{k},\quad e_{D_{k}}\geq F_{k}$ (1f)

in which $e_{t,k}$ is the SoC of EV $k$ during time step $t$.
$\underline{E}_{k}$ is the EV minimal SoC limit and $\overline{E}_{k}$ is the
EV maximum SoC. (1b) and (1c) model the EV power rating, (1d) models the SoC
evolution, (1e) models the upper and lower SoC limits, and (1f) models the
starting SoC and SoC charging target.

The total charging and discharge power is subjected to the station power limit

$\displaystyle\sum_{k\in\mathcal{K}}p_{t,k}$ $\displaystyle\leq L$ (1g)
$\displaystyle\sum_{k\in\mathcal{K}}b_{t,k}$ $\displaystyle\leq L$ (1h)

The control policy must be causal (non-anticipatory) [36] and only depends on
past and current information

$\displaystyle\\{b_{t,k},p_{t,k}|k\in\mathcal{K},t\in\mathcal{T}_{k}\\}\in\text{Causal
Control Policies}$ (1i)

###### Remark 1

Generalization to different charging scenarios. (1) provides a generalized
formulation in different charging scenarios. In V2G, both $B_{k}$ and $P_{k}$
are non-zero, while in single-directional smart charging, or V1G, the EV will
not inject power into the grid and thus $P_{k}$ is set to zero. In cases when
assuming a linear battery model, $P_{k}$, $B_{k}$, $c_{k}$ and $\eta_{k}$ are
constants, while in nonlinear battery models, these parameters are dependent
on the SoC.

##  IV Solution Method

We take a two-step approach to solving the EVCS control problem (1). First, we
formulate the charging session of a single EV as a price arbitrage problem,
which satisfies all constraints in (1) except the EVCS power limit constraints
(1f) and (1g), and solve it using an analytical SDP method [35]. Second, we
aggregate the control policies from the EVs in an active session at the
current time step and prioritize the distribution of control signals according
to an LLF approach [14] to ensure compliance with the EVCS power limit.

###  IV-A Decomposition to Arbitrage Problems

The core of the proposed control method is a single energy storage price
arbitrage problem formulated as an SDP [35].

###### Remark 2

V2G decomposition. Because we do not assume the EV charging actions would
impact market clearing prices, we can decompose (1) into parallel arbitrage
problems by relaxing the charging station power limit constraint (1f) and
(1g), which are the only coupling factors among all EVs. Each resulting sub-
problem becomes an arbitrage problem that maximizes the arbitrage profit (or,
equivalently, minimizes the charging session electricity cost) while meeting
the final SoC target. We will discuss in later sections how we aggregate
results from all EVs to incorporate the EVCS power limits.

To solve the arbitrage sub-problem under the causality policy constraint, we
adopt an SDP approach with the following formulation (for simplicity, we omit
the EV index $k$ in the following formulation, but the following formulation
is for a single EV):

$\displaystyle Q_{t-1}(e_{t-1}\,|\,\lambda_{t})$
$\displaystyle=\max_{b_{t},p_{t}}\lambda_{t}(p_{t}-b_{t})-c_{k}(e_{t-1})p_{t}+V_{t}(e_{t}\,|\,\lambda_{t})$
(2a) $\displaystyle V_{t}(e_{t}\,|\,\lambda_{t})$
$\displaystyle=\mathbb{E}_{\lambda_{t+1}}\Big{[}Q_{t}(e_{t}\,|\,\lambda_{t+1})\Big{|}\lambda_{t}\Big{]}$
(2b) subject to (1b)–(1e).

We model the time-varying price $\lambda_{t}$ as an order-1 Markov process, as
described in III-C, in which the price distribution over a time period $t+1$
depends on the realized price over $t$. $Q_{t-1}(e_{t-1}\,|\,\lambda_{t})$ is
the maximized current period profit given the SoC at the start of the time
period, while $V_{t}(e_{t}\,|\,\lambda_{t})$ is the expected value function
representing the opportunity value of energy stored in the battery at the end
of the time period. In the context of an EV charging session,
$V_{t}(e_{t}\,|\,\lambda_{t})$ represents the minimum expected cost of the
remainder of the session based on the current SoC, current RTP, and the RTP
uncertainty model.

###  IV-B Solving Nonlinear Battery Models

We extend the solution approach from [35] by incorporating the SoC-dependency
of battery behavior parameters into the first-order optimality condition
expression as

$\displaystyle q_{t-1,i}(e_{t-1})=$
$\displaystyle\pi_{t,i}\Big{(}\frac{\partial p_{t}}{\partial
e_{t-1}}-\frac{\partial b_{t}}{\partial
e_{t-1}}\Big{)}-c_{k}(e_{t-1})\frac{\partial p_{t}}{\partial
e_{t-1}}-\frac{\partial c_{k}(e_{t-1})}{\partial e_{t-1}}p_{t}$
$\displaystyle+v_{t,i}(e_{t})\frac{\partial e_{t}}{\partial e_{t-1}}=0$ (3)

where $q_{t-1,i}(e_{t-1})$ is the derivative of
$Q_{t-1}(e_{t-1}\,|\,\lambda_{t})$, or the storage device’s marginal
opportunity value. According to the Karush-Kuhn-Tucker (KKT) conditions and
(1d), we obtain the following:

$\displaystyle{\partial p_{t}}/{\partial e_{t-1}}=$ (4a)
$\displaystyle\begin{cases}\eta+(p_{t}/\eta){\cdot}({\partial\eta}/{\partial
e})\indent\indent&\quad\text{if \eqref{p1_c2} not binding}\\\ {\partial
P}/{\partial e}&\quad\text{if \eqref{p1_c2} binding}\end{cases}$
$\displaystyle{\partial b_{t}}/{\partial e_{t-1}}=$ (4b)
$\displaystyle\begin{cases}-1/\eta-(b_{t}/\eta){\cdot}({\partial\eta}/{\partial
e})\indent\indent&\text{if \eqref{p1_c1} not binding}\\\ {\partial
B}/{\partial e}&\text{if \eqref{p1_c1} binding}\end{cases}$
$\displaystyle{\partial e_{t}}/{\partial e_{t-1}}=$ (4c)
$\displaystyle\begin{cases}0\qquad\qquad\qquad\qquad\qquad\qquad\text{if
\eqref{p1_c1} or \eqref{p1_c2} not binding}\\\ 1-(1/\eta){\cdot}({\partial
P}/{\partial e})+(1/\eta)^{2}{\cdot}({\partial\eta}/{\partial e})P\\\
+\eta{\cdot}({\partial B}/{\partial e})+B{\cdot}{\partial\eta}/{\partial
e}\qquad\quad\text{if \eqref{p1_c1} and \eqref{p1_c2} binding}\end{cases}$

By replacing the partial derivative expressions given by (4) in (IV-B) for
full power rating (binding) and partial (non-binding) charging or discharging
cases, we obtain an analytical marginal opportunity value function update
expression. Note that the expressions in (4) involve the optimization
variables $p_{t}$ and $b_{t}$. The expressions are solved by approximating
$p_{t}$ and $b_{t}$ by the power ratings $B$ and $P$ corresponding to the
current SoC. The full formulation of this equation is deferred to Appendix -B.

###  IV-C Arbitrage Policy

We use the developed analytical SDP algorithm to calculate a marginal
opportunity value function for each charging session for each EV. At each time
step, the control decision $p_{t}$ and $b_{t}$ for each connected EV can be
determined by comparing the corresponding marginal opportunity value function
and the observed realized RTP $\lambda_{t}$. The value difference between the
EV battery’s marginal value and the marginal grid price will trigger a
charging, discharging, or idling control signal. The full marginal value
function and control policy calculation methods are deferred to Appendix -C
and Appendix -B.

###### Remark 3

Lagrangian relaxation of the final SoC constraint. We apply a Lagrangian
relaxation to incorporate the final SoC charging target constraint into the
SDP by assuming an arbitrarily large penalty ($1000/MWh) for not achieving the
charging target. This enables the marginal value function corresponding to the
EV departing time step to act as an inverse activation function, with a linear
penalty cost to the battery SoC until reaching the specified charging target.
This enforces the battery to charge regardless of the price when approaching
the end of the charging session in V1G and V2G cases to meet the charging
target.

When applying this control policy with a nonlinear system assumption is that
the storage device parameters $B$, $P$, $\eta$, and c used in the elaboration
of the control decision are approximations of the real storage device behavior
parameters. The control decision will be calculated using a trained model from
historical price data and executed in the testing environment. The testing
environment will limit the control inputs to a range within the true storage
behavior constraints.

###  IV-D EV Fleet Control Simulation

The EVCS control algorithm incorporates nonlinear battery parameter curve
approximations, data provided by the EV users (SoC charging targets and
session duration), the marginal value function calculation algorithm, an LLF
prioritization step, and the control policy outlined in IV-C. The algorithm is
executed in real-time as follows:

  1. 1.

Set $t\to t+1$.

  2. 2.

Calculate the marginal value function for the EVs that arrive to the EVCS at
time $t$. For the value function computation, set $t$ as the starting time,
the provided session duration as the time horizon $T$ and the provided SoC
session target as the final SoC $e^{f}$. The resulting marginal value function
will be used as the basis of the control policy for its corresponding EV for
the duration of the current session.

  3. 3.

Identify all connected EVs. Compute the ratio of time elapsed in the current
session to the total session duration for each connected EV. Based on the
calculated ratios, sort the EVs in descending order. This step is aligned with
an LLF [14] scheduling approach, which prioritizes EVs with the least time to
achieve their target.

  4. 4.

In the order defined by step 3), execute the control policy for each EV as
described in IV-C (i.e., comparing the EV marginal value to the grid marginal
price). If the facility power limit is reached, set the remaining EVs’ power
control signals to zero and go to step 1).

  5. 5.

Go to step 1) until reaching the target simulation time.

Note that the control signals in step 4) will be truncated by the testing
environment if they are outside the range of the actual battery behavior
model. The LLF sorting step provides a lightweight solution to aggregate the
individual value function results and comply with the facility’s power limit.
This enables the modular nature of the algorithm components and prevents
exponential computation time growth as the number of EVs in the EVCS increase.

##  V Case Study

###  V-A Data and Experiment Design

We test the proposed control algorithm using the 2019 New York Independent
System Operator (NYISO) price data. Price uncertainty is modeled using a 1st-
order Markov process trained with 2016-2018 NYISO price data. We include
prices from four zones to demonstrate performance results in different price
patterns: NYC, LONGIL, NORTH, and WEST.

A 101-sample resolution version of the SoC-dependent battery parameter curves
shown in Figure 1 is considered the ground truth and used as the testing
environment. We assume that the controller has access to a 10-sample
resolution version of the same battery nonlinear parameter curves. Although
the proposed method can handle different parameter curves for each EV, we
assume identical parameter curves in this case study for simplicity. With this
resolution gap between the environment and the valuation process, we
demonstrate the effectiveness of the proposed algorithm in providing efficient
control with a limited amount of data to approximate the battery model.

We consider six scenarios to test the proposed algorithm and establish
benchmarks for comparison. The following two scenarios are used as benchmarks:

  1. 1.

PF (perfect forecast). We perform a deterministic optimization using real-time
prices in Julia/Gurobi. The optimization problem setup can be found in
Appendix -D. This benchmark scenario represents the lowest possible EVCS
operating cost.

  2. 2.

UC (uncontrolled charging). EVs start charging as soon as they arrive at the
EVCS and charge with the maximum allowable power rating at all times until
reaching their charging target. Facility power limits are fairly distributed
among actively charging EVs. The control logic is implemented in Julia. This
is the second benchmark case and represents an EVCS without a control policy
or V2G capability.

and the remaining scenarios are solved with the proposed custom algorithm
implemented in Julia:

  1. 1.

NL-V2G. We perform SDP control assuming V2G capability and approximate EV
battery and charger behavior with a 10-sample version of the nonlinear
parameter curves.

  2. 2.

NL-V1G. Similar to NL-V2G, but assuming no V2G capability.

  3. 3.

L-V2G. Similar to NL-V2G, but approximating EV battery and nonlinear charger
parameters with constant values. Power ratings are set to the nominal charger
capability, one-way charging and discharging efficiencies to 95% and marginal
battery degradation cost to $15MWh

  4. 4.

L-V1G. Similar to L-V2G, but assuming no V2G capability.

The EVCS consists of 21 bi-directional (unless otherwise noted by the scenario
being tested) 17.2 kW level-2 chargers and has a 150 kW power limit, leading
to an over-subscription ratio of 2.4. We assume 75 users have access to the
EVCS and all users own a 100 kWh EV. Users’ energy requested, and arrival and
departure times are obtained from the Caltech ACN dataset [37], specifically
using the 2019 JPL data with energy requests greater than 5 kWh. Figure 2
shows the distribution of arrival and departure times of the dataset that will
be used for the simulation. A starting SoC of 10% is assumed for all arrivals.
SDP control is performed in all scenarios using the described 1st-order Markov
process price prediction.

Figure 2: Histogram of EVCS users’ arrival and departure times from the ACN
(JPL 2019) Dataset.

###### Remark 4

Charging target compliance. For a charging session to be successful in any of
the scenarios, the control algorithm must achieve a final SoC within 5% of the
user’s SoC charging target. We define the charging target compliance
performance metric as the ratio between successful and total charging
sessions. Note that the user may input a target that is infeasible due to a
short session duration. The charging target compliance metric will exclude
these infeasible cases.

All computations were performed on a personal laptop with an Intel Core
i9-10885H 2.5GHz CPU and 32 GB memory. The benchmark (PF) using MILP is solved
using Gurobi [38], while the proposed algorithms and the EVCS simulation are
implemented in Julia.

###  V-B Cost Savings and Charging Target Compliance

Figure 3 shows a sample charging session using both NL-V2G and NL-V1G
scenarios. V2G achieves cost reduction by charging during low price periods
and capturing additional revenue through energy arbitrage, and V1G reduces
EVCS operating costs only through its smart charging capability. Note that
embedding the final SoC requirement in the value-to-go function calculation
results in successful charging sessions for all the shown cases while
minimizing total cost and EV battery cycling during the session.

(a) V2G (b) V1G

Figure 3: Performance of the proposed control policy in a 9-hour charging
session in both (a) NL-V2G and (b) NL-V1G scenarios Figure 4: EVCS operating
cost savings of all scenarios simulated with the proposed method for the
uncontrolled case. Perfect Forecast represents the highest possible EVCS
operating cost savings under the considered simulation conditions. Figure 5:
EVCS cumulative operating costs during 2019 for all simulated NYISO zones.
Perfect-forecast V2G/V1G, stochastic V2G/V1G and uncontrolled charging
scenarios are shown.

Figure 4 shows the EVCS operating cost savings achieved by the proposed
algorithm and the PF scenario as a percentage of the uncontrolled charging
(UC) scenarios as well as the charging compliance results across the four
considered NYISO zones. NL-V2G results in an average operating cost savings of
35% over UC, with savings reaching up to 56% in the WEST zone while
maintaining average charging compliance of 95%. The cost savings average drops
to 17% across zones if NL-V1G is used. This demonstrates the impact of bi-
directional charging capability on EVCS operating costs. Both L-V2G and L-V1G
result in increased average cost savings (37% and 18% respectively) compared
to NL-V2G and NL-V1G, but at the expense of charging compliance (83% for L-V2G
and 73% for L-V1G). L-V2G and L-V1G increased cost savings come from
estimating the nonlinear power ratings as constant, which causes mismatches
between the control signal and the actual power rating capability at a given
time step. This mismatch leads to EVs not achieving their session charging
target, which results in lower power purchased from the grid and reduced
charging target compliance. Additionally, Figure 5 shows the cumulative costs
of the V2G/V1G PF, NL-V2G, NL-V1G, and UC cases for the simulated year.

###  V-C V2G Energy Equivalent Mileage

Most EVs have a battery warranty based on the production time and the drive
mileage. To this end, V2G puts on additional discharges to the battery and may
accelerate the expiration of the manufacturer warranty. In this section, we
study how much additional energy is discharged in V2G and the equivalent
mileage consumption to understand how much stress V2G would put on battery
warranties.

Table II compares the total EVCS energy input (charged) and output
(discharged) in NL-V2G and NL-V1G scenarios across all NYISO zones. On
average, 7.3% of the total energy charged in the NL-V2G scenario is used for
discharging to arbitrage electricity prices. From these results, an equivalent
mileage value for an EV participating in a V2G EVCS can be estimated by
translating the energy output from the station to mileage through an EPA EV
range estimate [39]. This value becomes relevant when calculating the impact
of V2G on EV warranty, which is regulated to be the first of 8 years or
100,000 mi (15 years or 150,000 mi proposed in California) [40]. As a case
study, we use the Tesla Model X EPA estimated range of 348 mi for a full
charge (100 kWh) [41]. Subsequently, a fraction of yearly energy output from
the NYC EVCS (9.6 MWh), proportional to the number of charging sessions
corresponding to a particular EV (139 sessions out of a total of 2967), is
equivalent to a mileage value of 1565 mi. This corresponds to a 12% increase
in the average mileage driven per year (USDOT average is 13500 mi [42]), which
would lead to passing 100,000 mi approximately nine months earlier than the
baseline average mileage. Additionally, using the EVCS total energy output, an
EPA estimated range of 348 mi and the cost savings achieved by V2G result in
an incremental EVCS operating cost savings benefit of $0.125/kWh and
$0.036/mi. Using our proposed method to estimate V2G equivalent mileage under
different control policies would increase accuracy when performing long-term
EVCS TEA studies.

TABLE II: Energy Balances (in MWh) for a 1-year EVCS simulation Zone | NYC | NORTH | WEST | LONGIL
---|---|---|---|---
NL-V1G Charged | 152.44 | 153.38 | 153.40 | 152.59
NL-V2G Charged | 163.85 | 161.08 | 177.48 | 168.65
NL-V2G Discharged | 9.63 | 5.84 | 20.25 | 13.84

###  V-D Computation Times

Table III shows computation times for 1-year simulations of three cases all
using nonlinear battery models: 1) PF-MILP: nonlinear V2G optimization with
perfect price forecast formulated using MILP and solved using Gurobi; 2) PF-
DP: nonlinear V2G optimization with perfect price forecast solved with the
proposed algorithm, note that in this case there is no uncertainty, so the
proposed algorithm is essentially dynamic programming; 3) NL-V2G: nonlinear
V2G optimization solved using the proposed stochastic dynamic programming
algorithm. Hence, PF-MILP and PF-DP are deterministic, while SDP is
stochastic.

The computation time result shows the computation tractability of the proposed
algorithm in both deterministic and stochastic optimization. The comparison
between PF-MILP and PF-DP is an apple-to-apple comparison as both algorithms
solve a deterministic EVCS problem. PF-DP, the deterministic version of our
proposed algorithm yields an average result within 1.5% of the solution times
provided by the MILP formulation, while the computation time is around 60x
faster. Our proposed SDP-based algorithm achieves a computation time 7.5x
faster than the MILP. Note that SDP is solving multi-stage stochastic
optimization while the MILP is solving deterministic optimization, while both
use nonlinear battery models. Therefore, our proposed algorithm can also be
considered a faster, more efficient, and open-source alternative for solving
deterministic smart charging control case studies.

TABLE III: Computation times (in seconds) for a 1-year V2G simulation Zone | NYC | NORTH | WEST | LONGIL
---|---|---|---|---
PF-MILP | 5567 | 6123 | 6117 | 6214
PF-DP | 97 | 93 | 90 | 97
NL-V2G | 830 | 825 | 829 | 814

##  VI Conclusion

We proposed and tested an EVCS controller based on a nonlinear analytical
stochastic dynamic programming algorithm and least-laxity first scheduling.
Using historical prices from New York State, our proposed V2G algorithm
achieved 24% to 56% of EVCS operating cost savings compared to uncontrolled
charging while maintaining a 95% charging target compliance and accounting for
EV battery nonlinear behavior and price uncertainty in real-time. Our study
covers smart charging in which EVs are not discharged to the grid. Still, our
approach provides on average 17% cost savings by responding to grid price
variations compared to uncontrolled charging. We also show the importance of
considering nonlinear battery models in V2G optimization, in which the battery
power rating and efficiencies are dependent on the SoC, which is critical to
ensure the EV meets its charging target while responding to time-varying
prices. Finally, the proposed algorithm is open-source and not requiring any
thrid party solvers, while the computation time surpasses commercial solvers.
Hence, our approach is suitable for real-world implementations and scale-up
for large-scale EV fleet management.

In the future, we plan to improve the proposed approach in several directions.
The first is to integrate the solution method with data-driven probability
price prediction methods. The current solution to the V2G problem using
stochastic dynamic programming still requires a Markov process to be trained
using historical price data. However, designing the Markov process can be
complicated and limited by the quantity of historical price data. Second, our
result still shows that V2G is more likely to miss charging targets due to the
feeder capacity constraints, we will investigate approaches to manage the
charging constraint better and improve charging compliance. Finally, we will
test our proposed algorithm using more sophisticated charging scenarios such
as using heterogeneous EV fleets and considering local renewable generation
and study the connection between driving patterns and zone prices with the
control policy performance.

## References

  * [1] IEA, “Net zero by 2050,” _IEA_ , 2021. [Online]. Available: https://www.iea.org/reports/net-zero-by-2050
  * [2] B. Bilgin, P. Magne, P. Malysz, Y. Yang, V. Pantelic, M. Preindl, A. Korobkine, W. Jiang, M. Lawford, and A. Emadi, “Making the case for electrified transportation,” _IEEE Transactions on Transportation Electrification_ , vol. 1, pp. 4–17, 2015.
  * [3] P. Jaramillo, S. K. Ribeiro, P. Newman, S. Dhar, O. Diemuodeke, T. Kajino, D. Lee, S. Nugroho, X. Ou, A. H. Strømman, and J. Whitehead, “Transport. in ipcc, 2022: Climate change 2022: Mitigation of climate change. contribution of working group iii to the sixth assessment report of the intergovernmental panel on climate change,” _IPCC_ , 2022.
  * [4] NREL, “Transportation and mobility research: Electric vehicle smart charging at scale,” _NREL_ , 2022. [Online]. Available: https://www.nrel.gov/transportation/managed-electric-vehicle-charging.html
  * [5] M. Nour, S. M. Said, A. Ali, and C. Farkas, “Smart charging of electric vehicles according to electricity price,” in _2019 International Conference on Innovative Trends in Computer Engineering (ITCE)_ , 2019, pp. 432–437.
  * [6] J. K. Szinai, C. J. Sheppard, N. Abhyankar, and A. R. Gopal, “Reduced grid operating costs and renewable energy curtailment with electric vehicle charge management,” _Energy Policy_ , vol. 136, p. 111051, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S030142151930638X
  * [7] O. Sadeghian, A. Oshnoei, B. Mohammadi-ivatloo, V. Vahidinasab, and A. Anvari-Moghaddam, “A comprehensive review on electric vehicles smart charging: Solutions, strategies, technologies, and challenges,” _Journal of Energy Storage_ , vol. 54, p. 105241, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2352152X22012403
  * [8] M. Jahnes, L. Zhou, M. Eull, W. Wang, and M. Preindl, “Design of a 22kw transformerless ev charger with v2g capabilities and peak 99.5% efficiency,” _IEEE Transactions on Industrial Electronics_ , 2022.
  * [9] T. Long, Q.-S. Jia, G. Wang, and Y. Yang, “Efficient real-time ev charging scheduling via ordinal optimization,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 5, pp. 4029–4038, 2021.
  * [10] Y. Cao, H. Wang, D. Li, and G. Zhang, “Smart online charging algorithm for electric vehicles via customized actor–critic learning,” _IEEE Internet of Things Journal_ , vol. 9, no. 1, pp. 684–694, 2022.
  * [11] G. R. C. Mouli, M. Kefayati, R. Baldick, and P. Bauer, “Integrated pv charging of ev fleet based on energy prices, v2g, and offer of reserves,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 2, pp. 1313–1325, 2017.
  * [12] M. Ebrahimi, M. Rastegar, M. Mohammadi, A. Palomino, and M. Parvania, “Stochastic charging optimization of v2g-capable pevs: A comprehensive model for battery aging and customer service quality,” _IEEE Transactions on Transportation Electrification_ , vol. 6, no. 3, pp. 1026–1034, 2020.
  * [13] A. Gonzalez-Castellanos, D. Pozo, and A. Bischi, “Detailed li-ion battery characterization model for economic operation,” _International Journal of Electrical Power & Energy Systems_, vol. 116, p. 105561, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0142061519315765
  * [14] Y. Nakahira, N. Chen, L. Chen, and S. H. Low, “Smoothed least-laxity-first algorithm for ev charging,” in _Proceedings of the Eighth International Conference on Future Energy Systems_ , ser. e-Energy ’17. New York, NY, USA: Association for Computing Machinery, 2017, p. 242–251. [Online]. Available: https://doi.org/10.1145/3077839.3077864
  * [15] R. Fachrizal, M. Shepero, D. van der Meer, J. Munkhammar, and J. Widén, “Smart charging of electric vehicles considering photovoltaic power production and electricity consumption: A review,” _eTransportation_ , vol. 4, p. 100056, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2590116820300138
  * [16] H. M. Abdullah, A. Gastli, and L. Ben-Brahim, “Reinforcement learning based ev charging management systems–a review,” _IEEE Access_ , vol. 9, pp. 41 506–41 531, 2021.
  * [17] J. Liu, G. Lin, S. Huang, Y. Zhou, Y. Li, and C. Rehtanz, “Optimal ev charging scheduling by considering the limited number of chargers,” _IEEE Transactions on Transportation Electrification_ , vol. 7, no. 3, pp. 1112–1122, 2021.
  * [18] X. Wu, W. Wang, and J. Du, “Effect of charge rate on capacity degradation of lifepo4 power battery at low temperature,” _International Journal of Energy Research_ , vol. 44, no. 3, pp. 1775–1788, 2020. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/er.5022
  * [19] R. Zhang, B. Xia, B. Li, L. Cao, Y. Lai, W. Zheng, H. Wang, W. Wang, and M. Wang, “A study on the open circuit voltage and state of charge characterization of high capacity lithium-ion battery under different temperature,” _Energies_ , vol. 11, no. 9, 2018. [Online]. Available: https://www.mdpi.com/1996-1073/11/9/2408
  * [20] Y. Preger, H. M. Barkholtz, A. Fresquez, D. L. Campbell, B. W. Juba, J. Romàn-Kustas, S. R. Ferreira, and B. Chalamala, “Degradation of commercial lithium-ion cells as a function of chemistry and cycling conditions,” _Journal of The Electrochemical Society_ , vol. 167, no. 12, p. 120532, 2020.
  * [21] L. Zhou, M. Jahnes, M. Eull, W. Wang, G. Cen, and M. Preindl, “Robust control design for ride-through/trip of transformerless onboard bidirectional ev charger with variable-frequency critical-soft-switching,” _IEEE Transactions on Industry Applications_ , 2022.
  * [22] B. Tar and A. Fayed, “An overview of the fundamentals of battery chargers,” in _2016 IEEE 59th International Midwest Symposium on Circuits and Systems (MWSCAS)_ , 2016, pp. 1–4.
  * [23] S. Habib, M. M. Khan, F. Abbas, L. Sang, M. U. Shahid, and H. Tang, “A comprehensive study of implemented international standards, technical challenges, impacts and prospects for electric vehicles,” _IEEE Access_ , vol. 6, pp. 13 866–13 890, 2018.
  * [24] M. Eull, L. Zhou, M. Jahnes, and M. Preindl, “Bidirectional non-isolated fast charger integrated in the electric vehicle traction drivetrain,” _IEEE Transactions on Transportation Electrification_ , vol. 8, pp. 180–195, 2021.
  * [25] H. Pandžić and V. Bobanac, “An accurate charging model of battery energy storage,” _IEEE Transactions on Power Systems_ , vol. 34, no. 2, pp. 1416–1426, 2019.
  * [26] M. Farag, M. Fleckenstein, and S. Habibi, “Continuous piecewise-linear, reduced-order electrochemical model for lithium-ion batteries in real-time applications,” _Journal of Power Sources_ , vol. 342, pp. 351–362, 2017\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0378775316317396
  * [27] A. Sakti, K. G. Gallagher, N. Sepulveda, C. Uckun, C. Vergara, F. J. de Sisternes, D. W. Dees, and A. Botterud, “Enhanced representations of lithium-ion batteries in power systems models and their effect on the valuation of energy arbitrage applications,” _Journal of Power Sources_ , vol. 342, pp. 279–291, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S037877531631758X
  * [28] G. Rancilio, A. Lucas, E. Kotsakis, G. Fulli, M. Merlo, M. Delfanti, and M. Masera, “Modeling a large-scale battery energy storage system for power grid application analysis,” _Energies_ , vol. 12, no. 17, 2019. [Online]. Available: https://www.mdpi.com/1996-1073/12/17/3312
  * [29] K. Schwenk, S. Meisenbacher, B. Briegel, T. Harr, V. Hagenmeyer, and R. Mikut, “Integrating battery aging in the optimization for bidirectional charging of electric vehicles,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 6, pp. 5135–5145, 2021.
  * [30] Z. J. Lee, G. Lee, T. Lee, C. Jin, R. Lee, Z. Low, D. Chang, C. Ortega, and S. H. Low, “Adaptive charging networks: A framework for smart electric vehicle charging,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 5, pp. 4339–4350, 2021.
  * [31] O. Frendo, N. Gaertner, and H. Stuckenschmidt, “Real-time smart charging based on precomputed schedules,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 6, pp. 6921–6932, 2019.
  * [32] F. Ahmad, M. S. Alam, S. M. Shariff, and M. Krishnamurthy, “A cost-efficient approach to ev charging station integrated community microgrid: A case study of indian power market,” _IEEE Transactions on Transportation Electrification_ , vol. 5, no. 1, pp. 200–214, 2019.
  * [33] F. Zhang, Q. Yang, and D. An, “Cddpg: A deep-reinforcement-learning-based approach for electric vehicle charging control,” _IEEE Internet of Things Journal_ , vol. 8, no. 5, pp. 3075–3087, 2021.
  * [34] InsideEVs, “Tesla model s plaid fast charging results,” 2021. [Online]. Available: https://insideevs.com/news/515641/tesla-models-plaid-charging-analysis/
  * [35] N. Zheng, J. J. Jaworski, and B. Xu, “Arbitraging variable efficiency energy storage using analytical stochastic dynamic programming,” _IEEE Transactions on Power Systems_ , pp. 1–1, 2022.
  * [36] A. Shapiro, D. Dentcheva, and A. Ruszczynski, _Lectures on stochastic programming: modeling and theory_. SIAM, 2021.
  * [37] Z. J. Lee, T. Li, and S. H. Low, “ACN-Data: Analysis and Applications of an Open EV Charging Dataset,” in _Proceedings of the Tenth International Conference on Future Energy Systems_ , ser. e-Energy ’19, 2019.
  * [38] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2022. [Online]. Available: https://www.gurobi.com
  * [39] EPA, “Fuel economy - all-electric vehicles,” 2022. [Online]. Available: https://www.fueleconomy.gov/feg/evtech.shtml
  * [40] California Air Resources Board, “California code of regulations, section 1962.4: Zero-emission vehicle standards for 2026 and subsequent model year passenger cars and light-duty trucks,” 2021. [Online]. Available: https://ww2.arb.ca.gov/sites/default/files/2021-12/draft%20zero%20emission%20vehicle%20regulation%201962.4%20posted.pdf
  * [41] Tesla, “Tesla model x,” 2022. [Online]. Available: https://www.tesla.com/modelx
  * [42] USDOT Federal Highway Administration, “Average annual miles per driver by age group,” 2022. [Online]. Available: https://www.fhwa.dot.gov/ohim/onh00/bar8.htm
  * [43] B. Xu, M. Korpås, and A. Botterud, “Operational valuation of energy storage under multi-stage price uncertainties,” in _2020 59th IEEE Conference on Decision and Control (CDC)_ , 2020, pp. 55–60.

###  -A Stochastic Price Arbitrage using Markov Process

We discretize the stochastic real-time price as a first-order Markov Process
with N nodes per time step ($\\{\pi_{t,i}|\,t\in\mathcal{T},\,i\in N\\}$), and
a stage transition probability $\rho_{i,j,t}$ indicating the probability of
transitioning from price node $i$ at time period $t$ to price node $j$ at time
period $t+1$. Thus, we reformulate the stochastic arbitrage problem in (2)
using discretized Markov price processes as

$\displaystyle\begin{split}Q_{t-1,i}(e_{t-1})&=\max_{b_{t},p_{t}}\pi_{t,i}\cdot(p_{t}-b_{t})-cp_{t}+V_{t,i}(e_{t})\end{split}$
(5a) $\displaystyle V_{t,i}(e_{t})$
$\displaystyle=\sum_{j\in\mathcal{N}}\rho_{i,j,t}\cdot Q_{t,j}(e_{t})$ (5b)
subject to the EV constraints (1b)–(1e).

###  -B Value Function Computation

Our proposed algorithm is based on the following result that updates
$q_{t-1,i}$ from $v_{t,i}$, where $q_{t-1,i}$ is the derivative of $Q_{t-1,i}$
and $v_{t,i}$ is the derivative of $V_{t,i}$ [43]

$\displaystyle q_{t-1,i}(e)=$ (6)
$\displaystyle\begin{cases}[1+\eta{\cdot}{\partial B}/{\partial
e}+B{\cdot}{\partial\eta}/{\partial
e}]{\cdot}v_{t,i}(e+B\eta)-\pi{\cdot}{\partial B}/{\partial e}&\\\
\indent\indent\text{if $\pi_{t,i}\leq v_{t,i}(e+B\eta){\cdot}\eta$}\\\
\pi_{t,i}{\cdot}[1/\eta+(B/\eta){\cdot}({\partial\eta}/{\partial e})]&\\\
\indent\indent\text{if $v_{t,i}(e+B\eta){\cdot}\eta<\pi_{t,i}\leq
v_{t,i}(e){\cdot}\eta$}\\\ v_{t,i}(e)&\\\ \indent\indent\text{if
$v_{t,i}(e){\cdot}\eta<\pi_{t,i}\leq[v_{t,i}(e)/\eta+c]^{+}$}\\\
(\pi_{t,i}-c){\cdot}[\eta+(P/\eta){\cdot}({\partial\eta}/{\partial
e})]-P{\cdot}{\partial c}/{\partial e}&\\\ \indent\indent\text{if
$[v_{t,i}(e)/\eta+c]^{+}<\pi_{t,i}$}\text{$\leq[v_{t,i}(e-P/\eta)/\eta+c]^{+}$}\\\
[1-(1/\eta){\cdot}({\partial P}/{\partial
e})+{P{\cdot}(1/\eta)}^{2}{\cdot}({\partial\eta}/{\partial
e})]{\cdot}v_{t,i}(e-P/\eta)\\\ +(\pi_{t,i}-c){\cdot}{\partial P}/{\partial
e}-P{\partial c}/{\partial e}&\\\ \indent\indent\text{if
$\pi_{t,i}>[v_{t,i}(e-P/\eta)/\eta+c]^{+}$}\end{cases}$

Note that $\eta$, B, P and c depend on the starting SoC e as formulated in
(2). The function form of these nonlinear parameters in (6) is omitted to
simplify the mathematical representation.

We restate the solution algorithm to compute the value function from [35]
which enforces a final SoC value higher than a given threshold $e^{f}$. In
this algorithm implementation we discretize $v_{t,i}$ into a set of
$\mathcal{M}$ segments with value and SoC pairs

$\displaystyle\hat{v}_{t,i}=\\{\nu_{t,i,m}|m\in\mathcal{M}\\}$ (7)

associated with equally spaced SoC segments $e_{t,m}$. In our implementation,
we discretized the SoC into 1000 segments. The valuation algorithm is listed
as following

  1. 1.

Set $T\to t$ to start from the last time period; initialize the final value-
to-go function segments $\nu_{T,m}$ to zeros for $e_{T,m}>e^{f}$ and to a very
high value (we use $1000/MWh) for $e_{T,m}\leq e^{f}$. Note that the final
value function does not depends on price nodes.

  2. 2.

Go to the earlier time step by setting $t-1\to t$.

  3. 3.

During period $t$, go through each price node for $i\in\mathcal{N}$ and value
function segment $m\in\mathcal{M}$. Update the charge and discharge efficiency
corresponding to the SoC segment, compute (1) and store $q_{t-1,i}(e)$; note
that here $q_{t-1,i}(e_{t-1,m})$ is also discretized with the same granularity
as the value function.

  4. 4.

Calculate the value function of the previous time step as

$\displaystyle\nu_{t-1,i,m}=\textstyle\sum_{j\in\mathcal{N}}\rho_{i,j,t}q_{t-1,j}(e_{t,m})$
(8)

which is the derivative of (5b).

  5. 5.

Go to step 2) until reaching the first time step.

###  -C Control Policy for Single Storage Device

We restate the control policy from [35]. After the value function computation
step, control can be executed by responding to realized market prices and
looking for the closest price node $\pi_{t,i}$ such that
$\underline{\pi}_{t,i}\leq\lambda_{t}<\overline{\pi}_{t,i}$, then the storage
control decision is updated as

$\displaystyle p_{t}$ $\displaystyle=\min\\{\hat{p}_{t},e_{t-1}\eta\\}$ (9a)
$\displaystyle b_{t}$ $\displaystyle=\min\\{\hat{b}_{t},(E-e_{t-1})/\eta\\}$
(9b) where $\hat{p}_{t}$ and $\hat{b}_{t}$ are calculated as
$\displaystyle\\{\hat{p}_{t},\hat{b}_{t}\\}=$
$\displaystyle\begin{cases}\\{0,B\\}&\text{if $\lambda_{t}\leq
v_{t,i}(e+B\eta)\eta$}\\\ \\{0,\alpha\\}&\text{if
$v_{t,i}(e+B\eta)\eta<\lambda_{t}\leq v_{t,i}(e)\eta$}\\\ \\{0,0\\}&\text{if
$v_{t,i}(e)\eta<\lambda_{t}\leq[v_{t,i}(e)/\eta+c]^{+}$}\\\
\\{\beta,0\\}&\text{if $[v_{t,i}(e)/\eta+c]^{+}<\lambda_{t}$}\\\
&\quad\text{$\leq[v_{t,i}(e-P/\eta)/\eta+c]^{+}$}\\\ \\{P,0\\}&\text{if
$\lambda_{t}>[v_{t,i}(e-P/\eta)/\eta+c]^{+}$}\end{cases}$ (9c) in which
$\alpha$ and $\beta$ are given as follows $\displaystyle\alpha$
$\displaystyle=(v^{-1}_{t,i}(\lambda_{t}/\eta)-e_{t-1})/\eta$ (9d)
$\displaystyle\beta$
$\displaystyle=(e_{t-1}-v^{-1}_{t,i}((\lambda_{t}-c)\eta))/\eta$ (9e) where
$v^{-1}_{t,i}$ is the inverse function of $v_{t,i}$.

(9a) and (9b) enforce the battery SoC constraints over the discharge
$\hat{p}_{t}$ and charge $\hat{b}_{t}$ decisions. (9c) calculates control
decisions and following the same principle as to (2) but use the observed
price $\lambda_{t}$ instead of the price nodes $\pi_{t,i}$. The control policy
uses the approximated values of B, P, $\eta$ and c (all 4 could be nonlinear)
corresponding to the current SoC.

###  -D MILP Formulation

We start by showing below the multi-period arbitrage formulation, which is
equivalent to our proposed stochastic dynamic programming formulation if
assuming a deterministic price process $\pi_{t}$ and also add a penalty
function $C_{end}$ and a weighing parameter $\alpha$ to represent the cost
associated with missing a final SoC target (the choice of $\alpha$ must
achieve a comparable charging compliance target performance as the
deterministic dynamic programming solution):

$\displaystyle\max_{p_{t},b_{t}}\quad-\alpha
C_{end}+\sum^{T}_{t}\pi_{t}{\cdot}(p_{t}-b_{t})-cp_{t}\textbf{
s.t.}\quad\text{(1c),(1e),(1f)}$

We modify this model to a MILP model for the variable power rating, discharge
penalty and efficiency benchmark calculation with ten SoC-efficiency segment
pairs as

$\displaystyle\max_{p_{k,t},b_{k,t}}\quad-\alpha
C_{end}+\sum^{T}_{t}\sum_{k}^{K}\pi_{t}{\cdot}(p_{k,t}-b_{k,t})-c_{k}p_{k,t}$
(11a) $\displaystyle\textbf{s.t.}\quad 0\leq\sum_{k}^{K}b_{k,t}/B_{k}\leq
1,\;0\leq\sum_{k}^{K}p_{k,t}/P_{k}\leq 1$ (11b) $\displaystyle
e_{k,t}-e_{k,t-1}=-p_{k,t}/\eta_{k}+b_{k,t}\eta_{k}$ (11c) $\displaystyle
E_{1}u_{1,t}\leq e_{1,t}\leq E_{1}$ (11d) $\displaystyle E_{k}u_{k,t}\leq
e_{k,t}\leq E_{k}u_{k-1,t},\quad\forall k\in\\{2,...,K-1\\}$ (11e)
$\displaystyle 0\leq e_{K,t}\leq E_{K}u_{K-1,t}$ (11f) $\displaystyle
C_{end}=(\sum_{k}^{K}e_{k,T}-e^{f})^{2}$ (11g)

where $k$ is the index of the nonlinear power rating, efficiency and cycling
penalty approximation segments (10 segments in this case). (11a) is the
objective function which sums up all segments. (11b) and (11c) are the power
rating constraints and energy storage evolution constraints implement on all
segments. (11d)-(11f) model the piece-wise linear approximation to the battery
nonlinear parameter curves with binary variables $u_{k,t}$, which enforce the
lower SoC segment must be full before upper SoC segments can take on non-zero
values. For this paper’s simulations, an $\alpha$ value of 1000000 is used.

In the EVCS context, this optimization step fits in step 2) of the solution
algorithm proposed in 4.3 (it replaces the value function calculation step).
Every time an EV arrives at the EVCS, the MILP will provide an optimal
schedule for the EV to follow for the remainder of the session. If the EV is
deviated from the provided optimal schedule due to a low priority assigned
during the LLF scheduling step, a new optimal schedule for the remainder of
the session will be computed.


# Generation and propagation of solitary waves in nematic liquid crystals

Xingzhou Tang  Pritzker School of Molecular Engineering, University of
Chicago, Chicago, Illinois 60637, U.S.A.  Ali Mozaffari  Pritzker School of
Molecular Engineering, University of Chicago, Chicago, Illinois 60637, U.S.A.
Noe Atzin  Pritzker School of Molecular Engineering, University of Chicago,
Chicago, Illinois 60637, U.S.A.  Soumik Das  Smith School of Chemical and
Biomolecular Engineering, Cornell University, Ithaca, New York 14853, U.S.A.
Nicholas L. Abbott  Smith School of Chemical and Biomolecular Engineering,
Cornell University, Ithaca, New York 14853, U.S.A.  Juan J. de Pablo  Pritzker
School of Molecular Engineering, University of Chicago, Chicago, Illinois
60637, U.S.A.

###### Abstract

Solitons in nematic liquid crystals offer intriguing opportunities for
transport and sensing in microfluidic systems. Little is known about the
elementary conditions that are needed to create solitons in nematic materials.
In this work, theory, simulations and experiments are used to study the
generation and propagation of solitary waves (or ”solitons”) in nematic liquid
crystals upon the application of an alternating current (AC) electric field.
We find that these solitary waves exhibit ”butterfly”-like or ”bullet”-like
structures that travel in the direction perpendicular to the applied electric
field. Such structures propagate over long distances without losing their
initial shape. The theoretical model adopted here serves to identify some of
the key requirements that are needed to generate solitons in the absence of
electrostatic interactions. These include surface imperfections that introduce
a twist in the director, unequal elastic constants, and negative anisotropic
dielectric permittivity. The results of simulations are shown to be in good
agreement with our own experimental observations, serving to establish the
validity of the theoretical concepts advanced in this work.

##  I Introduction

Solitons are travelling wave packets that propagate at constant speed over
long distances without losing their shape. A non-linear feedback mechanism
during wave propagation serves to minimize their expansion (or dispersion)
Dauxois2006 ; Manton2004 . Solitons and solitary waves are abundant in nature;
examples include water solitons in narrow channels, or the collective motion
of proteins and DNA. They also arise in optical fibers, in magnets, and in
nuclear physics heimburg2005soliton ; appali2012comparison ;
andersen2009towards ; kartashov2019 . Several recent studies have reported the
formation of solitons in liquid crystals (LCs) when subjected to an applied AC
electric field Lam1992 . Such studies have also outlined potential
applications as long-distance information carriers for optical signals, or
particles and chemical species within nematic LCs.

Three-dimensional solitons can be generated in a nematic liquid crystal that
exhibits a negative dielectric permittivity anisotropy upon application of an
alternating electric field Li2018 ; Aya2020 . Such solitons are usually
observed in the vicinity of a surface imperfection, dust particles, or the
edges of the electrodes, and they tend to travel over long distances in the
direction perpendicular to the director field. Experiments indicate that, at
low frequencies, solitons start as butterfly-like structures with quadrupolar
symmetry. At higher applied fields and frequencies, the butterfly structure
can emit bullet-like structures - solitons - having a curved director field.
Such bullets can propagate at high speeds. In this work we introduce a model
for the tensorial order parameter that is shown to be capable of predicting
the formation of these solitary waves in good agreement with experimental
observations.

The results of simulations of such a model reveal that the key ingredients for
soliton generation are (i) a negative anisotropy of permittivity, (ii) an
irregularity that serves to nucleate solitons, and (iii) an AC field whose
intensity and frequency must fall within a narrow range of values
das2022programming . The results presented here are compared to experimental
data for 4’butyl-4-heptyl-bicyclohexyl-4-carbonitrile $(CCN-47)$ Zawadzki2012
; Lucchetti2021  \- a liquid crystal for which detailed experimental
observations have recently become available.

##  II Theory and simulation

We consider a simple 3D thermotropic nematic liquid crystal. At each point in
the material, the local orientational order is described in terms of the
director field $\hat{\mathbf{n}}(\mathbf{r},t)$ or the nematic order tensor
$Q_{ij}(\mathbf{r},t)$. The local fluid flow velocity is given by
$\mathbf{v}(\mathbf{r},t)$. To solve the underlying model, we follow the
method of Ref. Tang2017 ; Tang2019  in 3 dimensions. The magnitude and the
direction of the nematic order are allowed to vary, so that solitary waves and
topological defects are able to form and move freely. We quantify the nematic
order through a traceless, symmetric tensor of the form

$Q_{ij}=s\left(\frac{3}{2}n_{i}n_{j}-\frac{1}{2}\delta_{ij}\right).$ (1)

The free energy of the system can be expressed as

$\displaystyle F=\int d^{2}r\biggl{[}$
$\displaystyle-\frac{1}{2}a|Q|^{2}+\frac{1}{3}b|Q|^{3}+\frac{1}{4}c|Q|^{4}$
$\displaystyle+\frac{1}{2}\left(K_{11}-K_{24}\right)\left(Q_{il}\partial_{j}Q_{jl}\right)^{2}$
(2)
$\displaystyle+\frac{1}{2}\left(K_{22}-K_{24}\right)\left(\epsilon_{ijk}Q_{il}\partial_{j}Q_{kl}\right)^{2}$
$\displaystyle+\frac{1}{2}K_{33}\left(Q_{il}\partial_{i}Q_{kl}\right)^{2}+\frac{1}{2}K_{24}Tr\left(\bm{\Delta}^{2}\right)$
$\displaystyle-\Delta\epsilon_{0}\epsilon_{a}E_{i}Q_{ij}E_{j}-\frac{4}{3}\chi_{0}E_{i}\partial_{j}Q_{ij}$
$\displaystyle-\frac{4}{3}\chi_{+}E_{i}\partial_{j}\left(Q_{ik}Q_{jk}\right)-\frac{1}{3}\chi_{2}E_{k}\partial_{k}\left(Q_{ij}Q_{ij}\right)$
$\displaystyle-\frac{4}{9}\chi_{-}E_{i}\left(Q_{ik}\partial_{j}Q_{jk}-Q_{jk}\partial_{j}Q_{ik}\right)\biggr{]}.$
(3)

The first three terms, which consist of an expansion of the free energy in
powers of $Q$, determine the magnitude of $\xi\sim\sqrt{cK/b^{2}}$ in the
bulk, where gradients are small. The next four terms serve to quantify the
Frank elastic free energy corresponding to splay, twist, bend and biaxial
splay deformations of the material, respectively. In general, the full Oseen-
Frank free energy density is written as selinger2018

$\displaystyle f_{O-F}=$
$\displaystyle\frac{1}{2}K_{11}S^{2}+\frac{1}{2}K_{22}T^{2}+\frac{1}{2}K_{33}\bm{B}^{2}$
$\displaystyle-\frac{1}{2}K_{24}\left[\frac{1}{2}S^{2}+\frac{1}{2}T^{2}-Tr\left(\bm{\Delta}^{2}\right)\right],$
(4)

where the last term is the saddle-splay term. Note that, in the literature,
one finds several variations in the notation for the saddle-splay term
selinger2018 . In particular, instead of $K_{24}$, the corresponding
coefficient is sometimes written as $\frac{1}{2}K_{24}$ or as
$\left(K_{22}+K_{24}\right)$. Those variations are not important for the
arguments that follow. The often used approximation of equal elastic constants
corresponds to $K_{11}=K_{22}=K_{33}=2K_{24}\equiv K$.

The subsequent terms in the free energy expression of Equation (3) correspond
to the potential generated by the alternating current (AC) field. The first of
these terms is the dielectric energy contribution, in which $\epsilon_{0}$ is
the dielectric permittivity of vacuum and
$\epsilon_{a}=\epsilon_{\parallel}-\epsilon_{\perp}$ is the permittivity
anisotropy of the nematic material; $\epsilon_{\parallel}$ and
$\epsilon_{\perp}$ are the dielectric permittivity in the directions
perpendicular and parallel to the director, respectively. In materials with
$\epsilon_{a}>0$, the liquid crystal aligns with the electric field, whereas
in the negative case the alignment is transverse. A negative dielectric
permittivity is required in our model, so that the director field favors
alignment within the $x-y$ plane (the plane perpendicular to the electric
field), as opposed to the $z-$axis (the direction of the electric field). The
system experiences a sinusoidal electric field with frequency $\omega$ and
maximum voltage $E_{0}$ along the $z-$axis, $\bm{E}=E_{0}sin\left(\omega
t\right)\hat{\bm{z}}$.

The terms with the $\chi$ coefficients represent the contributions of
flexoelectricity Blow2013 . The $\chi_{0}$ and $\chi_{+}$ terms depend on
spatial variations of the nematic order parameter and the electric field. Both
terms are negligible due to the absence of spatial differences. The nematic
order parameter, denoted by $S$, decreases at the location of topological
defects or solitary waves. These regions are small (by a factor of about
$10^{3}$) when compared to the entire sample size. The nematic order parameter
can therefore be viewed as approximately uniform. The $\chi_{2}$ term vanishes
when we calculate the relaxational dynamics by taking a partial derivative of
the free energy. The $\chi_{-}$ term is a crucial component; the penalties for
splay and bend deformations of the director field help amplify any
irregularities created by surface imperfections, dust particles or
inhomogeneities, which give rise to the nucleation of solitary waves. We note
here that previous studies have also taken into account the $e_{11}$ and
$e_{33}$ terms representing the effects of flexoelectricity in the director
field lavrentovich2020 ; alexe1993flexoelectric , where the polarization is
given by
$\bm{P}=e_{11}\left(\hat{\mathbf{n}}\nabla\cdot\hat{\mathbf{n}}\right)-e_{33}\left(\nabla\times\hat{\mathbf{n}}\right)\times\hat{\mathbf{n}}$.
Here $e_{11}=2\chi_{0}+S\chi_{+}$, $e_{33}=S\chi_{-}$, where $S$ is the
nematic order parameter.

To model the time evolution, we use the following equations for dynamic
relaxation

$\displaystyle\frac{\partial Q_{ij}(\bm{r},t)}{\partial
t}=-\frac{1}{\gamma_{1}}\frac{\delta F}{\delta Q_{ij}(\bm{r},t)},$ (5)

where $\gamma_{1}$ is the rotational viscosity. Hydrodynamic effects are not
included in this first attempt to model nematic solitons.

##  III Formation and propagation of solitons

To facilitate comparison of the predictions of our model to experiments we
introduce several dimensionless quantities. In simulations, a characteristic
length is defined by $\xi\sim\sqrt{cK/b^{2}}$, based on our free energy (3). A
characteristic time $\tau$ is defined by the average time required by a bullet
to propagate by one unit length, which amounts to $10^{6}s$ steps. In
experiments, the characteristic length is the size of a bullet, approximately
$20\sim 50\mu m$. The characteristic time $\tau$ is the response time required
by the nematic to return to equilibrium upon cessation of the electric field,
$10^{-2}s$; that time is commensurate with the period of the AC field in our
experiments. The electric field and frequency can be rewritten as
$\widetilde{E}=\frac{\xi\chi_{-}}{K}E$ and $\widetilde{\omega}=\tau\omega$.

Both experiments and simulations reveal the existence of two types of solitary
waves lavrentovich2020 . The first type consists of ”butterfly” structures
that, as already pointed out, tend to form around surface irregularities at
low-to-intermediate applied fields and frequencies. The second consists of
”bullet” structures, which are elongated entities that move rapidly and are
observed at higher frequencies and voltages. Both structures bend along the
direction of electric field ($z-$axis), and show an asymmetry in the direction
of propagation of the solitons ($y-$direction).

Figure 1: Director field induced by a DC field. The initial state, shown in
(a), is completely uniform. (b) A small patch on the top surface is used to
create a local tilt in the director in that region. (c) and (d) show the
director patterns that arise above the surface patch as the direction of the
electric field is alternated.

The results of our simulations confirm that several conditions are necessary
for solitons to emerge in a nematic liquid crystal. These are (1) a negative
dielectric permittivity, (2) flexoelectricity, (3) a surface irregularity, (4)
distinct Frank elastic constants for different deformation modes, and, (5) a
periodic electric field. The dielectric permittivity, as mentioned the in last
section, must be negative in order for the liquid crystal to align in the
direction perpendicular to the electric field. At $45^{\circ}C$, the elastic
constants for splay, bend and twist deformations of $CCN-47$ are $K_{11}=7pN$,
$K_{22}=1pN$ and $K_{33}=5pN$ based on experimental measurements
das2022programming ; lucchetti2021optical . The $CCN-47$ liquid crystal offers
a large, negative dielectric permittivity $(-6)$ at $45^{\circ}C$ Ref.
Dhara2008 . The dominant term of the flexoelectricity, $\chi_{-}$, is taken as
$10^{-11}C/m$ Li2018 ; de1993physics . FIG. 1 shows the director pattern in
the center of the sample ($x-y$ plane at $z=0$ and $x-z$ plane at $y=0$); one
can appreciate the combined effects of flexoelectricity and surface
irregularities. In the figure, we use rods to represent the 3D character of
the director field. The color is light yellow when the director field aligns
along the x-y plane, and it is dark when the director field aligns along the
z-axis. The structures form in a short period of time after turning on the
electric field. In FIG. 1a, the system is initially uniform. After applying a
constant electric field, $\bm{E}=E_{0}\hat{\bm{z}}$, the flexoelectricity acts
to splay and bend the director field, and an asymmetry in the $z-$axis appears
in the lower region of FIG. 1a. FIG. 1 b shows that a small random tilt is
induced in the center of the surface, which breaks the symmetry along the
$y-$axis and generates a periodic pattern similar to that observed in
experiments at low-frequency. FIG. 1 c and d show the same system as in b, but
with the direction of the electric field reversed ($\bm{E}=E_{0}\hat{\bm{z}}$
and $\bm{E}=-E_{0}\hat{\bm{z}}$ respectively). The $x-z$ plane in FIGS. 1 a
and c are similar because the electric field is the same, but the structures
are quite different when compared to c and d. In these two cases, a different
direction of bending along the $z$-axis can be observed. One can appreciate
that the patterns in the x-y plane are splayed and bent towards the opposite
direction, but they are not complete mirror images of each other. A bias in
the configuration appears, which is caused by the difference of the Frank
constants. Under a one-constant assumption, that asymmetry would not arise.

Figure 2: Director field induced by an electric field as it reverses
direction. (a) Simulated system and (b) applied electric field. (c) Director
field within the dashed area in (a) on the central x-y plane (z = 0) at
different times.

FIG. 2 shows the director pattern when the electric field is reversed. In this
figure, the directors adopt a structure similar to that of a butterfly-like
soliton as a result of the surface pattern and the flexoelectricity. This
structure shrinks when the field is reduced, and it increases again when the
electric field is raised again. The soliton propagates along the direction
perpendicular to the alignment of the director field and the electric field.
The flexoelectricity bends and splays the director field and breaks the
symmetry in z, but the director field is mostly aligned along the plane
perpendicular to the electric field ($x-y$ plane) due to the negative
dielectric permittivity. The surface irregularity breaks the symmetry along y,
and creates a butterfly-like pattern around the patch. This structure can
expand into a periodic pattern of the director field at locations far from the
surface patch due to the effects of flexoelectric polarization. A slightly
different pattern can be seen if the direction of the electric field is
reversed because of the different Frank constants. As a result, an asymmetric
structure occurs under the influence of the AC field.

Figure 3: (a) State diagram showing the system’s morphology as a function of
voltage and frequency for a nematic liquid crystal under an applied AC field.
(b)$\sim$(g) Starting from a uniform state (b), as the voltage is increased
one observes the formation of ”butterflies” in (c), followed by the formation
of bullets emitted from the butterflies in (d). At higher voltages one enters
a chaotic regime in (f). (g) represents the R-state that arises at low
frequencies.

In simulations we consider a channel with periodic boundaries in the x and y
directions. The system size is $300\times 200\times 10$. A butterfly-like
irregularity is fixed on the top surface, and the director field everywhere
else is initially aligned along the $x-$axis. A sinusoidal electric field
$\bm{E}=E_{0}sin\left(\omega t\right)\hat{\bm{z}}$ is applied, with voltage
$\widetilde{E_{0}}=0\sim 100$, and frequency $\widetilde{\omega}=0\sim 800$.
Different states can be observed for different values and frequency of the
electric field. One can construct a phase diagram as a function of the
electric field (FIG. 3). Different states arise for different magnitudes of
the frequency and field. We use red lines to represent the director field in a
plane; the black color means the director field is uniform, and yellow and
green represent different orientations when the director field is aligned
along the z-axis lavrentovich2020 . Section $\left(I\right)$ and
$\left(V\right)$ correspond to cases with a small electric strength and a
small frequency, respectively. In Section $\left(I\right)$, the field is too
small to perturb the director field, which remains almost uniform. In Section
$\left(V\right)$, the electric field is high enough to alter the director
field, but the oscillation time of the field is longer than the response time
of the material and a periodic pattern appears. This pattern is referred to as
the R-state Aya2020  and it is similar to the structure observed under a DC
field. Section $\left(IV\right)$ shows that under high voltage, solitary waves
are rarely observed, and the pattern is chaotic. In Section $\left(II\right)$
we observe butterfly-like structures, which upon further strengthening of the
field start to emit bullets in the orange region between Section
$\left(II\right)$ and Section $\left(III\right)$. When the voltage is
increased even more, in Section $\left(III\right)$, a state with more bullets
and stripes but fewer butterflies can be observed.

Figure 4: (a) Free energy of the system in the state corresponding to bullets
emitted from a stable butterfly (orange band in FIG. 3). The inset shows a
zoomed-in portion of the free energy corresponding to the emission of an
individual bullet. The panels on the right show director field configurations
through cross polarizers corresponding to different times in the figure shown
in the inset. (b) Results of simulations for the threshold voltage for bullet
emission from a butterfly for different sample thicknesses. (c) Experimental
results for threshold voltage for different sample thicknesses.

It is of interest to examine the contributions of elastic energy $F_{e}$, the
free energy of dielectric polarization $F_{di}$, and the flexoelectric
polarization (equation (3)), as bullets are emitted from the wings of a stable
butterfly. We show this process for $\widetilde{E_{0}}=38$ and
$\widetilde{\omega}=400$ in FIG. 4(a). These energies are set to zero in the
uniform state (the director field is completely aligned along the $x-$axis).
The elastic energy $F_{e}$ includes a combination of splay, bend, and twist
contributions. Since the electric field is along the $z-$axis, the free energy
of dielectric polarization $F_{di}$ should increase when more regions of the
director field are aligned along the $z-$axis. From the discussion of the
flexoelectric free energy $F_{flex}$ in Section $\left(II\right)$, it is
apparent that splay and bend deformations are key contributors. In FIG. 4(a),
we show the data at the time when $E_{z}=E_{0}$, so that $F_{flex}$ is always
positive. In the circle and the zoomed in panel in the figure, we focus on the
peak of the free energy. The energies are small before $t=3.4$, where no
bullet is emitted from the butterfly. They then increase ($t=3.4$ to $t=3.7$)
as a new bullet is generated from the upper branch of the butterfly. The peak
reaches a maximum ($t=3.5$) when the bullet is about to be emitted from the
butterfly, and it decreases when the bullet moves away. There are multiply
peaks in the figure; each peak represents a new-born bullet emitted from a
stationary butterfly.

We also consider the effect of anchoring on soliton formation. FIG. 4(b) and
(c) show the state in which the butterfly structure starts to emit bullets for
channels of different thickness. We include results from simulations and
experiments. Both sets of results show that, in a narrower channel, where
anchoring has a stronger influence, a higher electric field is needed to emit
bullets. We find that a linear relationship exists between the frequency and
the electric field strength. When the voltage is large, the system can
rearrange its structure in a shorter time interval. If the frequency is low,
the material has sufficient time to undergo a complete structure change, and
the solitary structure is broken. As a result, a higher voltage requires that
a higher frequency be used if we want to generate a solitary structure. The
inverse relationship between thickness and voltage is only observed over a
small range of thickness. For higher thickness, the effect of the confining
surfaces is reduced, and solitons are no longer observed regardless of the
applied voltage is. This is also seen in experiments Li2018 .

Figure 5: Comparison of results from simulations (left) and experiments
(right) corresponding to the director field as observed through cross
polarizers for nematic liquid crystal under different applied fields. In (a),
several butterflies are observed. In (b), bullets are emitted from the
butterfly. In (c), rapidly moving bullets are observed.

FIG. 5 focuses on the three most relevant regions of the phase diagram. We
show results from simulations and experiments. When $\widetilde{E_{0}}=36$ and
$\widetilde{\omega}=400$ (FIG. 5 (a)), stable butterflies are observed and
they move at very low speed. When the electric field reaches
$\widetilde{E_{0}}=39$ at the same frequency (FIG. 5 (c)), the butterflies
shrink, and the system consists mostly of stripes and bullets. For electric
fields between these two cases, e.g. $\widetilde{E_{0}}=38$ (FIG. 5 (b)), we
can observe that the bullets emerge from the butterfly created by the
irregularity on the surface. The bullets may emanate from any wing of the
butterfly. The bullets travel over a long distance in a direction
perpendicular to the director. In experiments, the same structures are
observed for the same frequency but at slightly higher voltages
das2022programming . We attribute the discrepancy to the simplicity of our
model and the fact that hydrodynamic effects have not been included in our
model and calculations.

##  IV Conclusions

In this work we have investigated the formation and propagation of solitons in
a 3D nematic liquid crystal confined between to planar surfaces with planar
anchoring. A minimal model of the nematic LC has been used to show that
solitons are nucleated at a surface inhomogeneity with perpendicular
anchoring, and are generated upon exposure to an AC field. A negative
dielectric permittivity and flexoelectricity have been shown to be critical
elements of the model required for soliton generation. The AC field causes the
director field to adopt distinct states periodically and, for sufficiently
high frequencies, the material is unable to relax to its stationary structure,
thereby emitting solitons that travel at large speeds throughout the system.

The model has been used to construct a state diagram as a function of the
strength of the applied field and its frequency. Several regimes are predicted
in simulations, including a uniform state, a stripe state, a chaotic state,
and a soliton state, which only occurs over a narrow range of electric field
strength. Good agreement is obtained between the results of simulations and
our own experimental observations, serving to establish the validity of the
model. Importantly, the findings presented in this work provide foundational
knowledge with which to design new systems for the controlled production of
solitons, thereby opening new opportunities for their use for applications in
microfluidic transport, optical, and sensing technologies.

## References

  * (1) T. Dauxois and M. Peyrard, Physics of solitons.  Cambridge University Press, 2006.
  * (2) N. Manton and P. Sutcliffe, Topological solitons.  Cambridge University Press, 2004.
  * (3) T. Heimburg and A. D. Jackson, “On soliton propagation in biomembranes and nerves,” Proceedings of the National Academy of Sciences, vol. 102, no. 28, pp. 9790–9795, 2005.
  * (4) R. Appali, U. Van Rienen, and T. Heimburg, “A comparison of the hodgkin–huxley model and the soliton theory for the action potential in nerves,” in Advances in Planar Lipid Bilayers and Liposomes, vol. 16, pp. 275–299, Elsevier, 2012.
  * (5) S. S. Andersen, A. D. Jackson, and T. Heimburg, “Towards a thermodynamic theory of nerve pulse propagation,” Progress in neurobiology, vol. 88, no. 2, pp. 104–113, 2009.
  * (6) Y. V. Kartashov, G. E. Astrakharchik, B. A. Malomed, and L. Torner, “Frontiers in multidimensional self-trapping of nonlinear fields and matter,” Nature Reviews Physics, vol. 1, no. 3, pp. 185–197, 2019.
  * (7) L. Lam, “Solitons and field induced solitons in liquid crystals,” in Solitons in Liquid Crystals, pp. 9–50, Springer, 1992.
  * (8) B.-X. Li, V. Borshch, R.-L. Xiao, S. Paladugu, T. Turiv, S. V. Shiyanovskii, and O. D. Lavrentovich, “Electrically driven three-dimensional solitary waves as director bullets in nematic liquid crystals,” Nature communications, vol. 9, no. 1, pp. 1–10, 2018.
  * (9) S. Aya and F. Araoka, “Kinetics of motile solitons in nematic liquid crystals,” Nature communications, vol. 11, no. 1, pp. 1–10, 2020.
  * (10) S. Das, S. Roh, N. Atzin, A. Mozaffari, X. Tang, J. J. de Pablo, and N. L. Abbott, “Programming solitons in liquid crystals using surface chemistry,” Langmuir, vol. 38, no. 11, pp. 3575–3584, 2022.
  * (11) A. Zawadzki and H. Walton, “Measurements of the splay and bend elastic constants of 4’butyl-4-heptyl-bicyclohexyl-4-carbononitrile, ccn47,” Molecular Crystals and Liquid Crystals, vol. 569, no. 1, pp. 10–14, 2012.
  * (12) L. Lucchetti, G. Nava, R. Barboza, F. Ciciulla, and T. Bellini, “Optical force-based detection of splay and twist viscoelasticity of ccn47 across the nematic-to-smectic a transition,” Journal of Molecular Liquids, vol. 329, p. 115520, 2021.
  * (13) X. Tang and J. V. Selinger, “Orientation of topological defects in 2d nematic liquid crystals,” Soft Matter, vol. 13, no. 32, pp. 5481–5490, 2017.
  * (14) X. Tang and J. V. Selinger, “Theory of defect motion in 2d passive and active nematic liquid crystals,” Soft matter, vol. 15, no. 4, pp. 587–601, 2019\.
  * (15) J. V. Selinger, “Interpretation of saddle-splay and the oseen-frank free energy in liquid crystals,” Liquid Crystals Reviews, vol. 6, no. 2, pp. 129–142, 2018.
  * (16) M. L. Blow and M. T. da Gama, “Interfacial motion in flexo-and order-electric switching between nematic filled states,” Journal of Physics: Condensed Matter, vol. 25, no. 24, p. 245103, 2013.
  * (17) O. D. Lavrentovich, “Design of nematic liquid crystals to control microscale dynamics,” Liquid Crystals Reviews, vol. 8, no. 2, pp. 59–129, 2020.
  * (18) A.-L. Alexe-Ionescu, “Flexoelectric polarization and second order elasticity for nematic liquid crystals,” Physics Letters A, vol. 180, no. 6, pp. 456–460, 1993.
  * (19) L. Lucchetti, G. Nava, R. Barboza, F. Ciciulla, and T. Bellini, “Optical force-based detection of splay and twist viscoelasticity of ccn47 across the nematic-to-smectic a transition,” Journal of Molecular Liquids, vol. 329, p. 115520, 2021.
  * (20) S. Dhara and N. Madhusudana, “Physical characterisation of 4’butyl-4-heptyl-bicyclohexyl-4-carbonitrile,” Phase Transitions, vol. 81, no. 6, pp. 561–569, 2008.
  * (21) P.-G. De Gennes and J. Prost, The physics of liquid crystals.  No. 83, Oxford university press, 1993.