# Interaction-aware Traffic Prediction and Scenario-based Model Predictive
Control for Autonomous Vehicles on Highways*

Xiaorong Zhang1, Sahar Zeinali1, and Georg Schildbach1 *This work was
supported by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) – project number 460891204.1Institute for electrical engineering
in Medicine,University of luebeck, Luebeck, Germany {xiaorong.zhang,
sahar.zeinali<EMAIL_ADDRESS>
###### Abstract

This paper addresses the problem of traffic prediction and control of
autonomous vehicles on highways. A modified Interacting Multiple Model Kalman
filter algorithm is applied to predict the motion behavior of the traffic
participants by considering their interactions. A scenario generation
component is used to produce plausible scenarios of the vehicles based on the
predicted information. A novel integrated decision-making and control system
is proposed by applying a Scenario-based Model Predictive Control approach.
The designed controller considers safety, driving comfort, and traffic rules.
The recursive feasibility of the controller is guaranteed under the inclusion
of the ‘worst case’ as an additional scenario to obtain safe inputs. Finally,
the proposed scheme is evaluated using the HighD dataset. Simulation results
indicate that the vehicle performs safe maneuvers in different traffic
situations under the designed control framework.

##  I INTRODUCTION

###  I-A Motivation

During the past decades, the design of control systems for autonomous vehicles
on highways has been extensively studied. The primary purpose of these systems
is to safely control the ego vehicle (EV) by utilizing the predicted motion
states of the surrounding target vehicles (TVs) [1]. The predicted states are
usually uncertain, so generating safe, comfortable, energy-efficient, and
real-time capable control strategies is challenging.

###  I-B Literature Review

Physics-based [2], maneuver-based [3], and interaction-aware motion models [4]
are used for the state prediction of vehicles. Since the interconnections
between traffic participants are considered in the interaction-aware models,
they are a good choice for describing realistic scenarios. Specifically, the
mutual influence of the vehicles is usually expressed from a finite set of
trajectory clusters or Dynamic Bayesian Networks (DBNs) [5][6]. Moreover, a
novel interaction-aware traffic model is proposed by Lefkopoulos et al. [7],
combining the physics of the vehicles, the intention of the drivers, and a no-
collision assumption using an Interacting Multiple Model Kalman filter (IMM-
KF), which establishes a new scheme with improved computational efficiency.

Model Predictive Control (MPC) has been widely applied in designing a planner
and controller for the EV by considering several constraints, such as the
traffic rules, safety, and the comfort of driving [8, 9, 10]. The underlying
reason for the prevalence of MPC is its ability to handle explicit constraints
in an optimization problem with a moving horizon [11]. As a significant
variant of Stochastic MPC, Scenario-based MPC (SCMPC) generates corresponding
constraints in terms of the possible situational context[12]. It has been
successfully implemented under several highway traffic conditions [13] [14],
as it is easily compatible with the traffic prediction component and can
handle uncertainty using the information, contained in a few scenarios.

Safety is the most critical aspect of controlling the EV. This feature becomes
more challenging in emergency scenarios, e.g., unexpected deceleration of the
leading vehicle (LV) or a sudden cut-in of a TV. These circumstances are
identified as a safety-critical-event (SCE) [15] [16], where the EV applies
immediate braking against the crash, and it may lead to a deceleration of the
EV until a standstill. Adaptive Cruise Control (ACC) is a helpful tool for
dealing with SCEs, where the EV makes decisions based on the information about
the LV. However, a large time headway in this algorithm may lead to over-
conservative actions [17]. A safety controller is proposed based on ACC in
[18], where the EV uses a predefined deceleration profile when the LV performs
a maneuver in an SCE manner. Another representative solution to SCEs is the
rigorous formalizing mathematical model of Responsibility-Sensitive-Safety
(RSS)[19]. In this model, a safety distance is defined by assuming a ‘worst
case’ scenario, and the EV responds to an SCE by decelerating at a predefined
rate and without applying a full braking force. This approach might be
sensitive to the parameter design, and the subtle change of the parameter set
might lead to a totally different decision strategy[20].

###  I-C Contribution

This paper presents a method for controlling the EV on highways with a
comprehensive but simple structure, combining the interaction-aware motion
prediction of IMM-KF with the decision-making and control design based on
SCMPC. Furthermore, the safety of the EV under the control inputs is
guaranteed by theoretically proving the recursive feasibility of the SCMPC
under the consideration of the ‘worst case’ scenario. To the best of the
authors’ knowledge, no work handles all the mentioned problems with an
integrated SCMPC architecture and proves the feasibility of the algorithm at
the same time. The performance of the proposed approach is evaluated for
different HighD dataset scenarios, indicating that EV performs safe and
desirable maneuvers by applying the designed control architecture.

##  II CONTROL ARCHITECTURE

The proposed control architecture is shown in Fig. 1 and works as follows.
First, the mode states of the TVs are predicted based on the IMM-KF. This
information is sent to a scenario generation component to produce all possible
maneuvers of the TVs. Then, the most likely scenarios are filtered based on a
predefined probability threshold. Finally, an SCMPC-based control system
consisting of two control modes is established, which corresponds to
’following the current lane’ or ’changing the lane’. In addition to the
filtered scenarios, the ‘worst case’ scenario is also included in the design
of both controllers. A decision-making module chooses the desired maneuver of
the EV based on the cost function value of the two control modes.

Figure 1: Schematic diagram of the proposed control structure

##  III SCENARIO GENERATION

###  III-A Intention-based Policy Mode

In this subsection, the longitudinal and lateral policy modes of the vehicles
are described. In the longitudinal direction, the ’velocity tracking’ (VT) and
’distance keeping’ (DK) modes [7] are utilized, where the EV tracks a
reference velocity in the VT mode and keeps a safe distance from its LV in the
DK mode, respectively. In the lateral direction, three modes corresponding to
’lane 1’, ’lane 2’, and ’lane 3’ are applied to represent the target lane of
the vehicles, as shown in Fig. 2. For example, if EV is currently at lane 2,
it can change to the right lane, lane 1, or the left lane, lane 3, according
to the associated lateral mode. Thus, the total number of modes is 6, $M=6$.

Figure 2: Illustration of lanes and vehicles

The common state vector $x_{k}$ in all modes at time $k$ is

$x_{k}\triangleq\begin{bmatrix}\underbrace{\begin{matrix}p_{\text{lon},k}&v_{\text{lon},k}&a_{\text{lon},k}\end{matrix}}_{x_{\text{lon},k}}&\underbrace{\begin{matrix}p_{\text{lat},k}&v_{\text{lat},k}&a_{\text{lat},k}\end{matrix}}_{x_{\text{lat},k}}\end{bmatrix}^{\top},$
(1)

here $p_{*,k},v_{*,k},a_{*,k}$ are respectively the position, velocity, and
acceleration in the corresponding direction,
$*\in\\{\text{lon},\text{lat}\\}$. The unknown reference velocity
$v_{\text{ref},k}$ (VT mode) or the reference time gap $t_{\text{gap},k}$ (DK
mode) is also included in the longitudinal policy mode to be estimated.
Therefore, the full state vector $z_{k}$ in each policy mode at time $k$ is

$z_{k}\triangleq\begin{bmatrix}\underbrace{\begin{matrix}x_{\text{lon},k}&r_{\text{ref},k}\end{matrix}}_{z_{\text{lon},k}}&\underbrace{\begin{matrix}x_{\text{lat},k}\end{matrix}}_{z_{\text{lat},k}}\end{bmatrix}^{\top},$
(2)

where $r_{\text{ref},k}$ indicates the unknown reference velocity
$v_{\text{ref},k}$ or reference time gap $t_{\text{gap},k}$. We use $x_{k}$ or
$z_{k}$ to denote the common or full states of an arbitrary vehicle, and
$x_{k}^{(\circledast)}$ or $z_{k}^{(\circledast)}$ to clarify the common or
full state of a specific vehicle throughout the paper, where $\circledast\in$
{EV, TV, LV}. To make the policies closer to the driver’s real intention, two
linear quadratic regulator (LQR)-based feedback controllers are included in
longitudinal and lateral policy modes, the associated control gains are
defined as follows:

$\displaystyle K_{\text{lon}}^{(\diamond)}$
$\displaystyle=\begin{bmatrix}K_{\text{lon},1}^{(\diamond)}&K_{\text{lon},2}^{(\diamond)}&K_{\text{lon},3}^{(\diamond)}\end{bmatrix}^{\top},$
(3a) $\displaystyle K_{\text{lat}}$
$\displaystyle=\begin{bmatrix}K_{\text{lat},1}&K_{\text{lat},2}&K_{\text{lat},3}\end{bmatrix}^{\top}.$
(3b)

where $\diamond\in\\{\text{VT},\text{DK}\\}$. Denoting the sampling time as
$T$, the discrete form of each policy mode at time $k$ is

$\displaystyle\underbrace{\begin{bmatrix}z_{\text{lon},k+1}\\\
z_{\text{lat},k+1}\end{bmatrix}}_{z_{k+1}}$
$\displaystyle=\underbrace{\begin{bmatrix}F_{\text{lon},k}^{(\diamond)}&\mathbf{0}_{4\times
3}\\\ \mathbf{0}_{3\times
4}&F_{\text{lat},k}^{(\lambda)}\end{bmatrix}}_{F_{k}}\underbrace{\begin{bmatrix}z_{\text{lon},k}\\\
z_{\text{lat},k}\end{bmatrix}}_{z_{k}}+\underbrace{\begin{bmatrix}E_{\text{lon},k}^{(\diamond)}\\\
E_{\text{lat},k}^{(\lambda)}\end{bmatrix}}_{E_{k}}+\omega_{k},$ (4a)
$\displaystyle y_{k}$ $\displaystyle=\bm{I}_{7\times 7}z_{k}+\nu_{k},$ (4b)

where $y_{k}\in\mathbb{R}^{7}$ and $\bm{I}_{7\times 7}$ are the measurement
vector and observation matrix. $\omega_{k}$ and $\nu_{k}$ are the process
noise and measurement noise, which are assumed to be normally distributed,
with the covariance $Q$ and $R$, respectively. The system matrices for
longitudinal direction are

$F_{\text{lon},k}^{\text{(VT)}}=\begin{bmatrix}1&T&\frac{T^{2}}{2}&0\\\
0&1-\frac{K_{\text{lon},2}^{\text{(VT)}}T^{2}}{2}&T-\frac{K_{\text{lon},3}^{\text{(VT)}}T^{2}}{2}&\frac{K_{\text{lon},2}^{\text{(VT)}}T^{2}}{2}\\\
0&-K_{\text{lon},2}^{\text{(VT)}}T&1-K_{\text{lon},3}^{\text{(VT)}}T&K_{\text{lon},2}^{\text{(VT)}}T\\\
0&0&0&1\end{bmatrix},$ $\begin{split}&F_{\text{lon},k}^{\text{(DK)}}=\\\
&\scalebox{1.0}{$\begin{bmatrix}1-\frac{K_{\text{lon},1}^{\text{(DK)}}{T}^{3}}{6}&-\frac{K_{\text{lon},1}^{\text{(DK)}}{T}^{2}}{2}&-K_{\text{lon},1}^{\text{(DK)}}T&0\\\
T-\frac{K_{\text{lon},2}^{\text{(DK)}}{T}^{3}}{6}&1-\frac{K_{\text{lon},2}^{\text{(DK)}}{T}^{2}}{2}&-K_{\text{lon},2}^{\text{(DK)}}T&0\\\
\frac{{T}^{2}}{2}-\frac{K_{\text{lon},3}^{\text{(DK)}}{T}^{3}}{6}&T-\frac{K_{\text{lon},3}^{\text{(DK)}}{T}^{2}}{2}&1-K_{\text{lon},3}^{\text{(DK)}}T&0\\\
-\frac{K_{\text{lon},1}^{\text{(DK)}}v_{\text{lead},k}{T}^{3}}{6}&-\frac{K_{\text{lon},1}^{\text{(DK)}}v_{\text{lead},k}{T}^{2}}{2}&-K_{\text{lon},2}^{\text{(DK)}}v_{\text{lead},k}T&1\end{bmatrix}^{\top}$},\end{split}$

and the LQR-related input matrices are

$E_{\text{lon},k}^{\text{(VT)}}=\mathbf{0}_{4\times 1},$
$\begin{split}&\scalebox{1.0}{$E_{\text{lon},k}^{\text{(DK)}}$}=\\\
&\scalebox{1.0}{$\begin{bmatrix}\frac{K_{\text{lon},1}^{\text{(DK)}}{T}^{3}}{6}-1&\frac{K_{\text{lon},2}^{\text{(DK)}}{T}^{3}}{6}-T&\frac{K_{\text{lon},3}^{\text{(DK)}}{T}^{3}}{6}-\frac{T^{2}}{2}&0\\\
\frac{K_{\text{lon},1}^{\text{(DK)}}{T}^{2}}{2}&\frac{K_{\text{lon},2}^{\text{(DK)}}{T}^{2}}{2}-1&\frac{K_{\text{lon},3}^{\text{(DK)}}{T}^{2}}{2}-T&0\\\
K_{\text{lon},1}^{\text{(DK)}}T&K_{\text{lon},2}^{\text{(DK)}}T&K_{\text{lon},3}^{\text{(DK)}}T-1&0\\\
0&0&0&0\end{bmatrix}\begin{bmatrix}p_{\text{lon},k}^{(\text{LV})}\\\
v_{\text{lon},k}^{(\text{LV})}\\\ a_{\text{lon},k}^{(\text{LV})}\\\
0\end{bmatrix}$}.\end{split}$

Consider the position $p^{(\lambda)}$ of the center line of the target lane
$\lambda$, the system matrix and input matrix for the lateral direction and
for the target lane $\lambda$ ($\lambda=1,2,3$) are as follows:

$\displaystyle F_{\text{lat},k}^{(\lambda)}$
$\displaystyle=\begin{bmatrix}-\frac{K_{\text{lat},1}{T}^{3}}{6}&-\frac{K_{\text{lat},2}{T}^{3}}{6}&-\frac{K_{\text{lat},3}{T}^{3}}{6}\\\
-\frac{K_{\text{lat},1}{T}^{2}}{2}&-\frac{K_{\text{lat},2}{T}^{2}}{2}&-\frac{K_{\text{lat},3}{T}^{2}}{2}\\\
-K_{\text{lat},1}T&-K_{\text{lat},2}T&-K_{\text{lat},3}T-1\end{bmatrix},$ (5a)
$\displaystyle E_{\text{lat},k}^{(\lambda)}$
$\displaystyle=\begin{bmatrix}\frac{K_{\text{lat},1}{T}^{3}}{6}p^{(\lambda)}\\\
\frac{K_{\text{lat},1}{T}^{2}}{2}p^{(\lambda)}\\\
K_{\text{lat},1}Tp^{(\lambda)}\end{bmatrix}.$ (5b)

###  III-B Interaction-aware Estimation and Prediction

Since the motion behavior of each vehicle is influenced by its surrounding
vehicles, the state estimation and prediction of each vehicle are calculated
in descending priority order by considering the interaction between vehicles,
as proposed in [7]. In particular, the priority criteria are given by

  * (i)

If two vehicles are in the same lane, the preceding vehicle has higher
priority.

  * (ii)

If two vehicles are in different lanes, the vehicle with the higher
longitudinal progress over a specific horizon has higher priority.

For the sake of clarity, a policy mode corresponds to one model of the IMM-KF
in this work. In the IMM-KF [7], we consider the Markov jump linear system
(4), where the transition probability from mode $i$ to mode $j$ is denoted as
$\pi^{(i|j)}$, and $\pi^{(i|j)}\in[0,1]$, $i,j\in\\{1,2,...,6\\}$. Since the
reference parameter of the VT and DK mode is different, we mix individual
common estimates and initialize each mode in the first step as:

$\displaystyle c^{(i)}$
$\displaystyle=\sum_{j=1}^{M}\pi^{(j|i)}\mu_{k-1}^{(j)},$ (6a)
$\displaystyle\mu_{k-1}^{(j|i)-}$
$\displaystyle=\frac{\pi^{(j|i)}\mu_{k-1}^{(j)}}{c^{(i)}},$ (6b)
$\displaystyle\bar{x}_{k-1}^{(i)-}$
$\displaystyle=\sum_{j=1}^{M}\mu_{k-1}^{(j|i)-}\hat{x}_{k-1}^{(j)-},$ (6c)
$\displaystyle\begin{split}\bar{P}_{k-1}^{(i)-}&=\sum_{j=1}^{M}\mu_{k-1}^{(j|i)-}\left[{P}_{k-1}^{(j)-}\right.\\\
&\left.+\left(\bar{x}_{k-1}^{(i)-}-\hat{x}_{k-1}^{(j)-}\right)\left(\bar{x}_{k-1}^{(i)-}-\hat{x}_{k-1}^{(j)-}\right)^{\top}\right],\end{split}$
(6d)

where $\mu_{k-1}^{(j|i)-}$ is the mixing conditional mode probability,
$\hat{x}_{k-1}^{(j)-}$, ${P}_{k-1}^{(j)-}$ are the common state estimation and
its associated covariance of the mode $j$, which are part of the full state
estimation $\hat{z}_{k-1}^{(j)-}$, and its associated covariance
$\mathbb{P}_{k-1}^{(j)-}$. The fused common state estimation and its
associated covariance are $\bar{x}_{k-1}^{(i)-}$ and $\bar{P}_{k-1}^{(i)-}$.
The corresponding fused full state estimation and its associated covariance
are $\bar{z}_{k-1}^{(i)-}$ and $\bar{\mathbb{P}}_{k-1}^{(i)-}$. The second
step is to predict and update each policy mode as

$\displaystyle\hat{z}_{k-1}^{(i)+}$
$\displaystyle=F_{k-1}^{(i)}\bar{z}_{k-1}^{(i)-}+E_{k-1}^{(i)},$ (7a)
$\displaystyle\mathbb{P}_{k-1}^{(i)+}$
$\displaystyle=F_{k-1}^{(i)}\bar{\mathbb{P}}_{k-1}^{(i)-}F_{k-1}^{(i)\top}+Q_{k-1}^{(i)},$
(7b) $\displaystyle\tilde{y}_{k}^{(i)}$
$\displaystyle=y_{k}^{(i)}-\bm{I}_{7\times 7}\hat{z}_{k-1}^{(i)+},$ (7c)
$\displaystyle r_{k}^{(i)}$ $\displaystyle=\bm{I}_{7\times
7}\mathbb{P}_{k-1}^{(i)+}\bm{I}_{7\times 7}^{\top}+R_{k}^{(i)},$ (7d)
$\displaystyle L_{k}^{(i)}$
$\displaystyle=\mathbb{P}_{k-1}^{(i)+}\bm{I}_{7\times
7}^{\top}{r_{k}^{(i)}}^{-1},$ (7e) $\displaystyle\hat{z}_{k}^{(i)-}$
$\displaystyle=\hat{z}_{k-1}^{(i)+}+L_{k}^{(i)}\tilde{y}_{k}^{(i)},$ (7f)
$\displaystyle\mathbb{P}_{k}^{(i)-}$ $\displaystyle=(\bm{I}_{7\times
7}-L_{k}^{(i)}\bm{I}_{7\times 7})\mathbb{P}_{k-1}^{(i)+},$ (7g)

using the prior state estimate $\hat{z}_{k-1}^{(i)+}$ and its covariance
$\mathbb{P}_{k-1}^{(i)+}$, the innovation residual $\tilde{y}_{k}^{(i)}$ and
its covariance $r_{k}^{(i)}$, Kalman gain $L_{k}^{(i)}$, and posterior
predicted state estimate $\hat{z}_{k}^{(i)-}$ and covariance
$\mathbb{P}_{k}^{(i)-}$. Based on the estimation result of individual mode
$\hat{z}_{k}^{(i)-}$, the state prediction of each policy mode is

$\displaystyle\hat{z}_{t|k}=\phi(t,1)\hat{z}_{k}^{-}+\sum_{\delta=k+1}^{t}\phi(t,\delta)E_{\delta-1},$
(8a)
$\displaystyle\phi(t,\delta)=\begin{cases}(\Pi_{\eta=\delta}^{t-1}F_{\eta}^{\top})^{\top}&\text{if}\leavevmode\nobreak\
t>\delta\\\ \bm{I}_{7\times 7},&\text{if}\leavevmode\nobreak\
t=\delta\end{cases},$ (8b)

where $t=k+1,...,k+1+N$, and $N$ is the prediction horizon. In order to obtain
a collision-free prediction (’no-collision prediction’), a mixed integer
quadratic programming (MIQP) problem is formulated to get the modified state
estimation $\hat{z}_{k}^{(\text{proj)}-}$, where the safety constraints
between the studied vehicle and other vehicles which have higher priority are
considered. Note that the state estimation of each policy mode is still
$\hat{z}_{k}^{-}$, and only the state prediction is modified in terms of
$\hat{z}_{k}^{(\text{proj)}-}$ and (8). The state estimation error between
$\hat{z}_{k}^{-}$ and $\hat{z}_{k}^{(\text{proj)}-}$, and its covariance are
used to augment the innovation residual $\bar{y}_{k}^{{(i)}}$ and its
covariance ${{r}_{k}}^{{(i)}}$ as $\breve{y}_{k}^{{(i)}}$ and
${\tilde{r}_{k}}^{{(i)}}$. Then, the policy mode probability is updated based
on the augmented matrices

$\displaystyle\tilde{L}_{k}^{(i)}$
$\displaystyle=\frac{\text{exp}(-\frac{1}{2}\breve{y}_{k}^{{(i)}\top}{\tilde{r}_{k}}^{{(i)-1}}\breve{y}_{k}^{(i)})}{{\left\lvert
2\pi\tilde{r}_{k}^{(i)}\right\rvert}^{1/2}},$ (9a)
$\displaystyle\tilde{\mu}_{k}^{(i)}$
$\displaystyle=\frac{c^{(i)}\tilde{L}_{k}^{(i)}}{\sum_{j=1}^{M}c^{(j)}\tilde{L}_{k}^{(j)}}.$
(9b)

The final step is to mix state estimation and its covariance according to the
updated probability of the individual mode

$\displaystyle\hat{x}_{k}^{-}$
$\displaystyle=\sum_{i=1}^{M}\tilde{\mu}_{k}^{(i)}\hat{x}_{k}^{(i)-},$ (10a)
$\displaystyle P_{k}^{-}$
$\displaystyle=\sum_{i=1}^{M}\tilde{\mu}_{k}^{(i)}[P_{k}^{(i)-}+(\hat{x}_{k}^{-}-\hat{x}_{k}^{(i)-})(\hat{x}_{k}^{-}-\hat{x}_{k}^{(i)-})^{\top}].$
(10b)

The updated state estimation is also modified in terms of a MIQP problem to
guarantee safety over the whole prediction horizon. The readers are referred
to Lefkopoulos et al. [7] for more details about the mentioned method.

###  III-C Scenario Generation of TVs

A scenario is defined as a tuple of motion maneuvers for all TVs. Assume that
the number of investigated TVs is $V$, then a total of $M^{V}$ possible
scenarios can be generated. $\mu^{(n)}_{i}$ is the probability of TV $n$ with
the policy mode $i$, $i\in\\{1,2,...,6\\}$. Assuming statistical independence
of each vehicle’s no-collision prediction over the prediction horizon, then
the probability of the scenario $s$ is calculated by

$\text{Pr}(s)=\overset{V}{\underset{n=1}{\prod}}{\mu^{(n)}_{i}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ s=1,...,M^{V},$ (11)

where $\sum_{s=1}^{M^{V}}\text{Pr}{(s)}=1$. In order to have high-probability
scenarios, scenarios with a probability less than a predefined threshold
$\underline{P}$ are not considered. The probability of the remaining scenarios
is normalized by

$\overline{\text{Pr}}(s)=\frac{\text{Pr}(s)}{1-\sum_{\zeta=1}^{\theta}\text{Pr}(\zeta)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ s=1,...,M^{V}-\theta,$ (12)

where $\theta$ is the total number of scenarios with probability less than
$\underline{P}$ .

##  IV SCENARIO-BASED MODEL PREDICTIVE CONTROL

Based on the predicted scenarios of the TVs, a feasible trajectory for the EV
is calculated by solving a constrained finite-time optimal control problem
(CFTOCP) in a moving horizon fashion. The objective of the optimization
problem is to follow the planned reference trajectory with minimum effort and
with the consideration of safety constraints, traffic rules, and driving
comfort. The first computed control input of the CFTOCP is fed to the system
at each time step.

###  IV-A Vehicle Model

The Jerk Model [21] is applied to represent the dynamics of the EV:

$\underbrace{\begin{bmatrix}p_{*,k+1}^{\text{(EV)}}\\\
v_{*,k+1}^{\text{(EV)}}\\\
a_{*,k+1}^{\text{(EV)}}\end{bmatrix}}_{x_{*,k+1}^{\text{(EV)}}}=\underbrace{\begin{bmatrix}1&T_{\text{p}}&\frac{1}{2}{T_{\text{p}}}^{2}\\\
0&1&T_{\text{p}}\\\ 0&0&1\\\
\end{bmatrix}}_{A}\underbrace{\begin{bmatrix}p_{*,k}^{\text{(EV)}}\\\
v_{*,k}^{\text{(EV)}}\\\
a_{*,k}^{\text{(EV)}}\end{bmatrix}}_{x_{*,k}^{\text{(EV)}}}+\underbrace{\begin{bmatrix}\frac{1}{6}{T_{\text{p}}}^{3}\\\
\frac{1}{2}{T_{\text{p}}}^{2}\\\
T_{\text{p}}\end{bmatrix}}_{B}\underbrace{j_{*,k}^{\text{(EV)}}}_{u_{*,k}^{\text{(EV)}}},$
(13)

with the prediction time step ${T_{\text{p}}}$, the states
$x_{*,k}^{\text{(EV)}}$, the position $p_{*,k}^{\text{(EV)}}$, the velocity
$v_{*,k}^{\text{(EV)}}$, the acceleration $a_{*,k}^{\text{(EV)}}$, the control
inputs $u_{*,k}^{\text{(EV)}}$, and the jerk $j_{*,k}^{\text{(EV)}}$, where
$*\in\\{\text{lon},\text{lat}\\}$. (13) can be rewritten as follows:

$x_{k+1}^{\text{(EV)}}=\underbrace{\begin{bmatrix}A&0_{3\times 3}\\\
0_{3\times
3}&A\end{bmatrix}}_{\bar{A}}x_{k}^{\text{(EV)}}+\underbrace{\begin{bmatrix}B\\\
B\end{bmatrix}}_{\bar{B}}u_{k}^{\text{(EV)}},$ (14)

where $x_{k}^{\text{(EV)}}$ is defined in (1), and
$u_{k}^{\text{(EV)}}=\begin{bmatrix}u_{\text{lon},k}^{\text{(EV)}}&u_{\text{lat},k}^{\text{(EV)}}\end{bmatrix}^{\top}$.
Note that the prediction time step ${T_{\text{p}}}$ differs from the sampling
time step $T$ of the IMM-KF, which usually satisfies $T_{\text{p}}>T$.

###  IV-B Scenario-based Model Predictive Controller

Considering lane-keeping and lane-change as possible motion behaviors of the
EV, two control modes are proposed leading to different reference
trajectories. The first control mode aims to make the EV keep its velocity and
stay in the current lane, the second control mode aims to lead the EV to a
target lane while maintaining its speed. The decision-making module for
choosing the control input resulting from these controllers is integrated into
the controller design. The input corresponding to the minimal cost function
value is applied to the system.

In addition to the generated scenarios, a so-called ‘worst case’ scenario is
introduced to guarantee the recursive feasibility of CFTOCP. In this scenario,
the LV is assumed to be decelerating with its minimum acceleration over the
prediction horizon. We introduce two sequences of control inputs
$u_{0}^{\text{(EV)}},...,u_{N-1}^{\text{(EV)}}$ and
$\breve{u}_{0}^{\text{(EV)}},...,\breve{u}_{N-1}^{\text{(EV)}}$. The first
input sequence is calculated by avoiding collision between EV and the LV/TVs
under the generated scenarios, which is used to calculate the value of the
cost function with its associated states $x_{k}^{\text{(EV)}}$. The second
input sequence is obtained by considering the safety constraints under the
‘worst case’ scenario, its associated states are
$\breve{x}_{k}^{\text{(EV)}}$, and the terminal set of the states is
$\breve{\mathbb{X}}^{\text{(EV)}}_{f}$, which is described in the proof of
recursive feasibility. The first computed inputs $u_{0}^{\text{(EV)}}$ and
$\breve{u}_{0}^{\text{(EV)}}$ must be equal in order to guarantee the
recursive feasibility. The CFTOCP is formulated as

$\displaystyle J=$
$\displaystyle\min_{u_{k}^{\text{(EV)}},x_{k+1}^{\text{(EV)}}}\sum_{k=0}^{N-1}{\left\lVert
x_{k+1}^{\text{(EV)}}-x_{\text{ref},{k+1}}^{\text{(EV)}}\right\rVert}_{\bar{Q}}+{\left\lVert
u_{k}^{\text{(EV)}}\right\rVert}_{\bar{R}},$ (15a) $\displaystyle
s.t.\leavevmode\nobreak\ \leavevmode\nobreak\ $

$x_{k+1}^{\text{(EV)}}=f(x_{k}^{\text{(EV)}},u_{k}^{\text{(EV)}}),\leavevmode\nobreak\
\leavevmode\nobreak\ k=0,1,...,N-1,$

(15b) (15c) (15d) (15e) (15f) (15g) (15h)

Here $x_{\text{ref},{k+1}}^{\text{(EV)}}$ is the reference state of the
reference trajectory based on the relevant control mode. Note that we consider
the EV only changes one lane in the lane-change mode according to the real
traffic situation, so the number of reference trajectories in lane-change mode
depends on the current lane of the EV. $\bar{Q}\in\mathbb{R}^{6\times 6}$ and
$\bar{R}\in\mathbb{R}^{2\times 2}$ are positive definite weighting matrices
for tuning. The feasible state sets $\mathbb{X}^{\text{(EV)}}$ and
$\breve{\mathbb{X}}^{\text{(EV)}}$, and input set $\mathbb{U}^{\text{(EV)}}$
and $\breve{\mathbb{U}}^{\text{(EV)}}$ are limited by appropriate constraints,
as detailed below.

######  Remark 1.

If there is no LV in reality, it is assumed that there is an LV far away from
the EV.

######  Remark 2.

During the lane change of the EV, we call the lane-keeping control mode
deactivated when keeping the current lane is infeasible.

###  IV-C Constraints

####  IV-C1 State and input constraints

The traffic rules limit the velocity of the vehicles, and the acceleration and
jerk are bounded to have a comfortable driving feel. The lateral position of
EV is limited by the upper and lower bounds $[l_{\text{ub}},l_{\text{lb}}]$ of
the lane

$\displaystyle\noindent 0<v_{\text{lon},k}^{\text{(EV)}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ $\displaystyle l_{\text{lb}}\leq
p_{\text{lat},k}^{\text{(EV)}}\leq l_{\text{ub}}$ (16a)
$\displaystyle\noindent\underline{a}_{\text{lon}}\leq
a_{\text{lon},k}^{\text{(EV)}}\leq\overline{a}_{\text{lon}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle\underline{a}_{\text{lat}}\leq
a_{\text{lat},k}^{\text{(EV)}}\leq\overline{a}_{\text{lat}},$ (16b)
$\displaystyle\noindent\underline{j}_{\text{lon}}\leq
j_{\text{lon},k}^{\text{(EV)}}\leq\overline{j}_{\text{lon}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle\underline{j}_{\text{lat}}\leq
j_{\text{lat},k}^{\text{(EV)}}\leq\overline{j}_{\text{lat}},$ (16c)

where $\underline{\bullet}$ and $\overline{\bullet}$ denote the minimum and
maximum values of the associated variables.

####  IV-C2 Safety constraints

In highway traffic, drivers are required to maintain a safe distance from the
preceding vehicles in the same lane, which is translated into the constraint:

$d_{k}\geq\underline{d},$ (17)

and the safety distance $\underline{d}$ is computed by [22]

$\underline{d}=\tau v_{\text{lon},k}^{\text{(EV)}}+\bigtriangleup d,$ (18)

with the design parameters $\tau$ and $\bigtriangleup d$. If the reference
point of all vehicles is in their respective center, for example, choose
$\bigtriangleup d\geq\frac{l^{(\text{EV})}+l^{(\text{LV})}}{2}$, where
$l^{(\text{EV})}$ and $l^{(\text{LV})}$ are the length of EV and LV. During
the lane-change period, in addition to keeping a safe distance from the LVs in
both the current and target lane, the EV should also maintain a safe distance
with the TV behind it in the target lane [12]. The required safety distance
also satisfies (18). The safety constraint under the generated scenario is
based on (18), while the safety distance for considering the ‘worst case’
scenario collapses to $\bigtriangleup d$.

###  IV-D Recursive Feasibility of the SCMPC

The proof of the recursive feasibility is able to provide a mathematical
guarantee for the feasibility of the designed controller.

######  Definition 1.

(Recursive Feasibility) The SCMPC in lane-keeping mode is called recursively
feasible if a collision with the current LV is always avoidable, while for
SCMPC in lane-change mode, it means that no accident occurs between EV and the
other vehicles during the lane-change first, and then EV is always safe in the
target lane.

If the safety constraints in the ‘worst case’ scenario are satisfied, it
implies that the EV is able to handle all possible traffic circumstances under
the SCMPC. The corresponding capability of EV is represented by a parameter,
the minimal stopping horizon, which is defined as follows.

######  Definition 2.

(Minimal Stopping Horizon) Given the initial velocity
$v_{{\text{lon}},0}^{\text{(EV)}}$ and the minimal acceleration
$\underline{a}_{\text{lon}}^{\text{(EV)}}$ of EV, the minimal stopping horizon
$\underline{N}\in\mathbb{N}$ satisfies

$\underline{N}=\left\lceil{\frac{v_{{\text{lon}},0}^{\text{(EV)}}}{\lvert\underline{a}_{\text{lon}}^{\text{(EV)}}\rvert
T_{\text{p}}}}\right\rceil,$ (19)

where $\left\lceil\bullet\right\rceil$ is defined as the smallest integer that
is not smaller than a real number $\bullet$.

Considering the general traffic situation and rules, we make the following
assumption.

######  Assumption 1.

All vehicles only drive forwards, and the EV is only responsible for front
collisions.

######  Assumption 2.

$u_{k}^{\text{(EV)}}={\begin{bmatrix}0&0\end{bmatrix}}^{\top}$ is one element
of the feasible set $\breve{\mathbb{U}}$.

The recursive feasibility of the SCMPC is then proved below.

######  Theorem 1.

If SCMPC is initially feasible, and the prediction horizon
$N\geq\underline{N}$, then the controller is recursively feasible based on
Assumptions 1, 2.

###### Proof.

Let two initial control inputs of the generated normal scenarios and the
‘worst case’ scenario be
$\\{{{u}_{0|0}^{(\text{EV})}},{{u}_{1|0}^{(\text{EV})}},...,{{u}_{\underline{N}|0}^{(\text{EV})}},...,{{u}_{N|0}^{(\text{EV})}}\\}$
and
$\\{{\breve{u}_{0|0}^{(\text{EV})}},{\breve{u}_{1|0}^{(\text{EV})}},...,{\breve{u}_{\underline{N}|0}^{(\text{EV})}},...,{\breve{u}_{N|0}^{(\text{EV})}}\\}$.
Choose the second control sequences as initially feasible solution
$\\{{\breve{u}_{0|0}}^{{(\text{EV})}\star},{\breve{u}_{1|0}}^{{(\text{EV})}\star},...,{\breve{u}_{\underline{N}|0}}^{{(\text{EV})}\star},...,{\breve{u}_{N|0}}^{{(\text{EV})}\star}\\}$
in the proof, and its related state sequence is
$\\{{\breve{x}_{0|0}}^{{(\text{EV})}\star},{\breve{x}_{1|0}}^{{(\text{EV})}\star},...,{\breve{x}_{\underline{N}|0}}^{{(\text{EV})}\star},...,{\breve{x}_{N|0}}^{{(\text{EV})}\star}\\}$.
Due to the condition (19), the terminal set $\breve{\mathbb{X}}_{\text{f}}$
satisfies

${\breve{x}_{N|0}}^{{(\text{EV})}\star}=\begin{bmatrix}p_{{\text{lon}},\underline{N}|0}^{(\text{EV})}&0&0&p^{(\lambda)}&0&0\end{bmatrix}^{\top},$
(20)

where the stopping longitudinal position
$p_{{\text{lon}},\underline{N}|0}^{(\text{EV})}$ of EV is determined by its
initial position $p_{{\text{lon}},0}^{\text{(EV)}}$, initial velocity
$v_{{\text{lon}},0}^{\text{(EV)}}$ and minimal acceleration
${\underline{a}_{\text{lon}}^{\text{(EV)}}}$. Moreover, the terminal lateral
position of EV is the position of the center line of the target lane
$p^{(\lambda)}$ under the associated control mode. According to the system
dynamics (14), we apply ${\breve{u}_{0|0}}^{{(\text{EV})}\star}$ to the system
and obtain

$x_{1}^{{(\text{EV})}}=\bar{A}x_{0}^{{(\text{EV})}}+\bar{B}{\breve{u}_{0|0}}^{{(\text{EV})}\star}={\breve{x}_{1|0}}^{{(\text{EV})}\star}.$
(21)

Then the following is a feasible solution for the MPC problem initialized at
$x_{1}^{{(\text{EV})}}$:

$\begin{split}\\{{{u}_{1|1}},{{u}_{2|1}},...,{{u}_{\underline{N}|1}},...,{{u}_{N-1|1}},{{u}_{N|1}}\\}=\\\
\\{{\breve{u}_{1|0}}^{\star},{\breve{u}_{2|0}}^{\star},...,{\breve{u}_{\underline{N}|0}}^{\star},...,{\breve{u}_{N-1|0}}^{\star},{\begin{bmatrix}0&0\end{bmatrix}}^{\top}\\},\end{split}$
(22)

where ${\begin{bmatrix}0&0\end{bmatrix}}^{\top}\in\breve{\mathbb{U}}$. The
corresponding state sequence is

$\begin{split}\\{{{x}_{2|1}},{{x}_{3|1}},...,{{x}_{\underline{N}+1|1}},...,{{x}_{N|1}},{{x}_{N+1|1}}\\}=\\\
\\{{\breve{x}_{2|0}}^{\star},{\breve{x}_{3|0}}^{\star},...,{\breve{x}_{\underline{N}|0}}^{\star},...,{\breve{x}_{N|0}}^{\star},\bar{A}{\breve{x}_{N|0}}^{\star}+\bar{B}\cdot{\begin{bmatrix}0&0\end{bmatrix}}^{\top}\\},\end{split}$
(23)

where
$\bar{A}{\breve{x}_{N|0}}^{\star}+\bar{B}\cdot{\begin{bmatrix}0&0\end{bmatrix}}^{\top}={\breve{x}_{N|0}}^{\star}$.
Both sequences are feasible for the MPC problem because they satisfy the
dynamics and the constraints. ∎

##  V SIMULATION AND DISCUSSION

The presented approach is evaluated with one documentation (ID:$01$) of the
HighD Dataset [23], which records the motion states of $1047$ vehicles in
$900s$ with a sampling time of $0.04s$. We take two specific traffic
situations as case studies, including the EV’s initial motion states and the
TVs’ motion states during the simulation time. The designed control system
manipulates the motion behavior of the EV in the next time steps.

###  V-A Simulation Setup

Fig. 3 gives the initial traffic scenes of case studies. Table I shows the
parameters used in the two case studies. The sampling time of the simulation
is denoted as $T_{\text{s}}$. Table II shows the initial mode states of
vehicles with the proper units, and the length $l$ and width $w$ of vehicles
in Cases 1 and 2. Note that we only consider EV changes to lane 1 in the lane-
change mode of the simulation.

Figure 3: The initial traffic scenes for Case 1 (left) and Case 2 (right)
Table I: Design parameters used in the simulation study

Parameter $T_{\text{s}}(\text{s})$ $T(\text{s})$ $T_{\text{p}}(\text{s})$
Value $10.00$ $0.04$ $0.40$ Parameter $N$
$\underline{a}_{\text{lon}}^{(\text{EV})}(\text{m}/\text{s}^{2})$
$\overline{a}_{\text{lon}}^{(\text{EV})}(\text{m}/\text{s}^{2})$ Value $15$
$-4.00$ $4.00$ Parameter
$\underline{a}_{\text{lat}}^{(\text{EV})}(\text{m}/\text{s}^{2})$
$\overline{a}_{\text{lat}}^{(\text{EV})}(\text{m}/\text{s}^{2})$
$\underline{j}_{\text{lon}}^{(\text{EV})}(\text{m}/\text{s}^{2})$ Value
$-4.00$ $4.00$ $-5.00$ Parameter
$\overline{j}_{\text{lon}}^{(\text{EV})}(\text{m}/\text{s}^{2})$
$\underline{j}_{\text{lat}}^{(\text{EV})}(\text{m}/\text{s}^{2})$
$\overline{j}_{\text{lat}}^{(\text{EV})}(\text{m}/\text{s}^{2})$ Value $5.00$
$-4.00$ $4.00$ Parameter $\tau(\text{s})$ $l_{\text{ub}}(\text{m})$
$l_{\text{lb}}(\text{m})$ Value $0.40$ $33.80$ $21.00$ Parameter
$p_{c}^{(1),(2),(3)}(\text{m})$ $K_{\text{lat}}$
$\underline{a}_{\text{lon}}^{(\text{LV})}(\text{m}/\text{s}^{2})$ Value
$\\{22.98,26.88,31.3\\}$ $[1.15,3.39,3.58]^{\top}$ $-3.00$

Table II: States initial values and vehicle parameters

(a) Case 1 Parameter | $\\{p_{\text{lon},0},v_{\text{lon},0},a_{\text{lon},0},p_{\text{lat},0},v_{\text{lat},0},a_{\text{lat},0}\\}$ | $\\{l,w\\}$
---|---|---
EV | $\\{2.84,34.8,0.27,25.65,-0.09,-0.01\\}$ | $\\{4.85,2.02\\}$
TV1 | $\\{127.87,32.74,0.13,21.52,-0.21,0.2\\}$ | $\\{5.96,2.32\\}$
TV2 | $\\{141.19,23.04,0.05,25.42,0.1,-0.03\\}$ | $\\{14.35,2.5\\}$

(b) Case 2 Parameter | $\\{p_{\text{lon},0},v_{\text{lon},0},a_{\text{lon},0},p_{\text{lat},0},v_{\text{lat},0},a_{\text{lat},0}\\}$ | $\\{l,w\\}$
---|---|---
EV | $\\{181.5,25.07,-0.29,25.49,0.09,-0.01\\}$ | $\\{4.14,1.92\\}$
TV1 | $\\{151.54,32.59,0.28,22.21,-0.27,0.01\\}$ | $\\{4.75,2.02\\}$
TV2 | $\\{201.31,23.21,0.17,25.62,0.17,-0.02\\}$ | $\\{9.2,2.5\\}$

###  V-B Simulation Results

Figs. 4 and 9 show the open-loop trajectories of EV over time, the closed-loop
trajectory of EV, and TVs’ trajectory in Cases 1 and 2.

####  V-B1 Case 1

The longitudinal and lateral velocity profile of the TVs and the EV throughout
the entire sampling time are illustrated in Fig. 5. Fig. 6 displays the cost
function value of two controllers, where the associated cost of the lane-
keeping or lane-change controller is set as an arbitrary big value, like 5000,
when it is deactivated. As Fig. 6 shows, between $0s$ and $7.6s$, lane-keeping
is the desirable maneuver. Between $0s$ and $4s$, the EV keeps its velocity
and then decelerates until around $7.6s$ to have a safe distance with TV2.
After that, at $7.6s$, the cost of changing the lane becomes lower, prompting
the EV to initiate a lane-change maneuver. After approximately $1.5s$, the
lane-change maneuver is successfully executed. Subsequently, the EV maintains
a constant velocity in its new lane for the remaining time steps, while the
lane-change behavior is very costly and becomes impossible at around $11.8s$.
In particular, the increasing values between $9.2s$ and $11.8s$ lead to a
visual discontinuity in the lane-change cost-function curve, the largest of
which reaches a value of 56000.

Figure 4: The motion trajectory of vehicles in Case 1 Figure 5: The
longitudinal and lateral velocity profile in Case 1 Figure 6: The cost
function value of controllers in Case 1

The traffic prediction results for Case 1 at $5s$ and $9s$ are shown as two
examples in Fig. 7 and 8. There are three generated scenarios at 5s, which
includes TV1 performing a VT maneuver in lane 1, edge color in black, or in
lane 2, edge color in blue, and TV2 performing VT maneuvers in lane 1, edge
color in red, or in lane 2, edge color in blue. The probabilities of the
associated scenarios are
$\\{{\mu}^{(1)}_{1}{\mu}^{(2)}_{2},{\mu}^{(1)}_{2}{\mu}^{(2)}_{1},{\mu}^{(1)}_{2}{\mu}^{(2)}_{2}\\}={\\{0.241,0.186,0.573\\}}$.
At $9s$, TV1 might perform the VT maneuver in lane 2, edge color in black, or
in lane 3, edge color in magenta, and TV2 might perform a VT maneuver in lane
2, edge color in black. The associated probabilities of the generated
scenarios are
$\\{\mu^{(1)}_{2}\mu^{(2)}_{2},\mu^{(1)}_{3}\mu^{(2)}_{2}\\}={\\{0.913,0.087\\}}$.

Figure 7: Traffic prediction at $5s$ in Case 1 Figure 8: Traffic prediction at
$9s$ in Case 1

####  V-B2 Case 2

The vehicles’ motion trajectories, the velocity profile, and the cost function
value of two controllers are displayed in Figs. 9, 10 and 11. Between $0s$ and
$2s$, the EV maintains its speed under the control strategy while the cost
keeps increasing. The main reason which leads to this situation is that it is
impossible to maintain a safe distance from its LV, TV2, by consistently
keeping a higher velocity throughout the prediction horizon, so the
inconsistency between desired motion states and the calculated states from the
controller leads to an increased cost. Since then, the EV starts to reduce its
speed to the speed of TV2 until $4s$, which accordingly leads to a decrease in
the value of the cost function. Note that the lane-change becomes available
after about $5s$, before which we set its associated cost function value as
5000. After that, its cost remains higher than lane-keeping. Therefore, EV
stays in the current lane in the following time steps.

Figure 9: The motion trajectory of vehicles in Case 2 Figure 10: The
longitudinal and lateral velocity profile in Case 2 Figure 11: The cost
function value of controllers in Case 2

The traffic prediction results for Case 2 at $3s$ and $6s$ are shown in Figs.
12 and 13. There are two generated scenarios at $3s$, which includes TV1
performing a VT maneuver in lane 1, edge color in black, or in lane 2, edge
color in white, and TV2 performing VT maneuvers in lane 2, edge color in
black. The associated probabilities are
$\\{\mu^{(1)}_{1}\mu^{(2)}_{2},\mu^{(1)}_{2}\mu^{(2)}_{2}\\}={\\{0.860,0.140\\}}$.
At $6s$, there are three considered scenarios, where TV1 might perform a VT
maneuver in lane 1, edge color in black, and TV2 might perform a VT maneuver
in lane 1, in lane 2, or in lane 3, edge color in black, white, and blue. The
corresponding probabilities are
$\\{\mu^{(1)}_{1}\mu^{(2)}_{1},\mu^{(1)}_{1}\mu^{(2)}_{2},\mu^{(1)}_{1}\mu^{(2)}_{3}\\}={\\{0.130,0.727,0.143\\}}$.

Figure 12: Traffic prediction at $3s$ in Case 2 Figure 13: Traffic prediction
at $6s$ in Case 2

Simulation results of Cases 1 and 2 indicate that EV executes safe maneuvers
under the designed control structure while considering the interaction with
other vehicles.

##  VI CONCLUSIONS

In this paper, the interaction-aware estimation of motion states of the
vehicles has been studied using IMM-KF, and the associated state predictions
have been combined with the probability to represent the uncertain
environment. The generated scenarios, along with the ‘worst case’ scenario,
have been applied in formulating the safety constraints of the SCMPC. The
control system consists of lane-keeping and lane-change control modes, where
the control input with a lower cost function value is implemented into the
system. Moreover, the recursive feasibility of the method has been guaranteed
based on the no-collision between EV and LV before the minimal stopping
horizon of the EV. The proposed algorithm has been validated for two highway
scenes chosen from the HighD dataset. The simulation results demonstrate the
capability of the proposed control architecture to perform safe maneuvers.

## References

  * [1] C. Gkartzonikas and K. Gkritza, “What have we learned? a review of stated preference and choice studies on autonomous vehicles,” _Transportation Research Part C: Emerging Technologies_ , vol. 98, pp. 323–337, 2019.
  * [2] M. Brännström, E. Coelingh, and J. Sjöberg, “Model-based threat assessment for avoiding arbitrary vehicle collisions,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 11, no. 3, pp. 658–669, 2010.
  * [3] G. S. Aoude, V. R. Desaraju, L. H. Stephens, and J. P. How, “Driver behavior classification at intersections and validation on large naturalistic data set,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 13, no. 2, pp. 724–736, 2012.
  * [4] S. Lefèvre, D. Vasquez, and C. Laugier, “A survey on motion prediction and risk assessment for intelligent vehicles,” _ROBOMECH journal_ , vol. 1, no. 1, pp. 1–14, 2014.
  * [5] A. Lawitzky, D. Althoff, C. F. Passenberg, G. Tanzmeister, D. Wollherr, and M. Buss, “Interactive scene prediction for automotive applications,” in _2013 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2013, pp. 1028–1033.
  * [6] G. Agamennoni, J. I. Nieto, and E. M. Nebot, “A bayesian approach for driving behavior inference,” in _2011 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2011, pp. 595–600.
  * [7] V. Lefkopoulos, M. Menner, A. Domahidi, and M. N. Zeilinger, “Interaction-aware motion prediction for autonomous driving: A multiple model kalman filtering scheme,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 1, pp. 80–87, 2020.
  * [8] J. Kong, M. Pfeiffer, G. Schildbach, and F. Borrelli, “Autonomous driving using model predictive control and a kinematic bicycle vehicle model,” in _Intelligent Vehicles Symposium, Seoul, Korea_ , 2015.
  * [9] J. Zhou, B. Olofsson, and E. Frisk, “Interaction-aware moving target model predictive control for autonomous vehicles motion planning,” in _2022 European Control Conference (ECC)_ , 2022, pp. 154–161.
  * [10] Y. Chen, H. Yu, J. Zhang, and D. Cao, “Lane-exchanging driving strategy for autonomous vehicle via trajectory prediction and model predictive control,” _Chinese Journal of Mechanical Engineering_ , vol. 35, no. 1, p. 71, 2022.
  * [11] E. F. Camacho and C. B. Alba, _Model predictive control_. Springer science & business media, 2013.
  * [12] G. Schildbach and F. Borrelli, “Scenario model predictive control for lane change assistance on highways,” in _2015 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2015, pp. 611–616.
  * [13] C. H. Ulfsjöö and D. Axehill, “On integrating pomdp and scenario mpc for planning under uncertainty–with applications to highway driving,” in _2022 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2022, pp. 1152–1160.
  * [14] G. Cesari, G. Schildbach, A. Carvalho, and F. Borrelli, “Scenario model predictive control for lane change assistance and autonomous driving on highways,” _IEEE Intelligent transportation systems magazine_ , vol. 9, no. 3, pp. 23–35, 2017.
  * [15] T. A. Dingus, S. G. Klauer, V. L. Neale, A. Petersen, S. E. Lee, J. Sudweeks, M. A. Perez, J. Hankey, D. Ramsey, S. Gupta _et al._ , “The 100-car naturalistic driving study, phase ii-results of the 100-car field experiment,” United States. Department of Transportation. National Highway Traffic Safety Administration, Tech. Rep., 2006.
  * [16] M. A. Perez, J. D. Sudweeks, E. Sears, J. Antin, S. Lee, J. M. Hankey, and T. A. Dingus, “Performance of basic kinematic thresholds in the identification of crash and near-crash events within naturalistic driving data,” _Accident Analysis & Prevention_, vol. 103, pp. 10–19, 2017.
  * [17] Y. He, B. Ciuffo, Q. Zhou, M. Makridis, K. Mattas, J. Li, Z. Li, F. Yan, and H. Xu, “Adaptive cruise control strategies implemented on experimental vehicles: A review,” _IFAC-PapersOnLine_ , vol. 52, no. 5, pp. 21–27, 2019.
  * [18] S. Magdici and M. Althoff, “Adaptive cruise control with safety guarantees for autonomous vehicles,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 5774–5781, 2017.
  * [19] S. Shalev-Shwartz, S. Shammah, and A. Shashua, “On a formal model of safe and scalable self-driving cars,” _arXiv preprint arXiv:1708.06374_ , 2017.
  * [20] X. Xu, X. Wang, X. Wu, O. Hassanin, and C. Chai, “Calibration and evaluation of the responsibility-sensitive safety model of autonomous car-following maneuvers using naturalistic driving study data,” _Transportation research part C: emerging technologies_ , vol. 123, p. 102988, 2021.
  * [21] J. Sudano, “The alpha - beta - gamma tracking filter with a noisy jerk as the maneuver model,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 29, no. 2, pp. 578–580, 1993.
  * [22] R. Schmied, H. Waschl, and L. Del Re, “Extension and experimental validation of fuel efficient predictive adaptive cruise control,” in _2015 American Control Conference (ACC)_. IEEE, 2015, pp. 4753–4758.
  * [23] R. Krajewski, J. Bock, L. Kloeker, and L. Eckstein, “The highd dataset: A drone dataset of naturalistic vehicle trajectories on german highways for validation of highly automated driving systems,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_ , 2018, pp. 2118–2125.


# Latent Diffusion Model for DNA Sequence Generation

Zehui Li
Imperial College London
<EMAIL_ADDRESS>&Yuhao Ni11footnotemark: 1
Imperial College London
<EMAIL_ADDRESS>&Tim August B. Huygelen
University College London
<EMAIL_ADDRESS>&Akashaditya Das
Imperial College London
<EMAIL_ADDRESS>&Guoxuan Xia
Imperial College London
<EMAIL_ADDRESS>&Guy-Bart Stan
Imperial College London
<EMAIL_ADDRESS>&Yiren Zhao22footnotemark: 2
Imperial College London
<EMAIL_ADDRESS>Equal contribution.Correspondence should be addressed
to<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

The harnessing of machine learning, especially deep generative models, has
opened up promising avenues in the field of synthetic DNA sequence generation.
Whilst Generative Adversarial Networks (GANs) have gained traction for this
application, they often face issues such as limited sample diversity and mode
collapse. On the other hand, Diffusion Models are a promising new class of
generative models that are not burdened with these problems, enabling them to
reach the state-of-the-art in domains such as image generation. In light of
this, we propose a novel latent diffusion model, DiscDiff, tailored for
discrete DNA sequence generation. By simply embedding discrete DNA sequences
into a continuous latent space using an autoencoder, we are able to leverage
the powerful generative abilities of continuous diffusion models for the
generation of discrete data. Additionally, we introduce Fréchet Reconstruction
Distance (FReD) as a new metric to measure the sample quality of DNA sequence
generations. Our DiscDiff model demonstrates an ability to generate synthetic
DNA sequences that align closely with real DNA in terms of Motif Distribution,
Latent Embedding Distribution (FReD), and Chromatin Profiles. Additionally, we
contribute a comprehensive cross-species dataset of 150K unique promoter-gene
sequences from 15 species, enriching resources for future generative modelling
in genomics. We will make our code public upon publication.

##  1 Introduction

Designing synthetic DNA sequences for genetic modifications is traditionally
guided by organism-specific workflows derived from extensive laboratory
experiments. As the amount of data generates from these workflows expand, deep
generative models are well-positioned to enable a new frontier in synthetic
DNA sequence generation over a wide range of potential applications [12, 37].
Generative adversarial networks (GANs) [11] are a popular choice for the
generation of synthetic DNA sequences, as demonstrated by the studies of [18,
34, 37]. GANs can effectively generate sequences, however, it has been shown
that the generated samples lack diversity [8] and suffer from mode collapse at
training time [22].

Given the success of diffusion models in image generation [26], protein
synthesis [36], and circuit design [35]. The application of diffusion models
for DNA sequence generation may generate sequences of higher quality. In this
case, quality is a combination of sequence diversity and the ability to
capture the underlying distribution/motifs. In this paper, we propose a latent
diffusion model for discrete data generation and apply it to DNA sequence
generation. Our contributions are as follows:

#### A generalisable framework for discrete data generation:

We propose a latent diffusion model made from a transformation function and a
denoising model. It separates the learning of latent distribution and the
denoising model. A lightweight Variational Auto-encoder architecture is used
for the transformation function for 1-dimensional discrete sequences.

#### A generative model for DNA Sequence Generation:

We apply the proposed latent diffusion model to DNA sequence generation,
obtaining promoter-gene samples with the properties of realistic DNA sequences
across different species.

#### An evaluation metric for generated DNA sequences:

We introduce Fréchet Reconstruction Distance (FReD) as a numerical metric to
evaluate the quality of the generated DNA sequences. FReD is consistent with
motif distribution, the standard qualitative approach to evaluate the
generated DNA sequences.

#### A cross-species dataset for DNA sequence generation:

We curate and create a dataset containing 150K unique promoter-gene sequences
across 15 species. This dataset paves the way for the construction of a
generative model for promoter sequences

(a) Inference

(b) Training

Figure 1: Denoising View of Latent Diffusion Models for DNA sequences. (a) A
zoom-in view of the inference process. (b) A zoom-in view of the training
process. The blue square represents the noises attached to the latent
representation of the DNA sequences.

##  2 Background

#### Diffusion Models

Intuitively, Diffusion Models (DMs) incrementally construct realistic examples
from sampled noise. At each time step $t$, a denoising function is used to
predict a more realistic example $x_{t-1}$ from $x_{t}$. An alternative
perspective by Song et al. [31] uses stochastic differential equations (SDE)
to describe the diffusion models. This unifies previous scored-based models
[30] and DDPM [15], providing a powerful formulation for analyzing DMs.

In [31], the process of adding noises is modelled with an SDE with a Wiener
process:

$\mathop{}\\!\mathrm{d}{x}=f(x,t)\mathop{}\\!\mathrm{d}{t}+g(t)\mathop{}\\!\mathrm{d}{w}$
(1)

  * •

The function $f$ is called the drift coefficient. This indicates how the
process $x$ tends to evolve over time, ignoring any random fluctuations.
Multiplying by $\mathop{}\\!\mathrm{d}{t}$ scales the drift by the time
increment.

  * •

$g(t)$ is called the diffusion coefficient. It measures the volatility or
variability of the process. It only depends on time $t$.

  * •

$\mathop{}\\!\mathrm{d}{w}$ represents the increment of the standard Wiener
process (or Brownian motion) $w$ over an infinitesimally small time interval.

The reverse process, which generates the realistic samples, can be represented
by another SDE:

$\mathop{}\\!\mathrm{d}{x}=[f(x,t)-{g(t)}^{2}\nabla_{x}\log
p_{t}(x)]\mathop{}\\!\mathrm{d}{t}+g(t)\,\mathop{}\\!\mathrm{d}{\overline{w}}$
(2)

Where $\mathop{}\\!\mathrm{d}{\overline{w}}$ represents the increment of the
reverse Wiener process. $\nabla_{x}\log p_{t}(x)$ is the gradient of the log
Probability Density Function (PDF) of $x$ at time $t$. At the training time, a
score function $\mathbf{S_{\theta}}(x(t),t)$ is trained to approximate
$\nabla_{x}\log p_{t}(x)$; and then realistic samples can be generated by
walking through Equation 2.

#### Diffusion Models for Discrete Data

There are two reasons why standard DMs cannot deal with discrete data ($x$):
1) the diffusion process with the Wiener process defined by Equation 1 is
described in terms of continuous variables, it is undefined when $x$ is
discrete. 2) the solution led by the reverse process requires a differentiable
PDF $\log p_{t}(x)$ to exist, but discrete variable $x$ does not have a PDF.

To account for these challenges, two possible solutions can be used to develop
DMs for discrete data. They are 1) define a diffusion-like process to model
discrete data [3, 32, 1]. 2) map the discrete input into a continuous latent
space [9, 19, 13]. For more details, see Appendix A.

Differing from the prior work on latent DMs for discrete data [9, 33], which
is trained in an end-to-end manner, we aim to separate learning of the mapping
from discrete to the latent space and score function. This facilitates the
learning task by not learning the latent distribution and score-based
denoising model at the same time. It also removes the need to adjust the
weight of reconstruction and generation in the end-to-end training [26].

##  3 DiscDiff: A Latent Diffusion Model for Discrete Data

We introduce DiscDiff, a flexible latent diffusion model designed for discrete
data generation. This model is structured into two main components: a
transformation function and a denoising model. The transformation function is
implemented with a lightweight Variational-Auto-Encoder (VAE), with an encoder
$z=\mathbf{E}(x)$, which translates discrete input $x$ to a continuous latent
variable $z$, and a decoder $\tilde{x}=\mathbf{D}(z)$, which reverts $z$ back
to its discrete form $\tilde{x}$. The denoising model
$\mathbf{S_{\theta}}(z(t),t)$ is employed to learn the score function
$\nabla_{z}\log p_{t}(z)$ for the latent variable.

In the training process, the learning phases of the transformation function
and the denoising model are separated. The first phase focuses on learning the
transformation function $\tilde{x}=\mathbf{D}(\mathbf{E}(x))$, with a primary
goal of minimising the reconstruction loss for discrete variable $x$. The
second phase concentrates on training the score function estimator, aiming to
minimise the difference between the estimator and $\nabla_{z}\log p_{t}(z)$.

###  3.1 Transformation Function

Figure 2: VAE as a transformation function. (a) shows the components of the
encoder, consisting of conv1d and conv2d. It lifts the channel dimension and
reduces the length dimension of the input. (b) shows the overall view, the
latent variable is modelled with Diagonal Gaussian Distribution. The Decoder
is symmetric to the Encoder architecture, consisting of transConv and
upsampling. See Appendix D for the values of $c,w,h,\ldots$.

#### Architecture

In order to map discrete DNA data to continuous space, we employed a
lightweight Variational Autoencoder (VAE) that employs a transformation
function composed of Convolutional Neural Networks (CNNs). The structure of
the VAE is depicted in Figure 2, illustrating how the latent variable is
modelled with a multivariate Gaussian distribution. The discrete input passes
through the encoder, through multiple layers designed to increase the channel
dimension while reducing the length dimension. Next, a Conv2D block is used to
extract essential features from the channel-length surface, mapping the input
to a higher-dimensional space. In symmetry to the Encoder, the Decoder is
constructed with transConv and upsampling operations, effectively inverting
the encoding process.

The proposed VAE stands out for its lightweight nature. This is attributed to
the CNN-based architecture which enables it to efficiently extract information
from 1D discrete data and map it to a high-dimensional continuous space.
Additional details about the VAE architecture, including the incorporation of
a multi-kernel design [7] in the Conv1D block for feature extraction at
various scales, are provided in Appendix D.

#### Cross Entropy Loss

When training the VAE, we propose to use Cross Entropy (CE) as reconstruction
loss. The loss function is given by:

$\mathbf{L}_{\theta,\phi}=\underbrace{\mathbb{E}_{z\sim
q_{\phi}(z|x)}\left[-\sum_{t=1}^{k}\sum_{l=1}^{n}p(x_{t},l)\log
p_{\theta}(x_{t}|z)\right]}_{\text{Reconstruction Loss}}+\mathbb{E}_{x\sim
p(x)}\left[\underbrace{\mathrm{KL}(q_{\phi}(z|x)\,||\,\mathcal{N}(z;\mu,\Sigma))}_{\text{KL
Divergence}}\right]$ (3)

where $x\in\mathbb{R}^{k\times n}$ is onehot-encode data; $p_{\theta}(x|z)$ is
the probabilistic decoder output from $\mathbf{D_{\theta}}$; $q_{\phi}(z|x)$
is the probabilistic output from encoder $\mathbf{E_{\phi}}$ that represents
the approximate posterior of the latent variable $z$ given the input $x$;
$\mathcal{N}(z;\mu,\Sigma)$ is the prior on $z$, here we use diagonal Gaussian
Distribution.

###  3.2 Diffusion in the Latent Space

Once the transformation function with $\mathbf{D_{\theta}}$ and
$\mathbf{E_{\phi}}$ are trained in the first stage, the diffusion model is
applied to the latent representation $z=\mathbf{E_{\phi}}(x)$. With the goal
of minimizing the difference between the estimator and the score function in
the latent space.

$min(\mathbf{\hat{S}_{\theta}}(z(t),t)-\nabla_{z}\log p_{t}(z))^{2}$ (4)

It is trained by sampling $t$ from uniform distribution $U[1,T]$ and $z_{t}$
from $p_{t}(z_{t}|z_{0})$ with a known transition kernel defined by Equation
1. The training process can also be interpreted as a training a noise
predictor [26]:

$\mathbb{E}_{z,t\sim
U[1,T],\varepsilon\sim\mathcal{N}(0,1)}\left[\|\varepsilon-\mathbf{\varepsilon_{\theta}}(z(t),t)\|^{2}_{2}\right]$
(5)

We adopted the UNet backbone from LDM [27] to implement
$\varepsilon_{\theta}(z(t),t)$. The UNet can capture coarse-grained features
and generate samples with high qualities in various domains [16, 20]. Figure 1
illustrates the denoising view of training and inference processes for the
diffusion model. The details about the UNet configuration can be found in
Appendix E.

###  3.3 Justification for the Two-stage Training

The separation of VAE and the denoising model training can be justified by the
loss function of LSGM [33]. LSGM jointly trains a VAE and the denoising model
with a training loss consisting of a reconstruction term and KL divergence
between encoder distribution $p(z_{0}|x)$ and the prior distribution of the
latent variable $p(z_{0})$. The latter term can be expressed in terms of the
score function $\mathbf{S_{\theta}}(x(t),t)$. The optimisation of this loss is
challenging, requiring different weighting mechanisms for the training
objective and variance reduction techniques for stability. The separation of
training is equivalent to a special training schedule with the KL divergence
term being fixed in the first stage, and then training the denoising model in
the second stage. The performance gain of separating the training process has
been shown in the image generation domain [26, 25]. In the section below, we
provide the empirical evaluation of this method for DNA sequence generation.

Figure 3: Dataset and evaluation methods for DNA sequence generation.(a) shows
the format of the curated dataset from EPDnew (b) Motif Distribution compares
the distribution of subsequences in DNA and Latent Distance metrics compare
the distribution of latent embedding of DNA.

##  4 DNA Sequence Generation with DiscDiff

We explored the unconditional generation of both regulatory elements and
protein-encoding genes. Our problem setting is different from the previous
study DDSM [1], where the transcription profile is provided as an input at
both training and inference time. In this study, we focus on a scenario where
no transcription profile is given. The generative model needs to discover the
properties of regulatory elements and genes purely from DNA sequences.
Identifying and generating regulatory elements is essential for successful
genetic modifications [4]. To achieve this, we constructed a dataset, trained
the diffusion model for DNA generation, and provided an evaluation metric
called Fréchet Reconstruction Distance (FReD) to assess the quality of
generated samples. For an in-depth look into the training specifics, including
hyperparameters and training equipment, see Appendix B.

###  4.1 Data Construction and Representation

#### Dataset

We created a cross-species dataset for DNA generation from the refined
Eukaryotic Promoter Database (EPDnew) [10, 23], housing organism-specific
promoters annotated with downstream genes, transcription starting sites, and
expression levels. From EPDnew, we processed 160K unique sequences across
various species, each mapping to different cell types and expression levels,
detailed in Table 2 of Appendix C. This study emphasizes unconditional
generation using a subset of this data. We employed a 2048-base context window
centred on a gene’s transcription starting site. Each sequence split evenly
into an upstream promoter and a downstream DNA segment, spans 2048 bases in
total.

#### Representation of DNA Sequences

DNA sequences are composed of the nucleotide bases $A$, $T$, $G$, and $C$. For
a DNA sequence of length 2048, a one-hot encoding function, denoted as $f$, is
employed:

$f:\\{A,T,G,C\\}^{2048}\rightarrow\mathbb{R}^{4\times 2048}$

It is worth noting that DNA sequences sometimes contain the $N$ token,
indicating an unknown or ambiguous nucleotide. In the encoding scheme, such a
token is represented as $\\{0.25,0.25,0.25,0.25\\}$. This representation
allows generative models to directly model the aleatoric uncertainty arising
from the data collection process when using a cross-entropy reconstruction
loss for VAE training.

###  4.2 Evaluation Metrics

#### Motif Distribution

A motif in the context of molecular biology refers to a short sequence of DNA
or RNA that has a specific structure or function. Motif distribution is a
commonly used metric to evaluate the biological relevance of generated DNA
sequences, which is frequently employed in prior studies [1, 18, 34, 37]. DNA
motifs with biological functions can shed light on the functionality of
generated sequences. By comparing the distribution of motifs between real and
generated sequences, we can measure the ability of the generative model to
replicate the underlying patterning of natural DNA sequences. A high-quality
generative model should exhibit a motif distribution closely mirroring that of
genuine DNA sequences, ensuring that generated sequences retain crucial
biological signatures.

#### Fréchet Reconstruction Distance

We introduce the Fréchet Reconstruction Distance (FReD) to assess the quality
of generated DNA samples quantitatively. Unlike the Fréchet Inception Distance
used for image generation evaluation [14], FReD leverages the encoder of a
pre-trained Auto-Encoder (AE) to transform DNA sequences into embeddings. This
encoding process measures the disparity between generated sample distributions
and real-world sequences. For this purpose, we train an AE, following the
architecture outlined in Section 3, on a reference genome distinct from the
training data used for generation. The embeddings are subsequently derived
from this encoder.

#### Evaluating Generated Sequences with the Sei Framework

Sei [5] is a deep-learning framework dedicated to predicting the chromatin
profile of the human genome. We used Sei to evaluate the generated sequences
in two ways. First, we employ Sei to encode the sequences and subsequently
calculate the distribution distance. As highlighted in Section 5, there’s a
consistency between the results of FReD and Sei. Second, using Sei, we
obtained predictions for 21,907 distinct chromatin profiles
$Sei(x)\in\mathbb{R}^{50000\times 21907}$ for 50,000 generated sequences. For
each profile, we calculated the number of "hits" among our sequences. A "hit"
is defined as a sequence for which the predicted likelihood of aligning with a
profile exceeds 0.9.

##  5 Results

Figure 4: TATA-Box Distribution in Real and Generated DNA Sequences. The left
plot represents mammalian sequences, while the right showcases plant
sequences. Both generated samples closely mirror the distribution of their
respective real DNA sequences. Figure 5: TATA-Box Distribution Across Epochs
Relative to Real DNA. Subplots from left to right showcase the motif
distribution by DiscDiff at epochs 0, 200, 1000, and 3000. With increasing
epochs, the peak distribution around TSS converges towards that of the real
DNA, but the background distribution starts to diverge after 200 epochs.

#### Motif Distribution

To assess the quality of the generated samples, we generated 50,000 DNA
sequences for both mammalian and plant species using DiscDiff. Their motif
distributions are illustrated in Figure 4. The plots reveal consistency
between the TATA-box [17] distributions of real DNA sequences and our
generated promoters. Additionally, Figure 5 demonstrates the evolution of the
motif distribution throughout training. Noticeably, while the peak
distribution around the Transcription Starting Site (TSS) is converging to the
real DNA sequences, the background distribution seems to be diverging after
200 training epochs. This trend is also captured by FReD and Sei Distance.

Figure 6: Change of FReD and Sei Embedding Distribution Distance across epochs. Both metrics show a rapid decrease until epoch 200, then gradually rise until epoch 3000. This behaviour indicates the complexity of capturing DNA sequence quality with a single value. Despite improvements in modelling the TSS peak, divergences in background motifs appear over extended training. Model | FReD$\downarrow$ | Sei Embedding Distribution Distance$\downarrow$
---|---|---
Random | 251.7 | 250.7
Sample from Training Set | 3.38 | 0.10
VAE | 48.66 | 84.77
DiscDiff 200 | 22.10 | 45.23
DiscDiff 3000 | 29.25 | 57.28
Table 1: Comparison of FReD and Sei Embedding Distribution Distance across
different models.

#### Latent Distribution Distance

Figure 6 presents the change of FReD and Sei Embedding Distribution Distance
values relative to the training set across epochs. Notably, these metrics
exhibit strong correlation to the training set: a sharp decline is observed
from epoch 0 to 200, followed by a gradual increase up to epoch 3000. This
trend highlights the intricacy of measuring the quality of generated DNA
sequences using a singular numerical metric. We attribute the rise in these
metrics (epoch 200 to 3000) to divergences in the background motif
distribution. Even as the modelling of the TSS peak improves with prolonged
training. Embedding-based approaches tend to prioritize the holistic
representation of DNA sequences over specific details. However, latent
distribution distances remain crucial as they help distinguish genuine DNA
from random or subpar sequences. As per Table 1, when comparing VAE and
DiscDiff, VAE-generated examples fare less favourably in motif distribution
(Appendix F).

#### Chromatin Profiles of Generated Sequences

Figure 7 presents the chromatin profiles of the 50,000 generated and real DNA
sequences. The y-axis indicates the count of sequences corresponding to each
profile. Among these, we highlight the top 10 profiles with the highest
counts, omitting cell line names for clarity. There’s a striking resemblance
between the generated sequences (Figure 7(a)) and the training data (Figure
7(b)) in terms of distribution and top-ranking profiles. Notably, profiles
such as H3K4me3, H3K27me3 and H3K9me3 are predominant. H3Kxxme3 markers are
closely linked to promoter activity as they remodel chromatin to be more
accessible to transcription factors (essential proteins for promoter
regulation). This observation aligns with our expectations since our training
set is derived from the promoter database.

(a) Generated DNA

(b) Real DNA

Figure 7: Sei Predicted Chromatin Profiles Distribution of Generated and Real
DNA Sequences. Chromatin Profiles of (a) Generated DNA are aligned well with
(b) Real DNA. The y-axis indicates the count of sequences for each profile,
with the top 10 profiles highlighted. Predominant profiles like H3K4me3,
H3K27me3, and H3K9me3, are linked to promoter activity.

##  6 Discussion

We presented DiscDiff, a latent diffusion model for DNA sequence generation.
By connecting discrete DNA sequences to continuous spaces with an autoencoder,
DiscDiff taps into diffusion model capabilities. Our Fréchet Reconstruction
Distance (FReD) offers improved genomics evaluation, and our model’s alignment
with real DNA benchmarks showcases its efficacy.

Looking ahead, our approach raises important questions: How might we design
general transformation functions optimally for this framework? What properties
make a function particularly appealing for the diffusion model, ensuring the
production of superior-quality samples? It has been shown in image generation,
that the choice of auto-encoder will influence the quality of generated
samples [26]. On the practical side, exploring the conditional generation of
DNA sequences influenced by factors such as cell type, expression level, and
environment holds potential. We envision a scenario where one could generate a
functional gene or regulatory element based solely on specified criteria. In
addition, the scarcity of experimental data for certain types of DNA sequences
advocates the strategy of fine-tuning pre-trained diffusion models with
existing methods [2, 28]. Such advancements could redefine synthetic biology,
and we are optimistic about delving deeper into these avenues.

## References

  * Avdeyev et al. [2023] Pavel Avdeyev, Chenlai Shi, Yuhao Tan, Kseniia Dudnyk, and Jian Zhou.  Dirichlet diffusion score model for biological sequence generation.  arXiv preprint arXiv:2305.10699, 2023.
  * Blattmann et al. [2022] Andreas Blattmann, Robin Rombach, Kaan Oktay, Jonas Müller, and Björn Ommer.  Semi-parametric neural image synthesis.  Advances in Neural Information Processing Systems, 11, 2022.
  * Campbell et al. [2022] Andrew Campbell, Joe Benton, Valentin De Bortoli, Thomas Rainforth, George Deligiannidis, and Arnaud Doucet.  A continuous time framework for discrete denoising models.  Advances in Neural Information Processing Systems, 35:28266–28279, 2022.
  * Cazier and Blazeck [2021] Andrew P Cazier and John Blazeck.  Advances in promoter engineering: novel applications and predefined transcriptional control.  Biotechnology Journal, 16(10):2100239, 2021.
  * Chen et al. [2022] Kathleen M Chen, Aaron K Wong, Olga G Troyanskaya, and Jian Zhou.  A sequence-based global map of regulatory activity for deciphering human genetics.  Nature genetics, 54(7):940–949, 2022.
  * Chen et al. [2022] Ting Chen, Ruixiang Zhang, and Geoffrey Hinton.  Analog bits: Generating discrete data using diffusion models with self-conditioning.  arXiv preprint arXiv:2208.04202, 2022.
  * Cui et al. [2016] Zhicheng Cui, Wenlin Chen, and Yixin Chen.  Multi-scale convolutional neural networks for time series classification.  arXiv preprint arXiv:1603.06995, 2016.
  * Dhariwal and Nichol [2021] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  Advances in neural information processing systems, 34:8780–8794, 2021.
  * Dieleman et al. [2022] Sander Dieleman, Laurent Sartran, Arman Roshannai, Nikolay Savinov, Yaroslav Ganin, Pierre H Richemond, Arnaud Doucet, Robin Strudel, Chris Dyer, Conor Durkan, et al.  Continuous diffusion for categorical data.  arXiv preprint arXiv:2211.15089, 2022.
  * Dreos et al. [2015] René Dreos, Giovanna Ambrosini, Rouayda Cavin Périer, and Philipp Bucher.  The eukaryotic promoter database: expansion of epdnew and new promoter analysis tools.  Nucleic acids research, 43(D1):D92–D96, 2015.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
  * Gupta and Zou [2019] Anvita Gupta and James Zou.  Feedback gan for dna optimizes protein functions.  Nature Machine Intelligence, 1(2):105–111, 2019.
  * Han et al. [2022] Xiaochuang Han, Sachin Kumar, and Yulia Tsvetkov.  Ssd-lm: Semi-autoregressive simplex-based diffusion language model for text generation and modular control.  arXiv preprint arXiv:2210.17432, 2022.
  * Heusel et al. [2017] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  Advances in neural information processing systems, 30, 2017.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in neural information processing systems, 33:6840–6851, 2020.
  * [16] Jonathan Ho, Tim Salimans, Alexey Gritsenko, William Chan, Mohammad Norouzi, and David J Fleet.  Video diffusion models, 2022.  URL https://arxiv. org/abs/2204.03458.
  * Juo et al. [1996] Zong Sean Juo, Thang Kien Chiu, Paul M Leiberman, Igor Baikalov, Arnold J Berk, and Richard E Dickerson.  How proteins recognize the tata box.  Journal of molecular biology, 261(2):239–254, 1996.
  * Killoran et al. [2017] Nathan Killoran, Leo J Lee, Andrew Delong, David Duvenaud, and Brendan J Frey.  Generating and designing dna with deep generative models.  arXiv preprint arXiv:1712.06148, 2017.
  * Li et al. [2022] Xiang Li, John Thickstun, Ishaan Gulrajani, Percy S Liang, and Tatsunori B Hashimoto.  Diffusion-lm improves controllable text generation.  Advances in Neural Information Processing Systems, 35:4328–4343, 2022.
  * Liu et al. [2023] Haohe Liu, Zehua Chen, Yi Yuan, Xinhao Mei, Xubo Liu, Danilo Mandic, Wenwu Wang, and Mark D Plumbley.  Audioldm: Text-to-audio generation with latent diffusion models.  arXiv preprint arXiv:2301.12503, 2023.
  * Lu et al. [2022] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu.  Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps.  Advances in Neural Information Processing Systems, 35:5775–5787, 2022.
  * Metz et al. [2016] Luke Metz, Ben Poole, David Pfau, and Jascha Sohl-Dickstein.  Unrolled generative adversarial networks.  arXiv preprint arXiv:1611.02163, 2016.
  * Meylan et al. [2020] Patrick Meylan, René Dreos, Giovanna Ambrosini, Romain Groux, and Philipp Bucher.  Epd in 2020: enhanced data visualization and extension to ncrna promoters.  Nucleic Acids Research, 48(D1):D65–D69, 2020.
  * Mikolov et al. [2013] Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  Efficient estimation of word representations in vector space.  arXiv preprint arXiv:1301.3781, 2013.
  * Preechakul et al. [2022] Konpat Preechakul, Nattanat Chatthee, Suttisak Wizadwongsa, and Supasorn Suwajanakorn.  Diffusion autoencoders: Toward a meaningful and decodable representation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10619–10629, 2022.
  * Rombach et al. [2022] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 10684–10695, 2022.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
  * Ruiz et al. [2023] Nataniel Ruiz, Yuanzhen Li, Varun Jampani, Yael Pritch, Michael Rubinstein, and Kfir Aberman.  Dreambooth: Fine tuning text-to-image diffusion models for subject-driven generation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 22500–22510, 2023.
  * Song et al. [2020] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  arXiv preprint arXiv:2010.02502, 2020.
  * Song and Ermon [2019] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  Advances in neural information processing systems, 32, 2019.
  * Song et al. [2020] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * Sun et al. [2022] Haoran Sun, Lijun Yu, Bo Dai, Dale Schuurmans, and Hanjun Dai.  Score-based continuous-time discrete diffusion models.  arXiv preprint arXiv:2211.16750, 2022.
  * Vahdat et al. [2021] Arash Vahdat, Karsten Kreis, and Jan Kautz.  Score-based generative modeling in latent space.  Advances in Neural Information Processing Systems, 34:11287–11302, 2021.
  * Wang et al. [2020] Ye Wang, Haochen Wang, Lei Wei, Shuailin Li, Liyang Liu, and Xiaowo Wang.  Synthetic promoter design in escherichia coli based on a deep generative network.  Nucleic Acids Research, 48(12):6403–6412, 2020.
  * Wang et al. [2023] Zixiao Wang, Yunheng Shen, Wenqian Zhao, Yang Bai, Guojin Chen, Farzan Farnia, and Bei Yu.  Diffpattern: Layout pattern generation via discrete diffusion.  arXiv preprint arXiv:2303.13060, 2023.
  * Watson et al. [2023] Joseph L Watson, David Juergens, Nathaniel R Bennett, Brian L Trippe, Jason Yim, Helen E Eisenach, Woody Ahern, Andrew J Borst, Robert J Ragotte, Lukas F Milles, et al.  De novo design of protein structure and function with rfdiffusion.  Nature, pages 1–3, 2023.
  * Zrimec et al. [2022] Jan Zrimec, Xiaozhi Fu, Azam Sheikh Muhammad, Christos Skrekas, Vykintas Jauniskis, Nora K Speicher, Christoph S Börlin, Vilhelm Verendel, Morteza Haghir Chehreghani, Devdatt Dubhashi, et al.  Controlling gene expression with deep generative design of regulatory dna.  Nature communications, 13(1):5099, 2022.

##  Appendix A Extended Related Work of Discrete Diffusion Models

Two possible solutions has been used to develop DMs for discrete data. The
first type of method defines diffusion-like processes to model the discrete
data [3, 32]. The Dirichlet Diffusion Score Model proposed by [1] uses the
Dirichlet distribution as the stationary distribution in the diffusion process
and generates regulatory elements in DNA sequences.

The second approach is to map the discrete input into a continuous latent
space. This approach has been used by several works for language modelling [9,
19, 13], where text is mapped to a continuous space through word2vec [24] or
embedding obtained by neural networks. BitDiffusion [6] also falls into this
category, where discrete data is transformed into analogue bits and then real
numbers. There are several advantages of mapping discrete data into continuous
space. For example, effective sampling methods [21, 29] to accelerate the
generation of DMs are still applicable. Moreover, compression is allowed when
mapping the discrete data into the latent space, reducing the computational
overhead for both training and inference of DMs.

##  Appendix B Training Details of Latent Diffusion Model

The model was trained utilizing the NVIDIA RTX A6000 graphics card. During the
first stage, a Variational Autoencoder (VAE) was trained on a dataset
comprising 160,000 sequences spanning 15 different species. This training was
executed over 24 GPU hours with the primary goal of accurately reconstructing
the input sequences. The chosen learning rate was $0.0001$, paired with the
Adam Optimizer. A batch size of $128$ was set, and the KL divergence weight
was determined at $10^{-4}$. The resultant hidden representation, $z$, belongs
to $\mathbb{R}^{16\times 16\times 16}$, which is precisely half the
dimensionality of the original input $x\in\mathbb{R}^{4\times 2048}$.

In the second training phase, a denoising UNet was trained on three GPUs
across $3,000$ epochs. Our empirical observations highlighted that the UNet’s
performance improved with increased size. As such, given the GPU RAM capacity
of $49$ GB, we optimized the model by maximizing the number of residual blocks
and self-attention layers. The training also incorporated a cosine learning
scheduler, introducing a warm-up period for the initial $500$ epochs. The
learning rate was set at $0.00005$, and a batch size of $256$ was chosen for
the training. DDPM [15] is used for sampling.

For a focused evaluation of the model’s performance, particularly in mammalian
DNA sequence generation, we opted for a subset of the entire dataset. This
subset includes sequences from H. Sapiens (human), Rattus Norvegicus (rat),
Macaca mulatta, and Mus musculus (mouse), which collectively represent 50% of
the total dataset. Training on this subset not only streamlines the evaluation
process but also allows for a more precise assessment of the generative
algorithm’s efficacy and accuracy in producing mammalian sequences

##  Appendix C Dataset Details

Table 2: The statistics of EPDnew dataset.  | # Unique Values | Examples
---|---|---
genes | 130,014 | CAACAAGGCACAGCAC…
promoters | 159,125 | GGCCCGGTTTTTTAAT…
cell types | 2713 | [Esophageal Epithelial Cells, Hippocampus, …]
gene function description | 130,014 | [AHSP:alpha hemoglobin stabilizing protein, …]
species | 15 | [Homo sapiens (Human).,Zea mays (corn)., …]
expression levels | 29,349,475 | [250,264,338, …]
transcription starting site | 29,349,475 | [10,0,-1, …]

#### Source Description

The European Promoter Database (EPD) is a vital tool for biologists, with a
primary focus on gathering and categorising promoters of eukaryotes. The EPD
comprises promoter data sourced from numerous published articles. The EPDnew
database, or High-Throughput EPD (HT-EPD), provides an even more comprehensive
and rigorous dataset. The EPDnew has been enhanced by combining conventional
EPD promoters with in-house analyses of promoter-specific high-throughput
data, albeit for only selected organisms. This strategy enhances the precision
and commendable coverage of EPDnew. It’s worth noting that the evidence used
in EPDnew primarily arises from Transcription Start Site (TSS) mapping derived
from high-throughput experiments, notably CAGE and Oligocapping.

#### Data Selection

Data was obtained from the EPDnew database without any exclusions, allowing
for the incorporation of all available samples and associated promoters.
Sequences between -1024 and +1023 base pairs from the transcription start site
were chosen, ensuring a total length of 2048 base pairs - a convenient choice
due to its alignment to powers of 2. This range covers crucial promoter
elements like GC and TATA boxes, and also allows the algorithm to identify any
possible hidden patterns and interconnections preceding and succeeding the
transcription start site.

#### Time Frame

The temporal aspect of our data is not limited to a specific time frame or
year. The data was collected entirely as of August 2023.

#### Data collection and preprocessing

Data collection started by manually downloading data files through FTP,
followed by accessing fasta data via the website’s graphical interface.
Carefully extracted metadata from Python scripts was then integrated with
promoter sequences to form a comprehensive table.

Given the species-specific fragmentation of the acquired data, a iterative
method aided in bringing together diverse FASTA files into one consistent
table, with important attributes such as ’kingdom’ and ’species’ integrated
into the DataFrame at extraction. A crucial pre-processing step involved
categorising DNA sequences as ’upstream’ or ’downstream’ based on their
relative positioning. Normalisation and transformation techniques were
implemented where necessary to ensure data uniformity and integrity were
maintained.

Simultaneously, an additional dataset containing gene expression data
underwent similar aggregation processes. Promoter expression sample data was
synchronised with the main DataFrame. To optimize data storage and retrieval,
the processed data was archived as a Pickle file.

Data integration and duplicate removal processes were also undertaken.
Subsequently, the integrity and consistency of the data have been validated
through multiple checks, including the verification of value counts for
specific columns and a comparison of the data against EPDnew statistics.

#### Data Variables

  * •

Promoter Sequences: -1024 to 1023 bp from tss

  * •

PID: A unique identifier for each promoter in the EPDnew database.

  * •

Species: one of 15 species in our dataset

  * •

Average Expression Value: mean of expression values of every sample pertaining
to a specified EPDnew promoter

  * •

Sample Name: Encompasses details like cell line and environmental conditions.

  * •

Gene Description: description of gene function of the gene the promoter is
regulating

#### Data Size and Structure

As shown in Table 2, the data set we compiled from EPDnew comprises of 159,125
distinct promoter sequences derived from 15 eukaryotic species, which are
connected to 130,014 individual genes. Most of these genes have gene function
definitions associated with them in the gathered data. The promoters have
experimental expression levels and co-determined transcription start sites
associated with them. These are frequently recorded in various cell types per
promoter and might hold many experimental values for a single cell type. Our
methodically structured dataset is in tabular form, facilitating smooth and
comprehensive accessibility and manipulation for in-depth analysis.

##  Appendix D VAE Architecture

The convolutional VAE follows the encoder-decoder architecture. The encoder
learns how to encode the discrete DNA sequence data into a continuous latent
representation. Simultaneously, the decoder learns how to decode the latent
representations into DNA sequence data in continuous space before it is
quantised with ArgMax. A detailed illustration of the encoder architecture is
in figure LABEL:figure:appendix_vae1.

The decoder employs an architecture in symmetry with the encoder. In the
decoder the same network layers are performed in reverse order and inverse
operations are used in place of forward operations. Such inverse operation
pairs include transposed convolution and convolution operation, upsampling and
maxpooling operation.

Figure 8: VAE Encoder Architecture

Two sections can be distinguished from the architecture: 1) the discrete data
encoder (in red and yellow) and 2) the CNN Blocks (in purple)

The first section of the encoder encodes the one-hot encoded input into a
2-dimensional image-like representation. From the initial 4 channels due to
one-hot encoding, the channel dimension is first increased to 64, then it
increases further with steps of exponentially increasing number of kernels up
to 128. This is implemented using Multi-kernel Convolution Blocks with
perturbable kernel parameters combined with residual connections. By including
a max-pooling layer in each layer, its length shrinks by powers of 2 to 128.
This effectively expands our 4x2048x1 data into 128x128x1, an image-like
surface when viewed along the channel and length. This stage of the encoder
aims to capture high-level positional features along the length dimension via
pooling, for each channel dimension of the DNA sequence by transforming from
1-dimensional features into a 2-dimensional surface.

The second section of the encoder uses a CNN architecture for learning latent
representations. There are three layers of the CNN block, each performing 3x3
kernel on the channel-length surface followed by BatchNorm2D and LeakyReLU
activation. The CNN block in the encoder learns to produce different feature
maps of size 16x16 from the 128x128 surface, which are concatenated into a
3-dimensional latent representation of shape 16x16x32, effectively increasing
the width dimension. The latent representation can be partitioned into two
blocks of size 16x16x16, representing the mean and log-variance of the learned
features. Assuming priors of diagonal Gaussian distributions, we can take
samples from the standard Gaussian distribution and use the reparameterisation
trick to optimise for the variational learning.

The Multi-kernel Convolution Block is the elementary building block in the
first section of the VAE network. The block is structured to first perform
normalisation and activation first, then conv1D operations are performed on a
list of kernel sizes, where the results are then concatenated. An example
Multi-kernel Convolution Block with the list of kernel sizes [1, 3, 5] is
illustrated in figure 9.

Figure 9: An Example of Multi-kernel Convolution Block

##  Appendix E UNet

The UNet model used consists of four down and four up blocks. Each block
consists of eight sequential ResNet blocks, followed by a single upsampling or
downsampling block to change the number of channels. The third down block and
the second up block are enriched with cross-attention layers, one applied
after each ResNet block.

Each ResNet block transforms an input tensor through operational layers
including normalisation, swish non-linearity, upsampling, two convolutional
layers and time-embedding projection, culminating in an output tensor formed
by summing the input tensor and the output of the second convolutional layer.

Attention layers, integral to the refinement of the feature representation,
are incorporated in the third down and second up blocks, each consisting of
eight 64-dimensional attention heads, executed by a scaled dot product
attention mechanism without dropout and bias in the projection layers.

The UNet input and output contain 16 channels, each containing 16 x 16 arrays.
The down blocks have channel dimensions of [256,256,512,512] and the up blocks
mirror the encoder blocks in reverse order. In the forward process, we use N =
1000 steps.

##  Appendix F Extra Results

Figure 10: Compare the TATA-Box Distribution of VAE, DiscDiff, and Real DNA
Sequences