aainstitutetext: Faculty of Physics, University of Tabriz, Tabriz 51666-16471,
Iran

# The Cosmological Impact of Brane-Chern-Simons Massive Gravity

S. Kazempour,111Corresponding author. a and A. R. Akbarieh
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we present a novel extension of massive gravity theory; the
Brane-Chern-Simons massive gravity theory. We explore the cosmological
implications of this theory by deriving the background equations and
demonstrating the existence of self-accelerating solutions. Interestingly, our
theory suggests the existence of self-accelerating mechanisms that originate
from an effective cosmological constant, leading to intriguing possibilities
for understanding the nature of cosmic acceleration. Furthermore, we perform a
tensor perturbation analysis to investigate the propagation of gravitational
waves in this framework. We derive the dispersion relation for gravitational
waves and study their behavior in the Friedmann-Lemaître-Robertson-Walker
cosmology within the context of Brane-Chern-Simons massive gravity. Utilizing
the latest Union2 type Ia supernovae dataset comprising 557 SNIa events, we
provide observational support for our theoretical framework, indicating that
the Brane-Chern-Simons massive gravity theory is consistent with cosmological
observations.

##  1 Introduction

The late-time accelerated expansion of the universe remains one of the most
intriguing puzzles in modern cosmology, prompting extensive research into the
underlying mechanisms driving this phenomenon. While the inclusion of
cosmological constant or dark energy components within the framework of
general relativity provides a plausible explanation Weinberg:1988cp ;
Peebles:2002gy ; Copeland:2006wr ; Carroll:2003st ; Cai:2009zp ; Bamba:2012cp
, alternative approaches have also gained significant attention.

Modified gravity theories offer a compelling avenue to explore, encompassing a
wide range of models that modify the geometric sector of general relativity
Nojiri:2010wj ; Clifton:2011jh ; CANTATA:2021ktz ; Ishak:2018his ;
Abdalla:2022yfr . Within this broad category, curvature-based gravity
theories, such as $f(R)$ gravity Starobinsky:1980te , $f(G)$ gravity
Nojiri:2005jg , $f(P)$ gravity Erices:2019mkd , Lovelock gravity
Lovelock:1971yv , and Horndeski/Galileon scalar-tensor theories
Horndeski:1974wa ; Deffayet:2009wt , have been extensively studied. Another
branch of modified gravity involves torsion-based theories, including $f(T)$
gravity Bengochea:2008gz ; Cai:2015emx , $f(T,T_{G})$ gravity Kofinas:2014owa
, $f(T,B)$ gravity Bahamonde:2015zma , and scalar-torsion theories Geng:2011aj
.

In this work, we focus on a particular subclass of gravitational modification
known as massive gravity, where the graviton is endowed with a mass
deRham:2010ik ; deRham:2010kj ; Hinterbichler:2011tt ; Hassan:2011hr ;
Hassan:2011zd . The concept of massive gravity was first introduced by Fierz
and Pauli in 1939, who formulated a unique Lorentz-invariant linear theory
without ghosts in a flat spacetime Fierz:1939ix . However, subsequent work by
van Dam, Veltman, and Zakharov revealed the presence of a discontinuity (known
as the van Dam-Veltman-Zakharov discontinuity) when taking the massless limit
of the theory vanDam:1970vg ; Zakharov:1970cc . Vainshtein proposed extending
the theory to the nonlinear level to address this issue Vainshtein:1972sx ,
but this led to the discovery of a ghost instability, dubbed the Boulware-
Deser ghost Boulware:1972yco .

Significant progress was made by de Rham, Gabadadze, and Tolley (dRGT) who
developed a fully nonlinear massive gravity theory without the Boulware-Deser
ghost in a certain decoupling limit deRham:2010ik ; deRham:2010kj . While dRGT
massive gravity successfully explains the accelerated expansion of the
universe in an open Friedmann-Lemaître-Robertson-Walker (FLRW) geometry, it
faces challenges in accommodating a homogeneous and isotropic universe
DeFelice:2012mx . Additionally, the theory suffers from a strong coupling
problem and a nonlinear ghost instability, causing scalar and vector
perturbations to vanish Gumrukcuoglu:2011zh .

To address these issues, the quasi-dilaton massive gravity theory was
introduced DAmico:2012hia ; Gannouji:2013rwa , with further developments in
subsequent works DeFelice:2013dua ; Mukohyama:2014rca ; Tamanini:2013xia ;
EmirGumrukcuoglu:2014uog ; Cai:2013lqa ; Kahniashvili:2014wua ; Cai:2014upa ;
Gumrukcuoglu:2016hic ; Gumrukcuoglu:2017ioy ; Gumrukcuoglu:2020utx ;
DeFelice:2021trp ; Akbarieh:2021vhv ; Akbarieh:2022ovn ; Kazempour:2022let ;
Kazempour:2022giy ; Kazempour:2022xzy . However, even in this extended
framework, instabilities persist in the perturbation analysis
Gumrukcuoglu:2013nza ; Haghani:2013eya ; DAmico:2013saf .

In this paper, we propose a novel extension of massive gravity, dubbed Brane-
Chern-Simons massive gravity, which exhibits a self-accelerating solution and
instability-free perturbation. Brane-Chern-Simons gravity is a theoretical
framework that synergizes brane theory and Chern-Simons theory, offering a
unique perspective on the nature of gravity and its interplay with other
forces Izaurieta:2009hz ; Gomez:2020agz . By combining the concepts of branes
in higher-dimensional spaces and the dynamics of gauge fields in odd-
dimensional spaces, Brane-Chern-Simons gravity suggests a topological origin
of gravity and provides a rich framework for exploring new physics beyond the
standard model Izaurieta:2009hz ; Concha:2016kdz ; Gomez:2011zzd ;
Gomez:2020agz . The other research related to the Chern-Simons Gravity can be
found in these references Nojiri:2019nar ; Odintsov:2022hxu ; Nojiri:2020pqr ;
Concha:2014vka .

The main objectives of this work are to demonstrate the existence of a self-
accelerating solution in Brane-Chern-Simons massive gravity which is related
to the effective cosmological constant and to perform a perturbation analysis.
In particular, we will derive the modified dispersion relation for
gravitational waves within this theory. Furthermore, to constrain the
parameters of the theory, we employ observational data from Type Ia
supernovae. By utilizing Bayesian statistical methods, specifically focusing
on minimizing $X^{2}$, we can rigorously compare the predictions of the theory
with these observational datasets.

The paper is structured as follows: In Section 2, we introduce the Brane-
Chern-Simons massive gravity framework, derive the background equations of
motion, and extract self-accelerating solutions. Section 3, delves into the
cosmological perturbation analysis, focusing on tensor perturbations and the
dispersion relation of gravitational waves. Section 4, involves testing the
solution of Brane-Chern-Simons massive gravity against the latest Union2 type
Ia supernovae (SNIa) dataset. Finally, in Section 5, we summarize our key
findings and discuss future research directions.

##  2 Brane-Chern-Simons Massive Gravity

In this section, we introduce the Brane-Chern-Simons massive gravity action,
and we discuss the evolution of a cosmological background. The action includes
Planck mass $M_{\rm Pl}$, the Ricci scalar $R$, the constant $\varepsilon$,
the potential $V(\varphi)$, the gravitational constant $\kappa$, a dynamical
metric $g_{\mu\nu}$ and its determinant $\sqrt{-g}$. The action is given by

$\displaystyle S=\frac{M_{\rm Pl}^{2}}{2}\int
d^{4}x\Bigg{\\{}\sqrt{-g}\bigg{[}R+\varepsilon RV(\varphi)-2\kappa
V(\varphi)+2{m}_{g}^{2}U(\mathcal{K})\bigg{]}\Bigg{\\}},$

also, $m_{g}$ is the mass of graviton. In the following, we introduce the
$U({\cal K})$ of this action separately. It is clear that the mass of the
graviton comes up with the potential $U$ which consists of three parts.

$U(\mathcal{K})=U_{2}+\alpha_{3}U_{3}+\alpha_{4}U_{4},$ (2)

where $\alpha_{3}$ and $\alpha_{4}$ are dimensionless free parameters of the
theory. $U_{i}$ ($i=2,3,4$) is given by,

$\displaystyle U_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}[\mathcal{K}]^{2}-[\mathcal{K}^{2}]\big{)},$
$\displaystyle U_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\big{(}[\mathcal{K}]^{3}-3[\mathcal{K}][\mathcal{K}^{2}]+2[\mathcal{K}^{3}]\big{)},$
$\displaystyle U_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{24}\big{(}[\mathcal{K}]^{4}-6[\mathcal{K}]^{2}[\mathcal{K}^{2}]+8[\mathcal{K}][\mathcal{K}^{3}]+3[\mathcal{K}^{2}]^{2}-6[\mathcal{K}^{4}]\big{)},$
(3)

where the quantity "$[\cdot]$” is interpreted as the trace of the tensor
inside brackets. It is essential to mention that the building block tensor
$\mathcal{K}$ is defined as

$\mathcal{K}^{\mu}_{\nu}=\delta^{\mu}_{\nu}-\sqrt{g^{\mu\alpha}f_{\alpha\nu}},$
(4)

where $f_{\alpha\nu}$ is the fiducial metric, which is defined through

$f_{\alpha\nu}=\partial_{\alpha}\phi^{c}\partial_{\nu}\phi^{d}\eta_{cd}.$ (5)

Here $g^{\mu\nu}$ is the physical metric, $\eta_{cd}$ is the Minkowski metric
with $c,d=0,1,2,3$ and $\phi^{c}$ are the Stueckelberg fields which are
introduced to restore general covariance. According to our cosmological
application purpose, we adopt the Friedman-Lemaître-Robertson-Walker (FLRW)
Universe. So, the general expression of the corresponding dynamical and
fiducial metrics are given as follows,

$\displaystyle g_{\mu\nu}$ $\displaystyle={\rm
diag}\left[-N^{2},a^{2},a^{2},a^{2}\right],$ (6) $\displaystyle f_{\mu\nu}$
$\displaystyle={\rm diag}\left[-\dot{f}(t)^{2},1,1,1\right].$ (7)

Here it is worth pointing out that $N$ is the lapse function of the dynamical
metric, and it is similar to a gauge function. Also, it is clear that the
scale factor is represented by $a$, and $\dot{a}$ is the derivative with
respect to time. Furthermore, the lapse function relates the coordinate-time
$dt$ and the proper-time $d\tau$ via $d\tau=Ndt$ Scheel:1994yr ;
Christodoulakis:2013xha . Function $f(t)$ is the Stueckelberg scalar function
whereas $\phi^{0}=f(t)$ and $\frac{\partial\phi^{0}}{\partial t}=\dot{f}(t)$
Arkani-Hamed:2002bjr . Therefore, the point-like Lagrangian of the Brane-
Chern-Simons massive gravity in FLRW cosmology is

$\displaystyle\mathcal{L}=-\frac{3M_{\rm Pl}^{2}\big{[}1+\varepsilon
V(\varphi)\big{]}\dot{a}^{2}}{N^{2}a^{2}}-\kappa M_{\rm
Pl}^{2}V(\varphi)+\frac{m_{g}^{2}(a-1)M_{\rm
Pl}^{2}}{a^{3}N}\Bigg{\\{}\bigg{(}\alpha_{3}+\alpha_{4}+a\big{[}a(\alpha_{4}+4\alpha_{3}$
$\displaystyle+6)-2\alpha_{4}-5\alpha_{3}-3\big{]}\bigg{)}N-\bigg{(}a\big{[}a(\alpha_{4}+3\alpha_{3}+3)-2\alpha_{4}-3\alpha_{3}\big{]}+\alpha_{4}\bigg{)}\dot{f}(t)\Bigg{\\}}.$

In order to simplify expressions later, we define

$H\equiv\frac{\dot{a}}{Na}.$ (9)

###  2.1 Background Equations of Motion

In order to achieve a constraint equation we should take the unitary gauge
into consideration, which means that we choose $f(t)=t$. The significance of
the unitary gauge lies in the fact that on the classical level, the unphysical
fields could be eliminated from the Lagrangian with the use of gauge
transformations Grosse-Knetter:1992tbp .

In this procedure, a constraint equation can be derived by varying with
respect to $f$. So, the equation is given by

$\displaystyle\frac{\delta\mathcal{L}}{\delta f}=3m_{g}^{2}M_{\rm
Pl}^{2}\dot{a}\bigg{\\{}\alpha_{3}+\alpha_{4}-2a(1+2\alpha_{3}+\alpha_{4})+a^{2}(3+3\alpha_{3}+\alpha_{4})\bigg{\\}}=0.$
(10)

In this stage, the Friedman equation is achieved by varying with respect to
the lapse $N$,

$\displaystyle\frac{1}{M_{\rm Pl}^{2}a^{3}}\frac{\delta\mathcal{L}}{\delta
N}=3H^{2}(1+\varepsilon V(\varphi))-\kappa
V(\varphi)+\frac{m_{g}^{2}(a-1)}{a^{3}}\bigg{[}\alpha_{3}+\alpha_{4}+a\big{(}-3$
$\displaystyle-5\alpha_{3}-2\alpha_{4}+a(6+4\alpha_{3}+\alpha_{4})\big{)}\bigg{]}=0.$
(11)

The equation of motion for $a$ is

$\displaystyle\frac{1}{3M_{\rm Pl}^{2}a^{2}N}\frac{\delta\mathcal{L}}{\delta
a}=3H^{2}\big{[}1+\varepsilon V(\varphi)\big{]}-\kappa
V(\varphi)+\frac{2\dot{H}(1+\varepsilon
V(\varphi))}{N}+\frac{m_{g}}{r^{2}N^{3}}\Bigg{\\{}N\bigg{[}(1+2r)\big{(}1$
$\displaystyle+2\alpha_{3}+\alpha_{4}\big{)}+rN\big{(}rN(6+4\alpha_{3}+\alpha_{4})-(2+r)(3+3\alpha_{3}+\alpha_{4})\big{)}\bigg{]}-\alpha_{3}-\alpha_{4}\Bigg{\\}}=0,$

where

$\displaystyle r\equiv\frac{a}{N}.$ (13)

In the last part of this subsection, it should be noted that the Stuckelberg
field $f$ introduces time reparametrization invariance. So, there is a Bianchi
identity which relates the four equations of motion,

$\displaystyle\frac{\delta S}{\delta\sigma}\dot{\sigma}+\frac{\delta S}{\delta
f}\dot{f}-N\frac{d}{dt}\frac{\delta S}{\delta N}+\dot{a}\frac{\delta S}{\delta
a}=0.$ (14)

###  2.2 Self-Accelerating Background Solutions

In this step, we want to discuss solutions. It could be started with the
Stueckelberg constraint in Eq. (10), so, we have

$\displaystyle\alpha_{3}+\alpha_{4}-2a(1+2\alpha_{3}+\alpha_{4})+a^{2}(3+3\alpha_{3}+\alpha_{4})=0,$

it would be suitable to mention that the constant solutions of $a$ lead to the
effective energy density and behave similarly to a cosmological constant. As a
result, the two solutions of Eq. (2.2) are

$a_{\rm
SA}^{\pm}=\frac{1+2\alpha_{3}\pm\sqrt{1+\alpha_{3}+\alpha_{3}^{2}-\alpha_{4}}+\alpha_{4}}{3+3\alpha_{3}+\alpha_{4}}.$
(16)

The Friedman equation (2.1) could be written in a different form,

$\displaystyle 3H^{2}\bigg{(}1+\varepsilon V(\varphi)\bigg{)}=\Lambda_{\rm
SA}^{\pm}+\kappa V(\varphi).$ (17)

Considering self-accelerating solutions, a condition on the parameter
$\varepsilon V(\varphi)$ is provided by the Friedman equation (17). So, we
need to consider $\varepsilon V(\varphi)>-1$ to keep the left-hand side of Eq.
(17) positive. The importance of this the issue lies in the fact that when we
add ordinary matters to the right-hand side, throughout the matter-dominated
era, we will have the standard cosmology. It is worth mentioning that the
effective cosmological constant from the mass term is

$\displaystyle\Lambda_{\rm SA}^{\pm}\equiv\frac{m_{g}^{2}(1-a_{\rm
SA}^{\pm})}{a_{\rm SA}^{3\pm}}\bigg{[}\alpha_{3}+\alpha_{4}+a_{\rm
SA}^{\pm}\big{(}-3-5\alpha_{3}-2\alpha_{4}+a_{\rm
SA}^{\pm}(6+4\alpha_{3}+\alpha_{4})\big{)}\bigg{]}.$ (18)

According to Eq. (16), the above equation can be written as

$\displaystyle\Lambda_{\rm
SA}^{\pm}=\frac{m_{g}^{2}}{(\alpha_{3}+\alpha_{4})^{2}}\Bigg{\\{}2\sqrt{1+\alpha_{3}+\alpha_{3}^{2}-\alpha_{4}}\bigg{(}\pm
1\mp\alpha_{4}\pm\alpha_{3}\pm\alpha_{3}^{2}\bigg{)}+3\alpha_{4}(1+\alpha_{3})$
$\displaystyle-3\alpha_{3}(\alpha_{3}+1)-2\alpha_{3}^{3}-2\Bigg{\\}}.$ (19)

Therefore, from Eq. (2.1) we have,

$\displaystyle r_{\rm SA}=\frac{m_{g}^{2}N\bigg{[}a_{\rm
SA}^{2\pm}\big{(}\alpha_{3}+2\big{)}+N\big{(}2\alpha_{3}+3\big{)}-2a_{\rm
SA}^{\pm}N\big{(}\alpha_{3}+2\big{)}-\alpha_{3}-1\bigg{]}}{\big{(}a_{\rm
SA}^{\pm}-1\big{)}^{2}}\pm\Bigg{\\{}\frac{m_{g}^{2}N^{2}}{\big{(}a_{\rm
SA}^{\pm}-1\big{)}^{2}}\bigg{[}m_{g}^{2}\big{(}\alpha_{3}$
$\displaystyle+1\big{)}^{2}+\bigg{(}m_{g}^{2}\big{(}\alpha_{3}(\alpha_{3}+3)+3\big{)}+\kappa
V(\varphi)\big{(}\alpha_{3}+1\big{)}-3H^{2}(\alpha_{3}+1)\big{(}\varepsilon
V(\varphi)+1\big{)}\bigg{)}N^{2}-N\big{(}\alpha_{3}$
$\displaystyle+1\big{)}\bigg{(}m_{g}^{2}\big{(}2\alpha_{3}+3\big{)}+2\dot{H}\big{(}\varepsilon
V(\varphi)+1\big{)}\bigg{)}+a_{\rm
SA}^{2\pm}\bigg{(}m_{g}^{2}\big{(}\alpha_{3}+2\big{)}^{2}+\big{[}m_{g}^{2}(\alpha_{3}+3)-\kappa
V(\varphi)$ $\displaystyle+3H^{2}(\varepsilon
V(\varphi)+1)\big{]}\big{(}N(\alpha_{3}+2)-2\alpha_{3}-3\big{)}N-6\dot{H}+2\big{[}\varepsilon
V(\varphi)\big{(}(\alpha_{3}+2)N-2\alpha_{3}-3\big{)}$
$\displaystyle+N(\alpha_{3}+2)-2\alpha_{3}\big{]}\dot{H}\bigg{)}+2a_{\rm
SA}^{\pm}\bigg{(}(\alpha_{3}+1)\big{[}m_{g}^{2}\big{(}2\alpha_{3}+5\big{)}-\kappa
V(\varphi)+3H^{2}\big{(}\varepsilon V(\varphi)+1\big{)}\big{]}N$
$\displaystyle-
m_{g}^{2}N^{2}\big{(}\alpha_{3}+2\big{)}^{2}+(\alpha_{3}+1)\big{(}2\dot{H}(\varepsilon
V(\varphi)+1)-m_{g}^{2}(\alpha_{3}+2)\big{)}\bigg{)}\bigg{]}\Bigg{\\}}^{\frac{1}{2}}$
$\displaystyle\times\Bigg{\\{}\frac{m_{g}^{2}N^{2}\bigg{[}2a_{\rm
SA}^{\pm}\big{(}\alpha_{3}+2\big{)}-2\alpha_{3}-3\bigg{]}}{\big{(}a_{\rm
SA}^{\pm}-1\big{)}^{2}}+\bigg{(}\frac{m_{g}^{2}\big{[}3(\alpha_{3}+2)-2(2\alpha_{3}+5)a_{\rm
SA}^{\pm}+(\alpha_{3}+3)a_{\rm SA}^{2\pm}\big{]}}{\big{(}a_{\rm
SA}^{\pm}-1\big{)}^{2}}$ $\displaystyle-\kappa
V(\varphi)+3H^{2}\big{(}\varepsilon
V(\varphi)+1\big{)}\bigg{)}N^{3}+2N^{2}\dot{H}\big{(}\varepsilon
V(\varphi)+1\big{)}\Bigg{\\}}^{-1}.$

Actually, we have used the Stuckelberg equation (2.2) in order to eliminate
$\alpha_{4}$.

A key finding of this subsection is the identification of self-accelerating
solutions within the theory, characterized by an effective cosmological
constant, $\Lambda_{\rm SA}^{\pm}$. Importantly, these solutions do not
exhibit strong coupling issues, and they provide a well-behaved description of
the accelerated expansion of the universe. This suggests that the Brane-Chern-
Simons massive gravity theory offers a viable framework for explaining the
late-time cosmic acceleration without encountering the challenges typically
associated with strong coupling regimes.

##  3 Perturbation Analysis

In this section, we would like to analyze tensor perturbation in order to
calculate the mass of graviton for our theory which we introduced in the
previous section. Furthermore, we are trying to show the stability condition
of the system. In order to find the action for quadratic perturbation, the
physical metric is expanded in small fluctuation, $\delta g_{\mu\nu}$, around
a solution $g_{\mu\nu}^{(0)}$,

$g_{\mu\nu}=g_{\mu\nu}^{(0)}+\delta g_{\mu\nu}.$ (21)

In the following analysis, we keep terms to quadratic order in $\delta
g_{\mu\nu}$. As we demonstrate all analysis in the unitary gauge, there are
not any problems concerning the form of gauge invariant combinations.
Moreover, we write the actions expanded in the Fourier domain with plane
waves, i.e., $\vec{\nabla}^{2}\rightarrow-k^{2}$, $d^{3}x\rightarrow d^{3}k$.
We raise and lower the spatial indices on perturbations by $\delta^{ij}$ and
$\delta_{ij}$. We start by considering tensor perturbations around the
background,

$\delta g_{ij}=a^{2}h_{ij}^{\rm TT},$ (22)

where

$\partial^{i}h_{ij}=0\quad{\rm and}\quad g^{ij}h_{ij}=0.$ (23)

The tensor perturbed action in the second order can be calculated for each
part of the action separately. We write the Brane-Chern-Simons part of the
perturbed action in quadratic order

$\displaystyle S^{(2)}_{\rm B-C-S}=\frac{M_{\rm Pl}^{2}}{8}\int
d^{3}k\,dt\,a^{3}N\Bigg{\\{}\big{[}1+\varepsilon
V(\varphi)\big{]}\bigg{[}\frac{\dot{h}_{ij}\dot{h}^{ij}}{N^{2}}-(\frac{k^{2}}{a^{2}}+\frac{4\dot{H}}{N}+6H^{2})h^{ij}h_{ij}\bigg{]}$
$\displaystyle+2\kappa V(\varphi)h^{ij}h_{ij}\Bigg{\\}}.$

The second-order piece of the massive gravity sector of the perturbed action
can be written as

$\displaystyle S^{(2)}_{\rm massive}=\frac{M_{\rm Pl}^{2}}{8}\int
d^{3}k\,dt\,a^{3}Nm_{g}^{2}\Bigg{\\{}\bigg{[}\alpha_{3}+\alpha_{4}-3\big{(}1+2\alpha_{3}+\alpha_{4}\big{)}a+2\big{(}3+3\alpha_{3}+\alpha_{4}\big{)}a^{2}\bigg{]}\frac{\dot{f}(t)}{a^{2}N}$
$\displaystyle-\bigg{[}1+2\alpha_{3}+\alpha_{4}-3\big{(}3+3\alpha_{3}+\alpha_{4}\big{)}a+2\big{(}6+4\alpha_{3}+\alpha_{4}\big{)}a^{2}\bigg{]}\frac{1}{a^{2}}\Bigg{\\}}h^{ij}h_{ij}.$

Summing up the second order pieces of the perturbed actions $S^{(2)}_{\rm
B-C-S}$, and $S^{(2)}_{\rm massive}$ we obtain the total action in a second
order for tensor perturbations

$\displaystyle S^{(2)}_{\rm total}=\frac{M_{\rm Pl}^{2}}{8}\int
d^{3}k\,dt\,a^{3}N\bigg{\\{}\big{[}1+\varepsilon
V(\varphi)\big{]}\frac{\dot{h}^{ij}\dot{h}_{ij}}{N^{2}}-\Big{(}\big{[}1+\varepsilon
V(\varphi)\big{]}\frac{k^{2}}{a^{2}}+M_{\rm
GW}^{2}\Big{)}h^{ij}h_{ij}\bigg{\\}}.$

So, the dispersion relation of gravitational wave is

$\displaystyle M^{2}_{\rm
GW}=\big{(}4\frac{\dot{H}}{N}+6H^{2}\big{)}\big{[}1+\varepsilon
V(\varphi)\big{]}-2\kappa V(\varphi)+\Omega,$ (27)

where

$\displaystyle\Omega=\Bigg{\\{}2Na_{\rm
SA}^{3\pm}\big{[}\alpha_{3}+3\big{]}-N\big{(}\alpha_{3}+1\big{)}+a_{\rm
SA}^{2\pm}\bigg{[}\big{(}\alpha_{3}+2\big{)}\dot{f}(t)-2N\big{(}3\alpha_{3}+7\big{)}\bigg{]}$
$\displaystyle+a_{\rm
SA}^{\pm}\bigg{[}N\big{(}5\alpha_{3}+8\big{)}-\dot{f}(t)\big{(}\alpha_{3}+1\big{)}\bigg{]}\Bigg{\\}}\times\bigg{(}Na_{\rm
SA}^{2\pm}\big{[}a_{\rm SA}^{2\pm}-1\big{]}\bigg{)}^{-1}.$ (28)

It is important to note that we eliminate $\alpha_{4}$ using Eq. (16). The
sign of the mass square of gravitational waves plays a crucial role in
ensuring the stability of long-wavelength gravitational waves. When the mass
square is positive, it indicates a stable configuration. Conversely, a
negative mass square implies the presence of tachyons, leading to potential
instability. However, due to the characteristic timescale associated with the
tachyonic mass, which is on the order of the Hubble scale, the development of
any instability would take a duration comparable to the age of the universe,
thus having minimal impact on the stability of the system over observable
timescales.

##  4 Cosmological Tests

The exploration of Type Ia supernovae provided compelling evidence for the
accelerated expansion of the universe Copeland:2006wr ; Frieman:2008sn ;
Perlmutter:2003kf ; Yang:2019fjt . Here, we utilize the Union2 supernovae Ia
dataset, comprising 557 SNIa events Amanullah:2010vv , to investigate the
Brane-Chern-Simons Massive Gravity theory. The observational data from the
SNIa dataset will be presented in terms $\mu_{\rm obs}$, enabling a direct
comparison with the theoretical predictions of the model under scrutiny.

$\displaystyle\mu_{\rm th}(z_{i})=5\log_{10}D_{L}(z_{i})+\mu_{0},$ (29)

it is important to note that $\mu_{0}=42.38-5\log_{10}\mathcal{H}$, where
$\mathcal{H}$ represents the Hubble constant $H_{0}$ in units of $100\,{\rm
km/s/Mpc}$. This relationship is essential for establishing the connection
between the observed magnitude and the luminosity distance, which can be
expressed as follows:

$\displaystyle D_{L}(z)=(1+z)\int_{0}^{z}\frac{dx}{E(n;q)},$ (30)

furthermore, it is evident that $E=H/H_{0}$ and $q$ are the model parameters.
It is important to note that $X^{2}$ is defined as follows:

$X_{\mu}^{2}(q)=\sum_{i}\frac{[\mu_{\rm obs}(z_{i})-\mu_{\rm
th}(z_{i})]^{2}}{\sigma^{2}(z_{i})},$ (31)

here, $\sigma$ corresponds to the $1\sigma$ error, and the parameter $\mu_{0}$
is a nuisance parameter independent of the data points. To minimize
$X_{\mu}^{2}$ in equation (31), we expand it with respect $\mu_{0}$
Nesseris:2005ur ; DiPietro:2002cz .

$X_{\mu}^{2}(q)=\tilde{T}-2\mu_{0}\tilde{Y}+\mu_{0}^{2}\tilde{Z},$ (32)

where

$\displaystyle\tilde{T}(q)=\sum_{i}\frac{[\mu_{\rm obs}(z_{i})-\mu_{\rm
th}(z_{i};\mu_{0}=0,q)]^{2}}{\sigma_{\mu_{\rm obs}}^{2}(z_{i})},$
$\displaystyle\tilde{Y}(q)=\sum_{i}\frac{\mu_{\rm obs}(z_{i})-\mu_{\rm
th}(z_{i};\mu_{0}=0,q)}{\sigma_{\mu_{\rm obs}}^{2}(z_{i})},$
$\displaystyle\tilde{Z}=\sum_{i}\frac{1}{\sigma_{\mu_{\rm obs}}^{2}(z_{i})}.$
(33)

It should be explained that for $\mu_{0}=\frac{\tilde{Y}}{\tilde{Z}}$,
equation (32) has a minimum at

$\tilde{X}_{\mu}^{2}(q)=\tilde{T}(q)-\frac{\tilde{Y}^{2}(q)}{\tilde{Z}}.$ (34)

As it is clear that $X_{\mu,{\rm min}}^{2}=\tilde{X}_{\mu,{\rm min}}^{2}$ we
can consider minimizing $\tilde{X}_{\mu}^{2}$, which is independent of
$\mu_{0}$. It is important to note that the best-fit model parameter is
determined by minimizing $X^{2}=\tilde{X}_{\mu}^{2}$. At the same time, we
know that the corresponding $\mathcal{H}$ is dependent on
$\mu_{0}=\frac{\tilde{Y}}{\tilde{Z}}$ for the best-fit parameter. By
considering equation (17) and applying the change of variables
$a=\frac{1}{1+z}$ and $\frac{d}{dt}=-H(z+1)\frac{d}{dz}$, we can express the
dimensionless Hubble parameter for this case. Consequently, the solution to
the asymptotic state at small redshifts should be provided as follows:

$\displaystyle H(z)\sim
H_{0}+\frac{\Theta^{\pm}}{H_{0}^{3}}z-\frac{3\Theta^{\pm\rm
2}+\Theta^{\pm}H_{0}^{4}}{2H_{0}^{7}}z^{2}+.....,$

where $\Theta^{\pm}$ is defined as

$\displaystyle\Theta^{\pm}=-\frac{2}{3}\frac{\Lambda_{\rm SA}^{\pm}+\kappa
V(\varphi)}{\big{(}1+\varepsilon V(\varphi)\big{)}},$ (36)

it is important to note that $H_{0}$ represents the Hubble parameter at the
present time. Therefore, the dimensionless Hubble parameter can be expressed
as follows:

$\displaystyle E=\frac{H(z)}{H_{0}}=1+Fz-(\frac{3F^{2}+F}{2})z^{2}+......$

where

$\displaystyle F=\frac{\Theta^{\pm}}{H_{0}^{4}}.$ (38)

In this step, we present plots of $X^{2}$ and likelihoods as functions of the
parameter $F$. Based on our calculations, the best-fit yields a value of
$X_{\rm min}^{2}=544.044$, and the corresponding best-fit parameter is,

$\displaystyle F=$ $\displaystyle
0.5497^{+0.0251}_{-0.0255},\quad\quad\mbox{with 1$\sigma$ uncertainty},$ (39)
$\displaystyle F=$ $\displaystyle
0.5497^{+0.0500}_{-0.0515},\quad\quad\mbox{with 2$\sigma$ uncertainty}.$ (40)

Consequently, it is noteworthy to mention that the best-fit value for the
Brane-Chern-Simons massive gravity theory corresponds to $\mathcal{H}=0.6986$.

Figure 1: The distance modulus diagram for the best fit (red solid line) of
the parameter of the Brane-Chern-Simons massive gravity theory in comparison
with the 557 Union2 SNIa data points (black dots). Figure 2: The $X^{2}$
functions of parameter $F$ for the Brane-Chern-Simons massive gravity theory.
Figure 3: The likelihood as functions of parameter $F$ for the Brane-Chern-
Simons massive gravity theory.

According to Figures 1 and 3, the result of fitting the Brane-Chern-Simons
massive gravity theory to the cosmological data provides us with the optimal
value for the model parameter.

Figure 4: The best-fit diagram for the luminosity distance $D_{L}$ with
respect to velocity $V$ for the Brane-Chern-Simons massive gravity theory.

To evaluate the validity of this theory, we utilize the Union2 supernovae Ia
dataset, consisting of $557$ SNIa events. By employing the Bayesian statistics
method, we constrain the model parameters based on the SNIa data.
Specifically, we generate plots of $X^{2}$ and the likelihood function as
functions of the parameter $F$, allowing us to determine the best-fit value by
minimizing $X^{2}$.

Our analysis yields a best-fit value of $F=0.5497$ corresponding to a minimum
$X^{2}$ value of $544.044$. Using this optimal parameter value, we calculate
the best-fit value of $\mathcal{H}$ for the theory, obtaining
$\mathcal{H}=0.6986$.

Figure 1 illustrates the distance modulus diagram, depicting the relationship
between the distance modulus $\mu$ and the redshift $z$, along with the best-
fit curve for the Brane-Chern-Simons massive gravity theory based on
supernovae Ia data. Additionally, we provide a plot of velocity $V$ versus
luminosity distance $D_{L}$, demonstrating the agreement between theoretical
predictions and supernovae Ia observations. In Figure 4, we present a plot of
the velocity $V$ as a function of luminosity distance $D_{L}$ for the Brane-
Chern-Simons massive gravity theory. The velocity is calculated from the
redshift using the equation $V=c-\frac{2c}{(z+1)^{2}+1}$.

##  5 Conclusion

In this work, we introduce a novel extension of massive gravity theory, the
Brane-Chern-Simons massive gravity theory. We explore its potential to explain
the late-time acceleration of the universe by deriving the full set of
equations of motion for an FLRW background and analyzing self-accelerating
background solutions. A key finding is the identification of self-accelerating
solutions characterized by an effective cosmological constant, $\Lambda_{\rm
SA}^{\pm}$. These solutions exhibit well-behaved behavior without strong
coupling issues, offering a promising framework to explain the accelerated
expansion of the universe.

Furthermore, our analysis of tensor perturbations provides insights into the
nature of the graviton mass within this theory, leading to the derivation of
the dispersion relation for gravitational waves. We also examine the stability
of long-wavelength gravitational waves, where the sign of the mass square
plays a crucial role. A positive mass square indicates stability, while a
negative value suggests the presence of tachyons. However, the development of
any instability due to tachyonic behavior is mitigated by the long timescales
associated with the Hubble scale.

The validation of our model is demonstrated in the final section through a
comparison with the latest Union2 SNIa dataset, consisting of 557 supernovae
events. The Brane-Chern-Simons massive gravity theory exhibits an excellent
fit with the observational data, as illustrated in Fig. 1. This allows us to
constrain the model parameters and obtain a best-fit value of
$\mathcal{H}=0.6986$ for the Hubble parameter, reinforcing the theory’s
capability to describe the late-time acceleration of the universe. Future
studies, utilizing larger and more diverse datasets, will be crucial for
further refining and validating these findings, ultimately advancing our
understanding of the fundamental principles governing the universe’s behavior.

###### Acknowledgements.

This work is based upon research funded by the University of Tabriz, Iran
National Science Foundation (INSF), and Iran National Elites Foundation
(INEF), under project No. 4014244. The authors thank A. Emir Gumrukcuolu for
insightful comments. We also appreciate Nishant Agarwal for sharing notes and
computational codes related to tensor perturbations, which greatly facilitated
our research.

#### Note added.

This is also a good position for notes added after the paper has been written.

## References

  * (1) S. Weinberg, Rev. Mod. Phys. 61 (1989), 1-23 doi:10.1103/RevModPhys.61.1
  * (2) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75 (2003), 559-606 doi:10.1103/RevModPhys.75.559 [arXiv:astro-ph/0207347 [astro-ph]].
  * (3) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15 (2006), 1753-1936 doi:10.1142/S021827180600942X [arXiv:hep-th/0603057 [hep-th]].
  * (4) S. M. Carroll, M. Hoffman and M. Trodden, Phys. Rev. D 68 (2003), 023509 doi:10.1103/PhysRevD.68.023509 [arXiv:astro-ph/0301273 [astro-ph]].
  * (5) Y. F. Cai, E. N. Saridakis, M. R. Setare and J. Q. Xia, Phys. Rept. 493 (2010), 1-60 doi:10.1016/j.physrep.2010.04.001 [arXiv:0909.2776 [hep-th]].
  * (6) K. Bamba, S. Capozziello, S. Nojiri and S. D. Odintsov, Astrophys. Space Sci. 342 (2012), 155-228 doi:10.1007/s10509-012-1181-8 [arXiv:1205.3421 [gr-qc]].
  * (7) S. Nojiri and S. D. Odintsov, Phys. Rept. 505 (2011), 59-144 doi:10.1016/j.physrep.2011.04.001 [arXiv:1011.0544 [gr-qc]].
  * (8) T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, Phys. Rept. 513 (2012), 1-189 doi:10.1016/j.physrep.2012.01.001 [arXiv:1106.2476 [astro-ph.CO]].
  * (9) E. N. Saridakis et al. [CANTATA], Springer, 2021, ISBN 978-3-030-83714-3, 978-3-030-83717-4, 978-3-030-83715-0 doi:10.1007/978-3-030-83715-0 [arXiv:2105.12582 [gr-qc]].
  * (10) M. Ishak, Living Rev. Rel. 22 (2019) no.1, 1 doi:10.1007/s41114-018-0017-4 [arXiv:1806.10122 [astro-ph.CO]].
  * (11) E. Abdalla, G. Franco Abellán, A. Aboubrahim, A. Agnello, O. Akarsu, Y. Akrami, G. Alestas, D. Aloni, L. Amendola and L. A. Anchordoqui, et al. JHEAp 34 (2022), 49-211 doi:10.1016/j.jheap.2022.04.002 [arXiv:2203.06142 [astro-ph.CO]].
  * (12) A. A. Starobinsky, Phys. Lett. B 91 (1980), 99-102 doi:10.1016/0370-2693(80)90670-X
  * (13) S. Nojiri and S. D. Odintsov, Phys. Lett. B 631 (2005), 1-6 doi:10.1016/j.physletb.2005.10.010 [arXiv:hep-th/0508049 [hep-th]].
  * (14) C. Erices, E. Papantonopoulos and E. N. Saridakis, Phys. Rev. D 99 (2019) no.12, 123527 doi:10.1103/PhysRevD.99.123527 [arXiv:1903.11128 [gr-qc]].
  * (15) D. Lovelock, J. Math. Phys. 12 (1971), 498-501 doi:10.1063/1.1665613
  * (16) G. W. Horndeski, Int. J. Theor. Phys. 10 (1974), 363-384 doi:10.1007/BF01807638
  * (17) C. Deffayet, G. Esposito-Farese and A. Vikman, Phys. Rev. D 79 (2009), 084003 doi:10.1103/PhysRevD.79.084003 [arXiv:0901.1314 [hep-th]].
  * (18) G. R. Bengochea and R. Ferraro, Phys. Rev. D 79 (2009), 124019 doi:10.1103/PhysRevD.79.124019 [arXiv:0812.1205 [astro-ph]].
  * (19) Y. F. Cai, S. Capozziello, M. De Laurentis and E. N. Saridakis, Rept. Prog. Phys. 79 (2016) no.10, 106901 doi:10.1088/0034-4885/79/10/106901 [arXiv:1511.07586 [gr-qc]].
  * (20) G. Kofinas and E. N. Saridakis, Phys. Rev. D 90 (2014), 084044 doi:10.1103/PhysRevD.90.084044 [arXiv:1404.2249 [gr-qc]].
  * (21) S. Bahamonde, C. G. Böhmer and M. Wright, Phys. Rev. D 92 (2015) no.10, 104042 doi:10.1103/PhysRevD.92.104042 [arXiv:1508.05120 [gr-qc]].
  * (22) C. Q. Geng, C. C. Lee, E. N. Saridakis and Y. P. Wu, Phys. Lett. B 704 (2011), 384-387 doi:10.1016/j.physletb.2011.09.082 [arXiv:1109.1092 [hep-th]].
  * (23) K. Hinterbichler, Rev. Mod. Phys. 84 (2012), 671-710 doi:10.1103/RevModPhys.84.671 [arXiv:1105.3735 [hep-th]].
  * (24) S. F. Hassan and R. A. Rosen, Phys. Rev. Lett. 108 (2012), 041101 doi:10.1103/PhysRevLett.108.041101 [arXiv:1106.3344 [hep-th]].
  * (25) S. F. Hassan and R. A. Rosen, JHEP 02 (2012), 126 doi:10.1007/JHEP02(2012)126 [arXiv:1109.3515 [hep-th]].
  * (26) C. de Rham and G. Gabadadze, Phys. Rev. D 82 (2010), 044020 doi:10.1103/PhysRevD.82.044020 [arXiv:1007.0443 [hep-th]].
  * (27) C. de Rham, G. Gabadadze and A. J. Tolley, Phys. Rev. Lett. 106 (2011), 231101 doi:10.1103/PhysRevLett.106.231101 [arXiv:1011.1232 [hep-th]].
  * (28) M. Fierz and W. Pauli, Proc. Roy. Soc. Lond. A 173 (1939), 211-232 doi:10.1098/rspa.1939.0140
  * (29) H. van Dam and M. J. G. Veltman, Nucl. Phys. B 22 (1970), 397-411 doi:10.1016/0550-3213(70)90416-5
  * (30) V. I. Zakharov, JETP Lett. 12 (1970), 312
  * (31) A. I. Vainshtein, Phys. Lett. B 39 (1972), 393-394 doi:10.1016/0370-2693(72)90147-5
  * (32) D. G. Boulware and S. Deser, Phys. Rev. D 6 (1972), 3368-3382 doi:10.1103/PhysRevD.6.3368
  * (33) A. De Felice, A. E. Gumrukcuoglu and S. Mukohyama, Phys. Rev. Lett. 109 (2012), 171101 doi:10.1103/PhysRevLett.109.171101 [arXiv:1206.2080 [hep-th]].
  * (34) A. E. Gumrukcuoglu, C. Lin and S. Mukohyama, JCAP 03 (2012), 006 doi:10.1088/1475-7516/2012/03/006 [arXiv:1111.4107 [hep-th]].
  * (35) G. D’Amico, G. Gabadadze, L. Hui and D. Pirtskhalava, Phys. Rev. D 87 (2013), 064037 doi:10.1103/PhysRevD.87.064037 [arXiv:1206.4253 [hep-th]].
  * (36) R. Gannouji, M. W. Hossain, M. Sami and E. N. Saridakis, Phys. Rev. D 87 (2013), 123536 doi:10.1103/PhysRevD.87.123536 [arXiv:1304.5095 [gr-qc]].
  * (37) A. De Felice, A. Emir Gümrükçüoğlu and S. Mukohyama, Phys. Rev. D 88 (2013) no.12, 124006 doi:10.1103/PhysRevD.88.124006 [arXiv:1309.3162 [hep-th]].
  * (38) S. Mukohyama, JCAP 12 (2014), 011 doi:10.1088/1475-7516/2014/12/011 [arXiv:1410.1996 [hep-th]].
  * (39) N. Tamanini, E. N. Saridakis and T. S. Koivisto, JCAP 02 (2014), 015 doi:10.1088/1475-7516/2014/02/015 [arXiv:1307.5984 [hep-th]].
  * (40) A. Emir Gümrükçüoğlu, L. Heisenberg and S. Mukohyama, JCAP 02 (2015), 022 doi:10.1088/1475-7516/2015/02/022 [arXiv:1409.7260 [hep-th]].
  * (41) Y. F. Cai, F. Duplessis and E. N. Saridakis, Phys. Rev. D 90 (2014) no.6, 064051 doi:10.1103/PhysRevD.90.064051 [arXiv:1307.7150 [hep-th]].
  * (42) T. Kahniashvili, A. Kar, G. Lavrelashvili, N. Agarwal, L. Heisenberg and A. Kosowsky, Phys. Rev. D 91 (2015) no.4, 041301 [erratum: Phys. Rev. D 100 (2019) no.8, 089902] doi:10.1103/PhysRevD.91.041301 [arXiv:1412.4300 [astro-ph.CO]].
  * (43) Y. F. Cai and E. N. Saridakis, Phys. Rev. D 90 (2014) no.6, 063528 doi:10.1103/PhysRevD.90.063528 [arXiv:1401.4418 [astro-ph.CO]].
  * (44) A. E. Gumrukcuoglu, K. Koyama and S. Mukohyama, Phys. Rev. D 94 (2016) no.12, 123510 doi:10.1103/PhysRevD.94.123510 [arXiv:1610.03562 [hep-th]].
  * (45) A. E. Gumrukcuoglu, K. Koyama and S. Mukohyama, Phys. Rev. D 96 (2017) no.4, 044041 doi:10.1103/PhysRevD.96.044041 [arXiv:1707.02004 [hep-th]].
  * (46) A. E. Gumrukcuoglu, R. Kimura and K. Koyama, Phys. Rev. D 101 (2020) no.12, 124021 doi:10.1103/PhysRevD.101.124021 [arXiv:2003.11831 [gr-qc]].
  * (47) A. De Felice, S. Mukohyama and M. C. Pookkillath, JCAP 12 (2021) no.12, 011 doi:10.1088/1475-7516/2021/12/011 [arXiv:2110.01237 [astro-ph.CO]].
  * (48) A. R. Akbarieh, S. Kazempour and L. Shao, Phys. Rev. D 103 (2021), 123518 doi:10.1103/PhysRevD.103.123518 [arXiv:2105.03744 [gr-qc]].
  * (49) A. R. Akbarieh, S. Kazempour and L. Shao, Phys. Rev. D 105 (2022) no.2, 023501 doi:10.1103/PhysRevD.105.023501 [arXiv:2203.00901 [gr-qc]].
  * (50) S. Kazempour and A. R. Akbarieh, Phys. Rev. D 105 (2022) no.12, 123515 doi:10.1103/PhysRevD.105.123515 [arXiv:2204.05595 [gr-qc]].
  * (51) S. Kazempour, A. R. Akbarieh, H. Motavalli and L. Shao, Phys. Rev. D 106 (2022) no.2, 023508 doi:10.1103/PhysRevD.106.023508 [arXiv:2205.10863 [gr-qc]].
  * (52) S. Kazempour, A. R. Akbarieh and E. N. Saridakis, Phys. Rev. D 106 (2022) no.10, 103502 doi:10.1103/PhysRevD.106.103502 [arXiv:2207.13936 [gr-qc]].
  * (53) A. E. Gümrükçüoğlu, K. Hinterbichler, C. Lin, S. Mukohyama and M. Trodden, Phys. Rev. D 88 (2013) no.2, 024023 doi:10.1103/PhysRevD.88.024023 [arXiv:1304.0449 [hep-th]].
  * (54) Z. Haghani, H. R. Sepangi and S. Shahidi, Phys. Rev. D 87 (2013) no.12, 124014 doi:10.1103/PhysRevD.87.124014 [arXiv:1303.2843 [gr-qc]].
  * (55) G. D’Amico, G. Gabadadze, L. Hui and D. Pirtskhalava, Class. Quant. Grav. 30 (2013), 184005 doi:10.1088/0264-9381/30/18/184005 [arXiv:1304.0723 [hep-th]].
  * (56) F. Izaurieta, E. Rodriguez, P. Minning, P. Salgado and A. Perez, Phys. Lett. B 678 (2009), 213-217 doi:10.1016/j.physletb.2009.06.017 [arXiv:0905.2187 [hep-th]].
  * (57) F. Gómez, S. Lepe and P. Salgado, Eur. Phys. J. C 81 (2021) no.1, 9 doi:10.1140/epjc/s10052-020-08804-z [arXiv:2012.10578 [gr-qc]].
  * (58) P. K. Concha, R. Durka, C. Inostroza, N. Merino and E. K. Rodríguez, Phys. Rev. D 94 (2016) no.2, 024055 doi:10.1103/PhysRevD.94.024055 [arXiv:1603.09424 [hep-th]].
  * (59) F. Gomez, P. Minning and P. Salgado, Phys. Rev. D 84 (2011), 063506 doi:10.1103/PhysRevD.84.063506
  * (60) S. Nojiri, S. D. Odintsov, V. K. Oikonomou and A. A. Popov, Phys. Rev. D 100 (2019) no.8, 084009 doi:10.1103/PhysRevD.100.084009 [arXiv:1909.01324 [gr-qc]].
  * (61) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 105 (2022) no.10, 104054 doi:10.1103/PhysRevD.105.104054 [arXiv:2205.07304 [gr-qc]].
  * (62) S. Nojiri, S. D. Odintsov, V. K. Oikonomou and A. A. Popov, Phys. Dark Univ. 28 (2020), 100514 doi:10.1016/j.dark.2020.100514 [arXiv:2002.10402 [gr-qc]].
  * (63) P. K. Concha, D. M. Penafiel, E. K. Rodriguez and P. Salgado, Eur. Phys. J. C 74 (2014), 2741 doi:10.1140/epjc/s10052-014-2741-6 [arXiv:1402.0023 [hep-th]].
  * (64) M. A. Scheel, S. L. Shapiro and S. A. Teukolsky, Phys. Rev. D 51 (1995), 4208-4235 doi:10.1103/PhysRevD.51.4208 [arXiv:gr-qc/9411025 [gr-qc]].
  * (65) T. Christodoulakis, N. Dimakis and P. A. Terzis, J. Phys. A 47 (2014), 095202 doi:10.1088/1751-8113/47/9/095202 [arXiv:1304.4359 [gr-qc]].
  * (66) N. Arkani-Hamed, H. Georgi and M. D. Schwartz, Annals Phys. 305 (2003), 96-118 doi:10.1016/S0003-4916(03)00068-X [arXiv:hep-th/0210184 [hep-th]].
  * (67) C. Grosse-Knetter and R. Kogerler, Phys. Rev. D 48 (1993), 2865-2876 doi:10.1103/PhysRevD.48.2865 [arXiv:hep-ph/9212268 [hep-ph]].
  * (68) J. Frieman, M. Turner and D. Huterer, Ann. Rev. Astron. Astrophys. 46 (2008), 385-432 doi:10.1146/annurev.astro.46.060407.145243 [arXiv:0803.0982 [astro-ph]].
  * (69) S. Perlmutter and B. P. Schmidt, Lect. Notes Phys. 598 (2003), 195-217 doi:10.1007/3-540-45863-8_11 [arXiv:astro-ph/0303428 [astro-ph]].
  * (70) Y. Yang and Y. Gong, JCAP 06 (2020), 059 doi:10.1088/1475-7516/2020/06/059 [arXiv:1912.07375 [astro-ph.CO]].
  * (71) R. Amanullah, C. Lidman, D. Rubin, G. Aldering, P. Astier, K. Barbary, M. S. Burns, A. Conley, K. S. Dawson and S. E. Deustua, et al. Astrophys. J. 716 (2010), 712-738 doi:10.1088/0004-637X/716/1/712 [arXiv:1004.1711 [astro-ph.CO]].
  * (72) S. Nesseris and L. Perivolaropoulos, Phys. Rev. D 72 (2005), 123519 doi:10.1103/PhysRevD.72.123519 [arXiv:astro-ph/0511040 [astro-ph]].
  * (73) E. Di Pietro and J. F. Claeskens, Mon. Not. Roy. Astron. Soc. 341 (2003), 1299 doi:10.1046/j.1365-8711.2003.06508.x [arXiv:astro-ph/0207332 [astro-ph]].


Interior Point Differential Dynamic Programming

Andrei Pavlov, Student Member, IEEE,
Iman Shames, Member, IEEE,
and Chris Manzie, Senior Member, IEEE.
The authors are with the Department
of Electrical and Electronic Engineering, The University of Melbourne, Parkville,
Victoria, 3010 Australia, e-mail: {apavlov@student.,ishames<EMAIL_ADDRESS>
This paper introduces a novel Differential Dynamic Programming (DDP) algorithm for solving discrete-time finite-horizon optimal control problems with inequality constraints. Two variants, namely Feasible- and Infeasible-IPDDP algorithms, are developed using primal-dual interior-point methodology, and their local quadratic convergence properties are characterised. We show that the stationary points of the algorithms are the perturbed KKT points, and thus can be moved arbitrarily close to a locally optimal solution. Being free from the burden of the active-set methods, it can handle nonlinear state and input inequality constraints without a discernible increase in its computational complexity relative to the unconstrained case. The performance of the proposed algorithms is demonstrated using numerical experiments on three different problems: control-limited inverted pendulum, car-parking, and unicycle motion control and obstacle avoidance.

Finite Horizon Optimal Control, Differential Dynamic Programming, Interior Point Methods, Numerical Methods.

§ INTRODUCTION

Among optimisation algorithms that can effectively utilise the structure of the optimal control problems a distinct place is held by Differential Dynamic Programming (DDP) algorithm introduced by Mayne[1]. Its advantages include linear complexity in the length of prediction horizon and locally optimal feedback policies as an extra output. It has provable local quadratic convergence [2], while global convergence with inexact line-search can also be established[3]. DDP can address minimax-type problems [4] or be used for Model Predictive Control applications[5]. However, the DDP approach requires second-order derivatives of system dynamics, which makes it prohibitively expensive for solving large problems. As a possible remedy, one can omit the second-order information as in iterative-Linear-Quadratic Regulator (iLQR)[6], or estimate it using Quasi-Newton DDP[7], Sampled DDP[8] or Unscented DDP[9].

One major shortcoming of DDP is a lack of an elegant generalisation for the inequality-constrained problems. Methods reported in the literature usually fall into one of two distinct categories: penalty, barrier and Augmented Lagrangian (AL) methods; or active sets methods. The first family utilises penalty and barrier functions (or augments the Lagrangian function) to convert constrained problems into unconstrained, e.g., see [10][11]. Algorithms of the first family potentially suffer from major drawbacks, such as ill-conditioning, need for hand-tuning, slow convergence and/or appearance of saddle-points. The second family is based on active-set methodology, e.g., see [12, 13, 14, 15]. In contrast to the former family, these methods require the constraints to be explicitly dependent on the control variables, and include an extra routine aimed at identification of active/inactive constraints. This approach is known to have combinatorial complexity in the number of constraints in the worst case. To circumvent potential computational difficulties Tassa et. al.[16] propose Control-Limited Differential Dynamic Programming (CLDDP) algorithm where only box constraints on input are considered. A hybrid algorithm which uses ideas from both families was reported in [17]. To the best of our knowledge there are no rigorous proofs of convergence for constrained DDP algorithms for problems with nonlinear state and input constraints.

While general purpose primal-dual interior-point methods are particularly successful in practice (see [18] for details), there were no extension of DDP algorithm accommodating primal-dual interior-point techniques reported so far. Our aim is to fill the gap and bring the benefits of this approach. As demonstrated in the following, the primal-dual Interior-Point DDP algorithm seems to be one of the most natural extensions to DDP. It requires neither modifying the objective function nor identifying active/inactive constraints by a separate procedure. Importantly, it has provable local quadratic convergence for problems with nonlinear state and control constraints, which is a new result in the DDP framework.

The paper is organised as follows. In Section <ref>, we introduce the finite-time optimal control problem, and state a system of (perturbed) first-order conditions for local optimality. In Section <ref>, the main contribution of the paper, i.e., two extensions for the Differential Dynamic Programming (DDP) algorithm, namely Feasible- and Infeasible-Interior-Point DDP algorithms are described and their local quadratic convergence properties are established. In Section <ref> we consider numerical examples to demonstrate advantages in convergence relative to existing approaches (CLDDP and log-barrier DDP).

§ PRELIMINARIES

Consider a discrete-time system of the following form

\begin{equation*} \label{sys}
x_{t+1}= f(x_t,u_t),
\end{equation*}

subject to constraints $c(x_t,u_t)\leq 0$,
where $x_t\in \mathbb{R}^n$ and $u_t\in \mathbb{R}^m$ are state and control vectors at time $t$, functions $f:\mathbb{R}^n\times \mathbb{R}^m \rightarrow \mathbb{R}^n$ and $c:\mathbb{R}^n\times \mathbb{R}^m \rightarrow \mathbb{R}^l$ are twice continuously differentiable.

Consider a finite-time constrained optimal control problem for the above system at a known initial state $\bar x_0$:

\begin{equation} \label{CFTOCP}
\begin{aligned}
J^{\star}_N(\bar x_0)=&\min_{\bm x, \bm u} & &\sum_{t=0}^{N-1} q(x_t, u_t) + p(x_{N}) \\
&\text{\ s.t.} & &x_0=\bar x_0,\text{ and for }t\in\{0,\ldots,N-1\}:\\
& & &x_{t+1}=f(x_t, u_t),\ c(x_t, u_t)\leq 0,
\end{aligned}
\end{equation}

where $q:\mathbb{R}^n\times \mathbb{R}^m \rightarrow \mathbb{R}$ and $p:\mathbb{R}^n \rightarrow \mathbb{R}$ are twice continuously differentiable stage and terminal costs respectively,  $N\geq 1$ is the prediction horizon length, $\bm x$ and $\bm u$ are vectors of the corresponding decision variables, i.e., $\bm x=(x_0,\ldots,x_N)$ and $\bm u=(u_0,\ldots,u_{N-1})$.

The Lagrangian function for this optimisation problem is

\begin{multline*}\label{Lagrange}
L(\bm x, \bm u, \bm \lambda,  \bm s) = \sum_{t=0}^{N-1} q(x_t, u_t)+ s_t^T  c(x_t, u_t)+p(x_N)\\ +\lambda_0^T( \bar x_0-x_0)
+ \sum_{t=1}^{N}\lambda_{t}^T  (f(x_{t-1},u_{t-1})-x_{t}),
\end{multline*}

where $s_t\in \mathbb{R}^l$ and $\lambda_t\in \mathbb{R}^n$ are the dual variables.

Let $\bm c(\bm x, \bm u)=\big(c(x_0,u_0),\ldots,c(x_{N-1},u_{N-1})\big)$, $\bm \lambda=(\lambda_0,\ldots,\lambda_N)$, $\bm s=(s_0,\ldots,s_{N-1})$ and $\bm S=\operatorname{diag}[\bm s]$. The perturbed KKT system is defined as

\begin{equation}\label{KKT}
\begin{aligned}
\nabla_{\bm x} L(\bm x,\bm u,\bm \lambda, \bm s) &= 0,\\
\nabla_{\bm u} L(\bm x,\bm u,\bm \lambda, \bm s) &= 0,\\
\nabla_{\bm \lambda} L(\bm x,\bm u,\bm \lambda, \bm s) &= 0,\\
\bm S \bm c(\bm x, \bm u) +\bm \mu&=0,\\
\bm c(\bm x, \bm u)\leq0,\ \bm s&\geq 0,
\end{aligned}
\end{equation}

where $\nabla$ is the gradient operator, inequalities are understood to be element-wise, and $\bm \mu$ is a vector of perturbations $\mu>0$ of an appropriate dimension. The KKT system, i.e., system (<ref>) with $\mu =0$, defines the first-order necessary conditions for a local constrained minimiser.

The KKT system has a solution $(\bm x^{\star}, \bm u^{\star},\bm \lambda^{\star}, \bm s^{\star})$, which satisfies the following conditions:

* Strict complementary holds at the solution, i.e., $c(x_t^{\star}, u_t^{\star})<s^{\star}_t$ for $t\in\{0,\ldots,N-1\}$.

* The standard second-order constrained optimality conditions[Here we refer to the conditions of the type outlined in <cit.>]. hold at $(\bm x^{\star}, \bm u^{\star},\bm \lambda^{\star}, \bm s^{\star})$.

Under Assumption <ref> and some additional regularity requirements, the optimal solution of (<ref>) can be obtained as the limit point of the solutions to the perturbed KKT system for decaying perturbation, e.g., see [20].

For the clarity of notation, we drop decision variables' iteration indices. Instead we use the superscript ${}^+$ to denote the value of a variable at the next iteration. A variable's value at the current iteration is denoted by the variable name without this superscript.

§ PRIMAL-DUAL INTERIOR-POINT DIFFERENTIAL DYNAMIC PROGRAMMING

In this section we introduce the primal-dual Interior-Point Differential Dynamic Programming (IPDDP) algorithms for solving inequality constrained optimal control problems[As no time invariance properties are required, the proposed algorithms can be applied to the case where $f$, $c$ and $q$ functions are time-varying.]. Note that equality constraints can be addressed with the method of Lagrange multipliers by applying the minimax DDP technique[4] to the resulting Lagrangian function.

Starting from now we make no distinction between decision variable $x_0$ and initial state $\bar x_0$. This causes no ambiguity since we always operate with dynamically feasible state trajectories. To illustrate the underlying idea behind the algorithm we apply Bellman's principle of optimality to transform (<ref>) into

\begin{equation*}
\min_{\substack{u_0\ s.t.\\ c(x_0,u_0)\leq 0}}\Big[ q(x_0,u_0) + \min_{\substack{u_1\ s.t.\\ c(f(x_0,u_0),u_1)\leq 0}}\big[q(f(x_0, u_0), u_1) + \ldots \big]\Big]
\end{equation*}

Invoking the optimality principle, i.e., using $J_0^\star(x):= p(x)$ and

\begin{equation*}
J^{\star}_k(x)=\min_{\substack{u\ s.t.\\ c(x,u)\leq 0}}\Big[ q(x,u) + J^{\star}_{k-1}(f(x, u))\Big],
\end{equation*}

where $x \in \mathbb{R}^n$, $u \in \mathbb{R}^m$ and $k\in \{1,\ldots, N\}$, the optimisation problem (<ref>) is written as

\begin{equation}\label{nested}
J^{\star}_N(x_0)=\min_{\substack{u_0\ s.t.\\ c(x_0,u_0)\leq 0}}\Big[ q(x_0,u_0) + J^{\star}_{N-1}(f(x_0, u_0))\Big].
\end{equation}

This enables one to apply the dynamic programming method to solve the optimisation problem (<ref>) as a sequence of nested optimisation problems, which is the main idea behind DDP algorithm[1]. However, in this paper, instead of solving the nested problems (<ref>) sequentially as in [12, 13, 14, 15], we replace the constrained minimisation problems with their min-max primal-dual counterparts as follows

\begin{equation}\label{pdnested}
J^{\star}_k(x)=\min_{u}\max_{s\geq 0}\Big[ \ell(x, u, s)+ J^{\star}_{k-1}(f(x, u))\Big],
\end{equation}

where $\ell(x,u,s):=q(x,u) +  s^Tc(x,u)$.

§.§ Primal-Dual Feasible-Interior-Point DDP

Let $\mu$ be a strictly positive constant and assume that a current solution estimate for (<ref>) is given by tuple $(\bm x, \bm u, \bm s)$, such that $x_{t+1}=f(x_{t}, u_{t})$, $c(x_t, u_t)<0$ and $s_t>0$ for $t\in\{0,\ldots,N-1\}$. The solution approach aims to improve the solution estimate by calculating the control inputs updates that minimise quadratic models of (<ref>) in the vicinity of the current trajectory through a Backward Pass, and computing a new trajectory after a Forward Pass. The backward and forward passes are carried out repeatedly until a convergence criterion is satisfied.

§.§.§ Backward pass
Define quadratic functions

\begin{equation*}
V^{t}(x):= {\big(V_x^t\big)}^T(x-x_t) + \frac{1}{2}(x-x_t)^TV_{xx}^t(x-x_t),
\end{equation*}

where $x\in \mathbb{R}^n$ and $t\in\{0,\ldots,N\}$. Let $V_x^N=p_x(x_N)$ and $V_{xx}^N=p_{xx}(x_N)$ are the gradient and Hessian of function $p(x)$ at $x_N$ respectively, while the remaining coefficients $V_x^t$ and $V_{xx}^t$ are to be defined recursively later in the following.

Define $Q^t(x,u,s): = \ell(x, u,s)+V^{t+1}\big(f(x,u)\big)$ for $t\in\{0,\ldots,N-1\}$, and consider its second-order variation $\delta Q^t(\delta x, \delta u, \delta s)$ around $(x_t, u_t, s_t)$

\begin{multline*}
\delta Q^t(\delta x, \delta u, \delta s):={\begin{bmatrix}Q_x^t \\ Q_u^t \\ Q_s^t \end{bmatrix}}^T{\begin{bmatrix}\delta x \\ \delta u\\ \delta s \end{bmatrix}}
\\+ \frac{1}{2}
{\begin{bmatrix} \delta x \\ \delta u\\ \delta s \end{bmatrix}}^{T}
\begin{bmatrix}Q_{xx}^t&Q_{xu}^t&Q_{xs}^t\\
\end{bmatrix}
{\begin{bmatrix}\delta x \\ \delta u\\ \delta s \end{bmatrix}},
\end{multline*}

with the derivatives being evaluated at $(x_t,u_t,s_t)$:

\begin{align}
Q_s^t&=c(x_t,u_t),\ Q_{sx}^t=c_x,\ Q_{su}^t=c_u,\ Q^t_{ss}=0,\\
Q_x^t&=\ell_x+f_x^T V_x^{t+1},\ Q_u^t=\ell_u+f_u^T V_x^{t+1},\\
Q_{xx}^t&=\ell_{xx}+f_x^T V_{xx}^{t+1}f_x+V_{x}^{t+1}\cdot f_{xx},\\
Q_{uu}^t&=\ell_{uu}+f_u^T V_{xx}^{t+1}f_u+V_{x}^{t+1}\cdot f_{uu},\\
Q_{xu}^t&=\ell_{xu}+f_x^T V_{xx}^{t+1}f_u+V_{x}^{t+1}\cdot f_{xu},
\end{align}

where $\cdot$ denotes tensor contraction along an appropriate dimension.
The aim is to construct a quadratic model of $J^{\star}_{N-t}(x)$ in the vicinity of $x_t$ by providing suitable expressions for $V^t(x)$, and obtain local improvements to the decision variables. To this aim, we consider a solution to

\begin{equation}\label{local}
\min_{\delta u}\Big[\max_{\delta s} \delta Q^t(\delta x,\delta u,\delta s)\quad \text{s.t.}\quad s_t +\delta s\geq0 \Big],
\end{equation}

where this problem is a local approximation to
$$\min_{u}\max_{s\geq0} \Big[\ell(x, u,s)+V^{t+1}\big(f(x,u)\big)\Big].$$

Local minimisation of (<ref>) with respect to $\delta u$ results in

\begin{equation} \label{eq:deltau}
Q_u^t+Q_{ux}^t \delta x +Q_{uu}^t\delta u +Q_{us}^t \delta s =0.
\end{equation}

Additionally, for $\delta s$ to be a stationary point of the inner maximization problem in (<ref>), it must satisfy $s_t+\delta s \geq 0$ and

\begin{equation}\label{eq:deltas}
(s_t+\delta s)\odot(Q_s^t + Q_{sx}^t\delta x + Q_{su}^t \delta u)=0,
\end{equation}

where $\odot$ is the element-wise vector multiplication. The equation obtained from dropping all the second-order terms in (<ref>) and adding the perturbation vector $\bm \mu$ to its left-hand side, along with (<ref>)

yield the following parametric system

\begin{equation} \label{eqnsys}
\begin{bmatrix}
Q_{uu}^t & Q_{us}^t\\
S_t Q_{su}^t &  C_t
\end{bmatrix}\begin{bmatrix}
\delta u\\
\delta s
\end{bmatrix}= -\begin{bmatrix}
\end{bmatrix}- \begin{bmatrix}
S_t Q_{sx}
\end{bmatrix}\delta x,
\end{equation}

where $r_t:=S_tc(x_t,u_t)+\bm \mu$, $C_t:=\operatorname{diag}[c(x_t,u_t)]$ and $S_t:=\operatorname{diag}[s_t]$. The parametric solution of (<ref>) as a function of $\delta x$ is given by

\begin{equation}\label{updatelaws}
\begin{bmatrix}
\delta u\\
\delta s
\end{bmatrix}=\begin{bmatrix}
\alpha_t\\
\eta_t
\end{bmatrix}+\begin{bmatrix}
\beta_t\\
\theta_t
\end{bmatrix}\delta x,
\end{equation}

where the coefficients are established from

\begin{equation}
\label{coefficients}
\begin{bmatrix}
\alpha_t & \beta_t\\
\eta_t & \theta_t
\end{bmatrix}=-\begin{bmatrix}
Q_{uu}^t & Q_{us}^t \\
S_t Q_{su}^t & C_t
\end{bmatrix}^{-1}\begin{bmatrix}
Q_u^t & Q^t_{ux}\\
r_t & S_t Q^t_{sx}
\end{bmatrix}.
\end{equation}

Note that (<ref>) can be readily computed for $t=N-1$ as $V^N(x)$ is determined by function $p(x)$. Next, the expressions for the coefficients of $V^t(x)$ for $t\in\{0,\ldots,N-1\}$ are

\begin{equation*}\label{updateV}
\begin{aligned}
% V_0^t&=Q^t(x_t,u_t,s_t) + \alpha_t^T \hat Q_u^t + \frac{1}{2}\alpha_t^T \hat Q_{uu}^t \alpha_t=Q_0^t + \frac{1}{2}\alpha_t^T \hat Q_u^t,\\
V_x^t&=\hat Q_x^t+\beta_t^T\hat Q_u^t+(\hat Q_{ux}^t+ \hat Q_{uu}^t\beta_t)^T\alpha_t=\hat Q_x^t+\hat Q_{xu}^t\alpha_t,\\
V_{xx}^t&=\hat Q_{xx}^t+\hat Q_{xu}^t\beta_t +\beta_t^T \hat Q_{ux}^t+\beta_t^T \hat Q_{uu}^t\beta_t=\hat Q_{xx}^t+\hat Q_{xu}^t\beta_t,
\end{aligned}
\end{equation*}

\begin{align}
\hat Q_x^t&=Q_x^t - Q_{xs}^t C_t^{-1}r_t,\\
\hat Q_u^t&=Q_u^t - Q_{us}^t C_t^{-1}r_t,\\
\hat Q_{xx}^t&=Q_{xx}^t -Q_{xs}^t C_t^{-1}S_tQ_{sx}^t,\\
\hat Q_{xu}^t&=Q_{xu}^t -Q_{xs}^t C_t^{-1}S_tQ_{su}^t,\\
\hat Q_{uu}^t&=Q_{uu}^t -Q_{us}^t C_t^{-1}S_t Q_{su}^t. \label{hatQuu}
\end{align}

Expressions (<ref>) are established by replacing the standard DDP coefficients with the corresponding coefficients of the condensed system obtained from eliminating $\delta s$ in (<ref>).

§.§.§ Forward pass
Define the update functions

\begin{equation}\label{updateguess}
\begin{aligned}
\phi_t(x)&:=u_t+\alpha_t+\beta_t (x - x_t),\\
\psi_t(x)&:=s_t+\eta_t+\theta_t (x - x_t),
\end{aligned}
\end{equation}

denote a new iterate by $(\bm x^+, \bm u^+, \bm s^+)$, where $\bm x^+=(x_0^+,\ldots, x_N^+)$, $\bm u^+=(u_0^+,\ldots, u_{N-1}^+)$, $\bm s^+=(s_0^+,\ldots, s_{N-1}^+)$ and $x_0^+=x_0$. For $t\in\{0,\ldots,N-1\}$ we compute

\begin{align}
u_t^+&=\phi_t(x_t^+), \notag\\
s_t^+&=\psi_t(x_t^+), \label{newiter} \\
x_{t+1}^+&=f(x_t^+, u_t^+). \notag
\end{align}

Note that the strict constraints satisfaction for the new iterate should be ensured. In other words, the iterate should satisfy $c(x_t^+, u_t^+)<0$ and $s_t^+>0$ for $t\in\{0,\ldots,N-1\}$. A common approach is to employ a line-search procedure by including a step length $\gamma\in [0, 1]$ in the definition of functions as follows

\begin{equation*} \label{updateguess2}
\begin{aligned}
\phi_t(x, \gamma)&=u_t+\gamma \alpha_t+\beta_t (x - x_t),\\
\psi_t(x, \gamma)&=s_t+\gamma \eta_t+\theta_t (x - x_t).
\end{aligned}
\end{equation*}

Line-search procedure of this type is common for DDP methods, e.g., see  <cit.> or <cit.>.

§.§ Properties of Feasible-Interior-Point DDP

Due to the “direct shooting” nature of the proposed algorithm the state trajectory vectors are not independent variables, but are functions of the initial state and control inputs, i.e., $x_{t}=f(x_{t-1},u_{t-1})$ for $t\in\{0,\ldots,N\}$. We, however, continue to write $\bm x$ to shorthand the writing.

A tuple $w=(\bm x, \bm u, \bm s)$ is (strictly) primal-dual feasible, if it satisfies (with strict inequalities) $c(x_t, u_t)\leq0$, $s_t\geq 0$ and $x_{t+1}=f(x_t,u_t)$ for all $t\in\{0,\ldots,N-1\}$.

Matrices $\hat Q_{uu}^t$, defined as in (<ref>) (or later in (<ref>)), are positive definite for all $t\in\{0,\ldots,N-1\}$.

Under Assumption <ref>, if $w=(\bm x, \bm u, \bm s)$ is strictly primal-dual feasible then the linear operators

\begin{equation*}
Q_{uu}^t & Q_{us}^t\\
S_t Q_{su}^t & C_t
\end{bmatrix}^{-1},\ t\in\{0,\ldots,N-1\},
\end{equation*}

are continuous at $w$.

Matrices $\hat Q_{uu}^t=Q_{uu}^t - Q_{us}^t C_t^{-1}S_tQ_{su}^t$ are positive definite by assumption and $C_t$ are negative definite by strict primal-dual feasibility of $w$, thus invertible. This is a sufficient condition for existence of $P_t$, which can be found using a block-wise matrix inverse formula.

The matrix inverse (when exists) is continuous in $w$ when its matrix components are continuous function of $w$. As the continuity of $Q_{uu}^{N-1}$, $Q_{us}^{N-1}$ and $C_{N-1}$ follows directly from $q(x,u)$, $c(x,u)$ and $p(x)$ being twice continuously differentiable, we readily conclude continuity of $P_{N-1}$. Next, since $f(x,u)$ is twice continuously differentiable, we establish continuity of $Q_{uu}^{N-2}$ and $P_{N-2}$. Following the recursion we can conclude continuity of $P_t$ for $t=N-3$ to $t=0$.

Now consider a vector-valued function $F(w, \mu)$ defined as

\begin{equation*}
F(w, \mu):=\big(Q_u^0, \ldots, Q_u^{N-1}, r_0, \ldots, r_{N-1}\big),
\end{equation*}

where $Q_u^t$ are defined as in (<ref>) and $r_t=S_t c(x_t,u_t)+ \bm \mu$. Note that any primal-dual feasible tuple $w$, which is a solution of $F(w, \mu)=0$ for a given $\mu>0$, is strictly primal-dual feasible, as $r_t=0$ precludes zero components in $s_t$ and $c(x_t,u_t)$, and is a stationary point of the IPDDP algorithm, since all $\alpha_t$ and $\eta_t$ in (<ref>) are zero in this case. Moreover, such $w$ is a perturbed KKT point, i.e., a solution of the perturbed KKT system, as it is established in the following theorem, and thus can be moved arbitrary close to the (locally) optimal solution of (<ref>).

Let $w=(\bm x, \bm u, \bm s)$ be primal-dual feasible. If $w$ is a solution of $F(w, \mu)=0$ for a given $\mu>0$, then it is a perturbed KKT point, i.e., there exists $\bm \lambda = (\lambda_0,\ldots,\lambda_N)$ such that $(\bm x, \bm u, \bm \lambda, \bm s)$ satisfies (<ref>).

Note that satisfaction of $\nabla_{\lambda_t}L(\bm x, \bm u, \bm \lambda, \bm s)=f(x_{t-1},u_{t-1})-x_t=0$, $c(x_t, u_t)\leq 0$ and $s_t\geq0$ for all $t\in\{0,\ldots,N\}$ follows directly from the primal-dual feasibility hypothesis.

Now assume $w$ is a primal-dual solution of $F(w, \mu)=0$, meaning that $Q^t_u=0$ and $r_t=S_tc(x_t,u_t)+\bm \mu=0$ for all $t\in\{0, \ldots, N-1\}$, and note that it directly ensures satisfaction of $\bm S \bm c(\bm x, \bm u) +\bm \mu=0$ in (<ref>). Also note that
$\nabla_{x_N}L(\bm x, \bm u, \bm \lambda, \bm s)=0$ is satisfied when $\lambda_N=V_x^N$. For $t=N-1$ we have

\begin{align*}
%   Q^{N-1}_x&=\ell_x+f_x^TV_x^N=\ell_x+f_x^T\lambda_N,\\
%   Q^{N-1}_u&=\ell_u+f_u^TV_x^N=\ell_u+f_u^T\lambda_N,\\
\nabla_{x_{N-1}}L&=\ell_x+f_x^T \lambda_N -\lambda_{N-1}=Q^{N-1}_x-\lambda_{N-1}, \\
\nabla_{u_{N-1}}L&=\ell_u+f_u^T \lambda_N=Q^{N-1}_u=0,
\end{align*}

where the derivatives are evaluated at $(x_{N-1}, u_{N-1}, s_{N-1})$. Proceeding with $\lambda_t=V^{t}_x=Q^{t}_x$ we obtain $(\bm x, \bm u, \bm \lambda, \bm s)$, where $\bm \lambda =(\lambda_0,\ldots,\lambda_N)$, which satisfies system (<ref>).

Moreover, once $w$ is sufficiently close to a solution of $F(w, \mu)=0$, the IPDDP iterates converge to the solution at a quadratic rate. This is demonstrated in the following theorem.

Let $w=(\bm x, \bm u, \bm s)$ and $w^+=(\bm x^+, \bm u^+, \bm s^+)$, defined by (<ref>), are strictly primal-dual feasible, and $\mu>0$. There exist $M\geq 0$ and $\varepsilon>0$ such that if $\|w-w^\star\|\leq\varepsilon$, then

\begin{equation}\label{thm1:eq}
\begin{aligned}
\| w^+ - w^\star\| &\leq M \| w - w^\star\|^2,\\
\| w^+ - w^\star\| &< \| w - w^\star\|,
\end{aligned}
\end{equation}

where $w^\star=(\bm x^\star, \bm u^\star, \bm s^\star)$ is a primal-dual feasible solution of $F(w, \mu)=0$.

Consider functions defined by $Q_u^t$ and $r_t$, and note that they are differentiable, as all the terms involved are differentiable. Then, by the Taylor theorem there exist functions $h_t(w, w^\star)$ and $g_t(w, w^\star)$ such that

\begin{align}\label{taylor}
Q_u^t\big\rvert_{w^\star}= Q_u^t + Q_{ux}^t (x_t^\star-x_t) + Q_{uu}^t (u_t^\star - u_t) \\+ Q_{us}^t (s_t^\star - s_t)+ h_t(w,w^\star)&=0,\notag\\
r_t\big\rvert_{w^\star}= r_t+S_t Q_{sx}^t (x_t^\star - x_t) + S_t Q_{uu}^t (u_t^\star - u_t)  \notag\\+ C_t (s_t^\star - s_t)+ g_t(w,w^\star)&=0,\notag
\end{align}

where we used that $w^\star$ is a solution of $F(w,\mu)=0$, and the norms of the residual functions are bounded:

\begin{align*}
\|h_t(w,w^\star)\|&\leq H_t \|w-w^\star\|^2,\\
\|g_t(w,w^\star)\|&\leq G_t \|w-w^\star\|^2,
\end{align*}

where $H_t$ and $G_t$ for $t\in\{0,\ldots,N-1\}$ are some constants.
Denote $\Delta x_t = x_t^+ - x_t$, $\Delta u_t = u_t^+ - u_t$, $\Delta s_t = s_t^+ - s_t$ for $t\in\{0,\ldots,N-1\}$, and notice that $(\Delta x_t,\Delta u_t,\Delta u_t)$ is a solution of system (<ref>), as it belongs to parametric family (<ref>):

\begin{align*}
\Delta u_t&= u_t^+ - u_t = \phi (x_t^+) - u_t= \alpha_t + \beta_t \Delta x_t,\\
\Delta s_t&= s_t^+ - s_t = \psi (x_t^+) - s_t= \eta_t + \theta_t \Delta x_t.
\end{align*}

Now, by adding and subtracting $x_t^+$, $u_t^+$, $s_t^+$ in the appropriate parentheses of (<ref>), and using the fact that $(\Delta x_t,\Delta u_t,\Delta u_t)$ is a solution of (<ref>), we have

\begin{align*}\label{taylorsys}
Q_{ux}^t (x_t^\star-x_t^+) + Q_{uu}^t (u_t^\star -u_t^+)+ Q_{us}^t (s_t^\star -s_t^+)\\+ h_t(w,w^\star)&=0,\\
S_t Q_{sx}^t (x_t^\star -x_t^+) + S_t Q_{su}^t (u_t^\star -u_t^+) + C_t (s_t^\star -s_t^+)\\+ g_t(w,w^\star)&=0,
\end{align*}

which can be rewritten using operator $P_t$ as

\begin{equation*}
\begin{bmatrix}
u_t^\star -u_t^+\\
s_t^\star -s_t^+
\end{bmatrix}=-P_t\begin{bmatrix}
\end{bmatrix}-P_t \begin{bmatrix}
S_t Q^t_{sx}
\end{bmatrix} (x_t^\star-x_t^+).
\end{equation*}

Since $P_0$ is a continuous function of $w$, and $x_0^\star=x_0^+=x_0$, for $t=0$ this gives

\begin{align*}
\big\|[u_0^\star - u_0^+, s_0^\star - s_0^+]\big\| & \leq \big\|P_0\big\|\big\|(h_0(w,w^\star),g_0(w,w^\star))\big\|\\&\leq M_0\big\|w-w^\star\big\|^2,
\end{align*}

where $M_0=(H_0+G_0)\max_{w\in\Omega}\|P_0\|$, and the operator norm $\|P_0\|$ (induced by Euclidian vector norm) of continuous linear operator $P_0$ is bounded on a compact subset $\Omega$ of all strictly primal-dual feasible $w$, which can be chosen to be sufficiently big, i.e., such that it contains set of all strictly primal-dual feasible $w: \| w - w^\star\|\leq\varepsilon$ for any $\varepsilon \leq \bar \varepsilon$, where $\bar \varepsilon>0$ is a constant.
By using Taylor formula for $f(x_0^\star,u_0^\star)$, and bounding the residuals in $\|w-w^\star\|\leq \varepsilon$ we get:

\begin{align*}
\|x_1^\star - x_1^+\| &= \big\|f(x_0^\star,u_0^\star) - f(x_0^+,u_0^+)\big\| \\
&\leq \big\|f_x^T (x_0^\star-x_0^+) + f_u^T (u_0^\star-u_0^+)\big\| \\ &\quad + K_1 \big\| (x_0^\star-x_0^+, u_0^\star-u_0^+)\big\|^2\\
&\leq K_2 \big\| (x_0^\star-x_0^+, u_0^\star-u_0^+)\big\|\leq K_3 \| w-w^\star\| ^2,
\end{align*}

for some constants $K_1,K_2,K_3\geq0$. Continuing by induction, there are some constants $M_t\geq0$ for $t\in\{0,\ldots,N-1\}$, such that for a positive $\varepsilon\leq \bar \varepsilon$ and all strictly primal-dual feasible $w:\|w-w^\star\|\leq \varepsilon $ we have:

$$\big\|(x_t^\star - x_t^+, u_t^\star - u_t^+, s_t^\star - s_t^+)\big\| \leq M_t\big\|w-w^\star\big\|^2,$$

which proves the first inequality in (<ref>) by choosing $M=M_1+M_2+\ldots+M_N$:

\begin{align*}
\big\|w^\star - w^+\big\| & \leq\sum_{t=1}^N \big\|(x_t^\star - x_t^+, u_t^\star - u_t^+, s_t^\star - s_t^+)\big\|\\& \leq M \big\|w-w^\star\big\|^2.
\end{align*}

The second inequality in (<ref>) is obtained from the first one by pickings $\varepsilon$ sufficiently small, i.e., such that $\varepsilon<1/ M$ in addition to $\varepsilon\leq\bar \varepsilon$. Now, for all strictly primal-dual feasible $w$ such that $\|w-w^\star\|\leq \varepsilon$ we additionally have $
\big\|w^+-w^\star\big\| \leq \varepsilon M  \big\|w-w^\star\big\|<\big\|w-w^\star\big\|$.

§.§ Primal-Dual Infeasible-Interior-Point DDP

In the previous sections, we assumed that one has access to a strictly primal-dual feasible initial solution guess, which might be quite restrictive in practice. Here we propose an extension where this assumption is relaxed. In this section we explain what should be modified in the Feasible-IPDDP algorithm, and how the theoretical results translate to a new algorithm.

§.§.§ Backward pass

We use the same definitions for $V^{t}(x)$ and $Q^t(x)$, and follow the same steps to obtain the first equation of system (<ref>). Next, however, we introduce slack variables $y_t\geq 0$ for $t\in\{0,\ldots,N-1\}$ to transform the inequality constraints into equalities, i.e., $c(x_t, u_t) + y_t = 0$, and use the perturbed complementarity condition, formulated in terms of the dual and slack variables, i.e., $S_ty_t-\bm \mu=0$.
Finally, we proceed by setting the first order Taylor expansion for $c(x_t, u_t) + y_t$ and $S_ty_t-\bm \mu$ to zero. In this way, the parametric system of equations becomes

\begin{equation}\label{infeas:sys}
\begin{bmatrix} Q_{uu}^t & Q_{us}^t & 0\\
Q_{su}^t & 0    & I\\
0      & Y_t  & S_t
\end{bmatrix}
\begin{bmatrix} \delta u\\
\delta s\\
\delta y
\end{bmatrix}
=-\begin{bmatrix} Q_u^t\\
\end{bmatrix}-\begin{bmatrix} Q_{ux}^t\\
\end{bmatrix}\delta x,
\end{equation}

where $r_t^p:=c(x_t,u_t)+y_t$, $r_t^d:=S_ty_t - \bm \mu$ and $Y_t:=\operatorname{diag}[y_t]$.
Note that elimination of $\delta y = -r^p_t - Q_{su}\delta u - Q_{sx}\delta x$ from the system leads to a reduced system similar to (<ref>)

\begin{equation*}
\begin{bmatrix} Q_{uu}^t & Q_{us}^t\\
S_t Q_{su}^t & -Y_t\\
\end{bmatrix}
\begin{bmatrix} \delta u\\
\delta s
\end{bmatrix}
=-\begin{bmatrix} Q_u^t\\
\hat r_t
\end{bmatrix}-\begin{bmatrix} Q_{ux}^t\\
S_t Q_{sx}^t
\end{bmatrix}\delta x,
\end{equation*}

where $\hat r_t:=S_t r^p_t-r^d_t$.
Solving the parametric system

\begin{equation*}
\label{infeas:coefficients}
\begin{bmatrix}
\alpha_t & \beta_t\\
\eta_t & \theta_t\\
\chi_t & \zeta_t
\end{bmatrix}=-\begin{bmatrix} Q_{uu}^t & Q_{us}^t & 0\\
Q_{su}^t & 0    & I\\
0      & Y_t  & S_t
\end{bmatrix}^{-1}\begin{bmatrix}
Q_u^t & Q^t_{ux}\\
r_t^p & Q^t_{sx}\\
r_t^d & 0
\end{bmatrix},
\end{equation*}

yields the new coefficients

\begin{align}
\hat Q_x^t &=Q_x^t + Q_{xs}^t Y_t^{-1}\hat r_t,\\
\hat Q_u^t &=Q_u^t + Q_{us}^t Y_t^{-1}\hat r_t,\\
\hat Q_{xu}^t&=Q_{xu}^t +Q_{xs}^t Y_t^{-1}S_tQ_{su}^t,\\
\hat Q_{xx}^t&=Q_{xx}^t +Q_{xs}^t Y_t^{-1}S_tQ_{sx}^t,\\
\hat Q_{uu}^t&=Q_{uu}^t +Q_{us}^t Y_t^{-1}S_t Q_{su}^t. \label{hatQinfeasuu}
\end{align}

§.§.§ Forward pass

Define the update functions

\begin{align*}
\phi_t(x)&=u_t + \alpha_t + \beta_t (x - x_t),\\
\psi_t(x)&=s_t+\eta_t + \theta_t (x -x_t),\\
\xi_t(x)&=y_t+\chi_t + \zeta_t(x -x_t),
\end{align*}

and denote a new solution guess by $(\bm x^+, \bm u^+, \bm s^+, \bm y^+)$, where $\bm x^+=(x_0^+,\ldots, x_N^+)$, $\bm u^+=(u_0^+,\ldots, u_{N-1}^+)$, $\bm s^+=(s_0^+,\ldots, s_{N-1}^+)$, $\bm y^+=(y_0^+,\ldots, y_{N-1}^+)$ and $x_0^+=x_0$. For $t\in\{0,\ldots,N-1\}$ we compute

\begin{equation}\label{infeas:iter}
\begin{aligned}
x_{t+1}^+&=f(x_t^+, u_t^+).
\end{aligned}
\end{equation}

Note that now only positivity of $s_t^+$ and $y_t^+$ must be ensured. Satisfaction of $c(x,u)\leq 0$ is only guaranteed in the limit of the convergence, which is a common property for the so-called “infeasible” methods. This potentially facilitates algorithm implementation, but an appropriate line-search procedure is still required in practical implementations, see [19] for details.

§.§ Properties of Infeasible-Interior-Point DDP

Now we add slack variables $\bm y=(y_0,\ldots,y_{N-1})$ into the tuple $w = (\bm x, \bm u, \bm s, \bm y)$, and provide the definition of dual feasibility as follows.

A tuple $w=(\bm x, \bm u, \bm s, \bm y)$ is (strictly) dual feasible, if it satisfies (with strict inequalities) $s_t\geq 0$, $y_t\geq 0$ and $x_{t+1}=f(x_t,u_t)$ for all $t\in\{0,\ldots,N-1\}$.

Consider a vector-valued function $F(w, \mu)$ defined as

\begin{equation*}
F(w, \mu):=\big(Q_u^0, \ldots, Q_u^{N-1}, r_0^p, \ldots, r_{N-1}^p, r_0^d, \ldots, r_{N-1}^d\big),
\end{equation*}

where $r_t^p:=c(x_t,u_t) + y_t$ and $r_t^d:= S_ty_t -\bm \mu $.

Under Assumption <ref>, if $w = (\bm x,\bm u,\bm s,\bm y)$ is strictly dual feasible then the linear operators

\begin{equation*}
P_t=\begin{bmatrix} Q_{uu}^t & Q_{us}^t & 0\\
Q_{su}^t & 0    & I\\
0      & Y_t  & S_t
\end{bmatrix}^{-1},\ t\in\{0,\ldots,N-1\},
\end{equation*}

are continuous at $w$.

The proof is similar to that of Lemma <ref>.

Let $w=(\bm x, \bm u, \bm s, \bm y)$ be dual feasible. If $w$ is a solution of $F(w, \mu)=0$ for a given $\mu>0$, then it is a perturbed KKT point, i.e., there exists $\bm \lambda = (\lambda_0,\ldots,\lambda_N)$ such that $(\bm x, \bm u, \bm \lambda, \bm s)$ satisfies (<ref>).

Note that $r_t^d=c(x_t, u_t) +y_t=0$ for $t\in\{0,\ldots,N\}$ ensures primal-dual feasibility. With this in mind one can follow the same steps as in the proof of Theorem <ref>.

Let $w=(\bm x, \bm u, \bm s, \bm y)$ and $w^+=(\bm x^+, \bm u^+, \bm s^+, \bm y^+)$, defined by (<ref>), are strictly dual feasible, and $\mu>0$. There exist $M\geq 0$ and $\varepsilon>0$ such that if $\|w-w^\star\|\leq\varepsilon$, then

\begin{align}
\| w^+ - w^\star\| &\leq M \| w - w^\star\|^2, \\
\| w^+ - w^\star\| &< \| w - w^\star\|,
\end{align}

where $w^\star=(\bm x^\star, \bm u^\star, \bm s^\star, \bm y^\star)$ is a dual feasible solution of $F(w, \mu)=0$.

The proof follows the same steps as the proof of Theorem <ref>. We consider differentiable functions defined by $Q_u^t$, $r_t^p$ and $r_t^d$, and use its first-order Taylor expansion at $w^\star$:

\begin{align*}\label{inf:taylor}
Q_u^t\big\rvert_{w^\star}= Q_u^t + Q_{ux}^t (x_t^\star-x_t) + Q_{uu}^t (u_t^\star - u_t) \\+ Q_{us}^t (s_t^\star - s_t)+ h_t(w,w^\star)&=0,\\
r_t^p\big\rvert_{w^\star}= r_t^p + Q^t_{sx} (x_t^\star-x_t) + Q^t_{su}(u_t^\star-u_t) \\+ (y_t^\star-y_t) + g_t(w,w^\star)&=0,\\
r_t^d\big\rvert_{w^\star}= r_t^d+ Y_t (s_t^\star - s_t)+S_t (y_t^\star - y_t) \\+ k_t(w,w^\star)&=0,
\end{align*}

and the norms of the residuals are bounded

\begin{align*}
\|h_t(w,w^\star)\|&\leq H_t \|w-w^\star\|^2,\\ \|g_t(w,w^\star)\| &\leq G_t \|w-w^\star\|^2,\\
\|k_t(w,w^\star)\| & \leq K_t \|w-w^\star\|^2,
\end{align*}

where $H_t$, $G_t$ and $K_t$ for $t\in\{0,\ldots,N-1\}$ are some constants.
As before, we use $\Delta x_t = x_t^+ - x_t$, $\Delta u_t = u_t^+ - u_t$, $\Delta s_t = s_t^+ - s_t$, $\Delta y_t = y_t^+ - y_t$ for $t\in\{0,\ldots,N\}$, and note that $(\Delta x_t,\Delta u_t,\Delta u_t, \Delta y_t)$ is a solution of system (<ref>). Thus,

\begin{align*}
Q_{ux}^t (x_t^\star-x_t^+) + Q_{uu}^t (u_t^\star - u_t^+) + Q_{us}^t (s_t^\star - s_t^+)\\+ h_t(w,w^\star)&=0,\\
Q^t_{sx} (x_t^\star-x_t^+) + Q^t_{su}(u_t^\star-u_t^+) + (y_t^\star-y_t^+) \\+ g_t(w,w^\star)&=0,\\
Y_t (s_t^\star - s_t^+)+S_t (y_t^\star - y_t^+)  + k_t(w,w^\star)&=0,
\end{align*}

and using the new definition of operators $P_t$:

\begin{equation*}
\begin{bmatrix}
u_t^\star -u_t^+\\
s_t^\star -s_t^+\\
y_t^\star -y_t^+
\end{bmatrix}=- P_t\begin{bmatrix}
\end{bmatrix}- P_t \begin{bmatrix}
\end{bmatrix} (x_t^\star-x_t^+).
\end{equation*}

From here the remaining proof is straightforward: we start with $t=0$ and proceed iteratively to $t=N-1$ as it is described in the proof of Theorem <ref>.

§ NUMERICAL SIMULATIONS

Here we compare implementations of the proposed Interior-Point DDP (IPDDP) algorithms[ https://github.com/xapavlov/ipddp] with implementations of Control-Limited DDP (CLDDP)[21] and log-barrier DDP.

The wall-clock time required for solving problems using these methods depends mainly on the number of iterations and time spent per iteration. While we focus on the total number of iterations required for algorithms to converge to a (locally) optimal solution, we briefly note that the per-iteration complexity for each algorithm, due to their nature, is different. Specifically, CLDDP require solving $N$ box-QPs (of size $m$ with $l$ constraints) per iteration, IPDDP and log-barrier DDP require solving $N$ linear systems of equations (of “effective” size $m$) per iteration. Thus, one expects that each iteration of the CLDDP takes longer than that of its competitors. However, a fair wall-clock time comparison between methods is beyond the scope of this paper.

To ensure convergence of the proposed algorithms to the perturbed KKT points we rely on the line-search (see Remark <ref>), step filter[22] and regularisation of type $Q_{uu}^t+\gamma I$. The regularisation is essential when Assumption <ref> is violated, what commonly occurs at the early stages of the optimisation routine, but holds later on (in the considered examples). Perturbation $\mu$ is initialised as $\mu \leftarrow J(\bm x, \bm u)/(N l)$, where $J(\bm x, \bm u)$ is the value of the objective function. Perturbation is reduced as $\mu \leftarrow \min(\mu/\kappa, \mu^{1.2})$, where we pick the reduction factor $\kappa=5$ (unless stated otherwise),  every time $\| F(w,\mu)\|_{\infty}$ is sufficiently small[Since the infinity norms of the objective function gradient and constraints violations at a step are proportional with $\mu$, it can be used as a surrogate for the progress of the algorithm.], i.e., less than than $0.2\mu$ in our implementation.
In all trials we use random initial control inputs, which are sampled from interval $[-0.01, 0.01]$ at uniform, and conduct 40 trials for each experiment. We use the following definition for the logarithmic optimality error

$$E_J:= \log_{10}\big[J(\bm x, \bm u) - J(\bm x^{\star}, \bm u^{\star})\big],$$

where $J(\cdot, \cdot)$ is the objective function evaluated at a current solution guess $(\bm x, \bm u)$ or at a (numerically obtained) locally optimal solution $(\bm x^{\star}, \bm u^{\star})$.

§.§ Inverted pendulum

Consider a task of stabilising the inverted pendulum

\begin{equation*}
f(x,u) =\begin{bmatrix}
\varphi+h \omega\\
\omega+h\sin(\varphi)+h u
\end{bmatrix},
\end{equation*}

where state $x=[\varphi\ \omega]^T$ ($\varphi$ is angle and $\omega$ is angular velocity), $u$ is the control input, $h=0.05$; initial state is $x_0=[-\pi\ 0]^T$, control constraints are $-0.25\leq u\leq 0.25$, and $N=500$. We choose the quadratic stage and terminal costs as

\begin{equation*}
q(x,u) = 0.025(\varphi^2 + \omega^2 + u^2),\ p(x)= 5(\varphi^2 + \omega^2).
\end{equation*}

[ ][Feasible-IPDDP ]

[ ][Infeasible-IPDDP ]

[ ][Log-barrier DDP]

[ ][CLDDP]

Inverted pendulum: logarithmic optimality error $E_J$ on y-axis vs iteration number on x-axis.

Fig. <ref> shows that the log-barrier DDP, Feasible- and Infeasible-IPDDP consistently convergence to the optimal solution in about 150 iterations. On the other hand, most runs of Control-Limited DDP (CLDDP) take more than 200 iterations to converge. We believe that the bang-bang nature of the optimal control complicates the correct identification of active and inactive constraints for CLDPP algorithm. This potential issue is mitigated for the other algorithms by sufficiently perturbing the problem, and following the central path, i.e., a trace of stationary points defined for different values of perturbations, see Fig. <ref>.

Inverted pendulum: control inputs $u_t$ on y-axis vs time-index $t$ on x-axis.

Inverted pendulum: the smallest eigenvalue of $\hat Q_{uu}^t$ on y-axis vs iteration number on x-axis (10 trials of Feasible-IPDDP).

Inverted pendulum: the smallest eigenvalue of $\hat Q_{uu}^t$ on y-axis vs perturbation $\mu$ on x-axis.

Fig. <ref> demonstrates that Assumption <ref> is commonly violated at the early stage of the optimisation routine, i.e., $\hat Q_{uu}^t$ is not positive definite, and the regularisation has to applied. Later in this example, after about 70 iterations, Assumption <ref> is always satisfied, and the quadratic convergence rate steps in. This is because $\hat Q_{uu}^t\succ0$ in a neighbourhood (by continuity) of the central path: see Fig. <ref>, which depicts the smallest eigenvalues (over all $t\in\{0,\ldots,N\}$) of $\hat Q_{uu}^t$ evaluated at the perturbed KKT points along the central path.

§.§ Car parking problem

Next we consider a car parking problem as in [16]. The dynamics of a car with state $x=[r_x\ r_y\ \varphi\ v]^T$ ($r_x$ and $r_y$ are the x- and y-coordinates, $\varphi$ is the car's heading and $v$ is its velocity) and input $u=[w\ a]^T$ ($w$ and $a$ are the front wheels' steering angle and acceleration respectively) is

\begin{equation*}
\begin{bmatrix}r_x+b(v,w) \cos(\varphi)\\
r_y+b(v,w) \sin(\varphi)\\
\varphi +\sin^{-1}\big( \frac{hv}{d}\sin(w)\big)\\
v + ha
\end{bmatrix},
\end{equation*}

where $b(v,w) =d+h v\cos(\omega) -\sqrt{d^2-h^2 v^2\sin^2(w)}$, $h=0.03$ and $d=2$ is a distance between the front and back axles of a car; initial state is $x_0=[1\ 1\ 3\ \pi/2\ 0]^T$, control constraints are $-0.5\leq w\leq 0.5$ and $-2\leq a\leq 2$, and $N=500$. The cost functions are

\begin{align*}
q(x,u)= 0.01\big(H(r_x,0.1)+H(r_y,0.1) + w^2 + 0.01a^2\big),\\
p(x) = H(r_x,0.1)+H(r_y,0.1)+H(\varphi,0.01)+H(v,0.1),
\end{align*}

where $H(y, z) =\sqrt{y^2 + z^2} - z$.

[ ][Feasible-IPDDP ]

[ ][Infeasible-IPDDP ]

[ ][Log-barrier DDP]

[ ][CLDDP]

Car-parking problem: logarithmic optimality error $E_J$ on y-axis vs iteration number on x-axis.

In many cases CLDDP converges faster than its competitors (see Fig. <ref>), but presence of very distinct locally optimal solutions greatly influences its performance. IPDDP and log-barrier algorithms follow the central path, thus behave in a more predictable manner. The main advantage of the proposed algorithms compared to the log-barrier variant is their improved ability to follow the central path. To demonstrate this, we solve the perturbed KKT system for perturbation $\mu$ sampled from a certain interval, and investigate the maximum stepsizes that the  Feasible-IPDDP and log-barrier DDP methods can take immediately after the reduction of the perturbation paramter. As one can see from Fig. <ref>, the log-barrier DDP is more sensitive to the perturbation reduction when $\mu\in[0.0001, 0.0005]$, and takes vanishingly small steps in the case where $\mu\in[0.0006, 0.0010]$. Although such severe incidents might not commonly take place, it is known that the Newton's method applied to the log-barrier problems, in contrast with the primal-dual approach, generally suffers from the step-scaling issues occurring after the reduction of the perturbation[23].

[ ][$0.0001\leq\mu\leq0.0005$]

[ ][$0.0006\leq\mu\leq0.0010$]

Car-parking problem: the maximum accepted stepsize on y-axis vs the perturbation reduction factor $\kappa$ on x-axis for two different intervals of the central path.

§.§ Unicycle motion control

Here we consider a unicycle model

$$f(x,u)= \begin{bmatrix}r_x + h v \cos(\varphi)\\r_y + hv \sin(\varphi) \\ \varphi+h u \end{bmatrix},$$
where $x=[r_x\ r_y\ \varphi]^T$ is the state vector ($r_x$ and $r_y$ are the x- and y-coordinates, $\varphi$ is the heading) and $u$ is the control input, $h=0.1$ and $v=1.5$; the initial state is $x_0=[-10\ 0\ 0]^T$ and $N=600$. We use the quadratic stage and terminal costs

\begin{equation*}
q(x,u)=0.1(r_x^2 + r_y^2 +\varphi^2 + 0.1 u^2),\ p(x)=0.1(r_x^2 + r_y^2 +\varphi^2).
\end{equation*}

The input and state constraints are chosen as follows

\begin{align*}
-1.5\leq u \leq 1.5,\
-1\leq r_y \leq 1,\ \|r_x+5.5, r_y+1\|\geq 1,\\
\|[r_x+8,\ r_y-0.2]\|\geq 0.5,\
\|[r_x+2.5,\ r_y-1]\|\geq 1.5.
\end{align*}

Fig. <ref> depicts two distinct locally optimal trajectories; guessing a feasible solution for this problem is not trivial. Next we will compare the relaxed logarithmic barrier strategy[24] implemented as DDP algorithm and Infeasible-IPDDP.

Unicycle motion control: two distinct locally optimal trajectories (solid lines) with $r_y$ on y-axis and $r_x$ on x-axis; the boundary of obstacles depicted as dotted circles.

[ ][Infeasible-IPDDP]
[ ][Relaxed log-barrier DDP]
Unicycle motion control: perturbation $\mu$ on y-axis vs iteration number on x-axis.

For the relaxed log-barrier DDP method we use the following penalty on the inequality constraints components

\begin{equation*}
\beta_\delta(z) =\begin{cases}
-\log z  & z>\delta,\\
\frac{1}{2}\Big[\Big( \frac{z-2\delta}{\delta}\Big)^2-1\Big] - \log \delta & z\leq \delta,
\end{cases}
\end{equation*}

and update the relaxation parameter as $\delta\leftarrow\mu$. For this experiment we sample initial perturbation $\mu$ from interval $[0.5, 1]$ at uniform. Fig. <ref> illustrates the convergence of the methods, where IPDDP outperforms its competitor. Note that $\mu$ is decreased only when the corresponding perturbed (or barrier) problem is solved with sufficient accuracy, i.e., when the norm of $F(\bm w, \mu)$ for Infeasible-IPDDP or the gradient of the penalised objective function for the relaxed log-barrier DDP are less than $0.2 \mu$, thus acting as an indicator of progress.

§ CONCLUSIONS

Two variants of Primal-Dual Interior-Point DDP algorithms, namely Feasible- and Infeasible-IPDDP have been proposed in this paper. The former operates with strictly feasible iterations, while the latter allows to start with infeasible solution guess, and ensures constraint satisfaction upon convergence. We prove the local quadratic convergence of the algorithms to its stationary points, and show that the stationary points of both algorithms satisfy the (perturbed) system of first-order conditions for optimality. The numerical simulations on three examples demonstrate the ability of the proposed algorithms to converge to the locally optimal solutions of the discrete-time optimal control problems with nonlinear dynamics and constraints, often being superior to its closest competitors in the total number (and variability) of iterations required to converge.

A future research direction is to establish the global convergence of IPDDP methods to the local optima. Also, we believe that the numerical performance can be further improved by implementing adaptive schemes for selecting the perturbation parameter and/or predictor-corrector steps. Other promising extensions include generalisation to differentiable manifolds.

[1]
D. Mayne, “A second-order gradient method for determining optimal trajectories
of non-linear discrete-time systems,” International Journal of
Control, vol. 3, no. 1, pp. 85–95, 1966.

[2]
J. F. A. De O. Pantoja, “Differential dynamic programming and newton's
method,” International Journal of Control, vol. 47, no. 5, pp.
1539–1553, 1988.

[3]
L.-Z. Liao and C. A. Shoemaker, “Convergence in unconstrained discrete-time
differential dynamic programming,” IEEE Transactions on Automatic
Control, vol. 36, no. 6, pp. 692–706, 1991.

[4]
J. Morimoto, G. Zeglin, and C. G. Atkeson, “Minimax differential dynamic
programming: Application to a biped walking robot,” in Proceedings
2003 IEEE/RSJ International Conference on Intelligent Robots and Systems
(IROS 2003)(Cat. No. 03CH37453), vol. 2.1em plus 0.5em minus
0.4emIEEE, 2003, pp. 1927–1932.

[5]
Y. Tassa, T. Erez, and W. D. Smart, “Receding horizon differential dynamic
programming,” in Advances in neural information processing systems,
2008, pp. 1465–1472.

[6]
W. Li and E. Todorov, “Iterative linear quadratic regulator design for
nonlinear biological movement systems.” in ICINCO (1), 2004, pp.

[7]
S. Sen and S. J. Yakowitz, “A quasi-Newton differential dynamic programming
algorithm for discrete-time optimal control,” Automatica, vol. 23,
no. 6, pp. 749–752, 1987.

[8]
J. Rajamäki, K. Naderi, V. Kyrki, and P. Hämäläinen, “Sampled
differential dynamic programming,” in 2016 IEEE/RSJ International
Conference on Intelligent Robots and Systems (IROS).1em plus 0.5em
minus 0.4emIEEE, 2016, pp. 1402–1409.

[9]
Z. Manchester and S. Kuindersma, “Derivative-free trajectory optimization with
unscented dynamic programming,” in 2016 IEEE 55th Conference on
Decision and Control (CDC).1em plus 0.5em minus 0.4emIEEE,
2016, pp. 3642–3647.

[10]
E. Pellegrini and R. P. Russell, “Applications of the multiple-shooting
differential dynamic programming algorithm with path and terminal
constraints,” in AAS/AIAA Astrodynamics Specialist Conference,
Stevenson, WA, 2017.

[11]
B. Plancher, Z. Manchester, and S. Kuindersma, “Constrained unscented dynamic
programming,” in 2017 IEEE/RSJ International Conference on Intelligent
Robots and Systems (IROS).1em plus 0.5em minus 0.4emIEEE,
2017, pp. 5674–5680.

[12]
D. M. Murray and S. J. Yakowitz, “Constrained differential dynamic programming
and its application to multireservoir control,” Water Resources
Research, vol. 15, no. 5, pp. 1017–1027, 1979.

[13]
S. Yakowitz, “The stagewise kuhn-tucker condition and differential dynamic
programming,” IEEE transactions on automatic control, vol. 31, no. 1,
pp. 25–30, 1986.

[14]
T. Lin and J. Arora, “Differential dynamic programming technique for
constrained optimal control,” Computational Mechanics, vol. 9, no. 1,
pp. 27–40, 1991.

[15]
Z. Xie, C. K. Liu, and K. Hauser, “Differential dynamic programming with
nonlinear constraints,” in 2017 IEEE International Conference on
Robotics and Automation (ICRA).1em plus 0.5em minus 0.4emIEEE,
2017, pp. 695–702.

[16]
Y. Tassa, N. Mansard, and E. Todorov, “Control-limited differential dynamic
programming,” in 2014 IEEE International Conference on Robotics and
Automation (ICRA).1em plus 0.5em minus 0.4emIEEE, 2014, pp.

[17]
G. Lantoine and R. P. Russell, “A hybrid differential dynamic programming
algorithm for constrained optimal control problems. part 1: theory,”
Journal of Optimization Theory and Applications, vol. 154, no. 2, pp.
382–417, 2012.

[18]
J. Gondzio, “Interior point methods 25 years later,” European Journal
of Operational Research, vol. 218, no. 3, pp. 587–601, 2012.

[19]
J. Nocedal and S. Wright, Numerical optimization.1em plus 0.5em
minus 0.4emSpringer Science & Business Media, 2006.

[20]
A. El-Bakry, R. A. Tapia, T. Tsuchiya, and Y. Zhang, “On the formulation and
theory of the newton interior-point method for nonlinear programming,”
Journal of Optimization Theory and Applications, vol. 89, no. 3, pp.
507–541, 1996.

[21]
Y. Tassa, “iLQG/DDP trajectory optimization,” MATLAB Central File Exchange:

[22]
R. Fletcher and S. Leyffer, “Nonlinear programming without a penalty
function,” Mathematical programming, vol. 91, no. 2, pp. 239–269,

[23]
M. H. Wright, “Why a pure primal Newton barrier step may be infeasible,”
SIAM Journal on Optimization, vol. 5, no. 1, pp. 1–12, 1995.

[24]
J. Hauser and A. Saccon, “A barrier function method for the optimization of
trajectory functionals with constraints,” in Proceedings of the 45th
IEEE Conference on Decision and Control.1em plus 0.5em minus
0.4emIEEE, 2006, pp. 864–869.