# Learning to Optimize in Model Predictive Control

Jacob Sacks∗ and Byron Boots∗ ∗University of Washington, Seattle WA 98105,
USA. {jsacks6<EMAIL_ADDRESS>
###### Abstract

Sampling-based Model Predictive Control (MPC) is a flexible control framework
that can reason about non-smooth dynamics and cost functions. Recently,
significant work has focused on the use of machine learning to improve the
performance of MPC, often through learning or fine-tuning the dynamics or cost
function. In contrast, we focus on learning to _optimize_ more effectively. In
other words, to improve the update rule within MPC. We show that this can be
particularly useful in sampling-based MPC, where we often wish to minimize the
number of samples for computational reasons. Unfortunately, the cost of
computational efficiency is a reduction in performance; fewer samples results
in noisier updates. We show that we can contend with this noise by learning
how to update the control distribution more effectively and make better use of
the few samples that we have. Our learned controllers are trained via
imitation learning to mimic an expert which has access to substantially more
samples. We test the efficacy of our approach on multiple simulated robotics
tasks in sample-constrained regimes and demonstrate that our approach can
outperform a MPC controller with the same number of samples.

##  I Introduction

Model Predictive Control (MPC) is a powerful, practical tool for solving
sequential decision problems on real-world systems. MPC has been successfully
used in a variety of tasks including autonomous helicopter aerobatics [1],
aggressive off-road driving [2, 3, 4], manipulation [5, 6], and humanoid robot
locomotion [7]. Recent work [4] has shown that many popular MPC algorithms can
be unified through the generic framework of dynamic mirror descent (DMD) [8],
a first-order online learning algorithm. This perspective provides an
opportunity to improve performance of existing algorithms by drawing on
powerful optimization techniques. Most modern approaches to optimization use
fixed update rules tailored to specific classes of problems. Recently,
research has explored _learning_ to optimize [9], where the update rule is
specified by a function approximator, such as a neural network, that can
improve optimization performance with experience.

In this work, we leverage the optimization perspective of sampling-based MPC
and adopt the learning-to-optimize framework in order to improve the update
rule. This is in contrast to most existing learning-based approaches to MPC,
which either focus on learning a good dynamics model [10, 11, 12, 3, 13, 14,
15], introducing a learned cost-shaping term into the objective [16], coupling
MPC with a learned value function [17, 18, 19, 20], or learning a good warm-
start for MPC to refine [21]. Other approaches have focused on performing MPC
with a learned latent space of high-dimensional observations [22, 23, 24, 25,
26, 27, 15], low-level skills [28], or controls such that sampling because
more efficient [29]. Another promising avenue has explored differentiating
through optimal controllers or planners to learn components of the
optimization pipeline. These methods propose to learn or fine-tune the
dynamics and cost functions of the controllers [16, 30, 31, 32, 33, 34, 29] or
parameters of differentiable planners [35, 36, 37, 38] end-to-end. Compared
with these approaches, we fix the dynamics and cost function and instead focus
on improving the optimization process.

For many practical sampling-based MPC algorithms, the primary challenge is
finding a good trade-off between speed and accuracy. Sampling-based MPC
algorithms work by using simple policies to sample control sequences, which
are used to roll out the dynamics function and compute a sample-based
approximation of the gradient of the objective function. This approximate
gradient is then used to update the sampling policy. Using complex dynamics
and cost functions can make each rollout computationally expensive. To contend
with this problem, one could use fewer samples to decrease computation, but
this can increase the noise in the sample-based gradient, leading to poor
performance.

In this paper, our objective is to _learn_ how to more effectively update the
control distribution with a small number of samples. To this end, we employ
imitation learning to train fast, low-sample controllers to imitate an expert
which makes use of additional samples. The learned optimizer is better able to
integrate information in the sample-constrained regime. Our key contributions
are that we:

  1. 1.

Leverage the gradient-based interpretation of many sampling-based MPC
algorithms and show how to improve performance by learning a better update
rule.

  2. 2.

Propose to use structured sampling techniques to provide more information to
the learned update than is contained in the noisy gradient, which enables us
to make better use of fewer samples.

  3. 3.

Empirically evaluate our proposed approach on multiple simulated robotics
tasks, in which the learned optimizer has restricted access to samples.

Our experiments show that the learned controller is indeed able to make better
use of fewer samples while remaining competitive or outperforming the expert
with the same number of samples. This illustrates the utility of the learning-
to-optimize framework in the domain of control and indicates the potential to
improve the viability of sampling-based MPC controllers on embedded platforms.

##  II Background

###  II-A Model Predictive Control

We consider the problem of controlling a discrete-time stochastic dynamical
system with states $x_{t}\in\mathbb{R}^{N}$ and controls
$u_{t}\in\mathbb{R}^{M}$. The system chooses controls using a policy
$\pi_{\theta_{t}}$ with parameters $\theta_{t}\in\Theta$, where $\Theta$ is
the set of feasible parameters. After applying the control, the system incurs
the instantaneous cost $c(x_{t},u_{t})$ and transitions to the next state
$x_{t+1}$ according to the dynamics

$x_{t+1}\sim f(x_{t},u_{t}),$ (1)

where $f:\mathbb{R}^{N}\times\mathbb{R}^{M}\rightarrow\mathbb{R}^{N}$ is a
stochastic transition map. Over a time horizon $H$, we sample a control
sequence $U_{t}\triangleq(u_{t},u_{t+1},\dots,u_{t+H-1})$, which results in a
state trajectory $X_{t}\triangleq(x_{t},x_{t+1},\cdots,x_{t+H})$. The total
cost incurred is

$C(X_{t},U_{t})=\sum_{h=0}^{H-1}c(x_{t+h},u_{t+h})+c_{term}(x_{t+H}),\vspace{-0.5ex}$
(2)

where $c_{term}(\cdot)$ is a terminal cost function. The goal of MPC is to
find the optimal set of parameters
$\bm{\theta}_{t}\triangleq(\theta_{t},\theta_{t+1},\cdots,\theta_{t+H-1})$ for
the sequence of policies
$\bm{\pi}_{\bm{\theta}_{t}}\triangleq(\pi_{\theta_{t}},\pi_{\theta_{t+1}},\cdots,\pi_{\theta_{t+H-1}})$.
At each time step, we solve

$\vspace{-0.5ex}\bm{\theta}_{t}\leftarrow\operatorname*{arg\,min}_{\bm{\theta}\in\bm{\Theta}}J(\bm{\pi}_{\bm{\theta}};x_{t}),$
(3)

where $J(\cdot)$ is a statistic defined on cost $C(X_{t},U_{t})$ such that its
minimum occurs at the optimal $\bm{\theta}_{t}$. In general, we do not have
access to the true dynamics function $f$ and instead approximate it with the
model $\hat{f}$, corresponding to the surrogate statistic
$\hat{J}(\bm{\pi}_{\bm{\theta}};x_{t})$.

This optimization problem can only be approximated in practice due to real-
time constraints. One commonly applied heuristic is to bootstrap the previous
approximate solution as an initialization for the current problem. This is
effective because the optimization problems between two consecutive time steps
share all control variables except the first and last. If our solution from
the previous problem is $\bm{\theta}_{t-1}$, then our warm start for the
current problem is given by

$\tilde{\bm{\theta}}_{t}=\Phi(\bm{\theta}_{t-1}),$ (4)

where $\Phi(\cdot)$ is called the shift operator [4]. A common choice is
$\tilde{\bm{\theta}}_{t}=(\theta_{t+1},\theta_{t+2},\dots,\theta_{t+H-1},\bar{\theta})$,
where $\bar{\theta}$ is a new parameter which reflects the expected final
action.

Recent work by Wagener et al. [4] showed that many common MPC algorithms fall
under the framework of an online learning algorithm known as dynamic mirror
descent (DMD) [8]. Online learning involves interactions between a learner and
an environment over $T$ rounds. In our case, the learner is the MPC algorithm,
which in round $t$ plays the decision $\tilde{\bm{\theta}}_{t}\in\bm{\Theta}$,
the shifted policy parameter sequence, along with side information $u_{t-1}$,
the control applied to the real system. The per-round loss is defined as
$\ell_{t}(\cdot)=\hat{J}(\cdot;x_{t})$, which is selected by the environment
via the state transition. At round $t$, DMD updates the parameters by the
rule:

$\bm{\theta}_{t}\leftarrow\operatorname*{arg\,min}_{\bm{\theta}\in\bm{\Theta}}\langle\gamma_{t}g_{t},\bm{\theta}\rangle+D_{\psi}(\bm{\theta}||\tilde{\bm{\theta}}_{t}),$
(5)

where $g_{t}=\nabla\ell_{t}(\tilde{\bm{\theta}}_{t})$, $\gamma_{t}>0$ is the
step size, and
$D_{\psi}(\bm{\theta}||\bm{\theta^{\prime}})=\psi(\bm{\theta})-\psi(\bm{\theta}^{\prime})-\langle\nabla\psi(\bm{\theta}^{\prime}),\bm{\theta}-\bm{\theta}^{\prime}\rangle$
is the Bregman divergence generated by a strictly convex function $\psi$ on
$\bm{\theta}$. Solving MPC with DMD is known as DMD-MPC and includes a family
of common MPC algorithms as special cases [4].

###  II-B Sampling-Based Model Predictive Control

A popular, practical sampling-based MPC algorithm is Model Predictive Path
Integral (MPPI) control [2, 3], which is a special case of DMD-MPC under
certain choices of objective function, control distribution, and Bregman
divergence [4]. Specifically, we assume the policies are open loop and choose
the exponential utility for the objective:

$\ell_{t}(\bm{\theta})=-\log\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}\exp\Big{(}-\frac{1}{\lambda}C(X_{t},U_{t})\Big{)}\Big{]},$
(6)

where $\lambda>0$ is a scaling parameter, also known as the temperature. Since
we generally assume that the cost function is non-differentiable with respect
to $\bm{\theta}$, we instead compute the gradients via a likelihood-ratio
derivative:

$\small\nabla\ell_{t}(\bm{\theta})=-\frac{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}\exp\Big{(}-\frac{1}{\lambda}C(X_{t},U_{t})\Big{)}\nabla_{\bm{\theta}}\log\pi_{\bm{\theta}}(U_{t})\Big{]}}{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}\exp\Big{(}-\frac{1}{\lambda}C(X_{t},U_{t})\Big{)}\Big{]}},$
(7)

We approximate these expectations with Monte Carlo sampling, which results in
a convex combination of gradients:

$\nabla\ell_{t}(\bm{\theta})=-\sum_{i=1}^{N}w_{i}\nabla_{\bm{\theta}}\log\pi_{\bm{\theta}}(U_{t}),\vspace{-0.5ex}$
(8)

with weights $w_{i}$ defined by the softmax operation

$w_{i}=\frac{e^{-\frac{1}{\lambda}C(X_{t}^{(i)},U_{t}^{(i)})}}{\sum_{j=1}^{N}e^{-\frac{1}{\lambda}C(X_{t}^{(j)},U_{t}^{(j)})}}.\vspace{-0.5ex}$
(9)

Because of this estimation, the gradients will be noisy, with noise that
scales inversely with number of samples. The more samples we use, the more
exact our approximations will be.

Next, we choose the Bregman divergence to be the KL divergence and the
distribution to be a factorized Gaussian:

$\small\pi_{\bm{\theta}}(U_{t})=\prod_{h=0}^{H-1}\pi_{\theta_{h}}(u_{t+h})=\prod_{h=0}^{H-1}\mathcal{N}(u_{t+h};\mu_{t+h},\Sigma_{t+h}),$
(10)

for some mean vectors $\mu_{t+h}$ and covariance matrices $\Sigma_{t+h}$.
Under these assumptions, Wagener et al. [4] showed that the solution to
Equation 5 becomes

$\displaystyle\mu_{t+h}$
$\displaystyle=(1-\gamma_{t}^{\mu})\tilde{\mu}_{t+h}+\gamma_{t}^{\mu}\frac{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}e^{-\frac{1}{\lambda}C(X_{t},U_{t})}u_{t+h}\Big{]}}{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}e^{-\frac{1}{\lambda}C(X_{t},U_{t})}\Big{]}}$
(11) $\displaystyle\Sigma_{t+h}$
$\displaystyle=(1-\gamma_{t}^{\sigma})\tilde{\Sigma}_{t+h}+\gamma_{t}^{\sigma}\frac{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}e^{-\frac{1}{\lambda}C(X_{t},U_{t})}m_{t+h}m_{t+h}^{T}\Big{]}}{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{f}}\Big{[}e^{-\frac{1}{\lambda}C(X_{t},U_{t})}\Big{]}},$

where $m_{t+h}=u_{t+h}-\mu_{t+h}$ and $\gamma_{t}^{\mu}$ and
$\gamma_{t}^{\sigma}$ are step sizes for the mean and covariance,
respectively. When we use the Monte-Carlo approximations to the gradients in
Equation 8, this results in the following update:

$\displaystyle\mu_{t+h}$
$\displaystyle=(1-\gamma_{t}^{\mu})\tilde{\mu}_{t+h}+\gamma_{t}^{\mu}\sum_{i=1}^{N}w_{i}u_{t+h}^{(i)}$
(12) $\displaystyle\Sigma_{t+h}$
$\displaystyle=(1-\gamma_{t}^{\sigma})\tilde{\Sigma}_{t+h}+\gamma_{t}^{\sigma}\sum_{i=1}^{N}w_{i}m_{t+h}^{(i)}m_{t+h}^{{(i)}^{T}},$

This update equation is the same as the MPPI update rule when
$\gamma_{t}^{\mu}=1$ and we do not update the covariance.

###  II-C Learning to Optimize Framework

Rather than hand design an update rule tailored to a specific subclass of
problems, the learning-to-optimize approach aims to learn a sequential update
rule from experience. For a set of optimizee parameters $\theta\in\Theta$ and
objective function $\ell(\theta)$, we find the minimizer
$\theta^{*}=\operatorname*{arg\,min}_{\theta\in\Theta}\ell(\theta)$ with an
iterative algorithm that has the update rule

$\theta_{t+1}=m_{\phi}(\theta_{t},t),\vspace{-0.5ex}$ (13)

where $m$ is the learned optimizer, which can be of any parameterized function
class with parameters $\phi$.

The majority of approaches to learning-to-optimize differentiate through the
optimization process using gradient descent [39, 40, 41, 42, 43] or use
reinforcement learning [44, 45, 46, 47]. However, we do not assume that the
optimization process is end-to-end differentiable and wish to avoid the high
sample complexity that often hinders reinforcement learning. Instead, we opt
to use imitation learning to train the optimizer. Chen et al. [43] also make
use of imitation learning, in which the experts are common hand-designed
optimizers that have access to full gradient information. However, in our
case, the learned optimizers only have access to noisier gradients than the
expert demonstrator and therefore less information.

Another major difference from prior work is that most literature in this area
targets optimizing deep neural networks. As such, they must contend with the
large parameter space of these models. Andrychowicz et al. [39] proposed to
use a coordinate-wise optimizer, in which the parameters of the optimizer are
shared across updates for all optimizee parameters. A downside to this
approach is that it throws away potentially useful information for improving
the learned update. Instead, since we have a moderate number of parameters, we
can jointly optimize the entire planning horizon of control distribution
parameters in order to capture relationships between time steps.

##  III Learning to Optimize for Control

###  III-A Design of the Learnable Optimzier

As shown in the previous section, the MPPI update rule in Equation 12
corresponds to performing mirror descent with an approximate gradient computed
from $N$ samples. Fewer samples results in a worse approximation and,
therefore, a noisier update. One possible avenue for improving performance
would be to employ more advanced first-order methods [48, 49, 50, 51, 52, 53].
However, by adopting the learning-to-optimize framework and replacing the
update rule with a learned optimizer, we can potentially do better than a
manually specified update. A naive approach would be to follow Equation 13 and
use the noisy gradient as input to the learned optimizer to produce the
updated parameters. However, our objective is to learn how to mitigate the
effect of a low number of samples on gradient noise, and the computation of
the noisy gradient itself potentially throws away information that may be
useful for improving the update. For instance, looking at Equation 12, we are
simply computing a weighted sum of the samples. This collapses the information
in each trajectory sample and its corresponding cost into a single vector.
Therefore, we propose instead to use the individual components which form the
gradient directly.

From Equation 12, we can see that the update is a function of the current mean
$\tilde{\mu}_{t+h}$ and covariance $\tilde{\Sigma}_{t+h}$, sample weights
$w^{(1:N)}_{t}$, and control samples $u_{t+h}^{(1:N)}$. The sample weights
themselves are actually a function of the total trajectory costs
$C^{(1:N)}_{t}$, where $C_{t}=C(X_{t},U_{t})$. One potential choice would be
to make each of these terms an input to the learned update:

$\mu_{t+h},\Sigma_{t+h}=m_{\phi}\Big{(}\tilde{\mu}_{t+h},\tilde{\Sigma}_{t+h},C^{(1:N)}_{t},u_{t+h}^{(1:N)}\Big{)}.$
(14)

A limitation of this choice of parameterization is that it assumes
independence of the updates between time steps in the rollouts. While this is
the case for vanilla MPPI [2, 3], we could potentially learn a better update
by incorporating information across time steps. However, if we parameterize
the optimizer with a fully-connected or recurrent neural network architecture,
this would result in a large number of parameters to learn, making
optimization difficult. Instead, we alter the way in which we sample from the
Gaussian policies to remedy this explosion in the dimensionality.

As proposed by Bhardwaj et al. [6], we make use of low-discrepancy Halton
sequences [54] to generate samples from the Gaussian policies. Normal pseudo-
random sequences often result in clusters of sampled points, leaving many
regions of the parameter space untouched. Low-discrepancy sequences are a
deterministic alternative that alleviate this problem by correlating each
point. A $D$-dimensional Halton sequence $x_{1},x_{2},\dots,x_{N}$, in which
$x_{i}\in\mathbb{R}^{D}$ is generated by

$x_{i}=(\phi_{p_{1}}(i),\dots,\phi_{p_{D}}(i)),\
\phi_{p_{b}}(i)=\sum_{j=1}^{\infty}a_{j}(p_{b})p_{b}^{-j},\vspace{-0.5ex}$
(15)

where $p_{1},\dots,p_{D}$ are consecutive prime numbers and
$a_{j}(p_{b})\in\\{0,1,\dots,p_{b}-1\\}$ such that the condition
$i=\sum_{j=1}^{\infty}a_{j}(p_{b})p_{b}^{j-1}$ holds. The Halton sequence is
sampled once at the beginning of the rollout and then transformed using the
mean and covariance of the Gaussian policy. While this can improve the
performance of sampling-based MPC, the main benefit is that it makes all
sampled control sequences a deterministic function of the current mean and
covariance. Therefore, the sampled control sequences can be excluded from the
learned update without loss of information.

As such, rather than optimizing each time step independently, we leverage the
structured nature of these samples to learn an update that optimizes the
entire trajectory jointly using only cost information. The resulting update is
then

$\bm{\mu}_{t},\bm{\Sigma}_{t}=m_{\phi}\Big{(}\tilde{\bm{\mu}}_{t},\tilde{\bm{\Sigma}}_{t},C^{(1:N)}_{t}\Big{)},$
(16)

where $\bm{\mu}_{t}\triangleq(\mu_{t},\mu_{t+1},\dots,\mu_{t+H-1})$,
$\bm{\Sigma}_{t}\triangleq(\Sigma_{t},\Sigma_{t+1},\dots,\Sigma_{t+H-1})$, and
$\tilde{\bm{\mu}}_{t}$ and $\tilde{\bm{\Sigma}}_{t}$ are defined similarly. We
can think about the Halton sequence as giving us a sense of what the
environment and cost landscape is like around the current state. Since the
learned optimizer can potentially make better use of its inputs than the
expert, we may be able to more effectively use fewer samples while maintaining
similar performance.

Finally, we note that Equation 12 is a convex combination of the previous
control parameters and the weighted samples. We can actually think about MPC
as a form of recurrent network, with the warm-started control distribution as
our form of memory about the previous time steps. From this perspective, the
step size is acting as a gating term which modulates how much information we
preserve about our history. The hidden state update in a gated recurrent unit
(GRU) [55] is of the same form, except the multiplicative gating term is also
learned. Inspired by this similarity, we use the following update:

$\begin{gathered}g_{t}^{\mu},g_{t}^{\sigma},h_{t}^{\mu},h_{t}^{\sigma}=m_{\phi}(\tilde{\bm{\mu}}_{t},\tilde{\bm{\Sigma}}_{t},C_{t}^{(1:N)})\\\
\bm{\mu}_{t}=(1-g_{t}^{\mu})\odot\tilde{\bm{\mu}}_{t}+g_{t}^{\mu}\odot
h_{t}^{\mu}\\\
\bm{\Sigma}_{t}=(1-g_{t}^{\sigma})\odot\tilde{\bm{\Sigma}}_{t}+g_{t}^{\sigma}\odot
h_{t}^{\sigma},\end{gathered}$ (17)

where $\odot$ is the Hadamard product. Here, $g_{t}^{\mu},g_{t}^{\sigma}$ are
the learned gating terms, which are passed through a sigmoid to ensure they
are between zero and one. Meanwhile, $h_{t}^{\mu},h_{t}^{\sigma}$ can be
interpreted as the updates proposed by the network. In our experiments, this
choice of parameterization significantly outperformed a simple fully-connected
network.

###  III-B Imitation Learning for Training the Optimizer

Input: Initial bootstrapped dataset $\mathcal{D}$, initial policy
$\bm{\pi}_{{\tilde{\bm{\theta}}}_{1}}$, initial state distribution $\rho$

Parameters: Iterations $K$, probabilities $\\{\beta_{k}\\}_{k=1}^{K}$,
rollouts per iteration $R$

for  _$k=1,2,\dots,K$_ do

Initialize dataset $\mathcal{D}_{i}\leftarrow\emptyset$

for  _r = 1, 2, …, R_ do

Sample initial state $x_{1}\sim\rho$

$\tilde{\bm{\theta}}_{1:T}$, $C_{1:T}^{(1:M)}$, $\bm{\theta}_{1:T}^{expert}$
$\leftarrow$ Rollout($x_{1},\bm{\pi}_{\tilde{\bm{\theta}}_{1}},\beta_{k}$)

Append
$\mathcal{D}_{i}\leftarrow\mathcal{D}_{i}\cup\\{(\tilde{\bm{\theta}}_{1:T},C_{1:T}^{(1:M)},\bm{\theta}_{1:T}^{expert})\\}$

end for

Aggregate datasets $\mathcal{D}\leftarrow\mathcal{D}\cup\mathcal{D}_{i}$

Train optimizer parameters $\phi$ on $\mathcal{D}$

end for

Algorithm 1 DAgger Training Loop

Unlike prior work in learning-to-optimize, we cannot assume that the
optimization process itself is differentiable, as it occurs online via
interactions with the environment. Even in the simulated case, we do not want
to assume that everything has been implemented in a differentiable fashion. We
could use reinforcement learning (RL), although it generally has high sample
complexity and may be slow to learn. Since we have access to a tuned optimal
controller, imitation learning is a promising direction for training the
optimizer. Our expert is an MPPI controller with unrestricted access to
samples. That is, we provide the controller with as many samples as needed to
achieve good performance. The learner is also an MPPI controller, but it has
access to fewer samples, and the standard update is replaced with the learned
optimizer. In our preliminary experiments, we tried using standard behavioral
cloning, in which we collect a dataset of expert demonstrations and train a
policy offline via regression. However, this did not work well due to
covariate shift between the expert and learner distributions. Instead, we used
DAgger [56] to perform imitation learning, which is an interactive algorithm
that aims to combat issues of covariate shift. The algorithm queries an expert
online for corrective labels on learner visited states. We outline the main
loop of DAgger in Algorithm 1.

Input: State $x_{1}$, policy $\bm{\pi}_{\tilde{\bm{\theta}}_{1}}$, probability
$\beta_{k}$

Parameters: Rollout length $T$, expert samples $N$, learner samples $M$

Output: Shifted parameters $\tilde{\bm{\theta}}_{1:T}$, sample costs
$C_{1:T}^{(1:M)}$, expert decisions $\bm{\theta}_{1:T}^{expert}$

for  _$t=1,2,\dots,T$_ do

Sample controls from policy
$\\{U_{t}^{(i)}\\}_{i=1}^{N}\sim\bm{\bm{\pi}}_{\tilde{\bm{\theta}}_{t}}$

Sample $X_{t}^{(i)}$ from dynamics $\hat{f}$ using $U_{t}^{(i)},x_{t}$

Compute costs $C_{t}^{(i)}\leftarrow C(X_{t}^{(i)},U_{t}^{(i)})$

Compute sample weights with Equation 9

Update $\bm{\tilde{\theta}}_{t}$ to $\bm{\theta}_{t}^{expert}$ using Equation
12

Sample $b\sim U(0,1)$

if  _$b\leq\beta_{k}$_ then

Set $\bm{\theta}_{t}\leftarrow\bm{\theta}_{t}^{expert}$

else

Update $\bm{\tilde{\theta}}_{t}$ to $\bm{\theta}_{t}^{learn}$ using Equation
17

Set $\bm{\theta}_{t}\leftarrow\bm{\theta}_{t}^{learn}$

end if

Sample $u_{t}\sim\pi_{\theta_{t}}$ or use mean $u_{t}\leftarrow\mu_{t}$

Apply control to system $x_{t+1}\sim f(x_{t},u_{t})$

Shift parameters $\tilde{\bm{\theta}}_{t+1}=\Phi(\bm{\theta}_{t})$

end for

Algorithm 2 DAgger Rollout Function

First, we begin by collecting a bootstrap dataset in which only the expert is
run. Next, each iteration $k$ of DAgger, we run $R$ rollouts according to
Algorithm 2 by sampling some initial state $x_{1}$ from a known initial state
distribution $\rho$. During a rollout, at each time step, we apply the
controls from the expert with probability $\beta_{k}$ and the learner with
probability $1-\beta_{k}$. The expert is always run in order to provide a
corrective target for training the policy at the next iteration. Both the
expert and learner controllers use the same trajectory samples, although the
learner only receives a subset of them. Generally, the mixing probabilities
$\beta_{k}$ are set according to a schedule such that we run the learner more
often in later iterations. In our experiments, we set $\beta_{k}=p^{k}$ for
some $p\in(0,1)$. After running the rollouts, we collect the warm-started
control distribution parameters $\tilde{\bm{\theta}}_{1:T}$, the trajectory
costs $C_{1:T}^{(1:M)}$, and the updated expert control distribution
parameters $\bm{\theta}_{1:T}^{expert}$ into a dataset $\mathcal{D}_{i}$. We
only collect the $M\leq N$ samples used by the learner. This data is then
aggregated into our main dataset $\mathcal{D}$ to train the optimizer.

##  IV Experiments

Implementation Details. In all experiments, our MPPI implementation is a
modified version of the one developed by Bhardwaj et al. [6]. This
implementation uses Halton sequences for generating control sequence samples
and smooths the sampled trajectories with 3rd degree B-splines. We use a fixed
diagonal covariance for the sampling distribution and do not perform
covariance adaptation. All hyperparameters were tuned using a grid search, and
the optimal number of samples is what the expert controller has access to
during data generation and training. Now, the optimal choice of
hyperparameters may be different for a given number of samples. Therefore, for
a fair comparison, we tune the MPPI hyperparameters separately for each sample
count used in our evaluation. Both MPPI and the neural networks are
implemented in PyTorch [57].

Task Details. We evaluate on simulated tasks:

  1. 1.

Cartpole: The task is to slide a cart along a rail to swing up the pole
attached via an unactuated joint using only actuation from the cart. Both the
expert and learner are given access to the true analytical dynamics. The
initial position of the cart and pole are randomized at every episode, which
lasts 200 time steps. An episode is considered successful if the pole is swung
up with a linear and angular velocity near zero.

  2. 2.

Franka Reacher: A 7 degree-of-freedom (DOF) Franka Panda robot arm must reach
a target goal from a fixed starting pose. The goal is randomly selected at the
beginning of each episode. Both the expert and learner use the same kinematic
model described in Bhardwaj et al. [6], which is different from the true
dynamics of the simulator (Nvidia’s Isaac Gym [58]). Each episode lasts for
500 time steps and is considered successful if the end effector reaches the
target position.

  3. 3.

Franka Obstacles: This task is identical to Franka Reacher, except now there
are two spherical obstacles placed in the environment which the arm must
avoid. The obstacle and the goal positions are randomized at the beginning of
each episode, which lasts for 600 time steps. An episode is considered
successful if the end effector reaches the goal while avoiding collisions.

Evaluation. We evaluate the performance of the learned optimizer by varying
the number of samples, up to the amount used by the expert. For each sample
amount, we compare against a standard MPPI implementation with access to the
same number of samples as the learned controller over 30 test rollouts. All
test rollouts use a fixed set of start states, goals, and obstacle locations.
This is achieved by setting the random seed value to a pre-defined test seed.
Our primary metric for comparison is success rate, which is defined as the
percentage of times the task goal was achieved out of all trials. In the
Franka Obstacles task, the placement of obstacles is randomized according to a
pre-specified distribution. Therefore, this task allows us to evaluate the
generalization capability of the learned optimizer to new environments which
are drawn from a similar distribution. Additionally, we report statistics of
the end effector test trajectories. Specifically, we compute a relative
trajectory length and average jerk as the ratio between the statistics for the
learned optimizer and baseline MPPI controller. The trajectory length is
averaged over all test runs, while the jerk is averaged over only successful
test runs.

Training Details. Prior work in learning-to-optimize made use of recurrent
architectures which can account for the history of the gradients and
optimization process. While we could potentially benefit from such an
architecture, we found that a simple multi-layer perceptron (MLP) was
sufficient to learn powerful optimizers. As such, in all experiments, the
learned optimizer is represented with a two-layer MLP using ReLU activation
functions. We use 1024, 2048, and 4096 hidden units per layer for the
Cartpole, Franka Reacher, and Franka Obstacles tasks, respectively. To prevent
overfitting, all networks are regularized with dropout [59] using a dropout
probability of 0.1. We use the Adam optimizer [52] with a learning rate of
$10^{-3}$ for Cartpole and Franka Reacher and $10^{-4}$ for Franka Obstacles.
We normalize the total trajectory costs based on the mean and standard
deviation of the training dataset. For all tasks, we bootstrapped the dataset
with 1024 trajectories, in which only the expert’s action was applied to the
system. We ran DAgger for 20 iterations with 128 rollouts per iteration and a
mixing probability schedule $\beta_{k}=0.8^{k}$. For each iteration, we train
the networks on the aggregated dataset for 1000 epochs with a batch size of 8.
The best performing network evaluated on a held-out validation set is saved
and used in the next iteration of DAgger.

TABLE I: Success Rate Across All Tasks. | # Samples | MPPI | L2O-MPC
---|---|---|---
Cartpole | 8 | 100.0 | 100.0
4 | 60.00 | 96.67
2 | 10.00 | 90.00
Franka Reacher | 64 | 100.0 | 100.0
32 | 90.00 | 100.0
16 | 63.33 | 100.0
8 | 20.00 | 80.00
4 | 10.00 | 63.33
2 | 3.333 | 16.67
Franka Obstacles | 512 | 80.00 | 80.00
256 | 76.67 | 80.00
128 | 73.33 | 76.67
64 | 63.33 | 76.67
32 | 33.33 | 66.67
16 | 6.667 | 46.67

##  V Results

We report the success rate for all tasks in Table I and refer to standard MPPI
by MPPI and MPC with the learned optimizer by L2O-MPC. Success rate is
computed for each task based on the criteria discussed in Section IV. For
Franka Obstacles, not every randomly generated scenario is feasible. Hence,
there is an upper-bound of an 80% success rate for both MPPI and L2O-MPC. We
can see that the performance of MPPI quickly drops off as the number of
control sequence samples is reduced. For Cartpole, the performance of L2O-MPC
remains fairly consistent even with a lower number of samples. While the
performance drop is more pronounced in the Franka experiments, L2O-MPC still
consistently matches or outperforms MPPI at each sample amount. In Franka
Reacher, L2O-MPC is able to withstand a $4\times$ decrease in the number of
samples while still achieving an $100\%$ success rate. Similarly, in Franka
Obstacles, L2O-MPC only incurs a $4\%$ decrease in performance under an
$8\times$ decrease in number of samples. This illustrates that L2O-MPC is
successfully able to generalize to new environments similar to those on which
it was trained.

Figure 1: Example trajectories of the Franka Obstacles task, in which the
Franka arm end effector (green) is tasked with reaching the goal (red) while
avoiding obstacles. The top row is a novel environment with three obstacles,
while the bottom row is an environment from the test set.

Figure 2: Trajectories of the Franka arm end effector when controlled by MPPI with 512 samples (blue), MPPI with 16 samples (pink), and L2O-MPC with 16 samples (purple) to move from the starting position (red) to the goal (green) while avoiding obstacles (cyan). Plots in the same row are from the same environment but viewed from differing perspectives. TABLE II: Trajectory Statistics for Franka Obstacles. # Samples | Length | Avg. Jerk
---|---|---
512 | 1.091 | 1.008
256 | 1.030 | 1.023
128 | 0.975 | 1.049
64 | 0.916 | 1.077
32 | 0.858 | 1.104
16 | 0.854 | 1.527

Figure 1 illustrates two different example trajectories of the Franka arm
avoiding different amounts of obstacles. These same environments are depicted
in Figure 2, where we show end effector trajectories under MPPI with 512
samples, MPPI with 16 samples, and L2O-MPC with 16 samples. We can see that
MPPI with 16 samples quickly diverges and is unable to reach the goal. On the
other hand, L2O-MPC with 16 samples is able to better make use of the samples
and still reach the goal while avoiding all obstacles. Moreover, the L2O-MPC
controller was trained only on environments that contain two obstacles.
Therefore, these results also indicate that the learned controller may
generalize to novel environments on which it was not trained.

Qualitatively, the L2O-MPC trajectories appear to be slightly more jittery
than the MPPI expert. In Table II, we provide the average relative jerk
between L2O-MPC and MPPI for successful test runs at different sample counts.
Indeed, we see that the L2O-MPC trajectories are less smooth than those of
MPPI, and this effect is exacerbated at lower sample counts. Additionally, we
provide the average relative trajectory length across all test environments.
With more samples, L2O-MPC has slightly longer trajectories than MPPI,
indicating that it is slower at reaching the goal. However, when given access
to fewer samples, L2O-MPC consistently has shorter trajectories, as it more
often reaches the goal. Therefore, while L2O-MPC is often jerkier and
sometimes slower than the expert with full samples, it succeeds more often in
achieving the desired objective in a timely fashion than MPPI given the same
number of samples.

##  VI Conclusion

We presented a method for improving upon standard sampling-based MPC
algorithms by learning a better update rule. This provides a novel way to
incorporate learning into model-based control algorithms, which is orthogonal
to the standard approaches of learning or fine-tuning the dynamics model
and/or cost function. We contend with noisy gradients by learning how to more
effectively update the control distribution. By using structured sampling
strategies, we are able to provide more information to the learned update and
better utilize fewer samples. We show through empirical evaluations that our
learned controllers remain competitive or outperform a baseline MPPI
controller with access to the same number of samples. This demonstrates the
viability of the learning-to-optimize framework in the context of control,
opening the door for a variety of techniques to be applied to improving the
performance of optimization-based controllers and planners. While we leveraged
imitation learning to train the optimizers, this is just one possible option
and an interesting direction for future work is to explore using reinforcement
learning to see if it can outperform the expert and model-free methods. Since
performance of sampling-based methods relies so heavily on thorough
exploration of the sample space, another possible avenue is to learn how to
generate better samples in addition to better updates.

## References

  * [1] P. Abbeel, A. Coates, and A. Y. Ng, “Autonomous Helicopter Aerobatics through Apprenticeship Learning,” _The International Journal of Robotics Research_ , vol. 29, no. 13, pp. 1608–1639, 2010\.
  * [2] G. Williams, P. Drews, B. Goldfain, J. M. Rehg, and E. A. Theodorou, “Aggressive driving with model predictive path integral control,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 1433–1440.
  * [3] G. Williams, N. Wagener, B. Goldfain, P. Drews, J. M. Rehg, B. Boots, and E. A. Theodorou, “Information theoretic MPC for model-based reinforcement learning,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 1714–1721.
  * [4] N. Wagener, C.-A. Cheng, J. Sacks, and B. Boots, “An Online Learning Approach to Model Predictive Control,” _arXiv preprint arXiv:1902.08967_ , 2019.
  * [5] V. Kumar, E. Todorov, and S. Levine, “Optimal Control with Learned Local Models: Application to Dexterous Manipulation,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 378–383.
  * [6] M. Bhardwaj, B. Sundaralingam, A. Mousavian, N. Ratliff, D. Fox, F. Ramos, and B. Boots, “STORM: An Integrated Framework for Fast Joint-Space Model-Predictive Control for Reactive Manipulation,” _arXiv preprint arXiv:2104.13542_ , 2021.
  * [7] T. Erez, K. Lowrey, Y. Tassa, V. Kumar, S. Kolev, and E. Todorov, “An integrated system for real-time Model Predictive Control of humanoid robots,” in _2013 13th IEEE-RAS International Conference on Humanoid Robots (Humanoids)_. IEEE, 2013, pp. 292–299.
  * [8] E. Hall and R. Willett, “Dynamical Models and Tracking Regret in Online Convex Programming,” in _International Conference on Machine Learning (ICML)_. PMLR, 2013, pp. 579–587.
  * [9] T. Chen, X. Chen, W. Chen, H. Heaton, J. Liu, Z. Wang, and W. Yin, “Learning to Optimize: A Primer and A Benchmark,” _arXiv preprint arXiv:2103.12828_ , 2021.
  * [10] J. Kocijan, R. Murray-Smith, C. E. Rasmussen, and A. Girard, “Gaussian process model based predictive control,” in _Proceedings of the 2004 American Control Conference (ACC)_ , vol. 3. IEEE, 2004, pp. 2214–2219.
  * [11] I. Lenz, R. A. Knepper, and A. Saxena, “DeepMPC: Learning Deep Latent Features for Model Predictive Control,” in _Robotics: Science and Systems (RSS)_. Rome, Italy, 2015.
  * [12] J. Fu, S. Levine, and P. Abbeel, “One-Shot Learning of Manipulation Skills with Online Dynamics Adaptation and Neural Network Priors,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 4019–4026.
  * [13] A. Nagabandi, G. Kahn, R. S. Fearing, and S. Levine, “Neural Network Dynamics for Model-Based Deep Reinforcement Learning with Model-Free Fine-Tuning,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 7559–7566.
  * [14] K. Chua, R. Calandra, R. McAllister, and S. Levine, “Deep Reinforcement Learning in a Handful of Trials using Probabilistic Dynamics Models,” _arXiv preprint arXiv:1805.12114_ , 2018.
  * [15] D. Hafner, T. Lillicrap, I. Fischer, R. Villegas, D. Ha, H. Lee, and J. Davidson, “Learning Latent Dynamics for Planning from Pixels,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 2555–2565.
  * [16] A. Tamar, G. Thomas, T. Zhang, S. Levine, and P. Abbeel, “Learning from the Hindsight Plan — Episodic MPC Improvement,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 336–343.
  * [17] M. Zhong, M. Johnson, Y. Tassa, T. Erez, and E. Todorov, “Value Function Approximation and Model Predictive Control,” in _2013 IEEE Symposium on Adaptive Dynamic Programming and Reinforcement Learning (ADPRL)_. IEEE, 2013, pp. 100–107.
  * [18] U. Rosolia and F. Borrelli, “Learning Model Predictive Control for iterative tasks. A Data-Driven Control Framework,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 7, pp. 1883–1896, 2017.
  * [19] K. Lowrey, A. Rajeswaran, S. Kakade, E. Todorov, and I. Mordatch, “Plan Online, Learn Offline: Efficient Learning and Exploration via Model-Based Control,” _arXiv preprint arXiv:1811.01848_ , 2018.
  * [20] M. Bhardwaj, A. Handa, D. Fox, and B. Boots, “Information Theoretic Model Predictive Q-Learning,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 840–850.
  * [21] T. Wang and J. Ba, “Exploring Model-Based Planning with Policy Networks,” _arXiv preprint arXiv:1906.08649_ , 2019.
  * [22] N. Wahlström, T. B. Schön, and M. P. Deisenroth, “From Pixels to Torques: Policy Learning with Deep Dynamical Models,” _arXiv preprint arXiv:1502.02251_ , 2015.
  * [23] M. Watter, J. T. Springenberg, J. Boedecker, and M. Riedmiller, “Embed to Control: A Locally Linear Latent Dynamics Model for Control from Raw Images,” _arXiv preprint arXiv:1506.07365_ , 2015.
  * [24] C. Finn and S. Levine, “Deep Visual Foresight for Planning Robot Motion,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 2786–2793.
  * [25] E. Banijamali, R. Shu, H. Bui, and A. Ghodsi, “Robust Locally-Linear Controllable Embedding,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2018, pp. 1751–1759.
  * [26] F. Ebert, C. Finn, S. Dasari, A. Xie, A. Lee, and S. Levine, “Visual Foresight: Model-Based Deep Reinforcement Learning for Vision-Based Robotic Control,” _arXiv preprint arXiv:1812.00568_ , 2018.
  * [27] J.-S. Ha, Y.-J. Park, H.-J. Chae, S.-S. Park, and H.-L. Choi, “Adaptive Path-Integral Autoencoder: Representation Learning and Planning for Dynamical Systems,” _32nd Conference on Neural Information Processing Systems (NeurIPS)_ , 2018\.
  * [28] A. Sharma, S. Gu, S. Levine, V. Kumar, and K. Hausman, “Dynamics-Aware Unsupervised Discovery of Skills,” _arXiv preprint arXiv:1907.01657_ , 2019.
  * [29] B. Amos and D. Yarats, “The Differentiable Cross-Entropy Method,” in _International Conference on Machine Learning (ICML)_. PMLR, 2020, pp. 291–302.
  * [30] B. Amos, I. D. J. Rodriguez, J. Sacks, B. Boots, and J. Z. Kolter, “Differentiable MPC for End-to-end Planning and Control,” _arXiv preprint arXiv:1810.13400_ , 2018\.
  * [31] P. Karkus, D. Hsu, and W. S. Lee, “QMDP-Net: Deep Learning for Planning under Partial Observability,” _arXiv preprint arXiv:1703.06692_ , 2017.
  * [32] M. Okada, L. Rigazio, and T. Aoshima, “Path Integral Networks: End-to-End Differentiable Optimal Control,” _arXiv preprint arXiv:1706.09597_ , 2017.
  * [33] M. Okada and T. Taniguchi, “Acceleration of Gradient-based Path Integral Method for Efficient Optimal and Inverse Optimal Control,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 3013–3020.
  * [34] M. Pereira, D. D. Fan, G. N. An, and E. Theodorou, “MPC-Inspired Neural Network Policies for Sequential Decision Making,” _arXiv preprint arXiv:1802.05803_ , 2018.
  * [35] A. Tamar, Y. Wu, G. Thomas, S. Levine, and P. Abbeel, “Value Iteration Networks,” _arXiv preprint arXiv:1602.02867_ , 2016.
  * [36] A. Srinivas, A. Jabri, P. Abbeel, S. Levine, and C. Finn, “Universal Planning Networks: Learning Generalizable Representations for Visuomotor Control,” in _International Conference on Machine Learning (ICML)_. PMLR, 2018, pp. 4732–4741.
  * [37] T. Yu, G. Shevchuk, D. Sadigh, and C. Finn, “Unsupervised Visuomotor Control through Distributional Planning Networks,” _arXiv preprint arXiv:1902.05542_ , 2019.
  * [38] M. Bhardwaj, B. Boots, and M. Mukadam, “Differentiable Gaussian Process Motion Planning,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 10 598–10 604.
  * [39] M. Andrychowicz, M. Denil, S. Gomez, M. W. Hoffman, D. Pfau, T. Schaul, B. Shillingford, and N. De Freitas, “Learning to learn by gradient descent by gradient descent,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , 2016, pp. 3981–3989.
  * [40] S. Ravi and H. Larochelle, “Optimization as a Model for Few-Shot Learning,” _International Conference on Learning Representations (ICLR)_ , 2016.
  * [41] K. Lv, S. Jiang, and J. Li, “Learning Gradient Descent: Better Generalization and Longer Horizons,” in _International Conference on Machine Learning (ICML)_. PMLR, 2017, pp. 2247–2255.
  * [42] Y. Chen, M. W. Hoffman, S. G. Colmenarejo, M. Denil, T. P. Lillicrap, M. Botvinick, and N. Freitas, “Learning to Learn without Gradient Descent by Gradient Descent,” in _International Conference on Machine Learning (ICML)_. PMLR, 2017, pp. 748–756.
  * [43] T. Chen, W. Zhang, Z. Jingyang, S. Chang, S. Liu, L. Amini, and Z. Wang, “Training Stronger Baselines for Learning to Optimize,” _Advances in Neural Information Processing Systems_ , vol. 33, 2020.
  * [44] K. Li and J. Malik, “Learning to Optimize,” _arXiv preprint arXiv:1606.01885_ , 2016.
  * [45] J. X. Wang, Z. Kurth-Nelson, D. Tirumala, H. Soyer, J. Z. Leibo, R. Munos, C. Blundell, D. Kumaran, and M. Botvinick, “Learning to reinforcement learn,” _arXiv preprint arXiv:1611.05763_ , 2016.
  * [46] K. Li and J. Malik, “Learning to Optimize Neural Nets,” _arXiv preprint arXiv:1703.00441_ , 2017.
  * [47] I. Bello, B. Zoph, V. Vasudevan, and Q. V. Le, “Neural Optimizer Search with Reinforcement Learning,” in _International Conference on Machine Learning (ICML)_. PMLR, 2017, pp. 459–468.
  * [48] Y. Nesterov, “A method of solving a convex programming problem with convergence rate $\mathcal{O}(1/k^{2})$,” in _Sov. Math. Dokl_ , vol. 27.
  * [49] N. Qian, “On the momentum term in gradient descent learning algorithms,” _Neural Networks_ , vol. 12, no. 1, pp. 145–151, 1999.
  * [50] J. Duchi, E. Hazan, and Y. Singer, “Adaptive Subgradient Methods for Online Learning and Stochastic Optimization.” _Journal of Machine Learning Research (JMLR)_ , vol. 12, no. 7, 2011.
  * [51] M. D. Zeiler, “ADADELTA: An Adaptive Learning Rate Method,” _arXiv preprint arXiv:1212.5701_ , 2012.
  * [52] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [53] S. J. Reddi, S. Kale, and S. Kumar, “On the Convergence of Adam and Beyond,” _arXiv preprint arXiv:1904.09237_ , 2019.
  * [54] J. H. Halton, “Algorithm 247: Radical-inverse quasi-random point sequence,” _Communications of the ACM_ , vol. 7, no. 12, pp. 701–702, 1964.
  * [55] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation,” _arXiv preprint arXiv:1406.1078_ , 2014.
  * [56] S. Ross, G. Gordon, and D. Bagnell, “A Reduction of Imitation Learning and Structured Prediction to No-Regret Online Learning,” in _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics (AISTATS)_. JMLR Workshop and Conference Proceedings, 2011, pp. 627–635.
  * [57] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 32, pp. 8026–8037, 2019.
  * [58] V. Makoviychuk, L. Wawrzyniak, Y. Guo, M. Lu, K. Storey, M. Macklin, D. Hoeller, N. Rudin, A. Allshire, A. Handa, _et al._ , “Isaac Gym: High Performance GPU-Based Physics Simulation For Robot Learning,” _arXiv preprint arXiv:2108.10470_ , 2021.
  * [59] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: A Simple Way to Prevent Neural Networks from Overfitting,” _The Journal of Machine Learning Research (JMLR)_ , vol. 15, no. 1, pp. 1929–1958, 2014\.


# On Dark Matter As A Geometric Effect in the Galactic Halo

Ugur Camci  ucamci<EMAIL_ADDRESS>Department of Chemistry and
Physics, Roger Williams University, One Old Ferry Road, Bristol, RI 02809, USA

###### Abstract

We obtain more straightforwardly some features of dark matter distribution in
the halos of galaxies by considering the spherically symmetric space-time,
which satisfies the flat rotational curve condition, and the geometric
equation of state resulting from the modified gravity theory. In order to
measure the equation of state for dark matter in the galactic halo, we provide
a general formalism taking into account the modified $f(X)$ gravity theories.
Here, $f(X)$ is a general function of $X\in\\{R,\mathcal{G},T\\}$, where
$R,\mathcal{G}$ and $T$ are the Ricci scalar, the Gauss-Bonnet scalar and the
torsion scalar, respectively. These theories yield that the flat rotation
curves appear as a consequence of the additional geometric structure
accommodated by those of modified gravity theories. Constructing a geometric
equation of state $w_{{}_{X}}\equiv p_{{}_{X}}/\rho_{{}_{X}}$ and inspiring by
some values of the equation of state for the ordinary matter, we infer some
properties of dark matter in galactic halos of galaxies.

###### pacs:

98.35.Gi, 95.35.+d, 04.50.Kd

##  I Introduction

While the evolution of the universe is driven by the dark energy (DE), the
formation of galaxies should be mainly determined by the dark matter (DM)
which provides the potential wells where baryons collapse to originate the
visible component of galaxies. The direct evidence of DM and DE follows from
distance measurements of type Ia supernovae (SNe Ia) indicating that the
expansion of the universe is speeding up (Riess et al., 1998; Perlmutter et
al., 1999), which is the most important evidence of a breakdown of general
relativity (GR). The accelerated expansion of the universe can be explained by
different models of gravity than GR, where the Einstein-Hilbert Lagrangian is
generalized by some general functions such that $f(R),f(T),f(R,\mathcal{G}),$
etc. The cosmic microwave background (CMB) anisotropy measurements are equally
strong indirect evidence for that type of matter and energy. All of the
projects to probe DE have the common feature of surveying wide areas to
collect large samples of galaxies, clusters or supernovae. An additional
indirect evidence for DM comes from the connection between galaxies/clusters
and their halos, and for that type of structures, much higher mass estimates
than expected from the observed stars and gas indicate a breakdown of Newton’s
theory of gravity Wechsler and Tinker (2018). The modified gravity models,
which are commonly considering as gravitational alternative for DM, can
explain the extra matter content to be needed in addition to the observational
data coming from galaxies. So the gravitational properties of galaxies and
clusters have to be checked by considering the modified theories of gravity.

The equation of state (EoS) for the ordinary matter, $w=p/\rho$, can take an
arbitrary value between $-1$ and $1$, where $w=-1$ corresponds to a
cosmological constant, which is also matched to an exterior Schwarzschild
vacuum solution, $w=0$ for dust matter and $w=1$ for a stiff matter. The EoS
parameter $w=1/3$ corresponds to the relativistic matter (the radiation), like
photons and massless neutrinos. For quintessence models the EoS parameter
range is $-1<w<-1/3$, and the matter with the property $w<-1$ is dubbed the
phantom energy. Due to the apparent abundances of DM in most of the galactic
halos, we can suggest a geometric EoS of the form $w_{X}\equiv p_{X}/\rho_{X}$
to the spherically symmetric space-times for studying gravitational properties
of DM in the galactic halo. Then, using the field equations of modified
gravity model and ignoring the contribution of ordinary baryonic matter to
that field equations, one can write those of field equations in terms of
$\rho_{X}$ and $p_{X}$ which are $X\in\\{R,T,\mathcal{G},...\\}$ term
contributions to energy density and pressure, where $R,T$ and $\mathcal{G}$
are the Ricci scalar, the torsion scalar and the Gauss-Bonnet scalar which is
defined by
$\mathcal{G}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\alpha\beta\mu\nu}R^{\alpha\beta\mu\nu}$
with $R_{\mu\nu}$ and $R_{\alpha\beta\mu\nu}$ being the Ricci tensor and the
Riemann tensor, respectively. There are some studies of trying to measure the
EoS of DM in the galactic halo where it has been considered some halo models
of galaxies predicting a significant amount of pressure in the DM fluid such
that some variations of scalar field DM Bharadwaj and Kar (2003); Matos et al.
2000a ; Matos et al. 2000b ; Peeples (2000); Arbey et al. (2003), or string
fluid Soleng (1995).

The theoretical expectations for large amounts of luminous matter in galaxies
are not in agreement with observations. An insufficient amount of standard
matter in galaxies is the main problem, and the existence of some non-luminous
matters, the DM, should be necessary to explain both dynamics and structure
formation of galaxies. It is well known that the extended flat rotation curves
of spiral galaxies provide some of the strongest evidence for mass
discrepancies in galaxies, and are usually interpreted as evidence for DM
halos. The observed flat rotational curves in outskirts of galaxies are
usually interpreted as evidence for the DM halos, i.e. most of the DM lies in
the halos of galaxies in which the rotation curves are strictly flat Selwood
(1999). According to the current information for galaxy formation, every
galaxy forms and evolves within a DM halo, and it is expected that visible
galaxies are hosted by all halos. The aim of this paper is to point out the
main astrophysical consequences that arise by considering the modified
theories of gravity in terms of a different perspective of the DM included in
galactic halo. Specially, inspiring by some values for the EoS parameter which
is used to the ordinary matter, we will investigate the kind of DM in halos of
galaxies.

The outline of this paper is as follows. In the following section, we briefly
describe the motion of massive test particles under the conditions for
circular orbits. In Section III, we present the form of the modified gravity
theory $f(X)$ , where $X\in\\{R,\mathcal{G},T,...\\}$, as well as write down
its equations of motion for each of $f(R)$, $f(R,G)$ and $f(T)$ gravity
theories that yields constant rotational curves in the outskirts of galaxies.
Finally, we conclude with summary and discussion of the obtained results in
Section IV.

##  II The Motion of Test Particles

In order to obtain results which are relevant to the galactic dynamics, we
assume that the galactic halo has spherical symmetry and that dragging effects
on material particles (stars and dust) are inappreciable. Therefore, we
restrict our study to the static and spherically symmetric metric. The most
general static and spherically symmetric metric can be written as

${\rm d}s^{2}=-A(r){\rm d}t^{2}+B(r){\rm d}r^{2}+d\Omega^{2}\,,$ (1)

where $d\Omega^{2}=C(r)\left({\rm d}\theta^{2}+\sin^{2}\theta{\rm
d}\varphi^{2}\right)$ and ($r$,$\theta$,$\varphi$) are the spherical
coordinates. Then, the equations of motion for a test particle in the space-
time (1) can be derived from the Lagrangian

$2\mathcal{L}=-A(r)\dot{t}^{2}+B(r)\dot{r}^{2}+C(r)\left[\dot{\theta}^{2}+\sin\theta^{2}\dot{\varphi}^{2}\right],$
(2)

where a dot means derivative with respect to the proper time. From the above
Lagrangian (2), the generalized momenta become

$\displaystyle p_{t}=-E=-A\dot{t}\,,\,\,p_{r}=B\,\dot{r}\,,$ $\displaystyle
p_{\theta}=L_{\theta}=C\,\dot{\theta}\,,\,\,p_{\varphi}=L_{\varphi}=C\,\sin^{2}\theta\,\dot{\varphi},$
(3)

where $E$ is the total energy of a test particle and $L_{i}$ are the
components of its angular momentum. Using
$L^{2}=L_{\theta}^{2}+L_{\varphi}^{2}/\sin^{2}\theta$ which is the first
integral corresponding to the squared total angular momentum, the norm of the
four-velocity ($g_{\mu\nu}u^{\mu}u^{\nu}=-1$) yields

$E^{2}=V(r)+AB\,\dot{r}^{2}\,,$ (4)

where $V(r)\equiv A(1+L^{2}/C)$ is the effective potential. Thus, the
conditions for circular orbits $\dot{r}=0$ and $\partial_{r}V(r)=0$ lead to

$E^{2}=\frac{A^{2}\,C^{\prime}}{A\,C^{\prime}-C\,A^{\prime}},\qquad
L^{2}=\frac{C^{2}\,A^{\prime}}{A\,C^{\prime}-C\,A^{\prime}}\,.$ (5)

The definition of tangential velocity of a test particle is given by Böhmer et
al. (2008)

$v_{tg}^{2}=\frac{C}{A}\left(\frac{d\Omega}{dt}\right)^{2}=\frac{C}{A}\frac{\dot{\Omega}^{2}}{\dot{t}^{2}}=\frac{A}{C}\frac{L^{2}}{E^{2}}\,.$
(6)

Then, using the constants of motion (5), it follows that the tangential
velocity of this test particle is

$v_{tg}^{2}=\frac{A^{\prime}/A}{C^{\prime}/C}\,,$ (7)

which has the form

$v_{tg}^{2}=\frac{rA^{\prime}}{2A}\,,$ (8)

for $C(r)=r^{2}$, where a prime represents derivative with respect to $r$.
Since the motion of test particles is defined via the geodesic equations, the
relation (7) is independent of any theory of gravity. In the flat rotation
curves region, where $v_{tg}\approx{\rm constant}$, integration of Eq.(7)
gives

$A(r)=A_{0}C(r)^{v_{tg}^{2}}\,,$ (9)

which takes the form

$A(r)=A_{0}r^{\ell}\,,$ (10)

for $C(r)=r^{2}$, where $A_{0}$ is an integration constant and
$\ell=2v_{tg}^{2}$. Thus, the metric (1) with the relation (9) describes the
geometry of the space-time for the DM dominated regions where the test
particles move in constant rotational curves Matos et al. 2000b ; Böhmer et
al. (2008).

##  III Modified Gravity Models for Flat Rotation Curves Region

The generic action that we will consider for modified gravity theories reads

$\mathcal{S}=\int{d^{4}x\sqrt{-g}\left[\frac{1}{2\kappa^{2}}f(X)+\mathcal{L}_{\rm
m}\right]}\,,$ (11)

where $\kappa^{2}=8\pi G_{N}$ is the standard gravitational coupling,
$g=\det(g_{\mu\nu})$, $\mathcal{L}_{\rm m}$ is the ordinary matter Lagrangian
and the function $f$ depends on $X\in\\{R,\mathcal{G},T,...\\}$. We also point
out that the galaxy is embedded in a modified spherically symmetric geometry,
generated by the non-zero contributions of the modified gravitational action.
Hereafter, in order to obtain effects of the DM in the halos of galaxies, we
have assumed that there exists no the luminous matter, i.e. $\mathcal{L}_{\rm
m}=0$. Before obtaining the field equations by varying the Lagrangian
according to the coefficients of metric (1), we note the fact that the metric
variable $B$ does not contributes to dynamics by using a point-like Lagrangian
approach which implies that the Lagrangian does not depends on the generalized
velocity $B^{\prime}$, but the equation of motion due to the variation of $B$
has to be considered as a further constraint equation Capozziello et al. 2007a
. Under the latter fact, the action for the geometrical part has the form

$\mathcal{S}_{f(X)}=\int\mathcal{L}_{f(X)}(A,B,C,X,A^{\prime},C^{\prime},X^{\prime})dr\,.$
(12)

Then, taking variation of this action with respect to the metric coefficients
$A,B$ and $C$ of the static spherically symmetric space-time (1), one can
write the field equations as

$\displaystyle F_{1}(q^{i},q^{\prime
i},A^{\prime\prime},C^{\prime\prime})=\rho_{{}_{X}}\,,$ (13) $\displaystyle
F_{2}(q^{i},q^{\prime i},A^{\prime\prime},C^{\prime\prime})=p_{{}_{X}}\,,$
(14) $\displaystyle p^{\prime}_{{}_{X}}+F_{3}(q^{i},q^{\prime
i})(\rho_{{}_{X}}+p_{{}_{X}})=Q(q^{i},q^{\prime i},X,X^{\prime})\,,$ (15)

where $F_{i}$ are the functions of the metric coefficients $q^{i}=\\{A,B,C\\}$
and their first and/or second derivatives, $\rho_{{}_{X}}(q^{i},q^{\prime
i},X,X^{\prime},X^{\prime\prime})$ and $p_{{}_{X}}(q^{i},q^{\prime
i},X,X^{\prime},X^{\prime\prime})$ are the density and pressure for the
modified gravity theory, respectively. Then, inspiring by some values for the
EoS parameter used in the ordinary matter, the geometric EoS defined by
$w_{X}\equiv p_{X}/\rho_{X}$ has capacity to inform us about the kind of DM in
the galactic halos. For instance, if the modified gravity model gives rise to
the geometrical EoS parameter as $w_{X}=1/3$, then it may be corresponding to
the dark radiation Ackerman et al. (2009) inside the galactic halo. Now, we
will take three examples to explain the above theoretical construction of the
extended gravity model that yields constant rotational curves in the outskirts
of galaxies even though it can be given so much examples of the idea
represented here.

  * •

$f(R)$ gravity: Considering the background metric (1), the point-like
Lagrangian for $f(R)$ gravity is given by Capozziello et al. 2007a

$\displaystyle\mathcal{L}_{f(R)}=f_{R}\left(\frac{A^{\prime}C^{\prime}}{\sqrt{AB}}+\sqrt{\frac{A}{B}}\frac{C^{\prime
2}}{2C}\right)+f_{RR}R^{\prime}\left(\frac{CA^{\prime}}{\sqrt{AB}}+2\sqrt{\frac{A}{B}}C^{\prime}\right)+\sqrt{AB}\left[Cf+(2-CR)f_{R}\right].$
(16)

After varying (16) with respect to $A,B$ and $C$, the field equations have the
form of (13)-(15), which read

$\displaystyle\frac{1}{B}\left(\frac{B^{\prime}C^{\prime}}{BC}-\frac{2C^{\prime\prime}}{C}+\frac{C^{\prime
2}}{2C^{2}}+\frac{2B}{C}\right)=\rho_{{}_{R}},$ (17)
$\displaystyle\frac{1}{B}\left(\frac{B^{\prime}C^{\prime}}{BC}+\frac{C^{\prime
2}}{2C^{2}}-\frac{2B}{C}\right)=p_{{}_{R}},$ (18) $\displaystyle
f_{R}\Big{[}\frac{A^{\prime\prime}}{A}+\frac{C^{\prime\prime}}{C}-\frac{A^{\prime}B^{\prime}}{2AB}+\frac{C^{\prime}}{2C}\left(\frac{A^{\prime}}{A}-\frac{B^{\prime}}{B}\right)-\frac{A^{\prime
2}}{2A^{2}}-\frac{C^{\prime
2}}{2C^{2}}\Big{]}+f^{\prime}_{R}\left(\frac{A^{\prime}}{A}-\frac{B^{\prime}}{B}+\frac{C^{\prime}}{C}\right)$
$\displaystyle\qquad+2f^{\prime\prime}_{R}+B(Rf_{R}-f)=0,$ (19)

where $\rho_{{}_{R}}$ and $p_{{}_{R}}$ are defined by

$\displaystyle\rho_{{}_{R}}\equiv\frac{1}{f_{R}}\Big{(}Rf_{R}-f+\frac{1}{B}\left[f^{\prime}_{R}\left(\frac{2C^{\prime}}{C}-\frac{B^{\prime}}{B}\right)+2f^{\prime\prime}_{R}\right]\Big{)}\,,$
(20) $\displaystyle
p_{{}_{R}}\equiv\frac{1}{f_{R}}\left[f-Rf_{R}-\frac{f^{\prime}_{R}}{B}\left(\frac{A^{\prime}}{A}+\frac{2C^{\prime}}{C}\right)\right]\,,$
(21)

where $f^{\prime}_{R}=f_{RR}R^{\prime}$ and
$f^{\prime\prime}_{R}=f_{RR}R^{\prime\prime}+f_{RRR}R^{\prime 2}$. For an
exact solution of the power law gravity $f(R)=f_{0}R^{n}$, the metric
coefficients are given by

$\displaystyle A(r)=A_{0}r^{\ell}\,,\qquad
B(r)=\frac{(1+2n-2n^{2})(4n^{2}-10n+7)}{(n-2)^{2}}\,,\qquad C(r)=r^{2}\,,$
(22)

with $\ell=2(n-1)(2n-1)/(2-n)$. Therefore, the tangential velocity for this
case becomes

$v_{tg}^{2}=\frac{(n-1)(2n-1)}{(2-n)},$ (23)

due to the relation $\ell=2v_{tg}^{2}$. Here, the density and pressure of
$R^{n}$ gravity have the form

$\rho_{{}_{R}}(r)=\frac{\rho_{0}}{r^{2}}\,,\qquad
p_{{}_{R}}(r)=\frac{p_{0}}{r^{2}}\,,$ (24)

where $\rho_{0}$ and $p_{0}$ are found as

$\displaystyle\rho_{0}=\frac{2(n-1)(8n^{3}-20n^{2}+11n+3)}{(2n^{2}-2n-1)(4n^{2}-10n+7)}\,,$
(25) $\displaystyle
p_{0}=\frac{2(1-n)(8n^{3}-24n^{2}+21n-1)}{(2n^{2}-2n-1)(4n^{2}-10n+7)}.$ (26)

The relation (23) and the form of $\rho_{{}_{R}}$ and $p_{{}_{R}}$ explicitly
represent that the constant tangential velocity $v_{tg}$, the density
$\rho_{{}_{R}}$ and pressure $p_{{}_{R}}$ vanish as $n=1$, which is the GR
case, that is, it can be concluded that the GR does not give any information
about the DM dominated region of galaxies. Furthermore, the geometric EoS
$w_{{}_{R}}=p_{{}_{R}}/\rho_{{}_{R}}$ is constant, and found to be

$w_{{}_{R}}=-\frac{(8n^{3}-24n^{2}+21n-1)}{(8n^{3}-20n^{2}+11n+3)}\,.$ (27)

  * •

$f(R,\mathcal{G})$ gravity: This theory of gravity is a generalization of
$f(R)$ gravity containing a Gauss$-$Bonnet scalar $\mathcal{G}$ in the action
(11). For this case, the point-like Lagrangian for the static spherical
symmetric metric (1) takes the form Bahamonde et al. (2020)

$\displaystyle\mathcal{L}_{f(R,\mathcal{G})}=\mathcal{L}_{f(R)}+\frac{f_{\mathcal{G}\mathcal{G}}\mathcal{G}^{\prime}}{\sqrt{AB}}\left(4A^{\prime}-\frac{A^{\prime}C^{\prime
2}}{BC}\right)-\mathcal{G}\sqrt{AB}f_{\mathcal{G}}\,,$ (28)

where $\mathcal{L}_{f(R)}$ is as given by (16). Then, solving the field
equations obtained from the above Lagrangian for
$f(R,\mathcal{G})=f_{0}R+f_{1}\sqrt{\mathcal{G}}$ gravity, it is found the
following exact spherically symmetric solution Bahamonde et al. (2020) :

$A(r)=A_{0}r^{\ell},\qquad B(r)=\frac{4-\ell^{2}}{4},\quad C(r)=r^{2}\,,$ (29)

where $f_{1}=-2f_{0}\ell\sqrt{2\ell(\ell-2)}/(\ell^{2}-2\ell+8)$, and the
density and pressure are

$\rho_{{}_{R\mathcal{G}}}(r)=\frac{2\ell^{2}}{(\ell^{2}-4)\,r^{2}},\quad
p_{{}_{R\mathcal{G}}}(r)=-\frac{2\ell(\ell+4)}{(\ell^{2}-4)\,r^{2}}\,,$ (30)

with the geometric EoS $w_{{}_{R\mathcal{G}}}=-1-4/\ell$ which reduces to
$w_{{}_{R\mathcal{G}}}=-1-2v_{tg}^{-2}$ since the relation $\ell=2v_{tg}^{2}$.

  * •

$f(T)$ gravity:
In this gravity theory, the gravitational contributions to the metric tensor
become a source of torsion $T$ rather than curvature $R$, i.e. the connection
considered in this theory of gravity is distinct from the regular Levi-Civita
connection which is replaced with its Weitzenböck analog. In order to be
consistent with the study Bahamonde and Camci (2019), we will use the
signature $(+,-,-,-)$ for spherically symmetric metric which has the form
$ds^{2}=A(r)^{2}dt^{2}-B(r)^{2}dr^{2}-C(r)^{2}(d\theta^{2}+\sin\theta^{2}d\varphi^{2})$.
The Lagrangian $\mathcal{L}_{f(T)}$ for the latter metric is given by
Bahamonde and Camci (2019)

$\displaystyle\mathcal{L}_{f(T)}=AB\left[C^{2}(f-Tf_{T})+2f_{T}\right]+2f_{T}AC\Big{[}\frac{AC^{\prime}}{B}\left(\frac{2A^{\prime}}{A}+\frac{C^{\prime}}{C}\right)-2\left(\frac{A^{\prime}}{A}+\frac{C^{\prime}}{C}\right)\Big{]}\,.\qquad$
(31)

After varying the above Lagrangian with respect to $A,B$ and $C$, the field
equations for power-law $f(T)=f_{0}T^{n}$ gravity can be written in the form
of (13)-(15) as follows

$\displaystyle\frac{2}{B^{2}}\left(\frac{2B^{\prime}C^{\prime}}{BC}-\frac{2C^{\prime\prime}}{C}-\frac{C^{\prime
2}}{C^{2}}+\frac{B^{2}}{C^{2}}\right)=\rho_{{}_{T}},$ (32)
$\displaystyle\frac{2}{B^{2}}\left(\frac{2A^{\prime}C^{\prime}}{AC}+\frac{C^{\prime
2}}{C^{2}}-\frac{B^{2}}{C^{2}}\right)=p_{{}_{T}},$ (33) $\displaystyle
p^{\prime}_{{}_{T}}+\left(\frac{A^{\prime}}{A}+\frac{C^{\prime}}{C}\right)(\rho_{{}_{T}}+p_{{}_{T}})+4(n-1)\frac{A^{\prime}C^{\prime}T^{\prime}}{AB^{2}CT}=0\,,$
(34)

where $\rho_{{}_{T}}$ and $p_{{}_{T}}$ are defined as

$\displaystyle\rho_{{}_{T}}\equiv(n-1)\left[\frac{T}{n}+\frac{4}{CB^{2}}\frac{T^{\prime}}{T}\left(C^{\prime}-B\right)\right]\,,\qquad
p_{{}_{T}}\equiv\frac{(1-n)}{n}T\,,$ (35)

which are the torsion contributions to energy density and pressure. If $n=1$,
then the theory becomes the Teleparallel Equivalent of General Relativity
(TEGR), and the energy density and pressure given by (35) vanish, just as
expected. Using the field equations for $f(T)=f_{0}T^{n}$ gravity, where
$n\neq 0,1,\frac{1}{2},\frac{5}{6},\frac{5}{4},\frac{3}{2}$, it is reported
the following exact solution in Ref. Bahamonde and Camci (2019)

$\displaystyle A(r)=A_{0}r^{\ell}\,,$ $\displaystyle
B(r)=\frac{(2n-1)(4n-1)}{(4n^{2}-8n+5)}\,,$ (36) $\displaystyle C(r)=r\,,$

where $\ell=4n(n-1)(2n-3)/(4n^{2}-8n+5)$. Here it is found the relation
$\ell=v_{tg}^{2}$ because of the form of metric coefficients in this case,
which yields that the tangential velocity has the form

$v_{tg}^{2}=\frac{4n(n-1)(2n-3)}{4n^{2}-8n+5}.$ (37)

For this solution, the density and pressure of $T^{n}$ gravity become

$\rho_{{}_{T}}(r)=\frac{\rho_{1}}{r^{2}}\,,\qquad
p_{{}_{T}}(r)=\frac{p_{1}}{r^{2}}\,,$ (38)

where $\rho_{1}$ and $p_{1}$ are obtained as

$\displaystyle\rho_{1}=\frac{8n(n-1)(2n-3)(6n-5)}{(2n-1)^{2}(4n-5)^{2}},\qquad
p_{1}=\frac{8n(n-1)(2n-3)^{2}}{(4n-5)(2n-1)^{2}}\,,$ (39)

which gives a constant geometric EoS $w_{{}_{T}}=p_{1}/\rho_{1}$ as

$w_{{}_{T}}=\frac{(2n-3)(4n-5)}{6n-5}.$ (40)

—————————-

Figure 1:  Plots showing the tangential velocity $v_{tg}$ (left panel) and the
geometrical equation of state $w_{{}_{R}}$ (right panel) given in (27) vs the
power $n$ of $R^{n}$ gravity where
$n_{0}=-0.1967807237,\,n_{1}=(1-\sqrt{3})/2$ and $n_{2}=(1+\sqrt{3})/2$.

—————————-

##  IV Discussions and conclusions

For the three examples considered in the previous section, the metric
coefficients have the form

$A(r)=A_{0}r^{\ell},\quad B(r)=B_{0}={\rm constant}\,,$ (41)

and $d\Omega^{2}=r^{2}\left({\rm d}\theta^{2}+\sin^{2}\theta{\rm
d}\varphi^{2}\right)$. In all examples for the considered $f(X)$ theory of
gravity, where $X\in\\{R,\mathcal{G},T\\}$, the geometric density $\rho_{X}$
and pressure $p_{X}$ with a spheroidal profile, i.e. $C(r)=r^{2}$, are found
as isothermal, which means that they are proportional with distance as
$1/r^{2}$, and the geometrical EoS $w_{X}$ is constant. In the first and third
examples discussed above section, we have shown that there is a connection
between power of the modified gravity theory and the corresponding geometric
EoS parameter describing the intrinsic property of the DM halos of galaxies.
In the second example, we have taken an extension of GR in which an extra term
included proportional with the square root of Gauss-Bonnet scalar
$\mathcal{G}$, and found that the geometric EoS is directly related with flat
rotational velocity of galaxies.

In order to test the consistency of results, it has been considered the
Newtonian limit of the model by Böhmer et al. (2008), in which they have
restricted their analysis to the constant velocity region, and found that the
model has a well-defined Newtonian limit. In the Newtonian limit the $g_{\rm
tt}$ component of the metric tensor has the form $A\approx 1+2\Phi_{N}$, where
$\Phi_{N}$ is the Newtonian gravitational potential, and it is obtained for
the constant velocity regions around galaxies in the limit of large $r$ as
follows Böhmer et al. (2008) :

$\Phi_{N}(r)\approx v_{\rm tg}^{2}\ln\left(\frac{r}{r_{0}}\right)\,,$ (42)

where we take $A_{0}$ and $\ell$ as $A_{0}=r_{0}^{-\ell}$ and $\ell=2v_{\rm
tg}^{2}$ in (41). Thus, Eq. (42) represents that the correction term to the
Newtonian potential in the DM dominated region, where the rotation curves are
strictly flat, must have a logarithmic dependence on the radial coordinate
$r$.

All test particles in stable circular motion move at the speed of light when
$v_{tg}=1$, but this gives rise to a contradiction by observations at the
galactic scale. Furthermore, the tangential velocity $v_{tg}$ tends to zero in
the limit of large $r$. Thus, the tangential velocity $v_{tg}$ has to be at
the interval $0<v_{tg}<1$. During early epochs of the universe, the DM
velocity is not so small than the speed of light, for example the relative
velocity is $v_{rel}\propto 0.3c$ at freeze-out epoch Tulin and Yu (2018).
While at later times such as the DM halos today and during the recombination
epoch, the DM velocity is very smaller. For instance, the tangential
velocities of DM halos in dwarf galaxies, spiral galaxies and galaxy clusters
are approximately proportional with $10^{-5},10^{-3}$ and $10^{-2}$ in units
of $c=1$ at distances large enough from the galactic center, respectively. In
the first and third examples, it is not only found all the spectrum of flat
rotation curves but also the geometric EoS for the power law $R^{n}$ and
$T^{n}$ theories of gravity at the interval $0<v_{tg}<1$. Now we conclude our
findings in three examples as follows:

  * •

$f(R)=f_{0}R^{n}$ gravity: The relation between $v_{tg}$ and the power $n$,
the Eq. (23), is a second order algebraic equation for $n$, and it has a
solution as

$n=\frac{1}{4}\left(-v_{tg}^{2}+3\pm\sqrt{v_{tg}^{4}+10v_{tg}^{2}+1}\,\right)\,.$
(43)

So, we can exactly calculate $n$ for some specific tangential velocities. For
instance, if $v_{tg}=10^{-3}$, which is the rotational velocity of spiral
galaxies, then we find $n=1.000001$ or $n=0.4999985$. For $n=1.000001$, the
geometrical density and pressure coefficients $\rho_{0}$ and $p_{0}$ become
$\rho_{0}=-4\times 10^{-6}$ and $p_{0}=8\times 10^{-6}$, while for
$n=0.4999985$ we find $\rho_{0}=1.000001$ and $p_{0}=-0.999999$. The latter
result is most important in the sense that the galactic halos for spiral
galaxies prefer almost dark energy with a geometric EoS
$w_{{}_{R}}=-0.999998$, which indicates that $n=0.4999985$ for the flat
rotational velocity region. A general picture for the relation between
$v_{tg}$ and the power $n$ is shown in Fig. 1 which gives us some important
information on the DM halos of galaxies. It is seen from this figure that
$n\in[n_{1},0.5)\cup(1,n_{2}]$. For $n=-0.0578084147$ and $n=0.0504825997$ we
arrive the dust DM ($w_{{}_{R}}=0$) and stiff DM ($w_{{}_{R}}=1$),
respectively. Also, for $n=n_{1}$ and $n=n_{2}$ that gives $v_{tg}=1$, the
geometric EoS yields $w_{{}_{R}}=-3$ which is an EoS parameter for the phantom
energy ($w_{{}_{R}}<-1$). At the interval $n_{1}<n<0.1321342251$, one get the
quintessence EoS parameter values, which are in the range
$-1<w_{{}_{R}}<-1/3$. In this model, the dark radiation ($w_{{}_{R}}=1/3$)
appears if $n=0$, for which $v_{tg}=0.7071067812$. For
$n=n_{0}=-0.1967807237$, the tangential velocity is $v_{tg}=0.8713186295$ and
$w_{{}_{R}}$ tends to the infinity, which represents that there is a phase
transition at this evolution stage of galaxy. Using the observed best-fit
values of $v_{\rm flat}$ to the flat rotation curves of sample galaxies Frank
et al. (2016); Carignan and Purton (1998); Kamada et al. (2017), we present
the tangential velocity of test particles, the corresponding power $n$, and
the equation of state parameter $w_{R}$ for $R^{n}$ gravity in Table I.

Capozziello et al. (2007b) have investigated the possibility that the observed
flatness of the rotation curves of spiral galaxies for $R^{n}$ gravity is not
evidence for the existence of DM, in which they have obtained the corrected
gravitational potential as

$\Phi(r)=-\frac{Gm}{2\,r}\left[1+\left(\frac{r}{r_{c}}\right)^{\beta}\right]\,,$
(44)

where $m$ is a pointlike mass, $r_{c}$ is an arbitrary parameter, and $\beta$
is a parameter depending on the power $n$. The above potential is still
asymptotically decreasing, but the corrected rotation curve is higher than the
Newtonian one and not flat. Because of that the Newtonian potential for the DM
dominated region should have a logarithmic dependence such as Eq. (42), the
correction term to the Newtonian potential in the DM dominated region, where
the rotation curves are flat, does not appear in Eq. (44). The difference in
the Newtonian limit in our model results in different values of the parameter
$n$ for $R^{n}$ gravity.

Table 1:  Properties of sample galaxies in $R^{n}$ gravity. An explanation of the quantities appeared in the columns: $v_{\rm flat}$ \- galactic flat rotation velocity data, $v_{\rm tg}$-tangential velocity, $c$-velocity of light, $n$-power of $R^{n}$, $w_{R}$\- equation of state parameter.  Galaxy | $v_{\rm flat}({\rm km\,s^{-1}})$ |  $v_{\rm tg}=v_{\rm flat}/c$ ( $10^{-3}$ ) | $n$ | $w_{R}=p_{0}/\rho_{0}$
---|---|---|---|---
Milky Way | $200$ | $0.667$ | $0.4999993327$ | $-0.9999991107$
|  |  | $1.000000445$ | $-2.000001558$
NGC2841 | $310$ | $1.034$ | $0.4999983963$ | $-0.9999978618$
|  |  | $1.000001069$ | $-2.000003747$
NGC3521 | $231$ | $0.77$ | $0.4999991107$ | $-0.9999988142$
|  |  | $1.000000593$ | $-2.000002081$
NGC3198 | $152$ | $0.507$ | $0.4999996144$ | $-0.9999994856$
|  |  | $1.000000257$ | $-2.000000894$
NGC2403 | $97$ | $0.324$ | $0.4999998424$ | $-0.99999979$
|  |  | $1.000000105$ | $-2.000000368$
UGC128 | $140$ | $0.47$ | $0.4999996687$ | $-0.9999995582$
|  |  | $1.000000221$ | $-2.000000768$
UGC5721 | $80$ | $0.267$ | $0.4999998931$ | $-0.9999998578$
|  |  | $1.000000071$ | $-2.000000244$
UGCA442 | $57$ | $0.19$ | $0.4999999459$ | $-0.9999999282$
|  |  | $1.000000036$ | $-2.000000121$
DDO154 | $45$ | $0.15$ | $0.4999999663$ | $-0.9999999551$
|  |  | $1.000000022$ | $-2.000000072$
  * •

$f(R,\mathcal{G})=f_{1}R+f_{2}\sqrt{\mathcal{G}}$ gravity:

In Ref. Matos et al. 2000b , it is constructed an exact solution of Einstein’s
field equations sourced by a scalar field assuming the flat rotation curve
condition for the galactic halo. Taking into account this extended theory of
gravity, it is worth notice that we have derived the same solution of Matos et
al. 2000b , which is given by (29). Furthermore, we introduced a geometric EoS
of the form $w_{{}_{R\mathcal{G}}}=-1-2/v_{tg}^{2}$. This gives us that the
halo has always the property $w_{{}_{R\mathcal{G}}}<-1$, which means that the
galactic halo should only be filled with phantom energy, and therefore there
is not exist the DM in the galactic halo for this extended theory of gravity.

—————————-

Figure 2:  Plots showing the tangential velocity $v_{tg}$ (upper left panel),
the Eq.(37), and the geometrical equation of state $w_{{}_{T}}$ (upper right
panel) given in (40) vs the power $n$ of $T^{n}$ gravity, where
$n_{0}=(5-\sqrt{5})/4,\,n_{1}=5/6$ and $n_{2}=(5+\sqrt{5})/4$. Since the
shapes at the upper right panel are tiny in the ranges of $0<n<0.5$ and
$1.5<n<n_{2}$, they are redrawn in the bottom two figures.

—————————-

  * •

$f(T)=f_{0}T^{n}$ gravity:

For this theory of gravity, the relation between $v_{tg}$ and the power $n$,
the Eq. (37), is a cubic equation for $n$ as the form

$n^{3}-a_{1}n^{2}+a_{2}n+a_{3}=0,$ (45)

where
$a_{1}=\frac{1}{2}\left(v_{tg}^{2}+5\right),\,a_{2}=\frac{1}{2}\left(2v_{tg}^{2}+3\right)$
and $a_{3}=-\frac{5}{8}v_{tg}^{2}$. This cubic equation can be simplified by
making the substitution $n=x+a_{1}/3$. In terms of the new variable $x$, Eq.
(45) then becomes $x^{3}+3Px-2Q=0$, where $P=(3a_{2}-a_{1}^{2})/9$ and
$Q=(2a_{1}^{3}-9a_{1}a_{2}-27a_{3})/54$. Defining the polynomial discriminant
$D=Q^{2}+P^{3}$, we can solve algebraically the latter cubic equation. If
$D>0$, one of the roots is real and the other two roots are complex
conjugates. If $D<0$, all roots are real and unequal. In the latter case,
defining $y=\arccos\left(Q/\sqrt{-P^{3}}\right)$, then the real valued
solutions of (45) are of the form
$n_{k}=\frac{a_{1}}{3}+2\sqrt{-P}\cos\left(\frac{2\pi
k}{3}+\frac{y}{3}\right)$, where $k\in\\{0,1,2\\}$ and $P\leq 0$. In addition
to the restrictions on $n$ such that $n\neq
0,1,\frac{1}{2},\frac{5}{6},\frac{5}{4},\frac{3}{2}$, we have additional
property of $n$ due to $v_{tg}\in(0,1)$ as
$n\in(0,\frac{1}{2})\cup(n_{0},n_{1})\cup(n_{1},1)\cup(\frac{3}{2},n_{2})$,
where $n_{0},n_{1}$ and $n_{2}$ are given in the caption of Fig. 2. For the
rotational velocity $v_{tg}=10^{-3}$ of spiral galaxies, the Eq. (45) has
three real roots $n_{1}=4.169\times 10^{-7},\,n_{2}=0.99999975$ and
$n_{3}=1.500000333$, and the roots $n_{2}$ and $n_{3}$ have the property such
that $w_{{}_{T}}=1.000003$ for $n_{2}$, which is very close to unity, and
$w_{{}_{T}}=1.67\times 10^{-7}$ for $n_{3}$, which is very close to zero. If
$n=5/2$, then $w_{{}_{T}}=1$ (stiff dark matter EoS parameter) and
$v_{tg}=\sqrt{3}$ which means that it exceeds the speed of light. Further, one
get the dark radiation EoS parameter, $w_{{}_{T}}=1/3$, if
$n=\frac{3}{2}\pm\frac{1}{\sqrt{6}}$ where both values are not in the interval
of $n\in(0,\frac{1}{2})\cup(n_{0},n_{1})\cup(n_{1},1)\cup(\frac{3}{2},n_{2})$.
So we conclude that neither stiff DM nor dust DM or dark radiation exist at
halos of galaxies in $T^{n}$ gravity. In Fig. 2, it is seen that
$w_{{}_{T}}<-1$ for $n\in(0,\frac{1}{2})\cup(n_{0},n_{1})$ which gives the
phantom energy region of galactic halos, and $w_{{}_{T}}>0$ for
$n\in(n_{1},1)\cup(\frac{3}{2},n_{2})$. In Table 2, we present the tangential
velocity of test particles, the corresponding power $n$, and the equation of
state parameter $w_{T}$ for $T^{n}$ gravity by using the observed best-fit
values of $v_{\rm flat}$ to the flat rotation curves of sample galaxies Frank
et al. (2016); Carignan and Purton (1998); Kamada et al. (2017).

In Ref. Finch and Said (2018), a model with $f(T)=T+\alpha T^{n}$ gravity was
used to reproduce galactic rotation curves, and found that in the
$n\rightarrow 1$ limit the results were not approaching the case of GR. They
have pointed out that the most promising region for $n$ is $1<n<\frac{3}{2}$,
but this range of $n$ is excluded in our work (see Fig.2). As claimed in Ref.
Bahamonde et al. (2020), the work in Ref. Finch and Said (2018) is problematic
in the sense that they used an incorrect perturbed solution derived from an
earlier paper Ruggiero and Radicella (2015), where the authors used incorrect
field equations in $f(T)$ gravity. In our study we have completely neglected
the effect of the baryonic matter on the space-time geometry. So, it would be
an interesting issue to analyze what could happen with including the baryonic
matter to the field equations for $f(T)$ gravity.

Table 2:  Properties of sample galaxies in $T^{n}$ gravity. An explanation of the quantities appeared in the columns: $v_{\rm flat}$ \- galactic flat rotation velocity data, $v_{\rm tg}$-tangential velocity, $c$-velocity of light, $n$-power of $T^{n}$, $w_{T}$\- equation of state parameter.  Galaxy | $v_{\rm flat}({\rm km\,s^{-1}})$ |  $v_{\rm tg}=v_{\rm flat}/c$ ( $10^{-3}$ ) | $n$ | $w_{T}=p_{1}/\rho_{1}$
---|---|---|---|---
Milky Way | $200$ | $0.667$ | $0.999999888$ | $1.000001344$
|  |  | $1.500000148$ | $7.4\times 10^{-8}$
|  |  | $1.8537\times 10^{-7}$ | $-2.999999851$
NGC2841 | $310$ | $1.034$ | $0.9999997327$ | $1.000003208$
|  |  | $1.500000356$ | $1.78\times 10^{-7}$
|  |  | $4.4548\times 10^{-7}$ | $-2.999999644$
NGC3521 | $231$ | $0.77$ | $0.9999998614$ | $1.000001783$
|  |  | $1.5000001783$ | $9.9\times 10^{-8}$
|  |  | $2.4704\times 10^{-7}$ | $-2.999999803$
NGC3198 | $152$ | $0.507$ | $0.9999999357$ | $1.000000772$
|  |  | $1.500000086$ | $4.3\times 10^{-8}$
|  |  | $1.071\times 10^{-7}$ | $-2.999999916$
NGC2403 | $97$ | $0.324$ | $0.9999999738$ | $1.000000314$
|  |  | $1.500000035$ | $1.75\times 10^{-8}$
|  |  | $4.374\times 10^{-8}$ | $-2.999999965$
UGC128 | $140$ | $0.47$ | $0.9999999448$ | $1.000000662$
|  |  | $1.500000074$ | $3.7\times 10^{-8}$
|  |  | $9.204\times 10^{-8}$ | $-2.999999927$
UGC5721 | $80$ | $0.267$ | $0.9999999822$ | $1.000000214$
|  |  | $1.500000024$ | $1.2\times 10^{-8}$
|  |  | $2.9704\times 10^{-8}$ | $-2.999999977$
UGCA442 | $57$ | $0.19$ | $0.999999991$ | $1.000000108$
|  |  | $1.500000012$ | $6\times 10^{-9}$
|  |  | $1.5042\times 10^{-8}$ | $-2.999999988$
DDO154 | $45$ | $0.15$ | $0.9999999944$ | $1.000000067$
|  |  | $1.500000007$ | $4\times 10^{-9}$
|  |  | $9.375\times 10^{-9}$ | $-2.999999992$

The cold dark matter (CDM) paradigm has been extremely successful in
reproducing expansion history and large-scale structure of the universe as
well as the observed DM halos of galaxies. In a CDM model, the DM in the
universe is arranged in DM halos of galaxies. It would be mentioned a
possibility that the CDM paradigm may break down on galactic scales. Supposing
the DM particles are warm, instead of cold, then it gives rise that they were
quasi-relativistic during kinetic decoupling from the thermal bath in the
early universe Bode et al. (2001). The results of this work are valid for
either cold or warm DMs.

In our approach, we need to mention that the energy density $\rho_{X}$ and the
pressure $p_{X}$, where $X\in\\{R,\mathcal{G},T\\}$, are effective geometric
quantities. So, these quantities may not have an exact physical interpretation
as for the similar physical fields. In this work, the modified gravity
theories such as $R^{n},T^{n}$ and
$f(R,\mathcal{G})=f_{1}R+f_{2}\sqrt{\mathcal{G}}$ have been introduced as a
possible way to explain the observed flat rotational velocities of galaxies
without the need of any DM component. It is seen from the above three examples
that the $R^{n}$ and $T^{n}$ gravities have more rich structures than the
$f(R,\mathcal{G})=f_{1}R+f_{2}\sqrt{\mathcal{G}}$ gravity. For future aim at
improving the halo model discussed here, it would be complementary to take
into account the ordinary matter in addition to the DM contribution in the
halos of galaxies.

## References

  * Riess et al. (1998) Riess, A. G., et al. [Supernova Search Team Collaboration]: Observational evidence from supernovae for an accelerating universe and a cosmological constant. Astrophys. J. 116, 1009 (1998)
  * Perlmutter et al. (1999) Perlmutter, S., et al. [Supernova Cosmology Project Collaboration]: Measurements of $\Omega$ and $\Lambda$ from 42 high-redshift supernovae. Astrophys. J. 517, 565 (1999)
  * Wechsler and Tinker (2018) Wechsler, R. H., Tinker, J. L.: The Connection Between Galaxies and Their Dark Matter Halos. Annu. Rev. Astron. Astrophys. 56, 435 (2018)
  * Bharadwaj and Kar (2003) Bharadwaj, S., Kar, S.: Modeling galaxy halos using dark matter with pressure. Phys. Rev. D 68, 023516 (2003)
  * (5) Matos, T., Guzmán, F. S., Ureña-López, L. A.: Scalar field as dark matter in the universe. Class. Quantum Grav. 17, 1707 (2000a)
  * (6) Matos, T., Guzmán, F. S., Núñez, D.: Spherical scalar field halo in galaxies. Phys. Rev. D 62, 061301 (R) (2000b)
  * Peeples (2000) Peebles, P. J. E.: Dynamics of a dark matter field with a quartic self-interaction potential. Astrophys. J. 534, L127 (2000) Phys. Rev. D 62, 023502 (2000)
  * Arbey et al. (2003) Arbey, A., Lesgourgues, J., Salati, P.: Galactic halos of fluid dark matter.Phys. Rev. D 68, 023511 (2003)
  * Soleng (1995) Soleng, H. H.: Dark matter and non-newtonian gravity from general relativity coupled to a fluid of strings. Gen. Relativ. Gravit. 27, 367 (1995)
  * Selwood (1999) Sellwood, J. A. Galaxy Dynamics: A Rutgers Symposium, eds David R. Merritt, M. Valluri, J. A. Sellwood, ASP Conference Series Volume 182, San Francisco, 1999; arXiv:astro-ph/9903184.
  * Böhmer et al. (2008) Böhmer, C. G., Harko, T., Lobo, F. S. N.: Dark matter as a geometric effect in $f(R)$ gravity. Astropart. Phys. 29, 386 (2008)
  * (12) Capozziello, S., Stabile, A., Troisi, A: Spherically symmetric solutions in $f(R)$ gravity via the Noether symmetry approach. Class. Quantum Grav. 24, 2153 (2007)
  * Ackerman et al. (2009) Ackerman, L., Buckley, M. R., Carroll, S. M., Kamionkowski, M.: Dark matter and dark radiation. Phys. Rev. D 79, 023519 (2009)
  * Bahamonde et al. (2020) Bahamonde, S., Dialektopoulas, K., Camci, U.: Exact Spherically Symmetric Solutions in Modified Gauss-Bonnet gravity from Noether symmetry approach. Symmetry 12, 68 (2020)
  * Bahamonde and Camci (2019) Bahamonde, S., Camci, U.: Exact Spherically Symmetric Solutions in Modified Teleparallel gravity. Symmetry 11, 1462 (2019)
  * Tulin and Yu (2018) Tulin, S., Yu, H.: Dark Matter Self-interactions and small scale structure. Phys. Rep. 730, 1 (2018)
  * Bode et al. (2001) Bode, P., Ostriker, J. P., Turok, N.: Halo Formation in warm dark matter models. Astrophys. J. 556, 93 (2001)
  * Frank et al. (2016) Frank, B. S., De Blok, W. J. G., Walter, F., Leroy, A., Carignan, C.: The impact of moplecular gas on mass models of nearby galaxies. Astron J. 151, 94 (2016)
  * Carignan and Purton (1998) Carignan, C., Purton, C.: The total mass of DDO154. Astrophys. J. 506, 125 (1998)
  * Kamada et al. (2017) Kamada, A., Kaplinghat, M., Pace, A. B., Yu, H.-B.: Self-Interacting Dark Matter Can Explain Diverse Galactic Rotation Curves. Phys. Rev. Lett. 119, 111102 (2017)
  * (21) Capozziello, S., Cardone, V.F., Troisi, A.: Low surface brightness galaxy rotation curves in the low energy limit of $R^{n}$ gravity: no need for dark matter? Mon. Not. R. Astron. Soc. 375, 1423 (2007)
  * Finch and Said (2018) Finch, A., Said, J. L.: Galactic rotation dynamics in $f(T)$ gravity. Eur. Phys. J. C 78, 560 (2018)
  * Bahamonde et al. (2020) S. Bahamonde, S., Said, J.L., Zubair, M.: Solar system tests in modified teleparallel gravity. J. Cosmol. Astropart. Phys. 10, 024 (2020)
  * Ruggiero and Radicella (2015) Ruggiero, M.L., Radicella, N.: Weak-field spherically symmetric solutiuons in $f(T)$ gravity. Phys. Rev. D 91, 104014 (2015)