# Black hole solution and thermal properties in $4D$ $AdS$ Gauss-Bonnet
massive gravity

Sudhaker Upadhyay<EMAIL_ADDRESS>Department of Physics, K. L. S.
College, Nawada, Bihar 805110, India  Department of Physics, Magadh
University, Bodh Gaya, Bihar 824234, India  Inter-University Centre for
Astronomy and Astrophysics (IUCAA), Pune,
Maharashtra 411007, India  School of Physics, Damghan University, P.O. Box
3671641167,
Damghan, Iran  Dharm Veer Singh<EMAIL_ADDRESS>Department of Physics,
Institute of Applied Science and Humanities, GLA University, Mathura, 281406
India.

###### Abstract

We consider an Einstein-Gauss-Bonnet massive gravity model in $4D$ $AdS$
spacetime to obtain a possible black hole solution and discuss the horizon
structure of this black hole. The real roots of the vanishing metric function
lead to various types of horizons. Furthermore, we derive various
thermodynamic quantities, thus insuring the validity of the first-law of
thermodynamics and the Smarr relation. The thermodynamic quantities are
modified in the presence of the massive gravity parameter and also discuss the
stability of the system from the heat capacity. Black hole space times can not
only possess standard thermodynamics, but also possess phase structures when
the cosmological constant is treated as the thermodynamic pressure. The
effects of massive parameter and GB parameter on phase transition are also
discussed. We also examine the phase structure through Maxwell’s equal area
law. The first order phase transition occurs only when pressure is lower than
its critical value. At critical pressure, the first order phase transition
terminates, and the second order phase transition occurs. No phase transition
occurs when pressure is greater than its critical value.

##  I Introduction

The complete understanding of the universe from the principles of general
relativity (GR) is one of the challenging tasks for present day science as GR
is not theoretically complete. The explanation of dark matter and dark energy
present in the Universe is one of them. Various efforts by various people were
made to complete the GR by modifying the Einstein theory of gravity, but the
conclusive theory is still missing. To be more precise, some of the examples
of modified theories of gravity are scalar-tensor theories [1, 2, 3, 4, 5, 6,
7], Lovelock gravity [8, 9, 10] and brane world cosmology [11, 12, 13]. Here,
we place emphasis on the higher derivative of gravity known as the Lovelock
theory, which extends Einstein’s theory to higher dimensions.

The Lovelock’s theorem [8, 9] states that Einstein’s GR with the cosmological
constant is the unique theory of gravity with four-dimensional spacetime,
diffeomorphism invariance, metricity, and second-order equations of motion.
Lovelock’s gravity [8], among the other gravity theories with higher
derivative curvature terms, has some special features. Gauss-Bonnet (GB)
gravity is a particular case of Lovelock gravity in higher dimensions [14]. In
Ref. [15], it is shown that for the GB black holes in $AdS$ spacetime, there
exists a minimal horizon radius below which the Hawking-Page phase transition
will not occur. Recently, the absorption cross section of planar scalar
massless waves impinging on solutions of the $4D$ Einstein-Gauss-Bonnet (EGB)
theory of gravity [16]. The solution depends also on the mass of the black
hole and the GB constant coupling.

One of the possibilities to modify GR by considering a massive graviton also
attracts the attention of various people. For instance, a ghost-free theory
with massive gravitons is developed in Ref. [17]. In curved spacetime, this
causes to the presence of ghost instabilities [18]. The black hole solutions
in the presence of massive gravity with their effects on the geometrical
structure is also studied [19, 20]. Hendi et al. [21] investigated charged
black hole solutions in GB massive gravity and their thermodynamics in $d$
dimensions.

In fact, the GB term in $4D$ is a total derivative and thus does not
contribute to the dynamics. However, the GB term contributes to the dynamics
of the gravitational field only for $d>4$. When the GB term in $4D$ is coupled
to a matter field, it contributes to dynamics. In an interesting work, Glavan
and Lin recently found EGB gravity in four dimensions by rescaling the GB
coupling constant $\alpha\to\alpha/D-4$, where the GB term contributes to the
local dynamics [22]. Indeed, taking the $D\to 4$ limit either breaks part of
the diffeomorphism. However it is possible to construct a Lorents violating
theory by breaking the diffeomorphism and it is invariant only under $3D$
spatial diffeomorphism as the same number of degree of freedom as general
relativity [23, 24]. One can find a consistent $4D$ EGB gravity that is
covariant under spatial diffeomorphism only but not under time diffeomorphism.
The consistent theory validates some results of the solution proposed by
Glavin and Lin [22]. The generlization of other spherically symmetric black
hole solution has been discussed in [25, 26, 27, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 28]. So, there is enough motivation to generalize the
result by studying the black hole solution for the massive EGB gravity in four
dimensions. This is an objective of the present work.

In this paper, we intend to generalize EGB action by adding a massive term and
obtain an exact solution of the massive EGB gravity in $4D$ $AdS$ spacetime.
We will discuss horizon structure and obtain various types of horizons.
Furthermore, we analyse thermal properties and the validity of the first-law
of thermodynamics. We also discuss the role of graviton mass and the GB
parameter in the stability of this black hole solution. Finally, we shed light
on the graviton mass and GB parameter dependence of phase transition and
critical points.

The remaining parts of the paper are organized as follows. In Sec. II, we
obtain a black hole solution and discuss the horizons of the EGB massive
gravity in $4D$ $AdS$ space. In Sec. III, we provide the thermodynamics of the
system. The stability of the black hole is analyzed in section IV. The phase
transition is presented in section V. The paper is summarized with concluding
remarks in the last section.

##  II Black hole solution in $4D$ EGB massive gravity

Let us begin by writing the action for the massive Einstein-Gauss-Bonnet (EGB)
gravity as following:

$\displaystyle S$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
d^{d}x\sqrt{-g}\left[{\cal R}-2\Lambda+{\alpha}{\cal
L_{GB}}+m^{2}\sum_{i}c_{i}{\cal U}_{i}(g,h)\right],$ (1)

where ${\cal R}$ is the curvature scalar, $\alpha$ is a Gauss-Bonnet (GB)
coupling constant with dimension [length]2,
$\mathcal{L_{GB}}:=R_{\mu\nu\gamma\delta}R^{\mu\nu\gamma\delta}-4R_{\mu\nu}R^{\mu\nu}+R^{2}$
is the GB Lagrangian density. Additionally, $m$ represents massive gravity
parameter related to the graviton mass, $h$ refers to the fixed symmetric
tensor, $c_{i}$ are free constants and $\mathcal{U}_{i}(g,f)$ refer to the
symmetric polynomials of the eigenvalues of matrix
$\mathcal{K}_{\nu}^{\mu}=\sqrt{g^{\mu\alpha}h_{\alpha\nu}}$ [33, 32]. The
components of $\mathcal{U}_{i}(g,h)$ is expressed by

$\displaystyle\mathcal{U}_{1}$ $\displaystyle=$
$\displaystyle\left[\mathcal{K}\right],$ $\displaystyle\mathcal{U}_{2}$
$\displaystyle=$
$\displaystyle\left[\mathcal{K}\right]^{2}-\left[\mathcal{K}^{2}\right],$
$\displaystyle\mathcal{U}_{3}$ $\displaystyle=$
$\displaystyle\left[\mathcal{K}\right]^{3}-3\left[\mathcal{K}\right]\left[\mathcal{K}^{2}\right]+2\left[\mathcal{K}^{3}\right],$
$\displaystyle\mathcal{U}_{4}$ $\displaystyle=$
$\displaystyle\left[\mathcal{K}\right]^{4}-6\left[\mathcal{K}^{2}\right]\left[\mathcal{K}\right]^{2}+8\left[\mathcal{K}^{3}\right]\left[\mathcal{K}\right]+3\left[\mathcal{K}^{2}\right]^{2}-6\left[\mathcal{K}^{4}\right],$
(2)

where the $[\ ]$ represents the trace of matrix ${\cal K}_{\nu}^{\mu}$. By
varying the action (1) with respect to metric tensor $g_{\mu\nu}$, the
equation of motion is calculated as

$\displaystyle G_{\mu\nu}+\Lambda g_{\mu\nu}+H_{\mu\nu}+m^{2}\chi_{\mu\nu}=0,$
(3)

where Einstein tensor ($G_{\mu\nu}$), Lanczos tensor ($H_{\mu\nu}$) and
massive tensor ($\chi_{\mu\nu}$) have following explicit expressions,
respectively:

$\displaystyle G_{\mu\nu}$ $\displaystyle=$ $\displaystyle
R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R,$ (4) $\displaystyle H_{\mu\nu}$
$\displaystyle=$
$\displaystyle-\frac{\alpha}{2}\left[8R^{\rho\sigma}R_{\mu\rho\nu\sigma}-4R_{\mu}^{\rho\sigma\lambda}R_{\nu\rho\sigma\lambda}-4RR_{\mu\nu}+8R_{\mu\lambda}R_{\nu}^{\lambda}\right.$
(5) $\displaystyle+$
$\displaystyle\,\,\left.g_{\mu\nu}\left(R_{\mu\nu\gamma\delta}R^{\mu\nu\gamma\delta}-4R_{\mu\nu}R^{\mu\nu}+R^{2}\right)\right],$
$\displaystyle\chi_{\mu\nu}$ $\displaystyle=$
$\displaystyle-\frac{c_{1}}{2}\left(\mathcal{U}_{1}g_{\mu\nu}-\mathcal{K}_{\mu\nu}\right)-\frac{c_{2}}{2}\left(\mathcal{U}_{2}g_{\mu\nu}-2\mathcal{U}_{1}\mathcal{K}_{\mu\nu}+2\mathcal{K}_{\mu\nu}^{2}\right)$
(6) $\displaystyle-$
$\displaystyle\,\,\frac{c_{3}}{2}(\mathcal{U}_{3}g_{\mu\nu}-3\mathcal{U}_{2}\mathcal{K}_{\mu\nu}+6\mathcal{U}_{1}\mathcal{K}_{\mu\nu}^{2}-6\mathcal{K}_{\mu\nu}^{3})$
$\displaystyle-$
$\displaystyle\,\,\frac{c_{4}}{2}(\mathcal{U}_{4}g_{\mu\nu}-4\mathcal{U}_{3}\mathcal{K}_{\mu\nu}+12\mathcal{U}_{2}\mathcal{K}_{\mu\nu}^{2}-24\mathcal{U}_{1}\mathcal{K}_{\mu\nu}^{3}+24\mathcal{K}_{\mu\nu}^{4}).$

To have a static spherically symmetric black hole solution, The spherically
symmetric solution after re-scaling the coupling constant by $\alpha/(D-4)$,
in the limit $D\rightarrow 4$, takes the line element as follows:

$ds^{2}=-f(r)dt^{2}+\frac{1}{f(r)}dr^{2}+r^{2}d\Omega_{d-2},$ (7)

where $d\Omega_{D-2}$ is the metric of a $(d-2)$-dimensional constant
curvature space. Now, we make ansatz for reference metric as

$h_{\mu\nu}=\mbox{diag}(0,0,c^{2},c^{2}\sin^{2}\theta),$ (8)

where $c$ is a positive constant. Using the metric ansatz (8), the components
of $\mathcal{U}_{i}$ (2) identify to [21]

$\displaystyle\mathcal{U}_{1}=\frac{(d-2)c}{r},\qquad\qquad\qquad\mathcal{U}_{3}=\frac{(d-2)(d-3)(d-4)c^{3}}{r^{3}},$
(9)
$\displaystyle\mathcal{U}_{2}=\frac{(d-2)(d-3)c^{2}}{r^{2}}\qquad\qquad\mathcal{U}_{4}=\frac{(d-2)(d-3)(d-4)(d-5)c^{4}}{r^{4}}.$
(10)

The $(r,r)$ components of Einstein field equation Eq. (3) in the limit $d\to
4$ gives

$\displaystyle
r^{5}-2r^{3}\alpha(f-1)f^{\prime}+r^{4}(f-1)+r^{2}\alpha(f-1)^{2}-\Lambda
r^{2}-m^{2}(cc_{1}r+c^{2}c_{2})=0,$ (11)

where prime denotes derivative with respect to $r$. We obtain the following
solution by solving the Eq. (11)

$\displaystyle
f_{\pm}(r)=1+\frac{r^{2}}{2\alpha}\left(1\pm\sqrt{1+4\alpha\left(\frac{2M}{r^{3}}-\frac{1}{l^{2}}-\frac{m^{2}}{2r^{2}}\left(cc_{1}r+2c^{2}c_{2}\right)\right)}\,\right),$
(12)

where $M$ is an integration constant related to the total mass of the black
hole. Here cosmological constant is expressed in terms of scale length factor
$l$ which, in general, takes value $(-3/l^{2})$ for $AdS$ solutions. This is
an exact solution of the field equation (3). In the solution (12), $M$ refers
to the integration constant related to the mass of the black hole, $\alpha$
refers to the Gauss-Bonnet coupling, and dimensionless quantity $m$ is related
to the mass of graviton. For the further analysis, it is convenient to work in
dimensionless parameters:

${\tilde{r}}=\frac{r}{l},\quad{\tilde{M}}=\frac{M}{l},\quad{\tilde{\alpha}}=\frac{\alpha}{l^{2}},\quad\text{and}\quad{\tilde{c}_{1}}=\frac{c_{1}}{l}.$
(13)

The solution (12) can be expressed in terms of the above parameters as

$\displaystyle
f_{\pm}({\tilde{r}})=1+\frac{{\tilde{r}}^{2}}{2{\tilde{\alpha}}}\left(1\pm\sqrt{1+4{\tilde{\alpha}}\left(\frac{2{\tilde{M}}}{{\tilde{r}^{3}}}-\frac{{m}^{2}}{2{\tilde{r}}^{2}}\left(c\tilde{c}_{1}{\tilde{r}}+2c^{2}c_{2}\right)\right)-1}\,\right),$
(14)

The solution (14) behaves asymptotically as

$\displaystyle
f_{-}=1-\frac{2{\tilde{M}}}{\tilde{r}}+{\tilde{r}^{2}}+\frac{m^{2}}{2}(c\tilde{c}_{1}r+c^{2}c_{2})+\frac{\tilde{r}^{2}}{2\tilde{\alpha}}+{{\cal
O}\left(\frac{1}{r^{3}}\right)},$ (15) $\displaystyle
f_{+}=1+\frac{2\tilde{M}}{\tilde{r}}-{\tilde{r}^{2}}-\frac{m^{2}}{2}(c\tilde{c}_{1}{\tilde{r}}+c^{2}c_{2})+\frac{\tilde{r}^{2}}{2\tilde{\alpha}}+{{\cal
O}\left(\frac{1}{\tilde{r}^{3}}\right)}.$ (16)

The $-$ve branch corresponds to the $4D$ $AdS$ EGB massive black hole, whereas
the $+$ve branch does not lead to a physically meaningful solution because the
positive sign in the mass term indicates graviton instabilities, so we will
stick to the $-$ve branch of (14). Furthermore if we take the limit
$\tilde{r}\to\infty$ ($\tilde{M}=0$), the $-$ve branch of the solution (14) is
asymptotically flat whereas the $+$ve branch of the solution (14) is
asymptotically $dS$ ($AdS$) depending on the sign of $\tilde{\alpha}$.

In the massless graviton limit ($m=0$), this solution reduces to

$\displaystyle
f_{\pm}(\tilde{r})=1+\frac{\tilde{r}^{2}}{2\tilde{\alpha}}\left(1\pm\sqrt{1+\frac{8\tilde{M}\tilde{\alpha}}{\tilde{r}^{3}}-{4\tilde{\alpha}}}\,\right).$
(17)

This black hole solution (17) coincides with the solution given by Glavan and
Lin in Ref. [22] when $m=0$ and matches with $AdS$ Schwarzschild massive black
hole solution when $\alpha=0$.

In order to study the horizon structure, we plot these thermodynamics
quantities in Fig. 1 as the function of horizon radius with different value of
massive gravity parameter with specific GB parameter $\tilde{\alpha}=0.1$ and
$\tilde{\alpha}=0.2$. The Fig. 1 describes the horizon structure of the $4D$
EGB massive black hole. The numerical analysis of $f(r_{+})=0$ reveals that it
is possible to find non-vanishing value of $\alpha$, $m$ and cosmological
constant ($\Lambda$) for which metric function is minimum. The metric function
$f(r_{+})=0$ gives three real roots $r_{+}$ (Cauchy horizon), $r_{-}$ (event
horizon) and $r_{c}$ (cosmological horizon).

---
Figure 1: Metric function plot $f(r)$ as a function of horizon radius with the variation of massive gravity parameter $({m})$ for GB parameter ${\tilde{\alpha}}=0.1$ and ${\tilde{\alpha}}=0.2$. Here, $c=1,\tilde{c}_{1}=-1,c_{2}=1$.  | $\tilde{\alpha}=0.1$ |  |  |  |  | $\tilde{\alpha}=0.2$ |
---|---|---|---|---|---|---|---
$m$ | $r_{-}$ | $r_{+}$ | $\delta$ | $m$ | $r_{-}$ | $r_{+}$ | $\delta$
0 | 0.1344 |  1.934 | 1.799 | $0$ | 0.0963 |  1.865 | 1.768
1.50 | 0.3685 |  1.715 | 1.346 | $1.50$ | 0.3716 |  1.5 | 1.128
2.09 | 0.3539 |  1.429 | 0.5996 | 1.73 | 0.5023 |  1.328 | 0.8257
2.50 | 0.3539 |  1.429 | 0.5996 | 1.90 | 0.5023 |  1.328 | 0.8257
Table 1: Cauchy ($\tilde{r}_{-}$) and event ($\tilde{r}_{+}$) horizons, and
$\delta=\tilde{r}_{+}-\tilde{r}_{-}$ for the $4D$ AdS EGB massive black hole
for GB parameter ${\tilde{\alpha}}=0.1$ and ${\tilde{\alpha}}=0.2$. Here,
$c=1,\tilde{c}_{1}=-1,c_{2}=1$.

From Fig. 1, it is clear that the size of the black hole increases with a
decrease in the GB parameter. The black hole has three (Cauchy, event, and
cosmological) horizons at $m>2.09$ with $\tilde{\alpha}=0.1$. However, the
event and cosmological horizons are possible for $m<2.09$. The horizon is a
decreasing function of massive parameter ($m$) (Fig. 1 right panel), an
increasing function of massive gravity parameter and finally independent of GB
parameter (Fig. 1 right panel). It should be pointed out that interestingly,
increasing and decreasing mentioned parameters, will lead to formation of the
different sized black hole.

##  III Thermodynamics

In this section, we study the thermodynamic properties of $4D$ $AdS$ Gauss-
Bonnet massive black hole. In particular, these are mass $(\tilde{M}_{+})$,
temperature $(\tilde{T}_{+})$, entropy $(\tilde{S}_{+})$ and free energy
$(\tilde{G}_{+})$. The mass of the black hole is derived by setting
$f_{-}(\tilde{r})|_{\tilde{r}=\tilde{r}_{+}}=0$ as

$\displaystyle\tilde{M}_{+}=\frac{1}{2\tilde{r}_{+}}\left[\tilde{r}_{+}^{2}+\tilde{\alpha}+{\tilde{r}_{+}^{4}}+{m^{2}\tilde{r}_{+}^{2}}\left(\frac{c\tilde{c}_{1}\tilde{r}_{+}}{2}+2c^{2}c_{2}\right)\right].$
(18)

The Hawking temperature for this black hole solution can be estimated from the
well-known area-law

$\tilde{T}_{+}=\frac{1}{4\pi}f^{\prime}(\tilde{r})|_{\tilde{r}=\tilde{r}_{+}},$
(19)

as follows

$\displaystyle\tilde{T}_{+}=T_{+}l=\frac{1}{4\pi\tilde{r}_{+}(r_{+}^{2}+2\tilde{\alpha})}\left[\tilde{r}_{+}^{2}-\tilde{\alpha}+{3\tilde{r}_{+}^{4}}+m^{2}\tilde{r}_{+}^{2}(c^{2}c_{2}+c\tilde{c}_{1}\tilde{r})\right].$
(20)

---
Figure 2: Temperature plot $(\tilde{T}_{+})$ as a function of horizon radius
with the variation of massive gravity parameter $(m)$ for GB parameter
$\tilde{\alpha}=0.1$ and $\tilde{\alpha}=0.2$. Here,
$c=1,\tilde{c}_{1}=-1,c_{2}=1$.

It is evident Fig. 2 that the temperature is negative when the two (Cauchy and
event) horizons merge, therefore, in this region, the obtained solutions are
not physical.

The thermodynamical quantities must follow the first law of thermodynamics

$d\tilde{M}_{+}=\tilde{T}d\tilde{S}_{+}.$ (21)

This leads to explicit expression for entropy
$\tilde{S}_{+}=\int\frac{d\tilde{M}_{+}}{\tilde{T}_{+}}d\tilde{r}_{+}$ as [28,
45]

$\displaystyle\tilde{S}_{+}=\pi\tilde{r}_{+}^{2}+4\pi\tilde{\alpha}\log[\tilde{r}]+\tilde{S}_{0},$
(22)

where $\tilde{S}_{0}$ is an integration constant.

It is worth mentioning that the definition of first-law of thermodynamics may
have additional terms also in the extended phase space. In order to check the
validity of these terms one should check the first-law in the non-differential
form called as the Smarr relation

$\tilde{M}_{+}=2\tilde{T}_{+}\tilde{S}_{+}+2\tilde{V}_{+}\tilde{P}-\tilde{A}\tilde{\alpha}-{\cal\tilde{C}}\tilde{c}_{1},$
(23)

where

$\displaystyle\tilde{P}=-\frac{\Lambda}{8\pi},\qquad\qquad\qquad\qquad\quad\tilde{V}_{+}=\left(\frac{\partial\tilde{M}_{+}}{\partial\tilde{P}}\right)_{\tilde{S}_{+},\tilde{\alpha},\tilde{c}_{1}},$
(24)
$\displaystyle{\cal\tilde{C}}=\left(\frac{\partial\tilde{M}_{+}}{\partial\tilde{c}_{1}}\right)_{\tilde{S}_{+},\tilde{\alpha},\tilde{P}},\qquad\qquad\qquad\tilde{A}=\left(\frac{\partial\tilde{M}_{+}}{\partial\tilde{\alpha}}\right)_{\tilde{S}_{+},\tilde{P},\tilde{c}_{1}}.$
(25)

From all the above additional thermodynamic quantities, it is easy to verify
the complete form of the first-law of thermodynamics in an extended space

$d\tilde{M}_{+}=\tilde{T}_{+}d\tilde{S}_{+}-\tilde{V}_{+}d\tilde{P}+\tilde{A}d\tilde{\alpha}+{\cal\tilde{C}}dc_{1}.$
(26)

To understand the global stability of the black hole, the behavior of the
Gibbs free energy $\tilde{G}_{+}$ is important to analyzed. The Gibbs free
energy, from the definition
$\tilde{G}_{+}=\tilde{M}_{+}-\tilde{T}_{+}\tilde{S}_{+}$, is derived by

$\displaystyle\tilde{G}_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{2\tilde{r}_{+}}\left[\tilde{r}_{+}^{2}+2\tilde{\alpha}+{\tilde{r}_{+}^{3}}+m^{2}\left(\frac{c\tilde{c}_{1}\tilde{r}_{+}^{2}}{2}+c^{2}c_{2}\tilde{r}_{+}\right)\right]$
(27) $\displaystyle-$
$\displaystyle\frac{(\tilde{r}_{+}^{2}+4\alpha\log[\frac{\tilde{r}}{\tilde{r}_{0}}])}{4\tilde{r}_{+}(\tilde{r}_{+}^{2}+2\tilde{\alpha})}\left[\tilde{r}_{+}^{2}-\tilde{\alpha}+{3\tilde{r}_{+}^{4}}+m^{2}\tilde{r}_{+}^{2}(c^{2}c_{2}+c\tilde{c}_{1}\tilde{r})\right].$

---
Figure 3: Plot of free energy $(\tilde{G}_{+})$ as a function of horizon
radius with the variation of massive gravity parameter $(m)$ for GB parameter
$\tilde{\alpha}=0.1$ and $\tilde{\alpha}=0.2$. Here, $c=1,c_{1}=-1,c_{2}=1$.

To analyze Gibbs free energy, we plot figure 3 for different values of massive
gravity parameter and GB parameter. Here, we see that Gibbs free energy
increases with GB parameter. However, with increasing $m$ the Gibbs free
energy first increases and then starts falling after a critical point.

##  IV Heat Capacitity and thermodynamics stability

Now, we study the heat capacity of the obtained black hole solution (12) in
the context of canonical ensemble by calculating the heat capacity. The heat
capacity is given by the following relation:

$\tilde{C}_{+}=\left(\frac{\partial\tilde{M}_{+}}{\partial\tilde{T}_{+}}\right)_{\tilde{\alpha}}=\left(\frac{\tilde{\partial}M_{+}}{\tilde{\partial}r_{+}}\right)_{\tilde{\alpha}}\left(\frac{\partial\tilde{r}_{+}}{\partial\tilde{T}_{+}}\right)_{\tilde{\alpha}}.$
(28)

Substituting the value of mass from (18) and temperature from (20) in Eq.
(28), we get

$\tilde{C}_{+}=\frac{{2\pi\tilde{r}^{2}(\tilde{r}^{2}+2\tilde{\alpha})^{2}}\left(1+c^{2}c_{2}\tilde{m}^{2}+c\tilde{c}_{1}m^{2}r+{3\tilde{r}^{2}}-\frac{\tilde{\alpha}}{\tilde{r}^{2}}\right)}{3(\tilde{r}^{6}+6\tilde{r}^{4}\tilde{\alpha})+[-(1+c^{2}c_{2}m^{2})\tilde{r}^{4}+(5+2c^{2}c_{2}m^{2})\tilde{r}^{2}\tilde{\alpha}+4c\tilde{c}_{1}m^{2}\tilde{r}^{3}\tilde{\alpha}+2\tilde{\alpha}^{2}]}.$
(29)

---
Figure 4: Plot of heat capacity ($\tilde{C}_{+}$) as a function of horizon
radius with the variation of massive gravity parameter $(m)$ for GB parameter
$\tilde{\alpha}=0.2$ with fixed value of $c=1,\tilde{c}_{1}=-1,c_{2}=1$.

In order to study the thermodynamical behavior of the system stability and
phase transition points, we plot this thermodynamics quantity in Fig. 4. For
Eq. (29), the heat capacity will have divergences at extremum points.

We know that roots and divergences present in the heat capacity represent
phase transition points of the system. Additionally, we should mention that
unstable systems undergo a phase transition and becomes stable. Corresponding
to the root there exists a phase transition from unstable (non-physical)
system to stable (physical) one. For small divergence a phase transition from
large black holes to small ones takes place and for large divergence a phase
transition occurs from smaller black holes to larger ones.

##  V  Phase transition of 4D EGB Massive black hole

In the extended phase space, the mass of the black hole treated as the
enthalpy rather then the internal energy of the system. The cosmological
constant is treated as the thermodynamics pressure $(\Lambda=-8\pi P)$ [42,
43] and corresponding conjugate quantity is volume which is given by

$\tilde{V}_{+}=\frac{4}{3}\pi\tilde{r}_{+}^{3}.$ (30)

The equations of state (EOS) can obtained by relation (23) and Hawking
temperature as

$\tilde{P}_{+}=\frac{\tilde{T}_{+}}{2\tilde{r}_{+}}+\frac{\tilde{T}_{+}\alpha}{\tilde{r}_{+}^{3}}+\frac{\tilde{\alpha}}{8\pi\tilde{r}_{+}^{4}}-\frac{1}{8\pi\tilde{r}_{+}^{2}}-\frac{m^{2}}{8\pi\tilde{r}_{+}^{2}}\left(c\tilde{c}_{1}\tilde{r}+c^{2}c_{2}\right),\qquad\qquad\tilde{v}=2\tilde{r}_{+},$
(31)

where $\tilde{v}$ represents specific volume. Using the properties of the
inflection points we can encode the information about the EoS in the critical
points. The critical point occurs when $\tilde{P}_{+}$ has an inflection point

$\frac{\partial\tilde{P}_{+}}{\partial\tilde{v}}=\frac{\partial^{2}\tilde{P}_{+}}{\partial\tilde{v}^{2}}=0.$
(32)

The corresponding critical point is determined by the following equation:

$\tilde{r}_{+}^{4}-12\tilde{\alpha}(\tilde{r}_{+}^{2}+\tilde{\alpha})+m^{2}\tilde{r}_{+}^{2}(c^{2}c_{2}\tilde{r}_{+}^{2}-6\tilde{\alpha}(c^{2}c_{2}+c\tilde{c}_{1}\tilde{r}_{+}))=0.$
(33)

$m$ | $\tilde{r}_{C}$ | $\tilde{T}_{C}$ | $\tilde{P}_{C}$ | $\frac{\tilde{P}_{C}\,\tilde{v}_{C}}{\tilde{T}_{C}}$
---|---|---|---|---
0 | 1.1370 | 0.0808 | 0.0126 | 0.3546
1 | 0.8525 | 0.1327 | 0.0440 | 0.5653
2 | 0.6714 | 0.3267 | 0.1666 | 0.6847
3 | 0.6050 | 0.6770 | 0.3920 | 0.7006
4 | 0.5757 | 1.1598 | 0.7164 | 0.7112
5 | 0.5606 | 1.7913 | 1.1373 | 0.7118
$\tilde{\alpha}$ |  |  |  |
0.1 | 0.8525 | 0.1327 | 0.0440 | 0.5653
0.2 | 1.1403 | 0.0214 | 0.0251 | 0.7028
0.3 | 1.3397 | 0.0590 | 0.0184 | 0.8356
0.4 | 1.4951 | 0.0459 | 0.0150 | 0.9771
0.5 | 1.6231 | 0.0370 | 0.0129 | 1.1317
Table 2: The table for critical temperature $\tilde{T}_{C}$, critical pressure
$\tilde{P}_{C}$ and ${\tilde{P}_{C}\,\tilde{v}_{C}}/{\tilde{T}_{C}}$
corresponding different value of $m$ with fixed value of
$c=1,\tilde{c}_{1}=-0.75,c_{2}=0.75$ and $\tilde{\alpha}=0.1$.

We find that Eq. (33) can not be solve analytically. In order to estimate the
critical points, namely, critical horizon radius ($\tilde{r}_{C}$), critical
temperature ($\tilde{T}_{+}$) and critical pressure $(\tilde{P}_{+})$ the
numerical solutions of Eq. (33) with variation of massive parameter and and GB
coupling are tabulated in the Table (II) and Table (III). It is interesting to
see that the critical pressure and critical temperature increase with
increasing critical radius and massive parameter at $\tilde{\alpha}=0.1$.
However, the critical pressure and critical temperature decrease with
increasing critical radius and GB coupling parameter for $m=1$.

---
Figure 5: The plots of temperature vs horizon radius for $\tilde{\alpha}=0.1$
(left) and $\tilde{\alpha}=0.2$ (right) with fixed values of
$c=1,\tilde{c}_{1}=-1,c_{2}=1$ and $m=1$ .

We study the phase transition by analyzing the Gibbs free energy for different
values of massive gravity parameter $(m)$ and GB coupling $(\alpha)$ as
depicted from Fig. 6. The characteristic swallow tail confirms that the
obtained values are the critical ones in which the phase transition takes
place. The thermodynamically stable (unstable) state occurs corresponding to
the lowest (highest) Gibbs free energy. As a result, the triangular loop in
the $\tilde{G}_{+}-\tilde{T}_{+}$ plane represents an unstable state. The
straight line corresponds to the first order phase transition as the
temperature rises in the system, as shown in the $\tilde{T}_{+}-\tilde{r}_{+}$
plane (Fig. 7). The curved portion of the isobar indicates the unstable state
due to the system’s possessing lower Gibbs free energy. We note that the net
change in Gibbs free energy around the triangular loop is zero.

We can also see that swallow tail shape occurs in case of
$\tilde{P}_{+}<\tilde{P}_{c}$ for the first order phase transition. However,
in case of $\tilde{P}_{+}=\tilde{P}_{c}$, the second order phase transition
occurs. The intercepts of the axis in the $\tilde{G}_{+}-\tilde{T}_{+}$ plane
decrease as the GB coupling parameter and mass gravity parameter increase,
whereas the intercepts of the $\tilde{G}$ -axis increase as the GB coupling
parameter and mass gravity parameter increase.

---
Figure 6: The plots of Gibbs free energy ($\tilde{G}_{+}$) vs temperature
($\tilde{T}_{+}$) for $\alpha=0.1$ (upper left), $\alpha=0.2$ (upper right),
$m=1.0$ (lower left) and $m=2.0$ (lower right) with fixed values of
$c=1,\tilde{c}_{1}=-0.7,c_{2}=0.75$ .
---
Figure 7: The plots of temperature ($\tilde{T}_{+}$) vs horizon radius
$(\tilde{r}_{+})$ and Gibbs free energy ($\tilde{G}_{+}$) vs temperature
($\tilde{T}_{+}$) with fixed values of $\tilde{\alpha}=0.1$ and for $m=1.0$
and $m=2.0$ $c=1,\tilde{c}_{1}=-1,c_{2}=1$.

When the temperature is lower than the critical temperature in the
$\tilde{T}_{+}-\tilde{r}_{+}$ plot, we can see three black holes (small,
intermediate, and large) for the $4DAdS$ EGB massive black hole with the same
massive gravity parameter of $m$ for a specific temperature range. The small
and large black holes are stable, but the intermediate black hole is unstable,
since the heat capacity $\tilde{C}_{+}$ is negative (see Fig. 4). Because of
the small free energy, $\tilde{T}_{+}<T_{\star}$ corresponds to a small black
hole, and $\tilde{T}_{+}>T_{\star}$ corresponds to a large black hole, where
$T_{\star}$ is the transition temperature, which has values of $0.088$ and
$0.048$ for $\tilde{\alpha}=0.1$ and $\tilde{\alpha}=0.2$, respectively. Due
to the same free energy, black holes can transit from one phase to another
phase at a critical temperature. In Fig. 7, point 5 represents the coexistence
of small and large black holes. The small black holes are represented by the
solid lines $1-L5$ or $1-5$, while the large black holes are represented by
the solid lines $R5-4$ or $5-4$. We noticed that the point $L5$ and $R5$ share
the same free energy.

##  VI Curvature Singularity and Phase Transition new thermodynamic geometry

we have proposed a new formalism of the thermodynamic geometry (NTG) which
explains the one-to-one correspondence between phase transitions and
singularities of the scalar curvature [44, 45].The NTG geometry is defined as

$dl^{2}_{NTG}=\frac{1}{T}\left(n^{j}_{i}\frac{\partial^{2}\Xi}{\partial
X^{i}\partial X^{j}}dX^{i}dX^{j}\right),$ (34)

where $(n^{j}_{i}=\text{diag}(-1,1,1,\ldots 1)$ and $Xi$ is thermodynamic
potential [44] and the geometrothermodynamics (GTD) metric is conformally
related to NTG metric such that this conformal transformation is singular at
unphysical points were generated in GTD metric [44].

To study the curvature singularity and phase transition in NTG. We write the
free energy in term of specific volume and temperature

$\displaystyle\tilde{G}_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{\tilde{v}_{+}}\left[\frac{\tilde{v}_{+}^{2}}{4}+2\tilde{\alpha}+\frac{\tilde{v}_{+}^{3}}{8}+m^{2}\left(\frac{c\tilde{c}_{1}\tilde{v}_{+}^{2}}{8}+\frac{c^{2}c_{2}\tilde{v}_{+}}{2}\right)-(\frac{\tilde{v}_{+}^{2}}{4}+4\alpha\log[\frac{\tilde{r}}{2}])\tilde{T}\right].$
(35)

The differential form of free energy
$d\tilde{F}_{+}=-\tilde{S}_{+}d\tilde{T}_{+}-\tilde{P}d\tilde{V}_{+}+\tilde{A}d\tilde{\alpha}+\tilde{\cal
C}d\tilde{c}_{1}$ combining with Eq. (35), we deduce

$\displaystyle\tilde{S}_{+}=-\left(\frac{\partial\tilde{F}_{+}}{\partial\tilde{T}_{+}}\right)_{\tilde{V},\tilde{\alpha},{\cal\tilde{C}}}=\frac{\pi}{4}\left[\tilde{v}^{2}+16\tilde{\alpha}\log\left(\frac{\tilde{v}}{\tilde{v}_{0}}\right)\right],$
(36)
$\displaystyle\tilde{P}=-\left(\frac{\partial\tilde{F}_{+}}{\partial\tilde{V}_{+}}\right)_{\tilde{V},\tilde{\alpha},{\cal\tilde{C}}}=\frac{\tilde{T}_{+}}{\tilde{v}_{+}}+\frac{8\tilde{v}_{+}\alpha}{\tilde{r}_{+}^{3}}+\frac{2\tilde{\alpha}}{\pi\tilde{v}_{+}^{4}}-\frac{1}{2\pi\tilde{v}_{+}^{2}}-\frac{m^{2}}{2\pi\tilde{v}_{+}^{2}}\left(\frac{c\tilde{c}_{1}\tilde{v}_{+}}{2}+c^{2}c_{2}\right),$
(37)
$\displaystyle\tilde{A}=-\left(\frac{\partial\tilde{F}_{+}}{\partial\tilde{\alpha}}\right)_{\tilde{V},\tilde{\alpha},{\cal\tilde{C}}}=\frac{1}{\tilde{v}}-4\pi\tilde{T}_{+}\log\left(\frac{\tilde{v}}{\tilde{v}_{0}}\right),$
(38)
$\displaystyle{\cal\tilde{C}}=-\left(\frac{\partial\tilde{F}}{\partial\tilde{c}_{1}}\right)_{\tilde{V},\tilde{\alpha},{\cal\tilde{C}}}=\frac{m^{2}\tilde{c}_{1}\tilde{v}^{2}}{16}.$
(39)

The phase transition is appear when a heat capacity diverges, when
$\tilde{C}_{+}>0$ the system is stable and $\tilde{C}_{+}<0$ the system is
unstable. The heat capacity is written as

$\tilde{C}_{+}=\frac{{2\pi\tilde{r}^{2}(\tilde{r}^{2}+2\tilde{\alpha})^{2}}\left(1+c^{2}c_{2}m^{2}+c\tilde{c}_{1}m^{2}r+{3\tilde{r}^{2}}-\frac{\tilde{\alpha}}{\tilde{r}^{2}}\right)}{3(\tilde{r}^{6}+6\tilde{r}^{4}\tilde{\alpha})+[-(1+c^{2}c_{2}m^{2})\tilde{r}^{4}+(5+2c^{2}c_{2}m^{2})\tilde{r}^{2}\tilde{\alpha}+4c\tilde{c}_{1}m^{2}\tilde{r}^{3}\tilde{\alpha}+2\tilde{\alpha}^{2}]}.$
(40)

Now, we are able to implement the NTG geometry to analysis the phase
transition behaviour of heat capacity (40). By substituting the
$\Xi=\tilde{H}=\tilde{M}_{+}+\tilde{E}+\tilde{P}\tilde{V}_{+}$ with
$X^{i}=(\tilde{S}_{+},\tilde{P},\tilde{\alpha},{\cal\tilde{C}})$ into Eq.
(34), we have

$g_{H}^{NTG}=\begin{bmatrix}-\left(\frac{\partial\tilde{T}_{+}}{\partial\tilde{S}_{+}}\right)_{\tilde{P},\tilde{c}_{1},\tilde{\alpha}}&0&0&0\\\
0&\left(\frac{\partial\tilde{V}_{+}}{\partial\tilde{S}_{+}}\right)_{\tilde{S}_{+},\tilde{c}_{1},\tilde{\alpha}}&\left(\frac{\partial\tilde{V}_{+}}{\partial\tilde{c}_{1}}\right)_{\tilde{S}_{+},\tilde{P},\tilde{\alpha}}&\left(\frac{\partial\tilde{A}}{\partial\tilde{P}}\right)_{\tilde{S}_{+},\tilde{c}_{1},\alpha}\\\
0&\left(\frac{\partial\tilde{V}_{+}}{\partial\tilde{C}}\right)_{\tilde{S}_{+},\tilde{P},\tilde{\alpha}}&\left(\frac{\partial\tilde{C}}{\partial\tilde{c}_{1}}\right)_{\tilde{S}_{+},\tilde{P},\tilde{\alpha}}&\left(\frac{\partial\tilde{A}}{\partial\tilde{c}_{1}}\right)_{\tilde{S}_{+},\tilde{P},\tilde{\alpha}}\\\
0&\left(\frac{\partial\tilde{A}}{\partial\tilde{P}}\right)_{\tilde{S}_{+},\tilde{c}_{1},\tilde{\alpha}}&\left(\frac{\partial\tilde{A}}{\partial\tilde{c}_{1}}\right)_{\tilde{S}_{+},\tilde{P},\tilde{\alpha}}&\left(\frac{\partial\tilde{A}}{\partial\tilde{\alpha}}\right)_{\tilde{S}_{+},\tilde{c}_{1},\tilde{P}}\end{bmatrix}$

.

Since all thermodynamic parameters are written as a function of
$(\tilde{T}_{+},\tilde{v},\tilde{c}_{1},\tilde{\alpha})$, it is convenient to
recast metric elements from the coordinate
$X_{i}=(\tilde{S},\tilde{P},\tilde{c}_{1},\tilde{\alpha})$ to the coordinate
$(\tilde{T}_{+},\tilde{v},\tilde{c}_{1},\tilde{\alpha})$ by using the
following Jacobians

$J=\frac{\partial(\tilde{S}_{+},\tilde{P},{\cal\tilde{C}},\tilde{\alpha})}{(\tilde{T}_{+},\tilde{v},{\cal\tilde{C}},\tilde{\alpha})}=\begin{bmatrix}0&\pi\left(\frac{\tilde{v}}{2}+\frac{4\tilde{\alpha}}{\tilde{v}}\right)&0&{4\pi\log\left(\frac{\tilde{v}}{\tilde{v}_{0}}\right)}\\\
\frac{\tilde{v}^{2}+8\tilde{\alpha}}{\tilde{v}^{3}}&\frac{\tilde{v}^{2}-\pi\tilde{T}_{+}\tilde{v}^{3}-8\tilde{(}1+3\pi\tilde{T}\tilde{v})}{\pi\tilde{v}^{5}}+\frac{m^{2}(4c^{2}c_{2}+c\tilde{c}_{1}\tilde{v}))}{4\pi\tilde{v}^{3}}&-\frac{m^{2}}{2\pi\tilde{v}_{+}^{2}}\left(\frac{\tilde{c}_{1}\tilde{v}}{2}+2cc_{2}\right)&\frac{2(1+4\pi\tilde{T}\tilde{v})}{\pi\tilde{v}^{4}}\\\
0&0&1&0\\\ 0&0&0&1\end{bmatrix}$

Finally the metric element convert to

${\hat{g}}=J^{T}g_{H}^{NTG}J=\begin{bmatrix}0&0&0&-\frac{4\pi}{\tilde{T}_{+}}\log\left(\frac{\tilde{v}}{\tilde{v}_{0}}\right)\\\
0&\frac{\tilde{v}^{2}-\pi\tilde{T}_{+}\tilde{v}^{3}-8\tilde{(}1+3\pi\tilde{T}_{+}\tilde{v})}{\pi\tilde{T}_{+}\tilde{v}^{3}}+\frac{m^{2}(4c^{2}c_{2}+c\tilde{c}_{1}\tilde{v}))}{4\pi\tilde{T}_{+}\tilde{v}}&0&0\\\
0&0&\frac{2}{\tilde{T}_{+}\tilde{v}}&0\\\
-\frac{4\pi}{\tilde{T}_{+}}\log\left(\frac{\tilde{v}}{\tilde{v}_{0}}\right)&0&0&0\end{bmatrix}$

. The denominator of the scalar curvature is

$D(R^{NTG})=\pi\left(\frac{\tilde{v}^{2}-\pi\tilde{T}_{+}\tilde{v}^{3}-8(1+3\pi\tilde{T}_{+}\tilde{v})}{\pi\tilde{v}^{5}}+\frac{m^{2}(4c^{2}c_{2}+c\tilde{c}_{1}\tilde{v}))}{4\pi\tilde{v}^{3}}\right)^{2}\log\left(\frac{\tilde{v}}{\tilde{v}_{0}}\right)$
(41)

---
Figure 8: The plots of heat capacity and scalar curvature vs specific volume
for $\tilde{\alpha}=0.1$, $m=2.0$ with fixed values of
$c=1,\tilde{c}_{1}=-1.0,c_{2}=1.0$ for
$\tilde{T}_{+}<\tilde{T}_{c},\tilde{T}_{+}=\tilde{T}_{c}$ and
$\tilde{T}_{+}>\tilde{T}_{c}$, where $\tilde{T}_{c}=0.3267$.

The plot of scalar curvature and heat capacity with respect to specific volume
are depicted Fig. 8. We observe that from the Fig. 8, there is two divergence
for heat capacity at $\tilde{T}_{+}<\tilde{T}_{c}$ and there are three
possible phases i.e., the small black hole (SBH), intermediate black hole
(IBH) and large black hole (LBH). At $\tilde{T}_{+}=\tilde{T}_{c}$, these
divergent points get closer and coincide to form a single divergent point. We
noticed that at $\tilde{T}_{+}>\tilde{T}_{c}$ the heat capacity is positive
and there is no phase transition, it means that the black hole is stable. The
scalar curvature is positive in the range of $0<\tilde{v}<\tilde{v}_{0}=1.7$
that implies a repulsive interaction between the microscopic black hole
molecules.

##  VII Conclusions

We have found a singular solution to the massive EGB gravity in $4D$ $AdS$
spacetime. In fact, we obtain both positive and negative branches of solutions
in which one is not physically correct and therefore can be dropped. The
horizon structure of the black hole is also discussed. We have found that the
real roots of the physically correct metric lead to three horizons (Cauchy,
event, and cosmological horizons) below a certain graviton mass limit.
However, only two horizons (event and cosmological horizons) are possible
after this mass limit.

The thermal properties of this black hole solution are also discussed in a
standard way. We derived the mass of the black hole by setting the metric
function to zero. The Hawking temperature is calculated from the well-known
area-law. We have seen that the black hole solution is not physical when
Cauchy and event horizons merge. We have further checked the validity of the
first law of black hole chemistry. Gibbs energy is also derived. In order to
check the stability of the black hole, we have computed its heat capacity.
Here we have one root and two divergent points. Finally, we have described the
effects of massive graviton and GB parameters on the phase transition and
critical points of the system. We have observed that the critical pressure and
critical temperature increase with increasing critical radius and graviton
mass. In contrast, the critical pressure and critical temperature decrease as
the critical radius and GB parameter increase. In fact, black holes undergo a
phase transition at a critical temperature. In the Gibbs free energy plots we
have seen characteristic swallow tail which locates the critical point where
phase transition takes place.

## References

  * [1] C. Barrabes and G.F. Bressange, Class. Quant. Grav. 14 (1997) 805.
  * [2] R.-G. Cai and Y.S. Myung, Phys. Rev. D 56 (1997) 3466.
  * [3] S. Capozziello and A. Troisi, Phys. Rev. D 72 (2005) 044022.
  * [4] T.P. Sotiriou, Class. Quant. Grav. 23 (2006) 5117.
  * [5] J.W. Moffat, JCAP 03 (2006) 004.
  * [6] V. Faraoni, Phys. Rev. D 75 (2007) 067302.
  * [7] K.-i. Maeda and Y. Fujii, Phys. Rev. D 79 (2009) 084026.
  * [8] D. Lovelock, J. Math. Phys. 12 (1971) 498.
  * [9] D. Lovelock, J. Math. Phys. 13 (1972) 874.
  * [10] N. Deruelle and L. Farina-Busto, Phys. Rev. D 41 (1990) 3696.
  * [11] J.M. Cline and H. Firouzjahi, Phys. Rev. D 64 (2001) 023505.
  * [12] T. Nihei, N. Okada and O. Seto, Phys. Rev. D 71 (2005) 063535.
  * [13] M. Demetrian, Gen. Rel. Grav. 38 (2006) 953.
  * [14] C. Lanczos, Ann. Math. 39, 842 (1938).
  * [15] R. Cai, Phys. Rev. D 65, 084014 (2002).
  * [16] H. C.D. Lima Junior, C. L.Benone and L. C.B. Crispino, Phys. Lett. B 811 (2020) 135921.
  * [17] M. Fierz, Helv. Phys. Acta 12 (1939) 3.
  * [18] D. G. Boulware and S. Deser, Phys. Rev. D 6 (1972) 3368.
  * [19] S. H. Hendi, S. Panahiyan, S. Upadhyay and B. Eslam Panah, Phys. Rev. D 95, 084036 (2017).
  * [20] S. Upadhyay, S. H. Hendi, S. Panahiyan and B. Eslam Panah, Prog. Theor. Exp. Phys. 093E01, 1-20 (2018).
  * [21] S. H. Hendi, S. Panahiyan and B. Eslam Panah, JHEP 01, 129 (2016).
  * [22] D. Glavan and C. Lin, Phys. Rev. Lett. 124, 081301 (2020).
  * [23] K. Aolki, M. A. Gorji and S. Mukohyama, Phys. Lett. B 810 (2020) 135843.
  * [24] K. Aolki, M. A. Gorji and S. Mukohyama, JCAP 2009 (2020) 014.
  * [25] P. G. S. Fernandes, Phys. Lett. B 805 (2020), 135468.
  * [26] R. A. Hennigar, D. Kubizňák, R. B. Mann and C. Pollack, JHEP 07 (2020) 027.
  * [27] S. G. Ghosh, D. V. Singh, R. Kumar and S. D. Maharaj, Annals Phys. 424 (2021), 168347.
  * [28] S. W. Wei and Y. X. Liu, Phys. Rev. D 101 (2020) no.10, 104018.
  * [29] D. V. Singh, S. G. Ghosh and S. D. Maharaj, Phys. Dark Univ. 30 (2020) 100730.
  * [30] D. V. Singh and S. Siwach, Phys. Lett. B. 408 135658 (2020).
  * [31] D. V. Singh, B. K. Singh and S. Upadhyay, Annals Phys. 434 (2021), 168642.
  * [32] B. K. Singh, R. P. Singh and D. V. Singh, Eur. Phys. J. Plus 136 (2021) no.5, 575.
  * [33] B. K. Singh, R. P. Singh and D. V. Singh, Eur. Phys. J. Plus 135 (2020) no.10, 862.
  * [34] N. Godani, D. V. Singh and G. C. Samanta, Phys. Dark Univ. 35 (2022), 100952.
  * [35] S. H. Hendi, M. Momennia, J. High Energ. Phys. 10 (2019) 207.
  * [36] B. Eslam Panah, Kh. Jafarzade, S. H. Hendi, Nuclear Physics B 961 (2020) 115269.
  * [37] S. H. Hendi, M. Momennia, Physics Letters B 777 (2018) 222.
  * [38] S. H. Hendi, B. Eslam Panah, S. Panahiyan, M. S. Talezadeh, Eur. Phys. J. C 77 (2017) 133.
  * [39] S. H. Hendi, B. Eslam Panah, S. Panahiyan ,Class. Quant. Grav. 33 (2016) 235007\.
  * [40] S. H. Hendi, S. Panahiyan, Phys. Rev. D 90 (2014) 124008.
  * [41] S. H. Hendi, B. E. Panah and S. Panahiyan, Fortsch. Phys. 66 (2018) no.3, 180000.
  * [42] D. Kubiznak and R.B. Mann, JHEP 07 (2012) 033.
  * [43] S. Gunasekaran, R.B. Mann and D. Kubiznak, JHEP 11 (2012) 110.
  * [44] S. A. Hosseini Mansoori and B. Mirza, Phys. Lett. B 799 (2019), 135040
  * [45] S. A. Hosseini Mansoori, Phys. Dark Univ. 31 (2021), 100776


# Feedback Linearization for Uncertain Systems via Reinforcement Learning

Tyler Westenbroek*, David Fridovich-Keil*, Eric Mazumdar*, Shreyas Arora,
Valmik Prabhu,
S. Shankar Sastry, and Claire J. Tomlin  This work was supported by HICON-
LEARN (design of HIgh CONfidence LEARNing-enabled systems), Defense Advanced
Research Projects Agency award number FA8750-18-C-0101, and Provable High
Confidence Human Robot Interactions, Office of Naval Research award number
N00014-19-1-2066. EECS, UC Berkeley<EMAIL_ADDRESS>
∗ indicates equal contribution.

###### Abstract

We present a novel approach to control design for nonlinear systems which
leverages model-free policy optimization techniques to learn a linearizing
controller for a physical plant with unknown dynamics. Feedback linearization
is a technique from nonlinear control which renders the input-output dynamics
of a nonlinear plant _linear_ under application of an appropriate feedback
controller. Once a linearizing controller has been constructed, desired output
trajectories for the nonlinear plant can be tracked using a variety of linear
control techniques. However, the calculation of a linearizing controller
requires a precise dynamics model for the system. As a result, model-based
approaches for learning exact linearizing controllers generally require a
simple, highly structured model of the system with easily identifiable
parameters. In contrast, the model-free approach presented in this paper is
able to approximate the linearizing controller for the plant using general
function approximation architectures. Specifically, we formulate a continuous-
time optimization problem over the parameters of a learned linearizing
controller whose optima are the set of parameters which best linearize the
plant. We derive conditions under which the learning problem is (strongly)
convex and provide guarantees which ensure the true linearizing controller for
the plant is recovered. We then discuss how model-free policy optimization
algorithms can be used to solve a discrete-time approximation to the problem
using data collected from the real-world plant. The utility of the framework
is demonstrated in simulation and on a real-world robotic platform.

##  I Introduction

Geometric nonlinear control theory has developed a powerful set of feedback
architectures which exploit the underlying structure of a control system to
simplify downstream tasks such as trajectory generation and tracking [1, 2].
However, geometric controllers often require an accurate model for the system
or a simple parameterization of the dynamics which can be readily identified.
Meanwhile, the model-free reinforcement learning literature [3, 4, 5] has
sought to automatically compute optimal feedback controllers for unknown
systems without relying on structural assumptions about the dynamics. However,
despite a recent resurgence of research into these methods [6, 7, 8], their
poor sample complexity has thus-far limited their applicability to many
robotics applications. In this paper we unify these disparate approaches by
using model-free reinforcement learning algorithms to optimize the performance
of a class of geometric controllers for a plant with unknown dynamics,
observing improvement in the performance of the trained controllers with
practical amounts of data collects from the plant.

Specifically, this paper focuses on a geometric technique called feedback
linearization, which renders the input-output dynamics of a nonlinear system
_linear_ under the application of an appropriately chosen feedback controller.
Once a linearizing controller has been constructed, it can be used in
conjunction with efficient tools from linear systems theory to generate and
track desired output trajectories for the full nonlinear plant [9, 10, 11].
Despite its widespread use throughout robotics [12, 13, 14, 1, 2] the primary
drawback of feedback linearization is that exact cancellation of the system’s
nonlinearities requires a precise model of the dynamics. Complex phenomena
such as friction, aerodynamic drag, and internal actuator dynamics yield
nonlinearities which may be challenging to model or identify. While there have
been extensive efforts to construct exact linearizing controllers for unknown
plants using extensions of adaptive linear control theory [15, 16, 17, 18, 19,
20, 21, 22], these methods require a highly structured representation of the
system’s nonlinearites and employ complicated parameter update schemes to
avoid singularities in the learned controller. A more detailed discussion of
these methods is provided in Section II-C.

Figure 1: Schematic diagram of our framework. By learning an appropriate
feedback linearizing controller, we render an initially unknown nonlinear
system (with state $x$, output $y$, and input $u$) _linear_ in an auxiliary
input $v$. The framework can make use of a nominal dynamics model (if
available).

In sharp contrast to the standard methods discussed above, we propose a
general framework for learning a linearizing controller for a plant with
unknown dynamics using model-free policy optimization techniques. Our approach
requires only that the order of the relationship between the inputs and
outputs is known, can utilize arbitrary function approximation schemes, and
remains singularity free during training. Concretely, our approach begins by
constructing a linearizing controller for an approximate dynamics model of the
plant (if available), and then augments this fixed model-based controller with
a learned affine correction term. We then define a continuous-time
optimization problem over the learned parameters which selects parameters
which better linearize the unknown plant. We provide conditions on the
structure of the learned component which ensure that the learning problem is
(strongly) convex. Our analysis draws connections between the familiar
persistency-of-excitation conditions from the adaptive control literature and
our formulation. Finally, we discuss how off-the-shelf reinforcement learning
algorithms can be used to solve discrete-time approximations to the
optimization problem, and demonstrate the utility of our framework by learning
linearizing controllers for highly uncertain robotic systems in both
simulation and on real-world hardware.

The rest of the paper is structured as follows. Section II introduces feedback
linearization and discusses prior approaches for learning linearizing
controllers. Section III details our approach, provides our theoretical
results, and discusses practical algorithms for solving the learning problem.
Our simulated and real-world robotic examples are presented in Section IV, and
Section V provides closing remarks and discusses future work.

##  II Feedback Linearization

This section outlines how to compute an input-output linearizing controller
for a known dynamics model and discusses prior data-driven methods for
learning a linearizing controller for an unknown plant. Due to space
constraints, we refer the interested reader to [1, Chapter 9] for a more
complete introduction to feedback linearization. In this paper, we consider
square control-affine systems of the form

$\displaystyle\dot{x}$ $\displaystyle=f(x)+g(x)u$ (1) $\displaystyle y$
$\displaystyle=h(x),$

with sate $x\in\mathbb{R}^{n}$, input $u\in\mathbb{R}^{q}$ and output
$y\in\mathbb{R}^{q}$. The mappings $f\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$,
$g\colon\mathbb{R}^{n}\to\mathbb{R}^{n\times q}$ and
$h\colon\mathbb{R}^{n}\to\mathbb{R}^{q}$ are each assumed to be smooth. We
restrict our attention to a compact subset $D\subset\mathbb{R}^{n}$ of the
state-space.

###  II-A Single-input single-output systems

We begin by introducing feedback linearization for single-input, single-output
(SISO) systems (i.e., $q=1$). We begin by examining the first time derivative
of the output:

$\displaystyle\dot{y}$
$\displaystyle=\frac{dh}{dx}(x)\cdot\Big{(}f(x)+g(x)u\Big{)}$
$\displaystyle=\underbrace{\frac{dh}{dx}(x)\cdot
f(x)}_{L_{f}h(x)}+\underbrace{\frac{dh}{dx}(x)\cdot g(x)}_{L_{g}h(x)}u$

Here the terms $L_{f}h(x)$ and $L_{g}h(x)$ are known as _Lie derivatives_ [1],
and capture the rate of change of $y=h(x)$ along the vector fields $f$ and
$g$, respectively. In the case that $L_{g}h(x)\neq 0$ for each $x\in D$, we
can apply the control law

$u(x,v)=\frac{1}{L_{g}h(x)}(-L_{f}h(x)+v)~{},$ (2)

which exactly ‘cancels out’ the nonlinearities of the system and enforces the
linear relationship $\dot{y}=v$. After application of this linearizing control
law, we can now control the output trajectory $y(t)$ through its first
derivative. However if the input does not affect the first time derivative of
the output (that is, if $L_{g}h\equiv 0$) then the control law (2) will be
undefined. In general, we can differentiate $y$ multiple times until the input
appears. Assuming that the input does not appear the first $\gamma-1$ times we
differentiate the output, the $\gamma$-th time derivative of $y$ will be of
the form

$y^{(\gamma)}=L_{f}^{\gamma}h(x)+L_{g}L_{f}^{\gamma-1}h(x)u.$

Here, $L_{f}^{\gamma}h(x)$ and $L_{g}L_{f}^{\gamma-1}h(x)$ are higher order
Lie derivatives (see [1, Chapter 9] for more details). If
$L_{g}L_{f}^{\gamma-1}h(x)\neq 0$ for each $x\in D$ then the control law

$u(x,v)=\frac{1}{L_{g}L_{f}^{\gamma-1}h(x)}\big{(}-L_{f}^{\gamma}h(x)+v\big{)}$

enforces the trivial linear relationship $y^{(\gamma)}=v$. We refer to
$\gamma$ as the _relative degree_ of the nonlinear system, which is simply the
order of its input-output relationship.

###  II-B Multiple-input multiple-output systems

Next, we consider (square) multiple-input, multiple-output (MIMO) systems
where $q>1$. As in the SISO case, we differentiate each of the output channels
until at least one input appears. Let $\gamma_{j}$ be the number of times we
need to differentiate $y_{j}$ (the $j$-th entry of $y$) for at least one input
to appear. Combining the resulting expressions for each of the outputs yields
an input-output relationship of the form

$[y_{1}^{(\gamma_{1})},\dots,y_{q}^{(\gamma_{q})}]^{T}=b(x)+A(x)u.$ (3)

Here, the matrix $A(x)\in\mathbb{R}^{q\times q}$ is known as the _decoupling
matrix_ and the vector $b(x)\in\mathbb{R}^{q}$ is known as the _drift term_.
If $A(x)$ is bounded away from singularity for each $x\in D$ then we observe
that the control law

$u(x,v)=A^{-1}(x)(-b(x)+v)$ (4)

where $v\in\mathbb{R}^{q}$ yields the decoupled linear system

$[y_{1}^{(\gamma_{1})},y_{2}^{(\gamma_{2})},\dots,y_{q}^{(\gamma_{q})}]^{T}=[v_{1},v_{2},\dots,v_{q}]^{T},$
(5)

where $v_{j}$ is the $j$-th entry of $v$ and $y_{j}^{(\gamma_{j})}$ is the
$\gamma_{j}$-th derivative of the $j$-th output. We refer to
$\gamma=(\gamma_{1},\gamma_{2},\dots,\gamma_{q})$ as the _vector relative
degree_ of the system. The decoupled dynamics (5) are LTI and can be compactly
represented with the _reference model_

$\dot{\xi}_{r}=A\xi_{r}+Bv_{r},$ (6)

where
$\xi_{r}=(y_{1},\dot{y}_{1},\dots,\dots,y_{1}^{\gamma_{1}-1},\dots,y_{q},\dots,y_{q}^{\gamma_{q}-1})$
and $A\in\mathbb{R}^{\|\gamma\|\times\ |\gamma\|}$ and
$B\in\mathbb{R}^{|\gamma|\times q}$ are dynamics matrices with appropriate
entries. The reference model can then be used to construct a desired
trajectory $\xi_{r}(\cdot)$ for the output of the nonlinear system (and its
derivatives), and to construct linear feedback controllers which can be used
in conjunction with (4) to track the reference (see [1, Theorem 9.14] for
details).

###  II-C Constructing linearizing controllers with data

The fundamental challenge model-based methods face when constructing a
linearizing controller for an unknown plant is that it is very difficult to
identify an estimate for the system’s decoupling matrix (or possibly its
inverse) which is guaranteed to be invertible. The predominant approaches for
learning linearizing controllers are founded on the linear model reference
adaptive control (MRAC) literature [15], and seek to recursively improve an
estimate for the true linearizing controller using data collected from the
plant [15, 16, 17, 18, 19, 23, 24, 25]. These methods update the parameters
for estimates of the decoupling matrix (or its inverse) online, but require
that the estimated matrix remain invertible at all times. Typically,
projection-based update rules are used to keep the parameters of the estimated
matrix in a region which is singularity-free. However, the construction of
these bounds requires a highly accurate yet simple parameterization of the
system, and assumes that the true parameters of the system lie within some
nominal set.

One alternative approach [26, 27, 28] is to perturb the estimated linearizing
control law to avoid singularities. These methods enable the use of more
general function approximation schemes but sacrifice some tracking
performance. Recently, non-parametric function approximators have been been
used to learn a linearizing controller [20, 29], but these methods still
require structural assumptions to avoid singularities. In the following
section we use model-free policy optimization algorithms to update the
paramters of a learned linearizing controller while avoiding singularities.

##  III Directly Learning a Linearizing Controller

Our goal is to learn a linearizing controller for the plant

$\displaystyle\dot{x}_{p}$ $\displaystyle=f_{p}(x_{p})+g_{p}(x_{p})u_{p}$ (7)
$\displaystyle y_{p}$ $\displaystyle=h_{p}(x_{p})$

which is unknown. Our approach can incorporate nominal model

$\displaystyle\dot{x}_{m}$ $\displaystyle=f_{m}(x_{m})+g_{m}(x_{m})u_{m}$ (8)
$\displaystyle y_{m}$ $\displaystyle=h_{m}(x_{m})$

which represents our "best guess" for the true dynamics of the plant. We will
assume that the plant and model have well-defined relative degrees
$(\gamma_{1}^{p},\dots,\gamma_{q}^{p})$ and
$(\gamma_{1}^{m},\dots,\gamma_{q}^{m})$, respectively, on some chosen set
$D\subset\mathbb{R}^{n}$. We make the following assumption about our plant and
model:

###### Assumption 1

The model system (8) and plant (7) have the same relative degree on $D$ in the
sense that
$(\gamma_{1}^{p},\dots,\gamma_{q}^{p})=(\gamma_{1}^{m},\dots,\gamma_{q}^{m})$.

This is a rather mild assumption, as the order of the relationship between the
inputs and outputs of the plant can usually be inferred from first principles,
even without a perfect model for the dynamics of the system. For example, for
Lagrangian systems, such as the manipulator arms we consider in Section IV, we
know the torques applied by an actuator produce an acceleration in the joints,
even if we don’t know the exact relationship between the two quantities.
Moreover, with this assumption in place, we know there are linearizing
controllers of the form

$\displaystyle u_{m}(x,v)$ $\displaystyle=\beta_{m}(x)+\alpha_{m}(x)v$
$\displaystyle u_{p}(x,v)$ $\displaystyle=\beta_{p}(x)+\alpha_{p}(x)v$

for the model and plant, respectively, which match the input-output dynamics
of both systems to a common reference model (6). Here,
$\beta_{p}(x),\beta_{m}(x)\in\mathbb{R}^{q}$ and
$\alpha_{p}(x),\alpha_{m}(x)\in\mathbb{R}^{q\times q}$. While we do not know
$u_{p}$ _a priori_ , we do know that

$\displaystyle\beta_{p}(x)$ $\displaystyle=\beta_{m}(x)+\Delta\beta(x)$
$\displaystyle\alpha_{p}(x)$ $\displaystyle=\alpha_{m}(x)+\Delta\alpha(x)$

for some continuous functions $\Delta\beta$ and $\Delta\alpha$. We construct
the following parameterized estimates for these functions:

$\Delta\beta(x)\approx\beta_{\theta_{1}}(x)\ \ \ \ \ \
~{}\Delta\alpha(x)\approx\alpha_{\theta_{2}}(x)$

Here, $\theta_{1}\in\Theta_{1}\subset\mathbb{R}^{K_{1}}_{1}$ and
$\theta_{2}\in\Theta_{2}\subset\mathbb{R}^{K_{2}}$ are parameters to be
trained by running experiments on the plant. We will assume that $\Theta_{1}$
and $\Theta_{2}$ are convex sets, and we will frequently abbreviate
$\theta=(\theta_{1},\theta_{2})\in\Theta_{1}\times\Theta_{2}:=\Theta$. We
assume that $\beta_{\theta_{1}}$ and $\alpha_{\theta_{2}}$ are continuous in
$x$ and continuously differentiable in $\theta_{1}$ and $\theta_{2}$,
respectively.

Figure 2: (a) Depiction of the tracking task for the double pendulum. The blue
curves represent the desired joint angles over time. The orange and green
curves represent trajectories generated by the ideal feedback linearizing
controller (perfect model information) and our learned policy. Both
trajectories begin with an initial tracking error, but quickly converge to the
desired motion. (b) Total $\ell_{2}$ error on a set-point tracking task for
the Baxter robot after 104 minutes of training for the desired linear system
is (dotted, green), the nominal feedback linearizing controller (red), and our
learned controller (blue)

Altogether, for a given $\theta=(\theta_{1},\theta_{2})\in\Theta$ our estimate
for the controller which exactly linearizes the plant is given by

$\hat{u}_{\theta}(x,v)=[\beta_{m}(x)+\beta_{\theta_{1}}(x)]+[\alpha_{m}(x)+\alpha_{\theta_{2}}(x)]v$
(9)

When no prior information about the dynamics of the plant is available (other
than its vector relative degree), we simply set $\beta_{m}\equiv 0$ and
$\alpha_{m}\equiv 0$ in the above expression. Next we define an optimization
problem which selects the parameters for the learned controller which best
linearize the plant.

###  III-A Continuous-time optimization problem

From Section II we know that the input-output dynamics of the plant are of the
form

$y^{(\gamma)}=b_{p}(x)+A_{p}(x)u$ (10)

where the terms $b_{p}$ and $A_{p}$ are unknown to us, and we have written the
highest order derivatives of the outputs as
$y^{(\gamma)}=(y_{1}^{(\gamma_{1})}\dots,y_{q}^{(\gamma_{q})})^{T}$ to
simplify notation. Under application of $\hat{u}_{\theta}$ the dynamics are
given by:

$y^{(\gamma)}=\underbrace{b_{p}(x)+A_{p}(x)\hat{u}_{\theta}(x,v)}_{W_{\theta}(x,v)}$
(11)

Since our goal is to find parameters which linearize the plant, we want to
find $\theta^{*}\in\Theta$ such that $W_{\theta^{*}}(x,v)\approx v$ for each
$x\in D$ and $v\in\mathbb{R}^{q}$. Thus, we define the point-wise loss
$\ell\colon\mathbb{R}^{n}\times\mathbb{R}^{q}\times\mathbb{R}^{K_{1}+K_{2}}\to\mathbb{R}$
by

$\ell(x,v,\theta)=\mathopen{}\mathclose{{}\left\lVert
v-W_{\theta}(x,v)}\right\rVert_{2}^{2},$ (12)

which provides a measure of how well the learned controller $\hat{u}_{\theta}$
linearizes the plant at the state $x$ when the virtual input $v$ is applied to
the linear reference model. Next, we specify a probability distribution $X$
over $\mathbb{R}^{n}$ with support $D\subset\mathbb{R}^{n}$ and let $V$ be the
uniform distribution over the set
$\mathopen{}\mathclose{{}\left\\{v\in\mathbb{R}^{q}\colon\mathopen{}\mathclose{{}\left\lVert
v}\right\rVert\leq 1}\right\\}$. We then define the weighted loss

$L(\theta)=\operatorname{\mathbb{E}}_{x\sim X,v\sim V}\ell(x,v,\theta)$ (13)

and select our optimal choice for the parameters of the learned controller via
the following optimization:

$(\textbf{P})\colon\ \min_{\theta\in\Theta}L(\theta)$ (14)

Here, the distribution $X$ models our preference for having an accurate
linearizing controller at different points in the state-space, and the
uniformity of $V$ ensures that optimal solutions of P accurately linearize the
plant for all possible choices of the virtual input to the reference model.
The primary challenge in solving P is that we do not know the terms in
$W_{\theta}(x,v)$, since we do not know $A_{p}(x)$ and $b_{p}(x)$. However,
using the relationship (11), we can query $W_{\theta}(x,v)$ by measuring
$y^{(\gamma)}$ when different inputs are applied to the plant. Thus, P can be
solved by running experiments on the plant and using any stochastic
optimization method which only requires access to the point-wise loss (13).
While we focus on the use of policy-gradient reinforcement learning algorithms
below, there are many possible model-free approaches for solving P which merit
further investigation. Importantly, since we have abstracted away the need for
a model when solving P, we can iteratively improve the performance of our
learned controller without requiring that it remain invertible at each stage
of the optimization. Next, we discuss when we can recover the true linearizing
controller for the plant by solving P:

###### Lemma 1

Suppose that there exists $\theta^{*}\in\Theta$ such that
$\hat{u}_{\theta^{*}}(x,v)=u_{p}(x,v)$ for each $x\in D$ and $v\in V$. Then
$\theta^{*}$ is a globally optimal solution of P.

###### Proof:

Note that if $u_{\theta^{*}}(x,v)=u_{p}(x,v)$ for each $x\in X$ and $v\in V$
then $L(\theta)=0$. Moreover, we clearly have $L(\theta)\geq 0$ for each
$\theta\in\Theta$. Thus, $\theta^{*}$ must be a global minimizer of the
optimization (14). ∎

However, P is generally non-convex meaning that in practice we can only hope
to find locally optimal solutions to the problem. Thus, we seek conditions
which simplify the structure of the optimization. The standard convergence
proofs in the adaptive control literature assume the learned controller is of
the form

$\beta_{\theta_{1}}(x)=\sum_{k=1}^{K_{1}}\theta_{k}^{1}\beta_{k}(x)\ \ \ \ \ \
\alpha_{\theta_{2}}(x)=\sum_{k=1}^{K_{2}}\theta_{k}^{2}\alpha_{k}(x)$ (15)

where $\mathopen{}\mathclose{{}\left\\{\beta_{k}}\right\\}_{k=1}^{K_{1}}$ and
$\mathopen{}\mathclose{{}\left\\{\alpha_{k}}\right\\}_{k=1}^{K_{2}}$ are
nonlinear continuous functions. When we adopt this structure our optimization
becomes convex, whose proof can be found in the Appendix:

###### Lemma 2

Assume that $\beta_{\theta_{1}}$ and $\alpha_{\theta_{2}}$ are of the form
(15). Then P is convex. Moreover, if
$\mathopen{}\mathclose{{}\left\\{\beta_{k}}\right\\}_{k=1}^{K_{1}}$ and
$\mathopen{}\mathclose{{}\left\\{\alpha_{k}}\right\\}_{k=1}^{K_{2}}$ are each
linearly independent then P is strongly convex.

Thus, when our learned controller takes the form (15) we can reliably use
iterative techniques to find its globally optimal solution. The proof of Lemma
2, which is given in the accompanying technical report, shows that $L$ is
actually quadratic in the parameters when the learned controller is of the
form (15). The key property being exploited here is that $W_{\theta}(x,v)$
becomes affines in $\theta$ when the controller is linear in the parameters.
Our conditions on $X$ and $V$ are analogous to the persistency of excitation
results from the adaptive control literature [17, Chapter 2], and ensure that
each component of the learned controller is excited while solving P. This is
the underlying reason why the optimization becomes strongly convex when the
components of our controller are linearly independent. Taken together, the
preceding lemmas imply our main theoretical result:

###### Theorem 1

Suppose that for some $\theta^{*}\in\Theta$ we have
$\hat{u}_{\theta^{*}}(x,v)=u_{p}(x,v)$ for each $x\in D$ and
$v\in\mathbb{R}^{q}$. Further assume that the learned controller is of the
form (15) and that the sets
$\mathopen{}\mathclose{{}\left\\{\beta_{k}}\right\\}_{k=1}^{K_{1}}$ and
$\mathopen{}\mathclose{{}\left\\{\alpha_{k}}\right\\}_{k=1}^{K_{2}}$ are
linearly independent. Then $\theta^{*}$ is the unique global (and local)
minimizer of P.

There are a number of well-studied bases, such as polynomials or radial basis
functions, which can be used to approximate any continuous function to a
desired degree of accuracy. Thus, by including enough terms in
$\mathopen{}\mathclose{{}\left\\{\beta_{k}}\right\\}_{k=1}^{K_{1}}$ and
$\mathopen{}\mathclose{{}\left\\{\alpha_{k}}\right\\}_{k=1}^{K_{2}}$ we can
theoretically recover $u_{p}$ to a pre-specified precision by solving P.
However, for high dimensional systems, the number of terms required in such
expansions can become prohibitively large. Other architectures, such as multi-
layer feed-forward neural networks, yield more compactly represented function
approximation schemes but complicate the analysis of P. However, since our
approach does not suffer from singularities during the optimization process,
in practice we can still use these powerful function approximators to find an
improved linearizing controller for high-dimensional systems with unknown
dynamics. We next demonstrate this point by discussing how policy optimization
algorithms can be used to solve a discrete-time approximation to P.

Figure 3: (a) The performance of our learned linearizing controller on two
high-performance quadrotor tracking tasks when the parameters of the inital
dynamics model are scaled by a factor of $0.6$. The first task is a figure-
eight, and the second is a corkscrew maneuver. In both maneuvers the quadrotor
also tracks an oscillating reference in the yaw angle. In both figures, the
green trajectory is the one taken by the learned controller and the blue
trajectory is that of the true linearizing controller for the system. In both
cases, the learned controller closely matches the desired behavior. (b)
Learning curves for learned linearizing controllers which are instantiated
with nominal dynamics models with three levels of scaling of the true dynamics
parameters.

###  III-B Discrete-time approximations and reinforcment learning

While the theory and optimization problem we developed in the previous section
are in continuous time, many real world plants have actuators which can only
be updated at a fixed sampling frequency. Thus in this section we formulate a
discrete time approximation to P which we cast as a canonical reinforcement
learning problem. Unfortunately, the analysis of linearizable systems in the
sampled-data setting becomes significantly more complex [30], which is why the
majority of the literature on adaptive control for linearizable systems is
formulated in continuous time.

Letting $\Delta t>0$ denote the sampling rate, we set $t_{k}=k\Delta t$ for
each $k\in\mathbb{N}$. Next, letting $x(\cdot)$ denote the state trajectory of
the plant, we then denote $x_{k}=x(t_{k})$ and let $u_{k}$ denote the control
applied on the interval $[t_{k},t_{k+1})$. This yields the following
difference equation for the dynamics of the plant:

$x_{k+1}=x_{k}+\underbrace{\int_{t_{k}}^{t_{k+1}}f_{p}(x(t))+g_{p}(x(t))u_{k}d\tau}_{F(x_{k},u_{k})}$

The map $F_{p}\colon\mathbb{R}^{n}\times\mathbb{R}^{q}\to\mathbb{R}^{n}$ will
generally no longer be affine in the control. Similarly, we let
$\xi(\cdot)=(y_{1}(\cdot),\dots,y_{1}^{(\gamma_{1}-1)}(\cdot),\dots,y_{q}(\cdot),\dots
y_{q}^{(\gamma_{q}-1)}(\cdot))$ denote the solution to the linearized portion
of the state, and set $\xi_{k}=\xi(t_{k})$ for each iterate. Now, if we
integrate (10) over the interval $[t_{k},t_{k+1})$ we obtain a difference
equation for the outputs of the form

$\xi_{k+1}=e^{A\Delta t}\xi_{k}+H(x_{k},u_{k})$ (16)

where again the mapping
$H\colon\mathbb{R}^{n}\times\mathbb{R}^{q}\to\mathbb{R}^{\gamma}$ will
generally no longer be affine in $u_{k}$. This is the primary hurdle to
applying the theory developed in the previous section to the current setting.
Nevertheless, we can still attempt to enforce an approximate linear
relationship between a virtual input $v_{k}$ and successive iterates of the
outputs. Letting $\bar{A}=e^{A\Delta t}$ and $\bar{B}=\int_{0}^{\Delta
t}e^{At}Bdt$, the sampled-data version of the reference model (6) becomes

$\xi_{k+1}=\bar{A}\xi_{k}+\bar{B}v_{k},$

which is the ideal linear behavior we would like to enforce by applying the
control $u_{k}=\hat{u}_{\theta}(x_{k},v_{k})$. To encourage actions which
better track the discrete-time reference model, we apply the point-wise loss
$\bar{\ell}\colon\mathbb{R}^{n}\times\mathbb{R}^{q}\times\mathbb{R}^{q}\to\mathbb{R}$
where

$\bar{\ell}(x_{k},u_{k},v_{k})=\|\bar{B}v_{k}-H(x_{k},u_{k})\|_{2}^{2}.$

We can calculate this quantity using (16) and by observing $\xi_{k}$ and
$\xi_{k+1}$, which can each be calculated by numerically differentiating the
outputs from the plant. This loss provides a measure of how well the control
$u_{k}$ enforces the desired change in the state of the linear system (as
specified by $v_{k}$) at the state $x_{k}$. We then define the following
reinforcement learning problem over the parameters of $\hat{u}_{\theta}$:

$\displaystyle\min_{\theta\in\Theta}\mathbb{E}_{x_{0}\sim X,v_{k}\sim
V,w_{k}\sim\mathcal{N}(0,\sigma_{w}^{2})}$
$\displaystyle\mathopen{}\mathclose{{}\left[\sum_{k=1}^{N}\bar{\ell}(x_{k},u_{k},v_{k})}\right]$
$\displaystyle\text{subject to: }\ \ \ x_{k+1}$
$\displaystyle=x_{k}+F(x_{k},u_{k}),\ \ \ x_{0}=x_{0}$ $\displaystyle u_{k}$
$\displaystyle=\hat{u}_{\theta}(x_{k},v_{k})+w_{k}$

Here, we sample initial conditions for the problem using our desired state
distribution $X$ and sample inputs to the linear system according to $V$ at
each time step. The zero mean added noise $w_{k}$ is used to encourage
exploration, and to make the effects of the policy random. Finally, $N$ is the
length of the training episodes. This problem can be solved to local
optimality using standard reinforcement learning algorithms [31, 32, 33, 34]
and by running experiments on the real-world hardware to evaluate
$\bar{l}(x_{k},u_{k},v_{k})$. While it is not immediately clear how to extend
the theoretical results from Section III-A to the present case, we intend to
address this issue in a forthcoming article.

##  IV Examples

We now use our approach to learn feedback linearizing policies for three
different systems and use the learned policies to construct tracking
controllers as in [1, Theorem 9.14]. In all cases, the input to the
parameterized policy replaces all angles with their sine and cosine.

###  IV-A Simulations

####  IV-A1 Double pendulum with polynomial policies

We first test our approach on a fully actuated double pendulum with state
$x=[\theta_{1},\theta_{2},\omega_{1},\omega_{2}]^{T}$, output
$y=[\theta_{1},\theta_{2}]^{T}$, where $\theta_{1}$ and $\theta_{2}$ represent
the angles of the two joints, with angular rates $\omega_{1}$ and $\omega_{2}$
respectively. The system has two inputs $u_{1}$ and $u_{2}$ that control the
torque at both joints. Although the system is relatively low dimensional and
fully actuated, it is highly nonlinear [35].

For the learning problem, we give the nominal model inaccurate estimates for
both the length and mass of each of the arms. Specifically, we scale each of
the parameters to $1/2$ their true values. The learned controller is comprised
of 150 radial basis functions, which are centered randomly throughout the
state-space. The learned controller is linear in its parameters, so as to
agree with Theorem 1. We use the REINFORCE algorithm [31], and baseline state-
value estimates with the average reward over all states. At each iteration (or
epoch) we collect $50$ rollouts of $0.25$ seconds each, and we train for 500
epochs. Figure 2 (a) presents the result of using our learned controller to
track a desired sinusoidal trajectory, and compares it to the performance of
the exact linearizing controller for the system. The difference between the
performance our learned controller and the theoretical ideal observed to be
very small. We do not plot the trajectory generated by the nominal model-based
controller, since it immediately diverges from the desired behavior. Similar
results were observed using the same learned controller to track other desired
trajectories for the system. The linear feedback matrix used in the tracking
controller was obtained by solving an infinite horizon LQR problem where the
state deviation was penalized 10 times more than the magnitude of the input.

####  IV-A2 14D quadrotor with neural network policies

Our second simulation environment uses the quadrotor model and feedback
linearization controller proposed in [36], which makes use of dynamic
extension [1]. In particular, the states for the model are
$(x,y,z,\psi,\theta,\varphi,\dot{x},\dot{y},\dot{z},p,q,r,\xi,\zeta)$ where
$x$, $y$ and $z$ are the Cartesian coordinates of the quadrotor, and $\psi$,
$\theta$ and $\varphi$ represent the roll, pitch and yaw of the quadrotor,
respectively. The next six states represent the time derivatives of these
state:
$\frac{d}{dt}(x,y,z,\psi,\theta,\varphi)=(\dot{x},\dot{y},\dot{z},p,q,r)$.
Finally, $\xi$ and $\zeta$ are the extra states obtained from the dynamic
extension procedure. The outputs for the model are the $x$, $y$, $z$ and
$\psi$ coordinates.

The inaccurate nominal dynamics model was constructed by multiplying the mass
and inertial constants of the true system by factors of $0.33$, $0.6$ and
$0.8$ for different experiments. The trained policies were feed-forward neural
networks with $\tanh$ activations with $2$ hidden layers of width $64$. For
each training epoch, 50 rollouts of length 25 were collected and the
parameters were updated using PPO. We trained both policies for 2500 epochs.
Figure 3 (b) illustrates how a better initial dynamics model leads to faster
learning of an accurate linearizing controller. In particular, the reward-per-
epoch is plotted for policies trained with the three scaling factors indicated
above. We observe that worse initial models result in worse policy
performance, given the same network architecture and training time. Figure 3
(a) demonstrates the ability of the learned controller to overcome significant
model mismatch (scaling factor of $0.6$) to match the desired linear tracking
behavior.

###  IV-B Robotic experiment: 7-DOF manipulator arm

We also evaluate our approach in hardware, on a 7-DOF Baxter robot arm. The
dynamics of this 14-dimensional system are extremely coupled and nonlinear.
Taking the 7 joint angles as output $y$, however, the system is input-output
linearizable with relative degree two. We use the system measurements (i.e.,
masses, link lengths, etc.) provided with Baxter’s pre-calibrated URDF [37]
and the OROCOS Kinematics and Dynamics Library (KDL) [38] to compute a nominal
feedback linearizing control law.

This nominal controller suffers from several inaccuracies. First, Baxter’s
actuators are series-elastic, meaning that each joint contains a torsion
spring [39] which is unmodeled, and the URDF itself may not be perfectly
accurate. Second, the OROCOS solver is numerical, which can lead to errors in
computing the decoupling matrix and drift term. Finally, our control
architecture is implemented in the Robot Operating System [40], which can lead
to minor timing inconsistency.

We use the PPO algorithm to tune the parameters of a $128\times 2$ neural
network with $\tanh$ activations for the learned component of the controller.
For each training epoch, 1250 rollouts of one timestep (0.05 s) each were
collected. We trained for 100 epochs, which took 104 minutes. Figure 2 (b)
summarizes typical results on tracking a square wave reference trajectory for
each joint angle with period 5s. The nominal feedback linearized model from
OROCOS has significant steady-state error. Our learned approach significantly
reduces, but does not eliminate, this error. We conjecture the remaining error
is a sign that the (relatively small) neural network may not be sufficiently
expressive.

##  V Discussion

While the methods we have presented avoid some of the challenges model-based
methods face when trying to construct a linearizing controller for an unknown
plant, we feel the techniques presented here should be used to complement
model-based design approaches rather than replace them. As we observed
empirically, a reasonable (though ultimately inaccurate) model can be used to
provide a better starting point for the learning process. Thus, we feel that
in practice our approach should be used primarily to overcome difficult to
model non-linearities, and should be used in combination with a simple nominal
dynamics models with parameters which are easily identified. Investigating the
performance of different learning algorithms in real-world settings merits
significant future research effort. On the theoretical side, we feel the
general way in which we constructed the optimization problem in III-A provides
a foundation for combining model-free policy optimization techniques with
geometric control architectures. Future work will investigate how to optimize
over different classes of controllers by extending the approach presented
here.

###  -A Proof of Lemma 2

First, we rearrange (10) into the form

$y_{p}^{\gamma}=b_{p}(x)+A_{p}(x)\beta_{m}(x)+A_{p}(x)\alpha_{m}(x)v\\\
+\sum_{i=1}^{K_{1}}\theta_{k}^{1}A_{p}(x)\beta_{k}(x)+\sum_{k=1}^{K_{2}}\theta_{k}^{2}A_{p}(x)\alpha_{k}(x)v$
(17)

to separate out the portions that depend on $\theta$. This can be further
condensed by putting $y^{\gamma}=\bar{W}(x,v)+\hat{W}(x,v)\theta$ where we set

$\bar{W}(x,v)=b_{p}(x)+A_{p}(x)\beta_{m}(x)+A_{p}(x)\alpha_{m}(x)v$
$\hat{W}(x,v)=A_{p}(x)[\beta_{1}(x),\dots,\beta_{K_{1}}(x),A_{1}(x)v,\\\
\dots,A_{K_{2}}(x)v]$

Letting $c(x,v)=(v-\bar{W}(x,v))$, we can rewrite

$\displaystyle\ell(x,v,\theta)$
$\displaystyle=\big{(}c(x,v)-\hat{W}(x,v)\theta\big{)}^{T}\big{(}c(x,v)-\hat{W}(x,v)\theta\big{)}$
$\displaystyle=\theta^{T}\hat{W}^{T}(x,v)\hat{W}(x,v)\theta-2\theta\hat{W}(x,v)c(x,v)$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +c(x,v)^{T}c(x,v)$

From here we observe that $L(\theta)=\theta^{T}W\theta+\theta^{T}F+d$ where
$W=\operatorname{\mathbb{E}}_{x\sim X,v\sim V}\hat{W}(x,v)^{T}\hat{W}(x,v)$ is
a positive semi-definite matrix, $F=\operatorname{\mathbb{E}}_{x\sim X,v\sim
V}\hat{W}(x,v)c(x,v)$ and $d=\operatorname{\mathbb{E}}_{x\sim X,v\sim
V}c(x,v)^{T}c(x,v)$. Thus, recalling that $\Theta$ is assumed to be a convex
set, we see that (14) is a convex optimization problem which will be strictly
convex if, and only if, $W$ is positive definite. Letting
$w^{1}_{k}\colon(x,v)\to A_{p}(x)\beta_{k}(x)$ and $w^{2}_{k}\colon(x,v)\to
A_{p}(x)\alpha_{k}(x)v$, we see that $W$ is nothing but the Grammian of the
set
$\Omega=\mathopen{}\mathclose{{}\left\\{w_{1}^{1},\dots,w_{K_{1}}^{1},\omega^{2}_{1},\dots,w_{K_{2}}^{2}}\right\\}$
on $C(D\times B_{v},\mathbb{R}^{q})$ with respect to an inner product which is
weighted by the distributions $X$ and $V$. Thus, $W$ will be positive definite
if and only if $\Omega$ is linearly independent on $C(D\times
B_{v},\mathbb{R}^{q})$. For the sake of contradiction assume that $\Omega$ is
not linearly independent. Then there exists scalars
$\mathopen{}\mathclose{{}\left\\{c_{k}^{1}}\right\\}_{k=1}^{2}$ and
$\mathopen{}\mathclose{{}\left\\{c_{k}^{2}}\right\\}_{k=1}^{K_{2}}$ such that
for each $x\in D$ and $v\in V$

$\sum_{k=1}^{K_{1}}c_{k}^{1}\omega_{k}^{1}(x,v)+\sum_{k=1}^{K_{2}}c_{k}^{2}\omega_{k}^{2}(x,v)=0$
(18)

Since we know that $A_{p}(x)$ is invertible for each $x\in D$, this statement
is equivalent to

$\sum_{k=1}^{K_{1}}c_{k}^{1}\beta_{k}(x)+\sum_{k=1}^{K_{2}}c_{k}^{2}\alpha_{k}(x)v=0.$
(19)

holding for each $x\in D$ and $v\in V$. However, it is not difficult to see
that this condition is ruled out in the case that
$\mathopen{}\mathclose{{}\left\\{\beta_{k}}\right\\}_{k=1}^{K_{1}}$ and
$\mathopen{}\mathclose{{}\left\\{\alpha_{k}}\right\\}_{k=1}^{K_{2}}$ are
linearly independent sets.

## References

  * [1] S. Sastry, _Nonlinear systems: analysis, stability, and control_. Springer Science & Business Media, 1999, vol. 10.
  * [2] A. Isidori, _Nonlinear control systems_. Springer Science & Business Media, 2013.
  * [3] D. P. Bertsekas and J. N. Tsitsiklis, _Neuro-dynamic programming_. Athena Scientific, 1996.
  * [4] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_ , 2018\.
  * [5] R. S. Sutton, D. A. McAllester, S. P. Singh, and Y. Mansour, “Policy gradient methods for reinforcement learning with function approximation,” in _Advances in neural information processing systems_ , 2000, pp. 1057–1063.
  * [6] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _International conference on machine learning_ , 2015, pp. 1889–1897.
  * [7] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [8] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [9] P. Martin, R. M. Murray, and P. Rouchon, “Flat systems, equivalence and trajectory generation,” 2003.
  * [10] R. E. Kalman _et al._ , “Contributions to the theory of optimal control,” _Bol. soc. mat. mexicana_ , vol. 5, no. 2, pp. 102–119, 1960.
  * [11] F. Borrelli, A. Bemporad, and M. Morari, _Predictive control for linear and hybrid systems_. Cambridge University Press, 2017.
  * [12] J. W. Grizzle, G. Abba, and F. Plestan, “Asymptotically stable walking for biped robots: Analysis via systems with impulse effects,” _IEEE Transactions on automatic control_ , vol. 46, no. 1, pp. 51–64, 2001.
  * [13] A. D. Ames, K. Galloway, K. Sreenath, and J. W. Grizzle, “Rapidly exponentially stabilizing control lyapunov functions and hybrid zero dynamics,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 4, pp. 876–891, 2014.
  * [14] D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in _2011 IEEE International Conference on Robotics and Automation_. IEEE, 2011, pp. 2520–2525.
  * [15] S. Sastry and M. Bodson, _Adaptive control: stability, convergence and robustness_. Courier Corporation, 1989\.
  * [16] J. J. Craig, P. Hsu, and S. S. Sastry, “Adaptive control of mechanical manipulators,” _The International Journal of Robotics Research_ , vol. 6, no. 2, pp. 16–28, 1987.
  * [17] S. S. Sastry and A. Isidori, “Adaptive control of linearizable systems,” _IEEE Transactions on Automatic Control_ , vol. 34, no. 11, pp. 1123–1131, 1989.
  * [18] K. Nam and A. Araposthathis, “A model reference adaptive control scheme for pure-feedback nonlinear systems,” _IEEE Transactions on Automatic Control_ , vol. 33, no. 9, pp. 803–811, 1988.
  * [19] I. Kanellakopoulos, P. V. Kokotovic, and A. S. Morse, “Systematic design of adaptive controllers for feedback linearizable systems,” in _1991 American Control Conference_. IEEE, 1991, pp. 649–654.
  * [20] J. Umlauft, T. Beckers, M. Kimmel, and S. Hirche, “Feedback linearization using gaussian processes,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 5249–5255.
  * [21] G. Chowdhary, H. A. Kingravi, J. P. How, and P. A. Vela, “Bayesian nonparametric adaptive control using gaussian processes,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 26, no. 3, pp. 537–550, 2014.
  * [22] ——, “Bayesian nonparametric adaptive control of time-varying systems using gaussian processes,” in _2013 American Control Conference_. IEEE, 2013, pp. 2655–2661.
  * [23] J. T. Spooner and K. M. Passino, “Stable adaptive control using fuzzy systems and neural networks,” _IEEE Transactions on Fuzzy Systems_ , vol. 4, no. 3, pp. 339–359, 1996.
  * [24] F.-C. Chen and H. K. Khalil, “Adaptive control of a class of nonlinear discrete-time systems using neural networks,” _IEEE Transactions on Automatic Control_ , vol. 40, no. 5, pp. 791–801, 1995.
  * [25] A. Yesildirek and F. L. Lewis, “Feedback linearization using neural networks,” in _Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94)_ , vol. 4. IEEE, 1994, pp. 2539–2544.
  * [26] E. B. Kosmatopoulos and P. A. Ioannou, “A switching adaptive controller for feedback linearizable systems,” _IEEE Transactions on automatic control_ , vol. 44, no. 4, pp. 742–750, 1999.
  * [27] ——, “Robust switching adaptive control of multi-input nonlinear systems,” _IEEE transactions on automatic control_ , vol. 47, no. 4, pp. 610–624, 2002\.
  * [28] C. P. Bechlioulis and G. A. Rovithakis, “Robust adaptive control of feedback linearizable mimo nonlinear systems with prescribed performance,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 9, pp. 2090–2099, 2008.
  * [29] J. Umlauft and S. Hirche, “Feedback linearization based on gaussian processes with event-triggered online learning,” _IEEE Transactions on Automatic Control_ , 2019.
  * [30] J. Grizzle and P. Kokotovic, “Feedback linearization of sampled-data systems,” _IEEE Transactions on Automatic Control_ , vol. 33, no. 9, pp. 857–859, 1988.
  * [31] R. S. Sutton and A. G. Barto, _Introduction to Reinforcement Learning_ , 1st ed. Cambridge, MA, USA: MIT Press, 1998\.
  * [32] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _International Conference on Machine Learning_ , 2015, pp. 1889–1897.
  * [33] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _CoRR_.
  * [34] D. Silver, G. Lever, N. Heess, T. Degris, D. Wierstra, and M. Riedmiller, “Deterministic policy gradient algorithms,” in _Proceedings of the 31st International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, 2014, pp. 387–395.
  * [35] T. Shinbrot, C. Grebogi, J. Wisdom, and J. A. Yorke, “Chaos in a double pendulum,” _American Journal of Physics_ , vol. 60, no. 6, pp. 491–499, 1992\.
  * [36] S. A. Al-Hiddabi, “Quadrotor control using feedback linearization with dynamic extension,” in _2009 6th International Symposium on Mechatronics and its Applications_. IEEE, 2009, pp. 1–3.
  * [37] R. Robotics, “Baxter,” _Retrieved Jan_ , vol. 10, p. 2014, 2013.
  * [38] H. Bruyninckx, “Open robot control software: the orocos project,” in _Proceedings 2001 ICRA. IEEE international conference on robotics and automation (Cat. No. 01CH37164)_ , vol. 3. IEEE, 2001, pp. 2523–2528.
  * [39] R. L. Williams II, “Baxter humanoid robot kinematics© 2017 dr. bob productions robert l. williams ii, ph. d., williar4@ ohio. edu mechanical engineering, ohio university, april 2017,” 2017.
  * [40] M. Quigley, K. Conley, B. P. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, and A. Y. Ng, “Ros: an open-source robot operating system,” in _ICRA Workshop on Open Source Software_ , 2009.